]>
Commit | Line | Data |
---|---|---|
ee299369 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
ee299369 | 17 | |
e5bbbc4e | 18 | //_________________________________________________________________________ |
19 | // Utility Class for handling Raw data | |
20 | // Does all transitions from Digits to Raw and vice versa, | |
21 | // for simu and reconstruction | |
22 | // | |
23 | // Note: the current version is still simplified. Only | |
24 | // one raw signal per digit is generated; either high-gain or low-gain | |
25 | // Need to add concurrent high and low-gain info in the future | |
26 | // No pedestal is added to the raw signal. | |
ee299369 | 27 | //*-- Author: Marco van Leeuwen (LBL) |
e5bbbc4e | 28 | |
ee299369 | 29 | #include "AliEMCALRawUtils.h" |
21cad85c | 30 | |
ee299369 | 31 | #include "TF1.h" |
32 | #include "TGraph.h" | |
e5bbbc4e | 33 | class TSystem; |
21cad85c | 34 | |
e5bbbc4e | 35 | class AliLog; |
72c58de0 | 36 | #include "AliRun.h" |
ee299369 | 37 | #include "AliRunLoader.h" |
e5bbbc4e | 38 | class AliCaloAltroMapping; |
ee299369 | 39 | #include "AliAltroBuffer.h" |
40 | #include "AliRawReader.h" | |
32cd4c24 | 41 | #include "AliCaloRawStreamV3.h" |
ee299369 | 42 | #include "AliDAQ.h" |
21cad85c | 43 | |
feedcab9 | 44 | #include "AliEMCALRecParam.h" |
ee299369 | 45 | #include "AliEMCALLoader.h" |
46 | #include "AliEMCALGeometry.h" | |
e5bbbc4e | 47 | class AliEMCALDigitizer; |
ee299369 | 48 | #include "AliEMCALDigit.h" |
20b636fc | 49 | #include "AliEMCAL.h" |
5e3106bc | 50 | #include "AliCaloCalibPedestal.h" |
21cad85c | 51 | |
ee299369 | 52 | ClassImp(AliEMCALRawUtils) |
21cad85c | 53 | |
ee299369 | 54 | // Signal shape parameters |
89d338a6 | 55 | Int_t AliEMCALRawUtils::fgTimeBins = 256; // number of sampling bins of the raw RO signal (we typically use 15-50; theoretical max is 1k+) |
e5bbbc4e | 56 | Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns |
09974781 | 57 | Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec |
ee299369 | 58 | |
59 | // some digitization constants | |
60 | Int_t AliEMCALRawUtils::fgThreshold = 1; | |
61 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule | |
e5bbbc4e | 62 | Int_t AliEMCALRawUtils::fgPedestalValue = 32; // pedestal value for digits2raw |
63 | Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled) | |
ee299369 | 64 | |
8cb998bd | 65 | AliEMCALRawUtils::AliEMCALRawUtils() |
b4133f05 | 66 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), |
67 | fNPedSamples(0), fGeom(0), fOption("") | |
8cb998bd | 68 | { |
b4133f05 | 69 | |
70 | //These are default parameters. | |
71 | //Can be re-set from without with setter functions | |
ee299369 | 72 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) |
b4133f05 | 73 | fOrder = 2; // order of gamma fn |
74 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
7643e728 | 75 | fNoiseThreshold = 3; // 3 ADC counts is approx. noise level |
76 | fNPedSamples = 4; // less than this value => likely pedestal samples | |
65bdc82f | 77 | |
78 | //Get Mapping RCU files from the AliEMCALRecParam | |
79 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
80 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
81 | ||
21cad85c | 82 | for(Int_t i = 0; i < 4; i++) { |
65bdc82f | 83 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
84 | } | |
85 | ||
72c58de0 | 86 | //To make sure we match with the geometry in a simulation file, |
87 | //let's try to get it first. If not, take the default geometry | |
33c3c91a | 88 | AliRunLoader *rl = AliRunLoader::Instance(); |
c61fe3b4 | 89 | if(!rl) AliError("Cannot find RunLoader!"); |
72c58de0 | 90 | if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) { |
91 | fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); | |
92 | } else { | |
93 | AliInfo(Form("Using default geometry in raw reco")); | |
937d0661 | 94 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); |
65bdc82f | 95 | } |
96 | ||
72c58de0 | 97 | if(!fGeom) AliFatal(Form("Could not get geometry!")); |
98 | ||
65bdc82f | 99 | } |
100 | ||
101 | //____________________________________________________________________________ | |
5544799a | 102 | AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry) |
103 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), | |
104 | fNPedSamples(0), fGeom(pGeometry), fOption("") | |
105 | { | |
106 | // | |
107 | // Initialize with the given geometry - constructor required by HLT | |
108 | // HLT does not use/support AliRunLoader(s) instances | |
109 | // This is a minimum intervention solution | |
110 | // Comment by MPloskon@lbl.gov | |
111 | // | |
112 | ||
113 | //These are default parameters. | |
114 | //Can be re-set from without with setter functions | |
115 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) | |
116 | fOrder = 2; // order of gamma fn | |
117 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
7643e728 | 118 | fNoiseThreshold = 3; // 3 ADC counts is approx. noise level |
119 | fNPedSamples = 4; // less than this value => likely pedestal samples | |
5544799a | 120 | |
121 | //Get Mapping RCU files from the AliEMCALRecParam | |
122 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
123 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
124 | ||
21cad85c | 125 | for(Int_t i = 0; i < 4; i++) { |
5544799a | 126 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
127 | } | |
128 | ||
129 | if(!fGeom) AliFatal(Form("Could not get geometry!")); | |
130 | ||
131 | } | |
132 | ||
133 | //____________________________________________________________________________ | |
65bdc82f | 134 | AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) |
135 | : TObject(), | |
136 | fHighLowGainFactor(rawU.fHighLowGainFactor), | |
b4133f05 | 137 | fOrder(rawU.fOrder), |
138 | fTau(rawU.fTau), | |
139 | fNoiseThreshold(rawU.fNoiseThreshold), | |
140 | fNPedSamples(rawU.fNPedSamples), | |
65bdc82f | 141 | fGeom(rawU.fGeom), |
142 | fOption(rawU.fOption) | |
143 | { | |
144 | //copy ctor | |
145 | fMapping[0] = rawU.fMapping[0]; | |
146 | fMapping[1] = rawU.fMapping[1]; | |
21cad85c | 147 | fMapping[2] = rawU.fMapping[2]; |
148 | fMapping[3] = rawU.fMapping[3]; | |
65bdc82f | 149 | } |
150 | ||
151 | //____________________________________________________________________________ | |
152 | AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) | |
153 | { | |
154 | //assignment operator | |
155 | ||
156 | if(this != &rawU) { | |
157 | fHighLowGainFactor = rawU.fHighLowGainFactor; | |
b4133f05 | 158 | fOrder = rawU.fOrder; |
159 | fTau = rawU.fTau; | |
160 | fNoiseThreshold = rawU.fNoiseThreshold; | |
161 | fNPedSamples = rawU.fNPedSamples; | |
65bdc82f | 162 | fGeom = rawU.fGeom; |
163 | fOption = rawU.fOption; | |
164 | fMapping[0] = rawU.fMapping[0]; | |
165 | fMapping[1] = rawU.fMapping[1]; | |
21cad85c | 166 | fMapping[2] = rawU.fMapping[2]; |
167 | fMapping[3] = rawU.fMapping[3]; | |
65bdc82f | 168 | } |
169 | ||
170 | return *this; | |
171 | ||
ee299369 | 172 | } |
65bdc82f | 173 | |
ee299369 | 174 | //____________________________________________________________________________ |
175 | AliEMCALRawUtils::~AliEMCALRawUtils() { | |
e5bbbc4e | 176 | //dtor |
65bdc82f | 177 | |
ee299369 | 178 | } |
65bdc82f | 179 | |
ee299369 | 180 | //____________________________________________________________________________ |
65bdc82f | 181 | void AliEMCALRawUtils::Digits2Raw() |
ee299369 | 182 | { |
183 | // convert digits of the current event to raw data | |
184 | ||
33c3c91a | 185 | AliRunLoader *rl = AliRunLoader::Instance(); |
ee299369 | 186 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); |
187 | ||
188 | // get the digits | |
189 | loader->LoadDigits("EMCAL"); | |
190 | loader->GetEvent(); | |
191 | TClonesArray* digits = loader->Digits() ; | |
192 | ||
193 | if (!digits) { | |
194 | Warning("Digits2Raw", "no digits found !"); | |
195 | return; | |
196 | } | |
65bdc82f | 197 | |
ee299369 | 198 | static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here |
199 | AliAltroBuffer* buffers[nDDL]; | |
200 | for (Int_t i=0; i < nDDL; i++) | |
201 | buffers[i] = 0; | |
202 | ||
e2c2134b | 203 | TArrayI adcValuesLow(fgTimeBins); |
204 | TArrayI adcValuesHigh(fgTimeBins); | |
ee299369 | 205 | |
ee299369 | 206 | // loop over digits (assume ordered digits) |
207 | for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { | |
208 | AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ; | |
209 | if (digit->GetAmp() < fgThreshold) | |
210 | continue; | |
211 | ||
212 | //get cell indices | |
213 | Int_t nSM = 0; | |
214 | Int_t nIphi = 0; | |
215 | Int_t nIeta = 0; | |
216 | Int_t iphi = 0; | |
217 | Int_t ieta = 0; | |
218 | Int_t nModule = 0; | |
65bdc82f | 219 | fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); |
220 | fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; | |
ee299369 | 221 | |
21cad85c | 222 | //Check which is the RCU, 0 or 1, of the cell. |
ee299369 | 223 | Int_t iRCU = -111; |
224 | //RCU0 | |
225 | if (0<=iphi&&iphi<8) iRCU=0; // first cable row | |
226 | else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; | |
227 | //second cable row | |
228 | //RCU1 | |
229 | else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; | |
230 | //second cable row | |
231 | else if(16<=iphi&&iphi<24) iRCU=1; // third cable row | |
21cad85c | 232 | |
233 | if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same | |
234 | ||
e36e3bcf | 235 | if (iRCU<0) |
236 | Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); | |
ee299369 | 237 | |
238 | //Which DDL? | |
239 | Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; | |
240 | if (iDDL >= nDDL) | |
241 | Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); | |
242 | ||
243 | if (buffers[iDDL] == 0) { | |
244 | // open new file and write dummy header | |
245 | TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); | |
21cad85c | 246 | //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C |
247 | Int_t iRCUside=iRCU+(nSM%2)*2; | |
248 | //iRCU=0 and even (0) SM -> RCU0A.data 0 | |
249 | //iRCU=1 and even (0) SM -> RCU1A.data 1 | |
250 | //iRCU=0 and odd (1) SM -> RCU0C.data 2 | |
251 | //iRCU=1 and odd (1) SM -> RCU1C.data 3 | |
252 | //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl; | |
253 | buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]); | |
ee299369 | 254 | buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; |
255 | } | |
256 | ||
257 | // out of time range signal (?) | |
258 | if (digit->GetTimeR() > GetRawFormatTimeMax() ) { | |
259 | AliInfo("Signal is out of time range.\n"); | |
260 | buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp()); | |
261 | buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin | |
262 | buffers[iDDL]->FillBuffer(3); // bunch length | |
263 | buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer | |
264 | // calculate the time response function | |
265 | } else { | |
e2c2134b | 266 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; |
ee299369 | 267 | if (lowgain) |
e2c2134b | 268 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold); |
ee299369 | 269 | else |
e2c2134b | 270 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold); |
ee299369 | 271 | } |
272 | } | |
273 | ||
274 | // write headers and close files | |
275 | for (Int_t i=0; i < nDDL; i++) { | |
276 | if (buffers[i]) { | |
277 | buffers[i]->Flush(); | |
278 | buffers[i]->WriteDataHeader(kFALSE, kFALSE); | |
279 | delete buffers[i]; | |
280 | } | |
281 | } | |
65bdc82f | 282 | |
ee299369 | 283 | loader->UnloadDigits(); |
284 | } | |
285 | ||
286 | //____________________________________________________________________________ | |
5e3106bc | 287 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, AliCaloCalibPedestal* pedbadmap) |
ee299369 | 288 | { |
65bdc82f | 289 | // convert raw data of the current event to digits |
ee299369 | 290 | |
c47157cd | 291 | digitsArr->Clear(); |
ee299369 | 292 | |
c47157cd | 293 | if (!digitsArr) { |
ee299369 | 294 | Error("Raw2Digits", "no digits found !"); |
295 | return; | |
296 | } | |
297 | if (!reader) { | |
298 | Error("Raw2Digits", "no raw reader found !"); | |
299 | return; | |
300 | } | |
301 | ||
32cd4c24 | 302 | AliCaloRawStreamV3 in(reader,"EMCAL",fMapping); |
ee299369 | 303 | // Select EMCAL DDL's; |
7643e728 | 304 | reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL |
feedcab9 | 305 | |
8cb998bd | 306 | //Updated fitting routine from 2007 beam test takes into account |
307 | //possibility of two peaks in data and selects first one for fitting | |
308 | //Also sets some of the starting parameters based on the shape of the | |
309 | //given raw signal being fit | |
310 | ||
311 | TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); | |
7643e728 | 312 | signalF->SetParameters(10.,5.,fTau,fOrder,0.); //set all defaults once, just to be safe |
8cb998bd | 313 | signalF->SetParNames("amp","t0","tau","N","ped"); |
7643e728 | 314 | signalF->FixParameter(2,fTau); // tau in units of time bin |
315 | signalF->FixParameter(3,fOrder); // order | |
48a56166 | 316 | |
ee299369 | 317 | Int_t id = -1; |
82cbdfca | 318 | Float_t time = 0. ; |
319 | Float_t amp = 0. ; | |
7643e728 | 320 | Float_t ped = 0. ; |
321 | Float_t ampEstimate = 0; | |
322 | Float_t timeEstimate = 0; | |
323 | Float_t pedEstimate = 0; | |
32cd4c24 | 324 | Int_t i = 0; |
325 | Int_t startBin = 0; | |
ee299369 | 326 | |
8cb998bd | 327 | //Graph to hold data we will fit (should be converted to an array |
328 | //later to speed up processing | |
329 | TGraph * gSig = new TGraph(GetRawFormatTimeBins()); | |
ee299369 | 330 | |
ee299369 | 331 | Int_t lowGain = 0; |
e5bbbc4e | 332 | Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref. |
ee299369 | 333 | |
32cd4c24 | 334 | // start loop over input stream |
335 | while (in.NextDDL()) { | |
336 | while (in.NextChannel()) { | |
7643e728 | 337 | |
338 | //Check if the signal is high or low gain and then do the fit, | |
339 | //if it is from TRU do not fit | |
340 | caloFlag = in.GetCaloFlag(); | |
341 | if (caloFlag != 0 && caloFlag != 1) continue; | |
342 | ||
5e3106bc | 343 | //Do not fit bad channels |
344 | if(pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) { | |
345 | //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow()); | |
346 | continue; | |
347 | } | |
348 | ||
32cd4c24 | 349 | // There can be zero-suppression in the raw data, |
350 | // so set up the TGraph in advance | |
351 | for (i=0; i < GetRawFormatTimeBins(); i++) { | |
7643e728 | 352 | gSig->SetPoint(i, i , -1); // init to out-of-range values |
32cd4c24 | 353 | } |
7643e728 | 354 | |
355 | Int_t maxTimeBin = 0; | |
356 | Int_t min = 0x3ff; // init to 10-bit max | |
357 | Int_t max = 0; // init to 10-bit min | |
32cd4c24 | 358 | while (in.NextBunch()) { |
7643e728 | 359 | |
32cd4c24 | 360 | const UShort_t *sig = in.GetSignals(); |
361 | startBin = in.GetStartTimeBin(); | |
7643e728 | 362 | if (maxTimeBin < startBin) { |
363 | maxTimeBin = startBin; // timebins come in reverse order | |
364 | } | |
365 | if (maxTimeBin < 0 || maxTimeBin >= GetRawFormatTimeBins()) { | |
366 | AliWarning(Form("Invalid time bin %d",maxTimeBin)); | |
367 | maxTimeBin = GetRawFormatTimeBins(); | |
32cd4c24 | 368 | } |
7643e728 | 369 | |
32cd4c24 | 370 | for (i = 0; i < in.GetBunchLength(); i++) { |
371 | time = startBin--; | |
7643e728 | 372 | gSig->SetPoint((Int_t)time, time, (Double_t) sig[i]) ; |
373 | if (max < sig[i]) max = sig[i]; | |
03a3ed3f | 374 | if (min > sig[i]) min = sig[i]; |
7643e728 | 375 | |
32cd4c24 | 376 | } |
377 | } // loop over bunches | |
56613d93 | 378 | |
7643e728 | 379 | gSig->Set(maxTimeBin+1); // set actual max size of TGraph |
380 | ||
381 | //Initialize the variables, do not keep previous values. | |
382 | // not really necessary to reset all of them (only amp and time at the moment), but better safe than sorry | |
383 | amp = -1 ; | |
384 | time = -1 ; | |
385 | ped = -1; | |
386 | ampEstimate = -1 ; | |
387 | timeEstimate = -1 ; | |
388 | pedEstimate = -1; | |
389 | if ( (max - min) > fNoiseThreshold) { | |
390 | FitRaw(gSig, signalF, maxTimeBin, amp, time, ped, | |
391 | ampEstimate, timeEstimate, pedEstimate); | |
392 | } | |
393 | ||
394 | if ( amp>0 && amp<2000 && time>0 && time<(maxTimeBin*GetRawFormatTimeBinWidth()) ) { //check both high and low end of amplitude result, and time | |
395 | //2000 is somewhat arbitrary - not nice with magic numbers in the code.. | |
396 | id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; | |
397 | lowGain = in.IsLowGain(); | |
398 | ||
399 | // check fit results: should be consistent with initial estimates | |
400 | // more magic numbers, but very loose cuts, for now.. | |
401 | // We have checked that amp an time values are positive so division for assymmetry | |
402 | // calculation should be OK/safe | |
403 | Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate); | |
404 | if ( (TMath::Abs(ampAsymm) > 0.1) || | |
405 | (TMath::Abs(time - timeEstimate) > 2*GetRawFormatTimeBinWidth()) ) { | |
5e3106bc | 406 | AliDebug(2,Form("Fit results ped %f amp %f time %f not consistent with expectations ped %f max-ped %f time %d", |
7643e728 | 407 | ped, amp, time, pedEstimate, ampEstimate, timeEstimate)); |
5e3106bc | 408 | |
7643e728 | 409 | // what should do we do then? skip this channel or assign the simple estimate? |
410 | // for now just overwrite the fit results with the simple estimate | |
411 | amp = ampEstimate; | |
412 | time = timeEstimate; | |
32cd4c24 | 413 | } |
7643e728 | 414 | |
415 | AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); | |
416 | // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp); | |
417 | // round off amplitude value to nearest integer | |
418 | AddDigit(digitsArr, id, lowGain, TMath::Nint(amp), time); | |
419 | } | |
420 | ||
32cd4c24 | 421 | // Reset graph |
422 | for (Int_t index = 0; index < gSig->GetN(); index++) { | |
7643e728 | 423 | gSig->SetPoint(index, index, -1) ; |
32cd4c24 | 424 | } |
425 | // Reset starting parameters for fit function | |
7643e728 | 426 | signalF->SetParameters(10.,5.,fTau,fOrder,0.); //reset all defaults just to be safe |
48a56166 | 427 | |
32cd4c24 | 428 | } // end while over channel |
429 | } //end while over DDL's, of input stream | |
ee299369 | 430 | |
431 | delete signalF ; | |
432 | delete gSig; | |
433 | ||
434 | return ; | |
435 | } | |
436 | ||
437 | //____________________________________________________________________________ | |
82cbdfca | 438 | void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) { |
439 | // | |
440 | // Add a new digit. | |
441 | // This routine checks whether a digit exists already for this tower | |
442 | // and then decides whether to use the high or low gain info | |
443 | // | |
444 | // Called by Raw2Digits | |
445 | ||
446 | AliEMCALDigit *digit = 0, *tmpdigit = 0; | |
82cbdfca | 447 | TIter nextdigit(digitsArr); |
448 | while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { | |
449 | if (tmpdigit->GetId() == id) | |
450 | digit = tmpdigit; | |
451 | } | |
452 | ||
453 | if (!digit) { // no digit existed for this tower; create one | |
a7ec7165 | 454 | if (lowGain && amp > fgkOverflowCut) |
82cbdfca | 455 | amp = Int_t(fHighLowGainFactor * amp); |
456 | Int_t idigit = digitsArr->GetEntries(); | |
457 | new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ; | |
458 | } | |
459 | else { // a digit already exists, check range | |
b4133f05 | 460 | // (use high gain if signal < cut value, otherwise low gain) |
82cbdfca | 461 | if (lowGain) { // new digit is low gain |
b4133f05 | 462 | if (digit->GetAmp() > fgkOverflowCut) { // use if stored digit is out of range |
82cbdfca | 463 | digit->SetAmp(Int_t(fHighLowGainFactor * amp)); |
464 | digit->SetTime(time); | |
465 | } | |
466 | } | |
b4133f05 | 467 | else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range |
82cbdfca | 468 | digit->SetAmp(amp); |
469 | digit->SetTime(time); | |
470 | } | |
471 | } | |
472 | } | |
473 | ||
474 | //____________________________________________________________________________ | |
7643e728 | 475 | void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & ped, Float_t & ampEstimate, Float_t & timeEstimate, Float_t & pedEstimate, const Float_t cut) const |
ee299369 | 476 | { |
477 | // Fits the raw signal time distribution; from AliEMCALGetter | |
7643e728 | 478 | // last argument: Float_t cut = 0.0; // indicating how much of amplitude w.r.t. max value fit should be above noise and pedestal |
479 | ||
480 | // initialize return values | |
481 | amp = 0; | |
482 | time = 0; | |
483 | ped = 0; | |
484 | ampEstimate = 0; | |
485 | timeEstimate = 0; | |
486 | pedEstimate = 0; | |
487 | ||
488 | // 0th step: remove plateau / overflow candidates | |
489 | // before trying to estimate amplitude, search for maxima etc. | |
490 | // | |
491 | Int_t nOrig = gSig->GetN(); // number of samples before we remove any overflows | |
492 | // Values for readback from input graph | |
493 | Double_t ttime = 0; | |
494 | Double_t signal = 0; | |
495 | ||
496 | /* | |
497 | // start: tmp dump of all values | |
498 | for (Int_t i=0; i<gSig->GetN(); i++) { | |
499 | gSig->GetPoint(i, ttime, signal) ; // get values | |
500 | printf("orig: i %d, time %f, signal %f\n",i, ttime, signal); | |
501 | } | |
502 | // end: tmp dump of all values | |
503 | */ | |
504 | ||
505 | // start from back of TGraph since RemovePoint will downshift indices | |
506 | for (Int_t i=nOrig-1; i>=0; i--) { | |
507 | gSig->GetPoint(i, ttime, signal) ; // get values | |
508 | if (signal >= (pedEstimate + fgkOverflowCut) ) { | |
509 | gSig->RemovePoint(i); | |
510 | } | |
511 | } | |
ee299369 | 512 | |
7643e728 | 513 | // 1st step: we try to estimate the pedestal value |
514 | Int_t nPed = 0; | |
515 | for (Int_t index = 0; index < gSig->GetN(); index++) { | |
82cbdfca | 516 | gSig->GetPoint(index, ttime, signal) ; |
7643e728 | 517 | // ttime < fNPedsamples used for pedestal estimate; |
518 | // ttime >= fNPedSamples used for signal checks | |
519 | if (signal >= 0 && ttime<fNPedSamples) { // valid value | |
520 | pedEstimate += signal; | |
82cbdfca | 521 | nPed++; |
ee299369 | 522 | } |
523 | } | |
e9dbb64a | 524 | |
82cbdfca | 525 | if (nPed > 0) |
7643e728 | 526 | pedEstimate /= nPed; |
e9dbb64a | 527 | else { |
7643e728 | 528 | //AliWarning("Could not determine pedestal"); |
529 | AliDebug(1,"Could not determine pedestal"); | |
530 | pedEstimate = 0; // good estimate for ZeroSupp data (non ZS data should have no problem with pedestal estimate) | |
e9dbb64a | 531 | } |
82cbdfca | 532 | |
7643e728 | 533 | // 2nd step: we look through the rest of the time-bins/ADC values and |
534 | // see if we have something that looks like a signal. | |
535 | // We look for a first local maxima, as well as for a global maxima | |
536 | Int_t locMaxFound = 0; | |
537 | Int_t locMaxId = 0; // time-bin index at first local max | |
538 | Float_t locMaxSig = -1; // actual local max value | |
539 | Int_t globMaxId = 0; // time-bin index at global max | |
540 | Float_t globMaxSig = -1; // actual global max value | |
541 | // We will also look for any values that look like they are in overflow region | |
542 | for (Int_t i=0; i<gSig->GetN(); i++) { | |
543 | gSig->GetPoint(i, ttime, signal) ; // get values | |
544 | ||
545 | // ttime < fNPedsamples used for pedestal estimate; | |
546 | // ttime >= fNPedSamples used for signal checks | |
547 | if (ttime >= fNPedSamples) { | |
548 | ||
549 | // look for first local maximum signal=ADC value | |
550 | if (!locMaxFound && signal > locMaxSig) { | |
551 | locMaxId = i; | |
552 | locMaxSig = signal; | |
e9dbb64a | 553 | } |
7643e728 | 554 | else if ( locMaxSig > (pedEstimate + fNoiseThreshold) ) { |
555 | // we enter this condition after signal<=max, but previous | |
556 | // max value was large enough. I.e. at least a significant local | |
557 | // maxima has been found (just before) | |
558 | locMaxFound = 1; | |
e9dbb64a | 559 | } |
7643e728 | 560 | |
561 | // also check for global maximum.. | |
562 | if (signal > globMaxSig) { | |
563 | globMaxId = i; | |
564 | globMaxSig = signal; | |
c61fe3b4 | 565 | } |
7643e728 | 566 | } // ttime check |
567 | } // end for-loop over samples after pedestal | |
568 | ||
569 | // OK, we have looked through the signal spectra, let's see if we should try to make the fit | |
570 | ampEstimate = locMaxSig - pedEstimate; // estimate using first local maxima | |
571 | if ( ampEstimate > fNoiseThreshold ) { // else it's just noise | |
572 | ||
573 | //Check that the local maximum we will use is not at the end or beginning of time sample range | |
574 | Double_t timeMax = -1; | |
575 | Int_t iMax = locMaxId; | |
576 | gSig->GetPoint(locMaxId, timeMax, signal) ; | |
577 | if (timeMax < 2 || timeMax > lastTimeBin-1) { // lastTimeBin is the lowest kept time-sample; current (Dec 2009) case | |
578 | // if (timeMax < 2 || timeMax > lastTimeBin-2) { // for when lastTimeBin is the lowest read-out time-sample, future (2010) case | |
579 | AliDebug(1,Form("Skip fit, maximum of the sample close to the edges : timeMax %3.2f, ampEstimate %3.2f",timeMax, ampEstimate)); | |
580 | return; | |
581 | } | |
582 | ||
583 | // Check if the local and global maximum disagree | |
584 | if (locMaxId != globMaxId) { | |
585 | AliDebug(1,Form("Warning, local first maximum %d does not agree with global maximum %d\n", locMaxId, globMaxId)); | |
586 | return; | |
587 | } | |
588 | ||
589 | // Get the maximum and find the lowest timebin (tailmin) where the ADC value is not | |
590 | // significantly different from the pedestal | |
591 | // first lower times edge a.k.a. tailmin | |
592 | Int_t tailMin = 0; | |
593 | Double_t tmptime = 0; | |
594 | for (Int_t i=iMax-1; i > 0; i--) { | |
595 | gSig->GetPoint(i, tmptime, signal) ; | |
596 | if((signal-pedEstimate) < fNoiseThreshold){ | |
597 | tailMin = i; | |
598 | break; | |
e9dbb64a | 599 | } |
5a056daa | 600 | } |
7643e728 | 601 | // then same exercise for the higher times edge a.k.a. tailmax |
602 | Int_t tailMax = lastTimeBin; | |
603 | for (Int_t i=iMax+1; i < gSig->GetN(); i++) { | |
604 | gSig->GetPoint(i, tmptime, signal) ; | |
605 | if ((signal-pedEstimate) <= (ampEstimate*cut + fNoiseThreshold)) { // stop fit at cut-fraction of amplitude above noise-threshold (cut>0 would mean avoid the pulse shape falling edge) | |
606 | tailMax = i; | |
607 | break; | |
608 | } | |
fb070798 | 609 | } |
fb070798 | 610 | |
7643e728 | 611 | // remove all points which are not in the distribution around maximum |
612 | // i.e. up to tailmin, and from tailmax | |
613 | if ( tailMax != (gSig->GetN()-1) ){ // else nothing to remove | |
614 | nOrig = gSig->GetN(); // can't use GetN call in for loop below since gSig size changes.. | |
615 | for(int j = tailMax; j < nOrig; j++) gSig->RemovePoint(tailMax); | |
fb070798 | 616 | } |
7643e728 | 617 | for(int j = 0; j<=tailMin; j++) gSig->RemovePoint(0); |
5a056daa | 618 | |
7643e728 | 619 | if(gSig->GetN() < 3) { |
620 | AliDebug(2,Form("Skip fit, number of entries in sample smaller than number of fitting parameters: in sample %d, fitting param 3", | |
621 | gSig->GetN() )); | |
622 | return; | |
623 | } | |
8cb998bd | 624 | |
7643e728 | 625 | timeEstimate = timeMax * GetRawFormatTimeBinWidth(); |
626 | ||
627 | // determine what the valid fit range is | |
628 | Double_t minFit = 9999; | |
629 | Double_t maxFit = 0; | |
630 | for (Int_t i=0; i < gSig->GetN(); i++) { | |
631 | gSig->GetPoint(i, ttime, signal); | |
632 | if (minFit > ttime) minFit=ttime; | |
633 | if (maxFit < ttime) maxFit=ttime; | |
634 | //debug: printf("no tail: i %d, time %f, signal %f\n",i, ttime, signal); | |
635 | } | |
636 | signalF->SetRange(minFit, maxFit); | |
637 | ||
638 | signalF->FixParameter(4, pedEstimate) ; | |
639 | signalF->SetParameter(1, timeMax); | |
640 | signalF->SetParameter(0, ampEstimate); | |
8cb998bd | 641 | |
642 | gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points | |
7643e728 | 643 | |
644 | // assign fit results | |
8cb998bd | 645 | amp = signalF->GetParameter(0); |
7643e728 | 646 | time = signalF->GetParameter(1) * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger? |
647 | ped = signalF->GetParameter(4); | |
648 | ||
5e3106bc | 649 | //BEG YS alternative methods to calculate the amplitude |
650 | Double_t * ymx = gSig->GetX() ; | |
651 | Double_t * ymy = gSig->GetY() ; | |
652 | const Int_t kN = 3 ; | |
653 | Double_t ymMaxX[kN] = {0., 0., 0.} ; | |
654 | Double_t ymMaxY[kN] = {0., 0., 0.} ; | |
655 | Double_t ymax = 0. ; | |
656 | // find the maximum amplitude | |
657 | Int_t ymiMax = 0 ; | |
658 | for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) { | |
659 | if (ymy[ymi] > ymMaxY[0] ) { | |
660 | ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude | |
661 | ymMaxX[0] = ymx[ymi] ; | |
662 | ymiMax = ymi ; | |
663 | } | |
664 | } | |
665 | // find the maximum by fitting a parabola through the max and the two adjacent samples | |
666 | if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) { | |
667 | ymMaxY[1] = ymy[ymiMax+1] ; | |
668 | ymMaxY[2] = ymy[ymiMax-1] ; | |
669 | ymMaxX[1] = ymx[ymiMax+1] ; | |
670 | ymMaxX[2] = ymx[ymiMax-1] ; | |
671 | if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) { | |
672 | //fit a parabola through the 3 points y= a+bx+x*x*x | |
673 | Double_t sy = 0 ; | |
674 | Double_t sx = 0 ; | |
675 | Double_t sx2 = 0 ; | |
676 | Double_t sx3 = 0 ; | |
677 | Double_t sx4 = 0 ; | |
678 | Double_t sxy = 0 ; | |
679 | Double_t sx2y = 0 ; | |
680 | for (Int_t i = 0; i < kN ; i++) { | |
681 | sy += ymMaxY[i] ; | |
682 | sx += ymMaxX[i] ; | |
683 | sx2 += ymMaxX[i]*ymMaxX[i] ; | |
684 | sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; | |
685 | sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; | |
686 | sxy += ymMaxX[i]*ymMaxY[i] ; | |
687 | sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ; | |
688 | } | |
689 | Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2); | |
690 | Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ; | |
691 | Double_t c = cN / cD ; | |
692 | Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ; | |
693 | Double_t a = (sy-b*sx-c*sx2)/kN ; | |
694 | Double_t xmax = -b/(2*c) ; | |
695 | ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude | |
696 | } | |
697 | } | |
698 | ||
699 | Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ; | |
700 | if (diff > 0.1) | |
701 | amp = ymMaxY[0] ; | |
702 | ||
703 | //END YS | |
704 | ||
7643e728 | 705 | } // ampEstimate > fNoiseThreshold |
ee299369 | 706 | return; |
707 | } | |
708 | //__________________________________________________________________ | |
709 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) | |
710 | { | |
8cb998bd | 711 | // Matches version used in 2007 beam test |
712 | // | |
ee299369 | 713 | // Shape of the electronics raw reponse: |
714 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
715 | // | |
7643e728 | 716 | // xx = (t - t0 + tau) / tau [xx is just a convenient help variable] |
717 | // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0 | |
718 | // F = 0 for xx < 0 | |
ee299369 | 719 | // |
720 | // parameters: | |
8cb998bd | 721 | // A: par[0] // Amplitude = peak value |
722 | // t0: par[1] | |
723 | // tau: par[2] | |
724 | // N: par[3] | |
725 | // ped: par[4] | |
ee299369 | 726 | // |
727 | Double_t signal ; | |
8cb998bd | 728 | Double_t tau =par[2]; |
e5bbbc4e | 729 | Double_t n =par[3]; |
8cb998bd | 730 | Double_t ped = par[4]; |
731 | Double_t xx = ( x[0] - par[1] + tau ) / tau ; | |
ee299369 | 732 | |
5a056daa | 733 | if (xx <= 0) |
8cb998bd | 734 | signal = ped ; |
ee299369 | 735 | else { |
e5bbbc4e | 736 | signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; |
ee299369 | 737 | } |
738 | return signal ; | |
739 | } | |
740 | ||
741 | //__________________________________________________________________ | |
742 | Bool_t AliEMCALRawUtils::RawSampledResponse( | |
743 | const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const | |
744 | { | |
745 | // for a start time dtime and an amplitude damp given by digit, | |
746 | // calculates the raw sampled response AliEMCAL::RawResponseFunction | |
747 | ||
ee299369 | 748 | Bool_t lowGain = kFALSE ; |
749 | ||
48a56166 | 750 | // A: par[0] // Amplitude = peak value |
751 | // t0: par[1] | |
752 | // tau: par[2] | |
753 | // N: par[3] | |
754 | // ped: par[4] | |
755 | ||
56e13066 | 756 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); |
48a56166 | 757 | signalF.SetParameter(0, damp) ; |
56e13066 | 758 | signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ; |
b4133f05 | 759 | signalF.SetParameter(2, fTau) ; |
760 | signalF.SetParameter(3, fOrder); | |
fe93d365 | 761 | signalF.SetParameter(4, fgPedestalValue); |
ee299369 | 762 | |
763 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { | |
56e13066 | 764 | Double_t signal = signalF.Eval(iTime) ; |
fe93d365 | 765 | |
7643e728 | 766 | // Next lines commeted for the moment but in principle it is not necessary to add |
ff10f540 | 767 | // extra noise since noise already added at the digits level. |
7643e728 | 768 | |
fe93d365 | 769 | //According to Terry Awes, 13-Apr-2008 |
770 | //add gaussian noise in quadrature to each sample | |
09974781 | 771 | //Double_t noise = gRandom->Gaus(0.,fgFEENoise); |
fe93d365 | 772 | //signal = sqrt(signal*signal + noise*noise); |
773 | ||
e2c2134b | 774 | // March 17,09 for fast fit simulations by Alexei Pavlinov. |
775 | // Get from PHOS analysis. In some sense it is open questions. | |
ff10f540 | 776 | //Double_t noise = gRandom->Gaus(0.,fgFEENoise); |
777 | //signal += noise; | |
7643e728 | 778 | |
ee299369 | 779 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; |
b4133f05 | 780 | if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits |
781 | adcH[iTime] = fgkRawSignalOverflow ; | |
ee299369 | 782 | lowGain = kTRUE ; |
783 | } | |
784 | ||
785 | signal /= fHighLowGainFactor; | |
786 | ||
787 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
b4133f05 | 788 | if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits |
789 | adcL[iTime] = fgkRawSignalOverflow ; | |
ee299369 | 790 | } |
791 | return lowGain ; | |
792 | } |