| 1 | #!/usr/bin/python
|
|---|
| 2 | #
|
|---|
| 3 | # Remo Dietlicher
|
|---|
| 4 | # ETH Zurich
|
|---|
| 5 | # Semesterarbeit
|
|---|
| 6 | #
|
|---|
| 7 | import pyfact
|
|---|
| 8 | from myhisto import *
|
|---|
| 9 | from hist import *
|
|---|
| 10 | import numpy as np
|
|---|
| 11 | import numpy.random as rnd
|
|---|
| 12 | from scipy import interpolate as ip
|
|---|
| 13 | from ROOT import *
|
|---|
| 14 | from time import time
|
|---|
| 15 | from optparse import OptionParser
|
|---|
| 16 |
|
|---|
| 17 | def ftcal(NChip, NEventsp, h,
|
|---|
| 18 | dfname = '/data00/fact-construction/raw/2011/11/24/20111124_014.fits.gz',
|
|---|
| 19 | calfname = '/data00/fact-construction/raw/2011/11/24/20111124_012.drs.fits.gz'):
|
|---|
| 20 |
|
|---|
| 21 | rd = pyfact.rawdata( dfname, calfname )
|
|---|
| 22 |
|
|---|
| 23 | fsampling = 1024./512. # sampling frequency
|
|---|
| 24 | freq = 250. # testfrequency
|
|---|
| 25 | P_nom = 1000./freq # nominal Period due to testfrequency
|
|---|
| 26 | DAMPING = 0.02 # Damping applied to correction
|
|---|
| 27 |
|
|---|
| 28 |
|
|---|
| 29 | # Functions needed for calibration
|
|---|
| 30 | # find negative-going crossings of mean value
|
|---|
| 31 | def Crossing(Data, Mean, t_nom, Start, Chip):
|
|---|
| 32 | TimeXing = np.zeros(rd.NROI)
|
|---|
| 33 | MeanXing = np.zeros(rd.NROI)
|
|---|
| 34 | NumXing = 0
|
|---|
| 35 | CellTime = CellTimef(np.roll(t_nom, -Start))
|
|---|
| 36 |
|
|---|
| 37 | for i in range(rd.NROI-1):
|
|---|
| 38 | if ((Data[i] > Mean) & (Data[i+1] < Mean)):
|
|---|
| 39 | MeanXing[NumXing] = i
|
|---|
| 40 |
|
|---|
| 41 | FirstCell = CellTime[i]
|
|---|
| 42 | SecondCell = CellTime[i+1]
|
|---|
| 43 |
|
|---|
| 44 | TimeXing[NumXing] = FirstCell+(SecondCell-FirstCell)/(1.-Data[i+1]/(Data[i]))*(1.-Mean/(Data[i]))
|
|---|
| 45 | NumXing += 1
|
|---|
| 46 |
|
|---|
| 47 | # calculate the frequency of the testwave
|
|---|
| 48 | freqp = fsampling*NumXing*1000./1024
|
|---|
| 49 | if(np.abs(freqp-freq) > 20):
|
|---|
| 50 | print "This is not a timecalibration run!!"
|
|---|
| 51 | raise SystemExit
|
|---|
| 52 |
|
|---|
| 53 | return MeanXing, TimeXing, NumXing
|
|---|
| 54 |
|
|---|
| 55 | # Determine time between crossings and apply correction
|
|---|
| 56 | def Corr(MeanXing, NumXing, TimeXing, t_nom, Start, Chip, Event):
|
|---|
| 57 |
|
|---|
| 58 | Ctot = np.zeros(NumXing)
|
|---|
| 59 |
|
|---|
| 60 | # loop over mean crossings
|
|---|
| 61 | for i in range(int(NumXing-1)):
|
|---|
| 62 | I1 = MeanXing[i]
|
|---|
| 63 | I2 = MeanXing[i+1]
|
|---|
| 64 | Period = TimeXing[i+1] - TimeXing[i]
|
|---|
| 65 |
|
|---|
| 66 |
|
|---|
| 67 | if(np.abs(1-Period/P_nom) > 0.2):
|
|---|
| 68 | continue
|
|---|
| 69 |
|
|---|
| 70 |
|
|---|
| 71 | if(Event == NEvents-1):
|
|---|
| 72 | h.list[Chip].dict["periods"].Fill(Period)
|
|---|
| 73 | if(Event == 0):
|
|---|
| 74 | h.list[Chip].dict["periods0"].Fill(Period)
|
|---|
| 75 |
|
|---|
| 76 |
|
|---|
| 77 |
|
|---|
| 78 | # calculate the frequency of the testwave
|
|---|
| 79 | freqp = 1000./Period
|
|---|
| 80 |
|
|---|
| 81 | C = (P_nom-Period)*DAMPING
|
|---|
| 82 | Ctot[i] = P_nom-Period
|
|---|
| 83 |
|
|---|
| 84 |
|
|---|
| 85 | numnom = Period*fsampling
|
|---|
| 86 | Correction = np.zeros(rd.NROI)
|
|---|
| 87 | Correction[I1+1:I2+1] += C/numnom
|
|---|
| 88 | Correction[:I1+1] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
|
|---|
| 89 | Correction[I2+1:] -= (I2-I1)*C/(numnom*(rd.NROI-(I2-I1)))
|
|---|
| 90 |
|
|---|
| 91 | Correction = np.roll(Correction, Start)
|
|---|
| 92 | t_nom += Correction
|
|---|
| 93 |
|
|---|
| 94 |
|
|---|
| 95 | if(Event > 200):
|
|---|
| 96 | h.list[Chip].dict["perioddevproj"].Fill(np.average(Ctot))
|
|---|
| 97 | h.list[Chip].dict["perioddev"].SetBinContent(Event+1, np.average(np.abs(Ctot)))
|
|---|
| 98 |
|
|---|
| 99 | return t_nom
|
|---|
| 100 |
|
|---|
| 101 |
|
|---|
| 102 | # Summing up the nominal times for each cell to the total cell time.
|
|---|
| 103 | def CellTimef(t_nom):
|
|---|
| 104 | CellTime = np.zeros(rd.NROI+1)
|
|---|
| 105 | for i in range(rd.NROI):
|
|---|
| 106 | CellTime[i+1] = CellTime[i] + t_nom[i]
|
|---|
| 107 |
|
|---|
| 108 | return CellTime
|
|---|
| 109 |
|
|---|
| 110 | # getting t_nom from CellTime
|
|---|
| 111 | def t_nomf(CellTime):
|
|---|
| 112 | t_nom = np.zeros(rd.NROI)
|
|---|
| 113 | for i in range(rd.NROI):
|
|---|
| 114 | t_nom[i] = CellTime[i+1]-CellTime[i]
|
|---|
| 115 |
|
|---|
| 116 | return t_nom
|
|---|
| 117 |
|
|---|
| 118 |
|
|---|
| 119 | rd.next()
|
|---|
| 120 | NPIX = NChip*9
|
|---|
| 121 | NEvents = NEventsp
|
|---|
| 122 |
|
|---|
| 123 | # create array to put the calibrated times for each cell in
|
|---|
| 124 | t_nom = np.ones([NPIX, rd.NROI])/fsampling
|
|---|
| 125 | DataTot = np.zeros([NEvents, NChip, rd.NROI])
|
|---|
| 126 |
|
|---|
| 127 | # loop over events
|
|---|
| 128 | for Event in range( NEvents ):
|
|---|
| 129 | print "Event ", Event
|
|---|
| 130 | Chip = 0
|
|---|
| 131 |
|
|---|
| 132 | # loop over every 9th channel in Chip
|
|---|
| 133 | for Pix in range(NPIX)[8::9]:
|
|---|
| 134 | print "Pix ", Pix
|
|---|
| 135 |
|
|---|
| 136 | # calculate t_nom for each pixel
|
|---|
| 137 | MeanXing, TimeXing, NumXing = Crossing(rd.acalData[Pix], np.average(rd.acalData[Pix]), t_nom[Pix], rd.startCells[Pix], Chip)
|
|---|
| 138 | t_nom[Pix] = Corr(MeanXing, NumXing, TimeXing, t_nom[Pix], rd.startCells[Pix], Chip, Event)
|
|---|
| 139 |
|
|---|
| 140 | DataTot[Event][Chip] = np.roll(rd.acalData[Pix], rd.startCells[Pix])
|
|---|
| 141 |
|
|---|
| 142 | Chip += 1
|
|---|
| 143 |
|
|---|
| 144 | rd.next()
|
|---|
| 145 |
|
|---|
| 146 | CellTime = np.zeros([NPIX, rd.NROI+1])
|
|---|
| 147 |
|
|---|
| 148 | DataMean = np.average(DataTot, 0)
|
|---|
| 149 |
|
|---|
| 150 | for Pix in range(NPIX):
|
|---|
| 151 | CellTime[Pix] = CellTimef(t_nom[Pix])
|
|---|
| 152 |
|
|---|
| 153 | t_nom = t_nom[8::9]
|
|---|
| 154 | CellTime = CellTime[8::9]
|
|---|
| 155 |
|
|---|
| 156 |
|
|---|
| 157 | return CellTime, DataMean
|
|---|
| 158 |
|
|---|
