]>
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" |
9f467289 | 51 | #include "AliCaloFastAltroFitv0.h" |
c8603a2b | 52 | #include "AliCaloNeuralFit.h" |
16605c06 | 53 | #include "AliCaloBunchInfo.h" |
54 | #include "AliCaloFitResults.h" | |
55 | #include "AliCaloRawAnalyzerLMS.h" | |
56 | #include "AliCaloRawAnalyzerPeakFinder.h" | |
57 | #include "AliCaloRawAnalyzerCrude.h" | |
9f467289 | 58 | |
ee299369 | 59 | ClassImp(AliEMCALRawUtils) |
21cad85c | 60 | |
ee299369 | 61 | // Signal shape parameters |
89d338a6 | 62 | Int_t AliEMCALRawUtils::fgTimeBins = 256; // number of sampling bins of the raw RO signal (we typically use 15-50; theoretical max is 1k+) |
e5bbbc4e | 63 | Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns |
09974781 | 64 | Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec |
ee299369 | 65 | |
66 | // some digitization constants | |
67 | Int_t AliEMCALRawUtils::fgThreshold = 1; | |
68 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule | |
e5bbbc4e | 69 | Int_t AliEMCALRawUtils::fgPedestalValue = 32; // pedestal value for digits2raw |
70 | Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled) | |
ee299369 | 71 | |
16605c06 | 72 | AliEMCALRawUtils::AliEMCALRawUtils(fitAlgorithm fitAlgo) |
b4133f05 | 73 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), |
9f467289 | 74 | fNPedSamples(0), fGeom(0), fOption(""), |
16605c06 | 75 | fRemoveBadChannels(kTRUE),fFittingAlgorithm(0),fRawAnalyzer(0) |
8cb998bd | 76 | { |
b4133f05 | 77 | |
78 | //These are default parameters. | |
79 | //Can be re-set from without with setter functions | |
9f467289 | 80 | //Already set in the OCDB and passed via setter in the AliEMCALReconstructor |
ee299369 | 81 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) |
b4133f05 | 82 | fOrder = 2; // order of gamma fn |
83 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
7643e728 | 84 | fNoiseThreshold = 3; // 3 ADC counts is approx. noise level |
85 | fNPedSamples = 4; // less than this value => likely pedestal samples | |
9f467289 | 86 | fRemoveBadChannels = kTRUE; //Remove bad channels before fitting |
16605c06 | 87 | fFittingAlgorithm = fitAlgo; |
88 | if (fitAlgo == kLMS) { | |
89 | fRawAnalyzer = new AliCaloRawAnalyzerLMS(); | |
90 | } | |
91 | else if (fitAlgo == kPeakFinder) { | |
92 | fRawAnalyzer = new AliCaloRawAnalyzerPeakFinder(); | |
93 | } | |
94 | else if (fitAlgo == kCrude) { | |
95 | fRawAnalyzer = new AliCaloRawAnalyzerCrude(); | |
96 | } | |
97 | else { | |
98 | fRawAnalyzer = new AliCaloRawAnalyzer(); | |
99 | } | |
100 | ||
65bdc82f | 101 | //Get Mapping RCU files from the AliEMCALRecParam |
102 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
103 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
104 | ||
21cad85c | 105 | for(Int_t i = 0; i < 4; i++) { |
65bdc82f | 106 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
107 | } | |
108 | ||
72c58de0 | 109 | //To make sure we match with the geometry in a simulation file, |
110 | //let's try to get it first. If not, take the default geometry | |
33c3c91a | 111 | AliRunLoader *rl = AliRunLoader::Instance(); |
c61fe3b4 | 112 | if(!rl) AliError("Cannot find RunLoader!"); |
72c58de0 | 113 | if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) { |
114 | fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); | |
115 | } else { | |
116 | AliInfo(Form("Using default geometry in raw reco")); | |
937d0661 | 117 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); |
65bdc82f | 118 | } |
119 | ||
72c58de0 | 120 | if(!fGeom) AliFatal(Form("Could not get geometry!")); |
121 | ||
65bdc82f | 122 | } |
123 | ||
124 | //____________________________________________________________________________ | |
16605c06 | 125 | AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry, fitAlgorithm fitAlgo) |
5544799a | 126 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), |
9f467289 | 127 | fNPedSamples(0), fGeom(pGeometry), fOption(""), |
16605c06 | 128 | fRemoveBadChannels(kTRUE),fFittingAlgorithm(0),fRawAnalyzer() |
5544799a | 129 | { |
130 | // | |
131 | // Initialize with the given geometry - constructor required by HLT | |
132 | // HLT does not use/support AliRunLoader(s) instances | |
133 | // This is a minimum intervention solution | |
134 | // Comment by MPloskon@lbl.gov | |
135 | // | |
136 | ||
137 | //These are default parameters. | |
138 | //Can be re-set from without with setter functions | |
9f467289 | 139 | //Already set in the OCDB and passed via setter in the AliEMCALReconstructor |
5544799a | 140 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) |
141 | fOrder = 2; // order of gamma fn | |
142 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
7643e728 | 143 | fNoiseThreshold = 3; // 3 ADC counts is approx. noise level |
144 | fNPedSamples = 4; // less than this value => likely pedestal samples | |
9f467289 | 145 | fRemoveBadChannels = kTRUE; //Remove bad channels before fitting |
16605c06 | 146 | fFittingAlgorithm = fitAlgo; |
147 | if (fitAlgo == kLMS) { | |
148 | fRawAnalyzer = new AliCaloRawAnalyzerLMS(); | |
149 | } | |
150 | else if (fitAlgo == kPeakFinder) { | |
151 | fRawAnalyzer = new AliCaloRawAnalyzerPeakFinder(); | |
152 | } | |
153 | else if (fitAlgo == kCrude) { | |
154 | fRawAnalyzer = new AliCaloRawAnalyzerCrude(); | |
155 | } | |
156 | else { | |
157 | fRawAnalyzer = new AliCaloRawAnalyzer(); | |
158 | } | |
159 | ||
5544799a | 160 | //Get Mapping RCU files from the AliEMCALRecParam |
161 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
162 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
163 | ||
21cad85c | 164 | for(Int_t i = 0; i < 4; i++) { |
5544799a | 165 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
166 | } | |
167 | ||
168 | if(!fGeom) AliFatal(Form("Could not get geometry!")); | |
169 | ||
170 | } | |
171 | ||
172 | //____________________________________________________________________________ | |
65bdc82f | 173 | AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) |
174 | : TObject(), | |
175 | fHighLowGainFactor(rawU.fHighLowGainFactor), | |
b4133f05 | 176 | fOrder(rawU.fOrder), |
177 | fTau(rawU.fTau), | |
178 | fNoiseThreshold(rawU.fNoiseThreshold), | |
179 | fNPedSamples(rawU.fNPedSamples), | |
65bdc82f | 180 | fGeom(rawU.fGeom), |
9f467289 | 181 | fOption(rawU.fOption), |
182 | fRemoveBadChannels(rawU.fRemoveBadChannels), | |
16605c06 | 183 | fFittingAlgorithm(rawU.fFittingAlgorithm), |
184 | fRawAnalyzer(rawU.fRawAnalyzer) | |
65bdc82f | 185 | { |
186 | //copy ctor | |
187 | fMapping[0] = rawU.fMapping[0]; | |
188 | fMapping[1] = rawU.fMapping[1]; | |
21cad85c | 189 | fMapping[2] = rawU.fMapping[2]; |
190 | fMapping[3] = rawU.fMapping[3]; | |
65bdc82f | 191 | } |
192 | ||
193 | //____________________________________________________________________________ | |
194 | AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) | |
195 | { | |
196 | //assignment operator | |
197 | ||
198 | if(this != &rawU) { | |
199 | fHighLowGainFactor = rawU.fHighLowGainFactor; | |
b4133f05 | 200 | fOrder = rawU.fOrder; |
201 | fTau = rawU.fTau; | |
202 | fNoiseThreshold = rawU.fNoiseThreshold; | |
203 | fNPedSamples = rawU.fNPedSamples; | |
65bdc82f | 204 | fGeom = rawU.fGeom; |
205 | fOption = rawU.fOption; | |
9f467289 | 206 | fRemoveBadChannels = rawU.fRemoveBadChannels; |
207 | fFittingAlgorithm = rawU.fFittingAlgorithm; | |
16605c06 | 208 | fRawAnalyzer = rawU.fRawAnalyzer; |
65bdc82f | 209 | fMapping[0] = rawU.fMapping[0]; |
210 | fMapping[1] = rawU.fMapping[1]; | |
21cad85c | 211 | fMapping[2] = rawU.fMapping[2]; |
212 | fMapping[3] = rawU.fMapping[3]; | |
65bdc82f | 213 | } |
214 | ||
215 | return *this; | |
216 | ||
ee299369 | 217 | } |
65bdc82f | 218 | |
ee299369 | 219 | //____________________________________________________________________________ |
220 | AliEMCALRawUtils::~AliEMCALRawUtils() { | |
e5bbbc4e | 221 | //dtor |
65bdc82f | 222 | |
ee299369 | 223 | } |
65bdc82f | 224 | |
ee299369 | 225 | //____________________________________________________________________________ |
65bdc82f | 226 | void AliEMCALRawUtils::Digits2Raw() |
ee299369 | 227 | { |
228 | // convert digits of the current event to raw data | |
229 | ||
33c3c91a | 230 | AliRunLoader *rl = AliRunLoader::Instance(); |
ee299369 | 231 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); |
232 | ||
233 | // get the digits | |
234 | loader->LoadDigits("EMCAL"); | |
235 | loader->GetEvent(); | |
236 | TClonesArray* digits = loader->Digits() ; | |
237 | ||
238 | if (!digits) { | |
239 | Warning("Digits2Raw", "no digits found !"); | |
240 | return; | |
241 | } | |
65bdc82f | 242 | |
ee299369 | 243 | static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here |
244 | AliAltroBuffer* buffers[nDDL]; | |
245 | for (Int_t i=0; i < nDDL; i++) | |
246 | buffers[i] = 0; | |
247 | ||
e2c2134b | 248 | TArrayI adcValuesLow(fgTimeBins); |
249 | TArrayI adcValuesHigh(fgTimeBins); | |
ee299369 | 250 | |
ee299369 | 251 | // loop over digits (assume ordered digits) |
252 | for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { | |
253 | AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ; | |
254 | if (digit->GetAmp() < fgThreshold) | |
255 | continue; | |
256 | ||
257 | //get cell indices | |
258 | Int_t nSM = 0; | |
259 | Int_t nIphi = 0; | |
260 | Int_t nIeta = 0; | |
261 | Int_t iphi = 0; | |
262 | Int_t ieta = 0; | |
263 | Int_t nModule = 0; | |
65bdc82f | 264 | fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); |
265 | fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; | |
ee299369 | 266 | |
21cad85c | 267 | //Check which is the RCU, 0 or 1, of the cell. |
ee299369 | 268 | Int_t iRCU = -111; |
269 | //RCU0 | |
270 | if (0<=iphi&&iphi<8) iRCU=0; // first cable row | |
271 | else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; | |
272 | //second cable row | |
273 | //RCU1 | |
274 | else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; | |
275 | //second cable row | |
276 | else if(16<=iphi&&iphi<24) iRCU=1; // third cable row | |
21cad85c | 277 | |
278 | if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same | |
279 | ||
e36e3bcf | 280 | if (iRCU<0) |
281 | Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); | |
ee299369 | 282 | |
283 | //Which DDL? | |
284 | Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; | |
285 | if (iDDL >= nDDL) | |
286 | Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); | |
287 | ||
288 | if (buffers[iDDL] == 0) { | |
289 | // open new file and write dummy header | |
290 | TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); | |
21cad85c | 291 | //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C |
292 | Int_t iRCUside=iRCU+(nSM%2)*2; | |
293 | //iRCU=0 and even (0) SM -> RCU0A.data 0 | |
294 | //iRCU=1 and even (0) SM -> RCU1A.data 1 | |
295 | //iRCU=0 and odd (1) SM -> RCU0C.data 2 | |
296 | //iRCU=1 and odd (1) SM -> RCU1C.data 3 | |
297 | //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl; | |
298 | buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]); | |
ee299369 | 299 | buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; |
300 | } | |
301 | ||
302 | // out of time range signal (?) | |
303 | if (digit->GetTimeR() > GetRawFormatTimeMax() ) { | |
304 | AliInfo("Signal is out of time range.\n"); | |
305 | buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp()); | |
306 | buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin | |
307 | buffers[iDDL]->FillBuffer(3); // bunch length | |
308 | buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer | |
309 | // calculate the time response function | |
310 | } else { | |
e2c2134b | 311 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; |
ee299369 | 312 | if (lowgain) |
e2c2134b | 313 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold); |
ee299369 | 314 | else |
e2c2134b | 315 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold); |
ee299369 | 316 | } |
317 | } | |
318 | ||
319 | // write headers and close files | |
320 | for (Int_t i=0; i < nDDL; i++) { | |
321 | if (buffers[i]) { | |
322 | buffers[i]->Flush(); | |
323 | buffers[i]->WriteDataHeader(kFALSE, kFALSE); | |
324 | delete buffers[i]; | |
325 | } | |
326 | } | |
65bdc82f | 327 | |
ee299369 | 328 | loader->UnloadDigits(); |
329 | } | |
330 | ||
331 | //____________________________________________________________________________ | |
16605c06 | 332 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap) |
ee299369 | 333 | { |
65bdc82f | 334 | // convert raw data of the current event to digits |
ee299369 | 335 | |
c47157cd | 336 | digitsArr->Clear(); |
ee299369 | 337 | |
c47157cd | 338 | if (!digitsArr) { |
ee299369 | 339 | Error("Raw2Digits", "no digits found !"); |
340 | return; | |
341 | } | |
342 | if (!reader) { | |
343 | Error("Raw2Digits", "no raw reader found !"); | |
344 | return; | |
345 | } | |
346 | ||
32cd4c24 | 347 | AliCaloRawStreamV3 in(reader,"EMCAL",fMapping); |
ee299369 | 348 | // Select EMCAL DDL's; |
7643e728 | 349 | reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL |
feedcab9 | 350 | |
16605c06 | 351 | // fRawAnalyzer setup |
352 | fRawAnalyzer->SetAmpCut(fNoiseThreshold); | |
353 | fRawAnalyzer->SetFitArrayCut(fNoiseThreshold); | |
354 | fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later | |
ee299369 | 355 | |
16605c06 | 356 | // channel info parameters |
ee299369 | 357 | Int_t lowGain = 0; |
e5bbbc4e | 358 | Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref. |
ee299369 | 359 | |
32cd4c24 | 360 | // start loop over input stream |
361 | while (in.NextDDL()) { | |
362 | while (in.NextChannel()) { | |
7643e728 | 363 | |
364 | //Check if the signal is high or low gain and then do the fit, | |
16605c06 | 365 | //if it is from TRU or LEDMon do not fit |
7643e728 | 366 | caloFlag = in.GetCaloFlag(); |
367 | if (caloFlag != 0 && caloFlag != 1) continue; | |
368 | ||
5e3106bc | 369 | //Do not fit bad channels |
9f467289 | 370 | if(fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) { |
5e3106bc | 371 | //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow()); |
372 | continue; | |
373 | } | |
374 | ||
16605c06 | 375 | vector<AliCaloBunchInfo> bunchlist; |
32cd4c24 | 376 | while (in.NextBunch()) { |
16605c06 | 377 | bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) ); |
378 | } // loop over bunches | |
7643e728 | 379 | |
16605c06 | 380 | Float_t time = 0; |
381 | Float_t amp = 0; | |
382 | ||
383 | if (fFittingAlgorithm == kLMS || fFittingAlgorithm == kPeakFinder || fFittingAlgorithm == kCrude) { | |
384 | // all functionality to determine amp and time etc is encapsulated inside the Evaluate call for these methods | |
385 | AliCaloFitResults fitResults = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2()); | |
386 | ||
387 | amp = fitResults.GetAmp(); | |
388 | time = fitResults.GetTof(); | |
389 | } | |
390 | else { // for the other methods we for now use the functionality of | |
391 | // AliCaloRawAnalyzer as well, to select samples and prepare for fits, | |
392 | // if it looks like there is something to fit | |
393 | ||
394 | // parameters init. | |
395 | Float_t ampEstimate = 0; | |
396 | short maxADC = 0; | |
397 | short timeEstimate = 0; | |
398 | Float_t pedEstimate = 0; | |
399 | Int_t first = 0; | |
400 | Int_t last = 0; | |
401 | Int_t bunchIndex = 0; | |
402 | // | |
403 | // The PreFitEvaluateSamples + later call to FitRaw will hopefully | |
404 | // be replaced by a single Evaluate call or so soon, like for the other | |
405 | // methods, but this should be good enough for evaluation of | |
406 | // the methods for now (Jan. 2010) | |
407 | // | |
408 | int nsamples = fRawAnalyzer->PreFitEvaluateSamples( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2(), bunchIndex, ampEstimate, maxADC, timeEstimate, pedEstimate, first, last); | |
7643e728 | 409 | |
16605c06 | 410 | if (ampEstimate > fNoiseThreshold) { // something worth looking at |
7643e728 | 411 | |
16605c06 | 412 | time = timeEstimate; |
413 | amp = ampEstimate; | |
414 | ||
415 | if ( nsamples > 1 ) { // possibly something to fit | |
416 | FitRaw(first, last, amp, time); | |
9f467289 | 417 | } |
16605c06 | 418 | |
419 | if ( amp>0 && time>0 ) { // brief sanity check of fit results | |
420 | ||
421 | // check fit results: should be consistent with initial estimates | |
422 | // more magic numbers, but very loose cuts, for now.. | |
423 | // We have checked that amp and ampEstimate values are positive so division for assymmetry | |
424 | // calculation should be OK/safe | |
425 | Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate); | |
426 | if ( (TMath::Abs(ampAsymm) > 0.1) ) { | |
427 | AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations ped %f max-ped %f time %d", | |
428 | amp, time, pedEstimate, ampEstimate, timeEstimate)); | |
429 | ||
430 | // what should do we do then? skip this channel or assign the simple estimate? | |
431 | // for now just overwrite the fit results with the simple estimate | |
432 | amp = ampEstimate; | |
433 | time = timeEstimate; | |
434 | } // asymm check | |
435 | } // amp & time check | |
436 | } // ampEstimate check | |
437 | } // method selection | |
438 | ||
439 | if (amp > fNoiseThreshold) { // something to be stored | |
440 | Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; | |
7643e728 | 441 | lowGain = in.IsLowGain(); |
442 | ||
16605c06 | 443 | // go from time-bin units to physical time fgtimetrigger |
444 | time = time * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger? | |
7643e728 | 445 | |
446 | AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); | |
447 | // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp); | |
448 | // round off amplitude value to nearest integer | |
449 | AddDigit(digitsArr, id, lowGain, TMath::Nint(amp), time); | |
450 | } | |
451 | ||
32cd4c24 | 452 | } // end while over channel |
453 | } //end while over DDL's, of input stream | |
16605c06 | 454 | |
ee299369 | 455 | return ; |
456 | } | |
457 | ||
458 | //____________________________________________________________________________ | |
82cbdfca | 459 | void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) { |
460 | // | |
461 | // Add a new digit. | |
462 | // This routine checks whether a digit exists already for this tower | |
463 | // and then decides whether to use the high or low gain info | |
464 | // | |
465 | // Called by Raw2Digits | |
466 | ||
467 | AliEMCALDigit *digit = 0, *tmpdigit = 0; | |
82cbdfca | 468 | TIter nextdigit(digitsArr); |
469 | while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { | |
470 | if (tmpdigit->GetId() == id) | |
471 | digit = tmpdigit; | |
472 | } | |
473 | ||
474 | if (!digit) { // no digit existed for this tower; create one | |
a7ec7165 | 475 | if (lowGain && amp > fgkOverflowCut) |
82cbdfca | 476 | amp = Int_t(fHighLowGainFactor * amp); |
477 | Int_t idigit = digitsArr->GetEntries(); | |
478 | new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ; | |
479 | } | |
480 | else { // a digit already exists, check range | |
b4133f05 | 481 | // (use high gain if signal < cut value, otherwise low gain) |
82cbdfca | 482 | if (lowGain) { // new digit is low gain |
b4133f05 | 483 | if (digit->GetAmp() > fgkOverflowCut) { // use if stored digit is out of range |
82cbdfca | 484 | digit->SetAmp(Int_t(fHighLowGainFactor * amp)); |
485 | digit->SetTime(time); | |
486 | } | |
487 | } | |
b4133f05 | 488 | else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range |
82cbdfca | 489 | digit->SetAmp(amp); |
490 | digit->SetTime(time); | |
491 | } | |
492 | } | |
493 | } | |
494 | ||
495 | //____________________________________________________________________________ | |
16605c06 | 496 | void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time) const |
497 | { // Fits the raw signal time distribution | |
498 | ||
499 | //-------------------------------------------------- | |
500 | //Do the fit, different fitting algorithms available | |
501 | //-------------------------------------------------- | |
502 | int nsamples = lastTimeBin - firstTimeBin + 1; | |
503 | TGraph *gSig = new TGraph( nsamples); | |
504 | for (int i=0; i<nsamples; i++) { | |
505 | Int_t timebin = firstTimeBin + i; | |
506 | gSig->SetPoint(timebin, timebin, fRawAnalyzer->GetReversed(timebin)); | |
7643e728 | 507 | } |
ee299369 | 508 | |
16605c06 | 509 | switch(fFittingAlgorithm) { |
510 | case kStandard: | |
511 | { | |
512 | if (gSig->GetN() < 3) { return; } // nothing much to fit | |
513 | //printf("Standard fitter \n"); | |
514 | // Create Graph to hold data we will fit | |
515 | TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); | |
516 | signalF->SetParameters(10.,5.,fTau,fOrder,0.); //set all defaults once, just to be safe | |
517 | signalF->SetParNames("amp","t0","tau","N","ped"); | |
518 | signalF->FixParameter(2,fTau); // tau in units of time bin | |
519 | signalF->FixParameter(3,fOrder); // order | |
520 | signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here | |
521 | signalF->SetParameter(1, time); | |
522 | signalF->SetParameter(0, amp); | |
523 | ||
524 | gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points | |
525 | ||
526 | // assign fit results | |
527 | amp = signalF->GetParameter(0); | |
528 | time = signalF->GetParameter(1); | |
e9dbb64a | 529 | |
16605c06 | 530 | delete signalF; |
531 | ||
532 | // cross-check with ParabolaFit to see if the results make sense | |
533 | FitParabola(gSig, amp); // amp is possibly updated | |
82cbdfca | 534 | |
16605c06 | 535 | //printf("Std : Amp %f, time %g\n",amp, time); |
536 | ||
537 | break; | |
538 | }//kStandard Fitter | |
539 | //---------------------------- | |
540 | case kFastFit: | |
541 | { | |
542 | //printf("FastFitter \n"); | |
543 | Double_t eSignal = 1; // nominal 1 ADC error | |
544 | Double_t dAmp = amp; | |
545 | Double_t eAmp = 0; | |
546 | Double_t dTime = time; | |
547 | Double_t eTime = 0; | |
548 | Double_t chi2 = 0; | |
549 | ||
550 | AliCaloFastAltroFitv0::FastFit(gSig->GetX(), gSig->GetY(), gSig->GetN(), | |
551 | eSignal, fTau, | |
552 | dAmp, eAmp, dTime, eTime, chi2); | |
553 | amp = dAmp; | |
554 | time = dTime * GetRawFormatTimeBinWidth(); | |
555 | //printf("FastFitter: Amp %f, time %g, eAmp %f, eTimeEstimate %g, chi2 %f\n",amp, time,eAmp,eTime,chi2); | |
556 | break; | |
557 | } //kFastFit | |
558 | //---------------------------- | |
559 | ||
560 | case kNeuralNet: | |
561 | { | |
562 | // Extracts the amplitude and the time information using a Neural Network approach | |
563 | // Author P. La Rocca | |
564 | //printf("NeuralNet \n"); | |
565 | Double_t input[5] = {0.}; | |
566 | Double_t maxgraph = 0.; | |
567 | Int_t maxindex = -10; | |
568 | // Values for readback from input graph | |
569 | Double_t ttime = 0; | |
570 | Double_t signal = 0; | |
571 | ||
572 | for (Int_t i=0; i < gSig->GetN(); i++) { | |
573 | gSig->GetPoint(i, ttime, signal); | |
574 | if(signal > maxgraph) { | |
575 | maxgraph = signal; | |
576 | maxindex = i; | |
577 | } | |
e9dbb64a | 578 | } |
16605c06 | 579 | if((maxindex-2) > -1 && (maxindex-2)< gSig->GetN()) { |
580 | gSig->GetPoint(maxindex-2, ttime, input[0]); | |
581 | input[0] /= maxgraph; | |
e9dbb64a | 582 | } |
16605c06 | 583 | if((maxindex-1) > -1 && (maxindex-1)< gSig->GetN()) { |
584 | gSig->GetPoint(maxindex-1, ttime, input[1]); | |
585 | input[1] /= maxgraph; | |
c61fe3b4 | 586 | } |
16605c06 | 587 | if((maxindex) > -1 && (maxindex)< gSig->GetN()) { |
588 | gSig->GetPoint(maxindex, ttime, input[2]); | |
589 | input[2] /= maxgraph; | |
e9dbb64a | 590 | } |
16605c06 | 591 | if((maxindex+1) > -1 && (maxindex+1)< gSig->GetN()) { |
592 | gSig->GetPoint(maxindex+1, ttime, input[3]); | |
593 | input[3] /= maxgraph; | |
7643e728 | 594 | } |
16605c06 | 595 | if((maxindex+2) > -1 && (maxindex+2)< gSig->GetN()) { |
596 | gSig->GetPoint(maxindex+2, ttime, input[4]); | |
597 | input[4] /= maxgraph; | |
598 | } | |
599 | ||
600 | AliCaloNeuralFit exportNN; | |
601 | amp = exportNN.Value(0,input[0],input[1],input[2],input[3],input[4])*maxgraph; | |
602 | time = exportNN.Value(1,input[0],input[1],input[2],input[3],input[4])+maxindex; | |
603 | ||
604 | break; | |
605 | } //kNeuralNet | |
606 | //---------------------------- | |
607 | }//switch fitting algorithms | |
fb070798 | 608 | |
16605c06 | 609 | delete gSig; // delete TGraph |
5a056daa | 610 | |
16605c06 | 611 | return; |
612 | } | |
8cb998bd | 613 | |
16605c06 | 614 | //__________________________________________________________________ |
615 | void AliEMCALRawUtils::FitParabola(const TGraph *gSig, Float_t & amp) const | |
616 | { | |
617 | //BEG YS alternative methods to calculate the amplitude | |
618 | Double_t * ymx = gSig->GetX() ; | |
619 | Double_t * ymy = gSig->GetY() ; | |
620 | const Int_t kN = 3 ; | |
621 | Double_t ymMaxX[kN] = {0., 0., 0.} ; | |
622 | Double_t ymMaxY[kN] = {0., 0., 0.} ; | |
623 | Double_t ymax = 0. ; | |
624 | // find the maximum amplitude | |
625 | Int_t ymiMax = 0 ; | |
626 | for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) { | |
627 | if (ymy[ymi] > ymMaxY[0] ) { | |
628 | ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude | |
629 | ymMaxX[0] = ymx[ymi] ; | |
630 | ymiMax = ymi ; | |
631 | } | |
632 | } | |
633 | // find the maximum by fitting a parabola through the max and the two adjacent samples | |
634 | if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) { | |
635 | ymMaxY[1] = ymy[ymiMax+1] ; | |
636 | ymMaxY[2] = ymy[ymiMax-1] ; | |
637 | ymMaxX[1] = ymx[ymiMax+1] ; | |
638 | ymMaxX[2] = ymx[ymiMax-1] ; | |
639 | if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) { | |
640 | //fit a parabola through the 3 points y= a+bx+x*x*x | |
641 | Double_t sy = 0 ; | |
642 | Double_t sx = 0 ; | |
643 | Double_t sx2 = 0 ; | |
644 | Double_t sx3 = 0 ; | |
645 | Double_t sx4 = 0 ; | |
646 | Double_t sxy = 0 ; | |
647 | Double_t sx2y = 0 ; | |
648 | for (Int_t i = 0; i < kN ; i++) { | |
649 | sy += ymMaxY[i] ; | |
650 | sx += ymMaxX[i] ; | |
651 | sx2 += ymMaxX[i]*ymMaxX[i] ; | |
652 | sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; | |
653 | sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; | |
654 | sxy += ymMaxX[i]*ymMaxY[i] ; | |
655 | sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ; | |
656 | } | |
657 | Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2); | |
658 | Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ; | |
659 | Double_t c = cN / cD ; | |
660 | Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ; | |
661 | Double_t a = (sy-b*sx-c*sx2)/kN ; | |
662 | Double_t xmax = -b/(2*c) ; | |
663 | ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude | |
664 | } | |
665 | } | |
666 | ||
667 | Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ; | |
668 | if (diff > 0.1) | |
669 | amp = ymMaxY[0] ; | |
670 | //printf("Yves : Amp %f, time %g\n",amp, time); | |
671 | //END YS | |
ee299369 | 672 | return; |
673 | } | |
16605c06 | 674 | |
ee299369 | 675 | //__________________________________________________________________ |
676 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) | |
677 | { | |
8cb998bd | 678 | // Matches version used in 2007 beam test |
679 | // | |
ee299369 | 680 | // Shape of the electronics raw reponse: |
681 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
682 | // | |
7643e728 | 683 | // xx = (t - t0 + tau) / tau [xx is just a convenient help variable] |
684 | // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0 | |
685 | // F = 0 for xx < 0 | |
ee299369 | 686 | // |
687 | // parameters: | |
8cb998bd | 688 | // A: par[0] // Amplitude = peak value |
689 | // t0: par[1] | |
690 | // tau: par[2] | |
691 | // N: par[3] | |
692 | // ped: par[4] | |
ee299369 | 693 | // |
694 | Double_t signal ; | |
8cb998bd | 695 | Double_t tau =par[2]; |
e5bbbc4e | 696 | Double_t n =par[3]; |
8cb998bd | 697 | Double_t ped = par[4]; |
698 | Double_t xx = ( x[0] - par[1] + tau ) / tau ; | |
ee299369 | 699 | |
5a056daa | 700 | if (xx <= 0) |
8cb998bd | 701 | signal = ped ; |
ee299369 | 702 | else { |
e5bbbc4e | 703 | signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; |
ee299369 | 704 | } |
705 | return signal ; | |
706 | } | |
707 | ||
708 | //__________________________________________________________________ | |
709 | Bool_t AliEMCALRawUtils::RawSampledResponse( | |
710 | const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const | |
711 | { | |
712 | // for a start time dtime and an amplitude damp given by digit, | |
713 | // calculates the raw sampled response AliEMCAL::RawResponseFunction | |
714 | ||
ee299369 | 715 | Bool_t lowGain = kFALSE ; |
716 | ||
48a56166 | 717 | // A: par[0] // Amplitude = peak value |
718 | // t0: par[1] | |
719 | // tau: par[2] | |
720 | // N: par[3] | |
721 | // ped: par[4] | |
722 | ||
56e13066 | 723 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); |
48a56166 | 724 | signalF.SetParameter(0, damp) ; |
56e13066 | 725 | signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ; |
b4133f05 | 726 | signalF.SetParameter(2, fTau) ; |
727 | signalF.SetParameter(3, fOrder); | |
fe93d365 | 728 | signalF.SetParameter(4, fgPedestalValue); |
ee299369 | 729 | |
730 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { | |
56e13066 | 731 | Double_t signal = signalF.Eval(iTime) ; |
fe93d365 | 732 | |
7643e728 | 733 | // Next lines commeted for the moment but in principle it is not necessary to add |
ff10f540 | 734 | // extra noise since noise already added at the digits level. |
7643e728 | 735 | |
fe93d365 | 736 | //According to Terry Awes, 13-Apr-2008 |
737 | //add gaussian noise in quadrature to each sample | |
09974781 | 738 | //Double_t noise = gRandom->Gaus(0.,fgFEENoise); |
fe93d365 | 739 | //signal = sqrt(signal*signal + noise*noise); |
740 | ||
e2c2134b | 741 | // March 17,09 for fast fit simulations by Alexei Pavlinov. |
742 | // Get from PHOS analysis. In some sense it is open questions. | |
ff10f540 | 743 | //Double_t noise = gRandom->Gaus(0.,fgFEENoise); |
744 | //signal += noise; | |
7643e728 | 745 | |
ee299369 | 746 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; |
b4133f05 | 747 | if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits |
748 | adcH[iTime] = fgkRawSignalOverflow ; | |
ee299369 | 749 | lowGain = kTRUE ; |
750 | } | |
751 | ||
752 | signal /= fHighLowGainFactor; | |
753 | ||
754 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
b4133f05 | 755 | if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits |
756 | adcL[iTime] = fgkRawSignalOverflow ; | |
ee299369 | 757 | } |
758 | return lowGain ; | |
759 | } |