# simulation of the real data for timecalibration
#
# Remo Dietlicher
# ETH Zürich / Semesterarbeit
#
#

import pyfact
from myhisto import *
from hist import *
import numpy as np
import numpy.random as rnd
from ROOT import *
from time import time
from optparse import OptionParser

NROI  = 1024 		# length of the DRS pipeline
NEvents  = 200	 	# Number of evaluated Event
fsampling = 1024./512.1 # sampling frequency
freq = 250.        	# testfrequency
P_nom = 1000./freq 	# nominal Period due to testfrequency
DAMPING = 0.02		# Damping applied to correction
NChip = 1		# Number of Chips

t_s = time()

parser = OptionParser()
parser.add_option('-e', '--extended', action = 'store_true', dest='extend', default = False )
(options, args) = parser.parse_args()

variante = "dat"
alg = "Remo"
save = "no"

if(options.extend):

	variante = raw_input("Welche Daten? (sim/dat): ")
	alg = raw_input("Welcher Algorithmus? (Remo/Oliver): ")
	NChip = int(raw_input("Wieviele Chips? ([1, ... ,160]): "))
	NEvents = int(raw_input("Wieviele Events? ([1, ... , 1000]) "))
	save = raw_input("Soll CellTime gespeichert werden? (yes/no) ")
	tag = raw_input("Soll ein spezieller tag an das File angehängt werden? ")
	
NPix = 9*NChip

h = hist_list(tcalsimHistograms, NChip, "Chip")
MeanXingHist = np.zeros(NROI)

def Datainit(variante, NEvents):

	if(variante == "sim"):
		y_off = 3.4
	
		# Start is the first bin being written on.
		Startp = rnd.random(NEvents)*(NROI - 2) 	# project [0, 1] onto int([0, lpipe-1])
		Start = range(NEvents)
		for i, val in enumerate(Startp):
			Start[i] = int(round(val))		# getting integer values needed in roll later
			
		ID = np.zeros(NROI)
		x = np.linspace(0., float(NROI), NROI+1)
		ID = np.sin(x*2*np.pi/NROI+np.pi)*5
		
		dt_true = np.zeros(NROI)
		for i in range(NROI):
			dt_true[i] = ID[i+1]-ID[i]
		
		versch = np.sum(dt_true)/1024
		dt_true = dt_true - versch
		print np.sum(dt_true)
	
		#create time-array, time in nanoseconds
		t_nom = np.ones(NROI)/fsampling
		t_true = t_nom + dt_true
	
		# simulation of the real pipeline
		rCellTime = np.zeros([NEvents, NROI])
		for Event in range(NEvents):
			for i in range(NROI-1):
				rCellTime[Event][i+1] = rCellTime[Event][i]+t_true[i]
			rCellTime[Event] = np.roll(rCellTime[Event],-Start[Event])
			rCellTime[-Start[Event]:] += NROI/fsampling
			
		print "t_true tot = ", np.sum(t_true)
		
		# calculation of the nominal pipeline
		def CellTimef(t_nom):
			CellTime = np.zeros(NROI+1)
			for i in range(NROI):
				CellTime[i+1] = CellTime[i] + t_nom[i]
		
			return CellTime	
			
		CellTime = CellTimef(t_nom)
		
		# Simulate Data
		Data = np.zeros([NEvents, NROI])
		for Event in range(NEvents):
			Data[Event] = np.sin(rCellTime[Event]*2.*np.pi/P_nom)+y_off
	
	if(variante == "dat"):
		t_true = 0
		Data = np.load("d1000.npy")
		Start = np.load("start_cells_1.npy")
		
		
	return Data, t_true, Start
	
Data, t_true, Start = Datainit(variante, NEvents)


# Functions needed for calibration
# calculation of the nominal pipeline
def CellTimef(t_nom):
	CellTime = np.zeros(NROI+1)
	for i in range(NROI):
		CellTime[i+1] = CellTime[i] + t_nom[i]
	
	return CellTime	
	
# getting t_nom from CellTime
def t_nomf(CellTime):
	t_nom = np.zeros(NROI)
	for i in range(NROI):
		t_nom[i] = CellTime[i+1]-CellTime[i]
		
	return t_nom	
	
		

# find negative-going crossings of mean value
def Crossing(Data, Mean, t_nom, Start, Chip):
	TimeXing = np.zeros(NROI)
	MeanXing = np.zeros(NROI)
	NumXing = 0
	CellTime = CellTimef(np.roll(t_nom, -Start))
	
	for i in range(NROI-1):
		if ((Data[i] > Mean) & (Data[i+1] < Mean)):
			MeanXing[NumXing] = i
			h.list[Chip].dict["lomeanxing"].Fill(np.mod(i+Start, NROI))
			
			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
	h.list[Chip].dict["meanfreq"].Fill(freqp)
	if(np.abs(freqp-freq) > 20):
		print "This is not a timing calibration run!!"
		raise SystemExit
						
	return MeanXing, TimeXing, NumXing
	
	
	
# find positive-going crossings of mean value
def CrossingPos(Data, Mean, t_nom, Start, Chip):
	TimeXing = np.zeros(NROI)
	MeanXing = np.zeros(NROI)
	NumXing = 0
	CellTime = CellTimef(np.roll(t_nom, -Start))
	
	for i in range(NROI-1):
		if ((Data[i] < Mean) & (Data[i+1] > Mean)):
			MeanXing[NumXing] = i
			h.list[Chip].dict["lomeanxing"].Fill(np.mod(i+Start, NROI))
			
			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
	h.list[Chip].dict["meanfreq"].Fill(freqp)
	if(np.abs(freqp-freq) > 20):
		print "This is not a timing calibration run!!"
		raise SystemExit
						
	return MeanXing, TimeXing, NumXing
	

	
def CrossingOliver(Data, Mean, t_nom, Start, Chip):

	TimeXing = np.zeros(NROI)
	MeanXing = np.zeros(NROI)
	NumXing = 0
	CellTime = CellTimef(t_nom)
	
	
	for i in range(NROI-1):
		if ((Data[i] > Mean) & (Data[i+1] < Mean)):
			MeanXing[NumXing] = np.mod(i+Start, NROI)
			
			FirstCell = CellTime[np.mod(i+Start,NROI)]
			SecondCell = CellTime[np.mod(i+1+Start, NROI)]
			
			if(SecondCell < FirstCell): SecondCell =+ NROI/fsampl
			
			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
	h.list[Chip].dict["meanfreq"].Fill(freqp)
	if(np.abs(freqp-freq) > 20):
		print "This is not a timing calibration 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(Event > 300):
			h.list[Chip].dict["period"].Fill(Period)
		
		
		if(np.abs(1-Period/P_nom) > 0.2):
			#print "Skipping zero crossing due to too large deviation of period."
			h.list[Chip].dict["loleftout"].Fill(MeanXing[i])
			h.list[Chip].dict["leftout"].Fill(np.mod(MeanXing[i]+Start, NROI))
			continue
			
		h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i])
		h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[i]+Start, NROI))
		MeanXingHist[np.mod(MeanXing[i]+Start, NROI)] += 1
		
		# calculate the frequency of the testwave
		freqp = 1000./Period
		h.list[Chip].dict["freq"].Fill(freqp)
		
		C = (P_nom-Period)*DAMPING
		Ctot[i] = P_nom-Period

		numnom = Period*fsampling
		Correction = np.zeros(NROI)
		Correction[I1+1:I2+1] += C/numnom
		Correction[:I1+1] -= (I2-I1)*C/(numnom*(NROI-(I2-I1)))
		Correction[I2+1:] -= (I2-I1)*C/(numnom*(NROI-(I2-I1)))
		"""
		Correction = np.zeros(NROI)
		Correction[I1+1:I2+1] += C/(I2-I1)
		Correction[:I1+1] -= C/(NROI-(I2-I1))
		Correction[I2+1:] -= C/(NROI-(I2-I1))

		numnom = Period*fsampling
		Correction = np.zeros(NROI)
		Correction[I1+1:I2+1] += C/numnom
		"""
			
		Correction = np.roll(Correction, Start)
		t_nom += Correction
	#t_nom = t_nom/np.sum(t_nom)*1024/fsampling
	
	if(Event > 200):
		h.list[Chip].dict["perioddevproj"].Fill(np.average(Ctot))
	h.list[Chip].dict["perioddev"].SetBinContent(Event+1, np.abs(np.average(Ctot)))	
	h.list[Chip].dict["lomeanxing"].Fill(MeanXing[NumXing])
	h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[NumXing]+Start, NROI))
	MeanXingHist[np.mod(MeanXing[NumXing]+Start, NROI)] += 1
		
			
	return t_nom

def CorrOliver(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
	
	CellTime = CellTimef(t_nom)
	for i in range(int(NumXing-1)):
		I1 = MeanXing[i]
		I2 = MeanXing[i+1]
		Period = TimeXing[i+1] - TimeXing[i]
		if(Period<0): Period += NROI/fsampling
		h.list[Chip].dict["period"].Fill(Period)
		
		if(np.abs(1-Period/P_nom) > 0.2):
			print "Skipping zero crossing due to too large deviation of period."
			h.list[Chip].dict["loleftout"].Fill(MeanXing[i])
			h.list[Chip].dict["leftout"].Fill(np.mod(MeanXing[i]+Start, NROI))
			continue
			
		# calculate the frequency of the testwave
		freqp = 1000./Period
		h.list[Chip].dict["freq"].Fill(freqp)
		
			
		h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i])
		h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[i]+Start, NROI))
		MeanXingHist[np.mod(MeanXing[i]+Start, NROI)] += 1
			
		Correction = (P_nom-Period)*DAMPING
		
		# Positive correction from MeanXing[i] to MeanXing[i+1]
		# and negative to all others
		Time = 0
		for j in range(NROI):
			diff = CellTime[j+1] - CellTime[j]
			if ((I2>I1 and j>I1 and j<=I2) or (I2<I1 and (j<=I2 or j>I1))):
				diff += Correction/np.mod(I2-I1+NROI, NROI)
			else:
				diff -= Correction/(NROI-np.mod(I2-I1+NROI, NROI))
			CellTime[j] = Time
			Time += diff
		CellTime[NROI] = Time

	t_nom = t_nomf(CellTime)
	
	h.list[Chip].dict["lomeanxing"].Fill(MeanXing[NumXing-1])	
	h.list[Chip].dict["meanxing"].Fill(np.mod(MeanXing[NumXing-1]+Start, NROI))
	
	return t_nom
	
	
# create array to put the calibrated times for each cell in
t_nom = np.ones(NROI)/fsampling
# 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
		if(alg == "Remo"):
			MeanXing, TimeXing, NumXing = Crossing(Data[Event], np.average(Data[Event]), t_nom, int(Start[Event]), Chip)
			h.list[Chip].dict["numxing"].Fill(NumXing)
			t_nom = Corr(MeanXing, NumXing, TimeXing, t_nom, int(Start[Event]), Chip, Event)
			h.list[Chip].dict["conv"].SetBinContent(Event+1, np.average(np.abs(t_nom-t_true)))
			

		if(alg == "Oliver"):
			MeanXing, TimeXing, NumXing = CrossingOliver(Data[Event], np.average(Data[Event]), t_nom, int(Start[Event]), Chip)
			h.list[Chip].dict["numxing"].Fill(NumXing)
			t_nom = CorrOliver(MeanXing, NumXing, TimeXing, t_nom, int(Start[Event]), Chip)
			h.list[Chip].dict["conv"].SetBinContent(Event+1, np.average(np.abs(t_nom-t_true)))
		
		Chip += 1

CellTime = CellTimef(t_nom)

print "t_tot = ", np.sum(t_nom)	

if(save == "yes"):
	np.save(str(alg)+"_"+str(variante)+"_"+str(NEvents)+"x"+str(NChip)+str(tag), CellTime)
	np.savetxt(str(alg)+"_"+str(variante)+"_"+str(NEvents)+"x"+str(NChip)+str(tag), CellTime[None], fmt="%7.4f")


for Chip in range(NChip):
	for i in range(NROI):
		h.list[Chip].dict["int_dev"].SetBinContent(i+1, CellTime[i]-i/fsampling)
		h.list[Chip].dict["diff_dev"].SetBinContent(i+1, t_nom[i]-1/fsampling)
		h.list[Chip].dict["data"].SetBinContent(i+1, Data[-1][i])
	
pyfact.SaveHistograms(h.list, "tcal_"+str(alg)+"_"+str(variante)+str(tag)+".root", "RECREATE")


t_f = time()
print "vergangene Zeit	= ",  (t_f - t_s)/60., " min."
	
print "Variante		= ", variante
print "Algorithmus		= ", alg
print "Speichern		= ", save
print "NChip			= ", NChip
print "NEvents			= ", NEvents
print "Histo saved as		= ", "tcal_"+str(alg)+"_"+str(variante)+str(tag)+".root"
if(save == "yes"):
	print "CellTime saved as	= ", str(alg)+"_"+str(variante)+"_"+str(NEvents)+"x"+str(NChip)+str(tag)+".npy"