#!/usr/bin/python
#
# Remo Dietlicher
# ETH Zurich
# Semesterarbeit
#
import pyfact
from myhisto import *
from hist import *
import numpy as np
import numpy.random as rnd
from scipy import interpolate as ip
from ROOT import *
from time import time
from optparse import OptionParser

def ftcal(NChip, NEventsp, h,
	      dfname = '/data00/fact-construction/raw/2011/11/24/20111124_014.fits.gz', 
		  calfname = '/data00/fact-construction/raw/2011/11/24/20111124_012.drs.fits.gz'):

	rd = pyfact.rawdata( dfname, calfname )

	fsampling = 1024./512.	# sampling frequency
	freq = 250.        		# testfrequency
	P_nom = 1000./freq 		# nominal Period due to testfrequency
	DAMPING = 0.02			# Damping applied to correction
	

	# Functions needed for calibration
	# find negative-going crossings of mean value
	def Crossing(Data, Mean, t_nom, Start, Chip):
		TimeXing = np.zeros(rd.NROI)
		MeanXing = np.zeros(rd.NROI)
		NumXing = 0
		CellTime = CellTimef(np.roll(t_nom, -Start))

		for i in range(rd.NROI-1):
			if ((Data[i] > Mean) & (Data[i+1] < Mean)):
				MeanXing[NumXing] = i

				FirstCell = CellTime[i]
				SecondCell = CellTime[i+1]

				TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
				NumXing += 1

		# calculate the frequency of the testwave
		freqp = fsampling*NumXing*1000./1024
		if(np.abs(freqp-freq) > 20):
			print "This is not a timecalibration run!!"
			raise SystemExit

		return MeanXing, TimeXing, NumXing

	# Determine time between crossings and apply correction
	def Corr(MeanXing, NumXing, TimeXing, t_nom, Start, Chip, Event):

		Ctot = np.zeros(NumXing)
		
		# loop over mean crossings
		for i in range(int(NumXing-1)):
			I1 = MeanXing[i]
			I2 = MeanXing[i+1]
			Period = TimeXing[i+1] - TimeXing[i]


			if(np.abs(1-Period/P_nom) > 0.2):
				continue
			

			if(Event == NEvents-1):
				h.list[Chip].dict["periods"].Fill(Period)
			if(Event == 0):
				h.list[Chip].dict["periods0"].Fill(Period)



			# calculate the frequency of the testwave
			freqp = 1000./Period

			C = (P_nom-Period)*DAMPING
			Ctot[i] = P_nom-Period
			
			
			numnom = Period*fsampling
			Correction = np.zeros(rd.NROI)
			Correction[I1+1:I2+1] += C/numnom
			Correction[:I1+1] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
			Correction[I2+1:] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))

			Correction = np.roll(Correction, Start)
			t_nom += Correction
			
			
			if(Event > 200):
				h.list[Chip].dict["perioddevproj"].Fill(np.average(Ctot))
			h.list[Chip].dict["perioddev"].SetBinContent(Event+1, np.average(np.abs(Ctot)))	

		return t_nom


	# Summing up the nominal times for each cell to the total cell time.
	def CellTimef(t_nom):
		CellTime = np.zeros(rd.NROI+1)
		for i in range(rd.NROI):
			CellTime[i+1] = CellTime[i] + t_nom[i]

		return CellTime

	# getting t_nom from CellTime
	def t_nomf(CellTime):
		t_nom = np.zeros(rd.NROI)
		for i in range(rd.NROI):
			t_nom[i] = CellTime[i+1]-CellTime[i]

		return t_nom


	rd.next()
	NPIX = NChip*9
	NEvents = NEventsp

	# create array to put the calibrated times for each cell in
	t_nom = np.ones([NPIX, rd.NROI])/fsampling
	DataTot = np.zeros([NEvents, NChip, rd.NROI])
	
	# loop over events
	for Event in range( NEvents ):
		print "Event ", Event
		Chip = 0

		# loop over every 9th channel in Chip
		for Pix in range(NPIX)[8::9]:
			print "Pix ", Pix

			# calculate t_nom for each pixel
			MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip)
			t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip, Event)
			
			DataTot[Event][Chip] = np.roll(rd.acalData[Pix], rd.startCells[Pix])
			
			Chip += 1

		rd.next()
	
	CellTime = np.zeros([NPIX, rd.NROI+1])
	
	DataMean = np.average(DataTot, 0)
	
	for Pix in range(NPIX):
		CellTime[Pix] = CellTimef(t_nom[Pix])

	t_nom = t_nom[8::9]
	CellTime = CellTime[8::9]


	return CellTime, DataMean