| 1 | #!/usr/bin/python
|
|---|
| 2 | #
|
|---|
| 3 | # Werner Lustermann
|
|---|
| 4 | # ETH Zurich
|
|---|
| 5 | #
|
|---|
| 6 | import pyfact
|
|---|
| 7 | from myhisto import *
|
|---|
| 8 | from hist import *
|
|---|
| 9 | import numpy as np
|
|---|
| 10 | import numpy.random as rnd
|
|---|
| 11 | from scipy import interpolate as ip
|
|---|
| 12 | from ROOT import *
|
|---|
| 13 | from time import time
|
|---|
| 14 | from optparse import OptionParser
|
|---|
| 15 |
|
|---|
| 16 | t_s = time()
|
|---|
| 17 |
|
|---|
| 18 | import os.path
|
|---|
| 19 | import sys
|
|---|
| 20 | #dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits'
|
|---|
| 21 | #calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits'
|
|---|
| 22 |
|
|---|
| 23 |
|
|---|
| 24 | #dfname = '/data00/fact-construction/raw/2011/11/24/20111124_113.fits.gz'
|
|---|
| 25 | #calfname = '/data00/fact-construction/raw/2011/11/24/20111124_111.drs.fits.gz'
|
|---|
| 26 |
|
|---|
| 27 | dfname = '/data00/fact-construction/raw/2012/02/23/20120223_205.fits.gz'
|
|---|
| 28 | calfname = '/data00/fact-construction/raw/2012/02/23/20120223_202.drs.fits.gz'
|
|---|
| 29 |
|
|---|
| 30 | if not os.path.isfile(dfname):
|
|---|
| 31 | print dfname, 'is not a file!!!!'
|
|---|
| 32 | sys.exit(1)
|
|---|
| 33 | if not os.path.isfile(calfname):
|
|---|
| 34 | print calfname, 'is not a file'
|
|---|
| 35 | sys.exit(2)
|
|---|
| 36 |
|
|---|
| 37 | rd = pyfact.rawdata( dfname, calfname )
|
|---|
| 38 |
|
|---|
| 39 | fsampling = 2. # sampling frequency
|
|---|
| 40 | freq = 250. # testfrequency
|
|---|
| 41 | P_nom = 1000./freq # nominal Period due to testfrequency
|
|---|
| 42 | DAMPING = 0.02 # Damping applied to correction
|
|---|
| 43 | NChip = 1 # Number of Chips
|
|---|
| 44 | NEvents = 300
|
|---|
| 45 |
|
|---|
| 46 |
|
|---|
| 47 | parser = OptionParser()
|
|---|
| 48 | parser.add_option('-e', '--extended', action = 'store_true', dest='extend', default = False )
|
|---|
| 49 | (options, args) = parser.parse_args()
|
|---|
| 50 |
|
|---|
| 51 | save = "yes"
|
|---|
| 52 |
|
|---|
| 53 | if(options.extend):
|
|---|
| 54 |
|
|---|
| 55 | NChip = int(raw_input("Wieviele Chips? ([1, ... , "+str(rd.NPIX/9)+"]): "))
|
|---|
| 56 | NEventsp = int(raw_input("Wieviele Events? ([1, ... , "+str(rd.NEvents)+"]) "))
|
|---|
| 57 | save = raw_input("Soll CellTime gespeichert werden? (yes/no) ")
|
|---|
| 58 |
|
|---|
| 59 |
|
|---|
| 60 | h = hist_list(tcalHistograms, NChip, "Chip")
|
|---|
| 61 |
|
|---|
| 62 |
|
|---|
| 63 | # Functions needed for calibration
|
|---|
| 64 |
|
|---|
| 65 | # find negative-going crossings of mean value
|
|---|
| 66 | def Crossing(Data, Mean, t_nom, Start, Chip):
|
|---|
| 67 | TimeXing = np.zeros(rd.NROI)
|
|---|
| 68 | MeanXing = np.zeros(rd.NROI)
|
|---|
| 69 | NumXing = 0
|
|---|
| 70 | CellTime = CellTimef(np.roll(t_nom, -Start))
|
|---|
| 71 |
|
|---|
| 72 | for i in range(rd.NROI-1):
|
|---|
| 73 | if ((Data[i] > Mean) & (Data[i+1] < Mean)):
|
|---|
| 74 | MeanXing[NumXing] = i
|
|---|
| 75 |
|
|---|
| 76 | FirstCell = CellTime[i]
|
|---|
| 77 | SecondCell = CellTime[i+1]
|
|---|
| 78 |
|
|---|
| 79 | TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
|
|---|
| 80 | NumXing += 1
|
|---|
| 81 |
|
|---|
| 82 | # calculate the frequency of the testwave
|
|---|
| 83 | freqp = fsampling*NumXing*1000./1024
|
|---|
| 84 | if(np.abs(freqp-freq) > 20):
|
|---|
| 85 | print "This is not a timecalibration run!!"
|
|---|
| 86 | raise SystemExit
|
|---|
| 87 |
|
|---|
| 88 | return MeanXing, TimeXing, NumXing
|
|---|
| 89 |
|
|---|
| 90 | def CrossingOliver(Data, Mean, t_nom, Start, Chip):
|
|---|
| 91 | TimeXing = np.zeros(lpipe)
|
|---|
| 92 | MeanXing = np.zeros(lpipe)
|
|---|
| 93 | NumXing = 0
|
|---|
| 94 | CellTime = CellTimef(t_nom)
|
|---|
| 95 |
|
|---|
| 96 |
|
|---|
| 97 | for i in range(lpipe-1):
|
|---|
| 98 | if ((Data[i] > Mean) & (Data[i+1] < Mean)):
|
|---|
| 99 | MeanXing[NumXing] = np.mod(i+Start, lpipe)
|
|---|
| 100 |
|
|---|
| 101 | FirstCell = CellTime[np.mod(i+Start,lpipe)]
|
|---|
| 102 | SecondCell = CellTime[np.mod(i+1+Start, lpipe)]
|
|---|
| 103 |
|
|---|
| 104 | if(SecondCell < FirstCell): SecondCell =+ lpipe/fsampl
|
|---|
| 105 |
|
|---|
| 106 | TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
|
|---|
| 107 | NumXing += 1
|
|---|
| 108 |
|
|---|
| 109 |
|
|---|
| 110 | return MeanXing, TimeXing, NumXing
|
|---|
| 111 |
|
|---|
| 112 | # Determine time between crossings and apply correction
|
|---|
| 113 | def Corr(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
|
|---|
| 114 |
|
|---|
| 115 | # loop over mean crossings
|
|---|
| 116 | for i in range(int(NumXing-1)):
|
|---|
| 117 | I1 = MeanXing[i]
|
|---|
| 118 | I2 = MeanXing[i+1]
|
|---|
| 119 | Period = TimeXing[i+1] - TimeXing[i]
|
|---|
| 120 |
|
|---|
| 121 |
|
|---|
| 122 | if(np.abs(1-Period/P_nom) > 0.2):
|
|---|
| 123 | #print "Skipping zero crossing due to too large deviation of period."
|
|---|
| 124 | h.list[Chip].dict["loleftout"].Fill(MeanXing[i])
|
|---|
| 125 | #h.list[Chip].dict["leftout"].Fill(np.mod((MeanXing[i]+Start), rd.NROI))
|
|---|
| 126 | continue
|
|---|
| 127 |
|
|---|
| 128 | #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[i])
|
|---|
| 129 | h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[i]+Start), rd.NROI))
|
|---|
| 130 |
|
|---|
| 131 | # calculate the frequency of the testwave
|
|---|
| 132 | freqp = 1000./Period
|
|---|
| 133 | h.list[Chip].dict["freq"].Fill(freqp)
|
|---|
| 134 |
|
|---|
| 135 | C = (P_nom-Period)*DAMPING
|
|---|
| 136 | numnom = Period*fsampling
|
|---|
| 137 | Correction = np.zeros(rd.NROI)
|
|---|
| 138 | Correction[I1+1:I2+1] += C/numnom
|
|---|
| 139 | Correction[:I1+1] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
|
|---|
| 140 | Correction[I2+1:] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
|
|---|
| 141 |
|
|---|
| 142 | Correction = np.roll(Correction, Start)
|
|---|
| 143 | t_nom += Correction
|
|---|
| 144 | #h.list[Chip].dict["lomeanxing"].Fill(MeanXing[-1])
|
|---|
| 145 | h.list[Chip].dict["meanxing"].Fill(np.mod((MeanXing[-1]+Start), rd.NROI))
|
|---|
| 146 |
|
|---|
| 147 | return t_nom
|
|---|
| 148 |
|
|---|
| 149 | def CorrOliver(MeanXing, NumXing, TimeXing, t_nom, Start, Chip):
|
|---|
| 150 |
|
|---|
| 151 | CellTime = CellTimef(t_nom)
|
|---|
| 152 | for i in range(int(NumXing-1)):
|
|---|
| 153 | I1 = MeanXing[i]
|
|---|
| 154 | I2 = MeanXing[i+1]
|
|---|
| 155 | Period = TimeXing[i+1] - TimeXing[i]
|
|---|
| 156 | if(Period<0): Period += rd.NROI/fsampl
|
|---|
| 157 |
|
|---|
| 158 | if(np.abs(1-Period/P_nom) > 0.2):
|
|---|
| 159 | #print "Skipping zero crossing due to too large deviation of period."
|
|---|
| 160 | continue
|
|---|
| 161 |
|
|---|
| 162 | Correction = (P_nom-Period)*DAMPING
|
|---|
| 163 |
|
|---|
| 164 | # Positive correction from MeanXing[i] to MeanXing[i+1]
|
|---|
| 165 | # and negative to all others
|
|---|
| 166 | Time = 0
|
|---|
| 167 | for j in range(rd.NROI):
|
|---|
| 168 | diff = CellTime[j+1] - CellTime[j]
|
|---|
| 169 | if ((I2>I1 and j>I1 and j<=I2) or (I2<I1 and (j<=I2 or j>I1))):
|
|---|
| 170 | diff += Correction/np.mod(I2-I1+rd.NROI, rd.NROI)
|
|---|
| 171 | else:
|
|---|
| 172 | diff -= Correction/(rd.NROI-np.mod(I2-I1+rd.NROI, rd.NROI))
|
|---|
| 173 | CellTime[j] = Time
|
|---|
| 174 | Time += diff
|
|---|
| 175 | CellTime[rd.NROI] = Time
|
|---|
| 176 |
|
|---|
| 177 | print CellTime[-1]
|
|---|
| 178 | t_nom = t_nomf(CellTime)
|
|---|
| 179 |
|
|---|
| 180 | return t_nom
|
|---|
| 181 |
|
|---|
| 182 |
|
|---|
| 183 |
|
|---|
| 184 | # Summing up the nominal times for each cell to the total cell time.
|
|---|
| 185 | def CellTimef(t_nom):
|
|---|
| 186 | CellTime = np.zeros(rd.NROI+1)
|
|---|
| 187 | for i in range(rd.NROI):
|
|---|
| 188 | CellTime[i+1] = CellTime[i] + t_nom[i]
|
|---|
| 189 |
|
|---|
| 190 | return CellTime
|
|---|
| 191 |
|
|---|
| 192 | # getting t_nom from CellTime
|
|---|
| 193 | def t_nomf(CellTime):
|
|---|
| 194 | t_nom = np.zeros(rd.NROI)
|
|---|
| 195 | for i in range(rd.NROI):
|
|---|
| 196 | t_nom[i] = CellTime[i+1]-CellTime[i]
|
|---|
| 197 |
|
|---|
| 198 | return t_nom
|
|---|
| 199 |
|
|---|
| 200 |
|
|---|
| 201 | rd.next()
|
|---|
| 202 | NPIX = NChip*9
|
|---|
| 203 | NEvents = NEventsp
|
|---|
| 204 |
|
|---|
| 205 | # create array to put the calibrated times for each cell in
|
|---|
| 206 | t_nom = np.ones([NPIX, rd.NROI])/fsampling
|
|---|
| 207 | # loop over events
|
|---|
| 208 | for Event in range( NEvents ):
|
|---|
| 209 | print "Event ", Event
|
|---|
| 210 | Chip = 0
|
|---|
| 211 |
|
|---|
| 212 | # loop over every 9th channel in Chip
|
|---|
| 213 | for Pix in range(NPIX)[8::9]:
|
|---|
| 214 | print "Pix ", Pix
|
|---|
| 215 |
|
|---|
| 216 | # calculate t_nom for each pixel
|
|---|
| 217 | MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip)
|
|---|
| 218 | h.list[Chip].dict["numxing"].Fill(NumXing)
|
|---|
| 219 | #h.list[Chip].dict["start"].Fill(rd.startCells[Pix])
|
|---|
| 220 | t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip)
|
|---|
| 221 |
|
|---|
| 222 | Chip += 1
|
|---|
| 223 |
|
|---|
| 224 | rd.next()
|
|---|
| 225 |
|
|---|
| 226 | t_nom = t_nom[8::9]
|
|---|
| 227 |
|
|---|
| 228 | if(save == "yes"):
|
|---|
| 229 | np.save(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
|
|---|
| 230 | np.savetxt(str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip), t_nom)
|
|---|
| 231 |
|
|---|
| 232 | for Chip in range(NChip):
|
|---|
| 233 | for i in range(rd.NROI):
|
|---|
| 234 | h.list[Chip].dict["int_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][:i])-i*1/fsampling)
|
|---|
| 235 | h.list[Chip].dict["diff_dev"].SetBinContent(i+1, np.sum(t_nom[Chip][i] - 1/fsampling))
|
|---|
| 236 |
|
|---|
| 237 | pyfact.SaveHistograms(h.list, str(dfname[-20:-8])+"_"+str(NEvents)+"x"+str(NChip)+".root", "RECREATE")
|
|---|
| 238 |
|
|---|
| 239 |
|
|---|
| 240 | t_f = time()
|
|---|
| 241 |
|
|---|
| 242 |
|
|---|
| 243 |
|
|---|
| 244 | print "time ellapsed = ", (t_f - t_s)/60., " min."
|
|---|
| 245 | print "Speichern = ", save
|
|---|
| 246 | print "NChip = ", NChip
|
|---|
| 247 | print "NEvents = ", NEventsp
|
|---|
| 248 | print "Histo saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".root"
|
|---|
| 249 | if(save == "yes"):
|
|---|
| 250 | print "CellTime saved as = ", str(dfname[-20:-8])+str(NEvents)+"x"+str(NChip)+".npy"
|
|---|
| 251 |
|
|---|
