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