]>
Commit | Line | Data |
---|---|---|
168c7b3c | 1 | // -*- mode: c++ -*- |
ee299369 | 2 | /************************************************************************** |
3 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
4 | * * | |
5 | * Author: The ALICE Off-line Project. * | |
6 | * Contributors are mentioned in the code where appropriate. * | |
7 | * * | |
8 | * Permission to use, copy, modify and distribute this software and its * | |
9 | * documentation strictly for non-commercial purposes is hereby granted * | |
10 | * without fee, provided that the above copyright notice appears in all * | |
11 | * copies and that both the copyright notice and this permission notice * | |
12 | * appear in the supporting documentation. The authors make no claims * | |
13 | * about the suitability of this software for any purpose. It is * | |
14 | * provided "as is" without express or implied warranty. * | |
15 | **************************************************************************/ | |
16 | ||
17 | /* $Id$ */ | |
ee299369 | 18 | |
de39a0ff | 19 | |
e5bbbc4e | 20 | //_________________________________________________________________________ |
21 | // Utility Class for handling Raw data | |
22 | // Does all transitions from Digits to Raw and vice versa, | |
23 | // for simu and reconstruction | |
24 | // | |
25 | // Note: the current version is still simplified. Only | |
26 | // one raw signal per digit is generated; either high-gain or low-gain | |
27 | // Need to add concurrent high and low-gain info in the future | |
28 | // No pedestal is added to the raw signal. | |
ee299369 | 29 | //*-- Author: Marco van Leeuwen (LBL) |
e5bbbc4e | 30 | |
ee299369 | 31 | #include "AliEMCALRawUtils.h" |
507751ce | 32 | #include <stdexcept> |
21cad85c | 33 | |
4fe71e02 | 34 | #include "TF1.h" |
35 | #include "TGraph.h" | |
6751f762 | 36 | #include <TRandom.h> |
e5bbbc4e | 37 | class TSystem; |
21cad85c | 38 | |
e5bbbc4e | 39 | class AliLog; |
72c58de0 | 40 | #include "AliRun.h" |
ee299369 | 41 | #include "AliRunLoader.h" |
e5bbbc4e | 42 | class AliCaloAltroMapping; |
ee299369 | 43 | #include "AliAltroBuffer.h" |
44 | #include "AliRawReader.h" | |
32cd4c24 | 45 | #include "AliCaloRawStreamV3.h" |
ee299369 | 46 | #include "AliDAQ.h" |
21cad85c | 47 | |
feedcab9 | 48 | #include "AliEMCALRecParam.h" |
ee299369 | 49 | #include "AliEMCALLoader.h" |
50 | #include "AliEMCALGeometry.h" | |
e5bbbc4e | 51 | class AliEMCALDigitizer; |
ee299369 | 52 | #include "AliEMCALDigit.h" |
916f1e76 | 53 | #include "AliEMCALRawDigit.h" |
20b636fc | 54 | #include "AliEMCAL.h" |
5e3106bc | 55 | #include "AliCaloCalibPedestal.h" |
9f467289 | 56 | #include "AliCaloFastAltroFitv0.h" |
c8603a2b | 57 | #include "AliCaloNeuralFit.h" |
16605c06 | 58 | #include "AliCaloBunchInfo.h" |
59 | #include "AliCaloFitResults.h" | |
7683df1d | 60 | #include "AliCaloRawAnalyzerFastFit.h" |
61 | #include "AliCaloRawAnalyzerNN.h" | |
16605c06 | 62 | #include "AliCaloRawAnalyzerLMS.h" |
63 | #include "AliCaloRawAnalyzerPeakFinder.h" | |
64 | #include "AliCaloRawAnalyzerCrude.h" | |
de39a0ff | 65 | #include "AliEMCALTriggerRawDigitMaker.h" |
66 | #include "AliEMCALTriggerSTURawStream.h" | |
67 | #include "AliEMCALTriggerData.h" | |
9f467289 | 68 | |
ee299369 | 69 | ClassImp(AliEMCALRawUtils) |
21cad85c | 70 | |
ee299369 | 71 | // Signal shape parameters |
a2339fea | 72 | Int_t AliEMCALRawUtils::fgTimeBins = 256; // number of sampling bins of the raw RO signal (we typically use 15-50; theoretical max is 1k+) |
73 | Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns | |
74 | Double_t AliEMCALRawUtils::fgTimeTrigger = 600E-9 ; // the time of the trigger as approximately seen in the data | |
ee299369 | 75 | |
76 | // some digitization constants | |
a2339fea | 77 | Int_t AliEMCALRawUtils::fgThreshold = 1; |
ee299369 | 78 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule |
a2339fea | 79 | Int_t AliEMCALRawUtils::fgPedestalValue = 0; // pedestal value for digits2raw, default generate ZS data |
80 | Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled) | |
ee299369 | 81 | |
168c7b3c | 82 | AliEMCALRawUtils::AliEMCALRawUtils( Algo::fitAlgorithm fitAlgo) |
b4133f05 | 83 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), |
9f467289 | 84 | fNPedSamples(0), fGeom(0), fOption(""), |
f4a4dc82 | 85 | fRemoveBadChannels(kTRUE),fFittingAlgorithm(0), |
86 | fTimeMin(-1.),fTimeMax(1.), | |
de39a0ff | 87 | fUseFALTRO(kFALSE),fRawAnalyzer(0), |
88 | fTriggerRawDigitMaker(0x0) | |
8cb998bd | 89 | { |
b4133f05 | 90 | |
91 | //These are default parameters. | |
92 | //Can be re-set from without with setter functions | |
9f467289 | 93 | //Already set in the OCDB and passed via setter in the AliEMCALReconstructor |
46f1d25f | 94 | fHighLowGainFactor = 16. ; // Adjusted for a low gain range of 82 GeV (10 bits) |
95 | fOrder = 2; // Order of gamma fn | |
96 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
97 | fNoiseThreshold = 3; // 3 ADC counts is approx. noise level | |
98 | fNPedSamples = 4; // Less than this value => likely pedestal samples | |
99 | fRemoveBadChannels = kFALSE; // Do not remove bad channels before fitting | |
100 | fUseFALTRO = kTRUE; // Get the trigger FALTRO information and pass it to digits. | |
4fe71e02 | 101 | SetFittingAlgorithm(fitAlgo); |
16605c06 | 102 | |
65bdc82f | 103 | //Get Mapping RCU files from the AliEMCALRecParam |
104 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
105 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
106 | ||
21cad85c | 107 | for(Int_t i = 0; i < 4; i++) { |
65bdc82f | 108 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
109 | } | |
110 | ||
72c58de0 | 111 | //To make sure we match with the geometry in a simulation file, |
112 | //let's try to get it first. If not, take the default geometry | |
33c3c91a | 113 | AliRunLoader *rl = AliRunLoader::Instance(); |
a2339fea | 114 | if (rl && rl->GetAliRun()) { |
115 | AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL")); | |
116 | if(emcal)fGeom = emcal->GetGeometry(); | |
117 | else { | |
118 | AliDebug(1, Form("Using default geometry in raw reco")); | |
119 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); | |
120 | } | |
121 | ||
72c58de0 | 122 | } else { |
e853f058 | 123 | AliDebug(1, Form("Using default geometry in raw reco")); |
937d0661 | 124 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); |
65bdc82f | 125 | } |
126 | ||
72c58de0 | 127 | if(!fGeom) AliFatal(Form("Could not get geometry!")); |
de39a0ff | 128 | |
129 | fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker(); | |
72c58de0 | 130 | |
65bdc82f | 131 | } |
132 | ||
133 | //____________________________________________________________________________ | |
168c7b3c | 134 | AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry, Algo::fitAlgorithm fitAlgo) |
5544799a | 135 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), |
9f467289 | 136 | fNPedSamples(0), fGeom(pGeometry), fOption(""), |
f4a4dc82 | 137 | fRemoveBadChannels(kTRUE),fFittingAlgorithm(0), |
138 | fTimeMin(-1.),fTimeMax(1.), | |
de39a0ff | 139 | fUseFALTRO(kFALSE),fRawAnalyzer(), |
140 | fTriggerRawDigitMaker(0x0) | |
5544799a | 141 | { |
142 | // | |
143 | // Initialize with the given geometry - constructor required by HLT | |
144 | // HLT does not use/support AliRunLoader(s) instances | |
145 | // This is a minimum intervention solution | |
146 | // Comment by MPloskon@lbl.gov | |
147 | // | |
148 | ||
149 | //These are default parameters. | |
150 | //Can be re-set from without with setter functions | |
9f467289 | 151 | //Already set in the OCDB and passed via setter in the AliEMCALReconstructor |
46f1d25f | 152 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) |
153 | fOrder = 2; // order of gamma fn | |
154 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
155 | fNoiseThreshold = 3; // 3 ADC counts is approx. noise level | |
156 | fNPedSamples = 4; // Less than this value => likely pedestal samples | |
157 | fRemoveBadChannels = kFALSE; // Do not remove bad channels before fitting | |
158 | fUseFALTRO = kTRUE; // Get the trigger FALTRO information and pass it to digits. | |
4fe71e02 | 159 | SetFittingAlgorithm(fitAlgo); |
46f1d25f | 160 | |
5544799a | 161 | //Get Mapping RCU files from the AliEMCALRecParam |
162 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
163 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
164 | ||
21cad85c | 165 | for(Int_t i = 0; i < 4; i++) { |
5544799a | 166 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
167 | } | |
168 | ||
169 | if(!fGeom) AliFatal(Form("Could not get geometry!")); | |
de39a0ff | 170 | |
171 | fTriggerRawDigitMaker = new AliEMCALTriggerRawDigitMaker(); | |
5544799a | 172 | } |
173 | ||
174 | //____________________________________________________________________________ | |
65bdc82f | 175 | AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) |
176 | : TObject(), | |
177 | fHighLowGainFactor(rawU.fHighLowGainFactor), | |
b4133f05 | 178 | fOrder(rawU.fOrder), |
179 | fTau(rawU.fTau), | |
180 | fNoiseThreshold(rawU.fNoiseThreshold), | |
181 | fNPedSamples(rawU.fNPedSamples), | |
65bdc82f | 182 | fGeom(rawU.fGeom), |
9f467289 | 183 | fOption(rawU.fOption), |
184 | fRemoveBadChannels(rawU.fRemoveBadChannels), | |
16605c06 | 185 | fFittingAlgorithm(rawU.fFittingAlgorithm), |
f4a4dc82 | 186 | fTimeMin(rawU.fTimeMin),fTimeMax(rawU.fTimeMax), |
de39a0ff | 187 | fUseFALTRO(rawU.fUseFALTRO), |
188 | fRawAnalyzer(rawU.fRawAnalyzer), | |
189 | fTriggerRawDigitMaker(rawU.fTriggerRawDigitMaker) | |
65bdc82f | 190 | { |
191 | //copy ctor | |
192 | fMapping[0] = rawU.fMapping[0]; | |
193 | fMapping[1] = rawU.fMapping[1]; | |
21cad85c | 194 | fMapping[2] = rawU.fMapping[2]; |
195 | fMapping[3] = rawU.fMapping[3]; | |
65bdc82f | 196 | } |
197 | ||
198 | //____________________________________________________________________________ | |
199 | AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) | |
200 | { | |
201 | //assignment operator | |
202 | ||
203 | if(this != &rawU) { | |
204 | fHighLowGainFactor = rawU.fHighLowGainFactor; | |
46f1d25f | 205 | fOrder = rawU.fOrder; |
206 | fTau = rawU.fTau; | |
207 | fNoiseThreshold = rawU.fNoiseThreshold; | |
208 | fNPedSamples = rawU.fNPedSamples; | |
209 | fGeom = rawU.fGeom; | |
210 | fOption = rawU.fOption; | |
9f467289 | 211 | fRemoveBadChannels = rawU.fRemoveBadChannels; |
212 | fFittingAlgorithm = rawU.fFittingAlgorithm; | |
de39a0ff | 213 | fTimeMin = rawU.fTimeMin; |
214 | fTimeMax = rawU.fTimeMax; | |
46f1d25f | 215 | fUseFALTRO = rawU.fUseFALTRO; |
216 | fRawAnalyzer = rawU.fRawAnalyzer; | |
217 | fMapping[0] = rawU.fMapping[0]; | |
218 | fMapping[1] = rawU.fMapping[1]; | |
219 | fMapping[2] = rawU.fMapping[2]; | |
220 | fMapping[3] = rawU.fMapping[3]; | |
de39a0ff | 221 | fTriggerRawDigitMaker = rawU.fTriggerRawDigitMaker; |
65bdc82f | 222 | } |
223 | ||
224 | return *this; | |
225 | ||
ee299369 | 226 | } |
65bdc82f | 227 | |
ee299369 | 228 | //____________________________________________________________________________ |
229 | AliEMCALRawUtils::~AliEMCALRawUtils() { | |
e5bbbc4e | 230 | //dtor |
65bdc82f | 231 | |
ee299369 | 232 | } |
65bdc82f | 233 | |
ee299369 | 234 | //____________________________________________________________________________ |
65bdc82f | 235 | void AliEMCALRawUtils::Digits2Raw() |
ee299369 | 236 | { |
237 | // convert digits of the current event to raw data | |
238 | ||
33c3c91a | 239 | AliRunLoader *rl = AliRunLoader::Instance(); |
ee299369 | 240 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); |
a2339fea | 241 | |
ee299369 | 242 | // get the digits |
243 | loader->LoadDigits("EMCAL"); | |
244 | loader->GetEvent(); | |
245 | TClonesArray* digits = loader->Digits() ; | |
246 | ||
247 | if (!digits) { | |
248 | Warning("Digits2Raw", "no digits found !"); | |
249 | return; | |
250 | } | |
a2339fea | 251 | |
ee299369 | 252 | static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here |
253 | AliAltroBuffer* buffers[nDDL]; | |
254 | for (Int_t i=0; i < nDDL; i++) | |
255 | buffers[i] = 0; | |
a2339fea | 256 | |
e2c2134b | 257 | TArrayI adcValuesLow(fgTimeBins); |
258 | TArrayI adcValuesHigh(fgTimeBins); | |
de39a0ff | 259 | |
ee299369 | 260 | // loop over digits (assume ordered digits) |
261 | for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { | |
262 | AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ; | |
de39a0ff | 263 | if(!digit){ |
264 | AliFatal("NULL Digit"); | |
ee299369 | 265 | } |
de39a0ff | 266 | else{ |
267 | if (digit->GetAmplitude() < fgThreshold) | |
a2339fea | 268 | continue; |
de39a0ff | 269 | |
270 | //get cell indices | |
271 | Int_t nSM = 0; | |
272 | Int_t nIphi = 0; | |
273 | Int_t nIeta = 0; | |
274 | Int_t iphi = 0; | |
275 | Int_t ieta = 0; | |
276 | Int_t nModule = 0; | |
277 | fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); | |
278 | fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; | |
279 | ||
280 | //Check which is the RCU, 0 or 1, of the cell. | |
281 | Int_t iRCU = -111; | |
282 | //RCU0 | |
283 | if (0<=iphi&&iphi<8) iRCU=0; // first cable row | |
284 | else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; | |
285 | //second cable row | |
286 | //RCU1 | |
287 | else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; | |
288 | //second cable row | |
289 | else if(16<=iphi&&iphi<24) iRCU=1; // third cable row | |
290 | ||
291 | if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same | |
292 | ||
293 | if (iRCU<0) | |
a2339fea | 294 | Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); |
de39a0ff | 295 | |
296 | //Which DDL? | |
297 | Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; | |
a2339fea | 298 | if (iDDL < 0 || iDDL >= nDDL){ |
299 | Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); | |
de39a0ff | 300 | } |
a2339fea | 301 | else{ |
302 | if (buffers[iDDL] == 0) { | |
303 | // open new file and write dummy header | |
304 | TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); | |
305 | //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C | |
306 | Int_t iRCUside=iRCU+(nSM%2)*2; | |
307 | //iRCU=0 and even (0) SM -> RCU0A.data 0 | |
308 | //iRCU=1 and even (0) SM -> RCU1A.data 1 | |
309 | //iRCU=0 and odd (1) SM -> RCU0C.data 2 | |
310 | //iRCU=1 and odd (1) SM -> RCU1C.data 3 | |
311 | //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl; | |
312 | buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]); | |
313 | buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; | |
314 | } | |
315 | ||
316 | // out of time range signal (?) | |
317 | if (digit->GetTimeR() > GetRawFormatTimeMax() ) { | |
318 | AliInfo("Signal is out of time range.\n"); | |
319 | buffers[iDDL]->FillBuffer((Int_t)digit->GetAmplitude()); | |
320 | buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin | |
321 | buffers[iDDL]->FillBuffer(3); // bunch length | |
322 | buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer | |
323 | // calculate the time response function | |
324 | } else { | |
325 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmplitude(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; | |
326 | if (lowgain) | |
327 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold); | |
328 | else | |
329 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold); | |
330 | } | |
331 | }// iDDL under the limits | |
de39a0ff | 332 | }//digit exists |
333 | }//Digit loop | |
a2339fea | 334 | |
ee299369 | 335 | // write headers and close files |
336 | for (Int_t i=0; i < nDDL; i++) { | |
337 | if (buffers[i]) { | |
338 | buffers[i]->Flush(); | |
339 | buffers[i]->WriteDataHeader(kFALSE, kFALSE); | |
340 | delete buffers[i]; | |
341 | } | |
342 | } | |
de39a0ff | 343 | |
ee299369 | 344 | loader->UnloadDigits(); |
345 | } | |
346 | ||
347 | //____________________________________________________________________________ | |
de39a0ff | 348 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap, TClonesArray *digitsTRG, AliEMCALTriggerData* trgData) |
ee299369 | 349 | { |
65bdc82f | 350 | // convert raw data of the current event to digits |
a2339fea | 351 | |
de39a0ff | 352 | if(digitsArr) digitsArr->Clear("C"); |
a2339fea | 353 | |
c47157cd | 354 | if (!digitsArr) { |
ee299369 | 355 | Error("Raw2Digits", "no digits found !"); |
356 | return; | |
357 | } | |
358 | if (!reader) { | |
359 | Error("Raw2Digits", "no raw reader found !"); | |
360 | return; | |
361 | } | |
de39a0ff | 362 | |
363 | AliEMCALTriggerSTURawStream inSTU(reader); | |
364 | ||
365 | AliCaloRawStreamV3 in(reader,"EMCAL",fMapping); | |
366 | ||
ee299369 | 367 | // Select EMCAL DDL's; |
7643e728 | 368 | reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL |
a2339fea | 369 | |
de39a0ff | 370 | fTriggerRawDigitMaker->Reset(); |
371 | fTriggerRawDigitMaker->SetIO(reader, in, inSTU, digitsTRG, trgData); | |
372 | ||
16605c06 | 373 | // fRawAnalyzer setup |
2cd0ffda | 374 | fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods |
375 | fRawAnalyzer->SetOverflowCut(fgkOverflowCut); | |
16605c06 | 376 | fRawAnalyzer->SetAmpCut(fNoiseThreshold); |
377 | fRawAnalyzer->SetFitArrayCut(fNoiseThreshold); | |
378 | fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later | |
a2339fea | 379 | |
16605c06 | 380 | // channel info parameters |
829ba234 | 381 | Int_t lowGain = 0; |
e5bbbc4e | 382 | Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref. |
a2339fea | 383 | |
32cd4c24 | 384 | // start loop over input stream |
385 | while (in.NextDDL()) { | |
916f1e76 | 386 | |
a2339fea | 387 | // if ( in.GetDDLNumber() != 0 && in.GetDDLNumber() != 2 ) continue; |
388 | ||
32cd4c24 | 389 | while (in.NextChannel()) { |
a2339fea | 390 | |
7643e728 | 391 | //Check if the signal is high or low gain and then do the fit, |
16605c06 | 392 | //if it is from TRU or LEDMon do not fit |
7643e728 | 393 | caloFlag = in.GetCaloFlag(); |
a2339fea | 394 | // if (caloFlag != 0 && caloFlag != 1) continue; |
395 | if (caloFlag > 2) continue; // Work with ALTRO and FALTRO | |
396 | ||
916f1e76 | 397 | //Do not fit bad channels of ALTRO |
398 | if(caloFlag < 2 && fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) { | |
a2339fea | 399 | //printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow()); |
400 | continue; | |
5e3106bc | 401 | } |
a2339fea | 402 | |
16605c06 | 403 | vector<AliCaloBunchInfo> bunchlist; |
32cd4c24 | 404 | while (in.NextBunch()) { |
a2339fea | 405 | bunchlist.push_back( AliCaloBunchInfo(in.GetStartTimeBin(), in.GetBunchLength(), in.GetSignals() ) ); |
16605c06 | 406 | } // loop over bunches |
0a851600 | 407 | if (bunchlist.size() == 0) continue; |
168c7b3c | 408 | if ( caloFlag < 2 ) |
409 | { // ALTRO | |
410 | Float_t time = 0; | |
411 | Float_t amp = 0; | |
412 | short timeEstimate = 0; | |
413 | Float_t ampEstimate = 0; | |
414 | Bool_t fitDone = kFALSE; | |
415 | Float_t chi2 = 0; | |
416 | Int_t ndf = 0; | |
417 | ||
418 | if ( fFittingAlgorithm == Algo::kFastFit || fFittingAlgorithm == Algo::kNeuralNet || | |
419 | fFittingAlgorithm == Algo::kLMS || fFittingAlgorithm == Algo::kPeakFinder || | |
420 | fFittingAlgorithm == Algo::kCrude) { | |
a2339fea | 421 | // all functionality to determine amp and time etc is encapsulated inside the Evaluate call for these methods |
422 | AliCaloFitResults fitResults = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2()); | |
423 | ||
424 | amp = fitResults.GetAmp(); | |
425 | time = fitResults.GetTime(); | |
426 | timeEstimate = fitResults.GetMaxTimebin(); | |
427 | ampEstimate = fitResults.GetMaxSig(); | |
428 | chi2 = fitResults.GetChi2(); | |
429 | ndf = fitResults.GetNdf(); | |
168c7b3c | 430 | if (fitResults.GetStatus() == Ret::kFitPar) { |
a2339fea | 431 | fitDone = kTRUE; |
432 | } | |
433 | } | |
434 | else { // for the other methods we for now use the functionality of | |
435 | // AliCaloRawAnalyzer as well, to select samples and prepare for fits, | |
436 | // if it looks like there is something to fit | |
437 | ||
438 | // parameters init. | |
439 | Float_t pedEstimate = 0; | |
440 | short maxADC = 0; | |
441 | Int_t first = 0; | |
442 | Int_t last = 0; | |
443 | Int_t bunchIndex = 0; | |
444 | // | |
445 | // The PreFitEvaluateSamples + later call to FitRaw will hopefully | |
446 | // be replaced by a single Evaluate call or so soon, like for the other | |
447 | // methods, but this should be good enough for evaluation of | |
448 | // the methods for now (Jan. 2010) | |
449 | // | |
450 | int nsamples = fRawAnalyzer->PreFitEvaluateSamples( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2(), bunchIndex, ampEstimate, maxADC, timeEstimate, pedEstimate, first, last); | |
451 | ||
452 | if (ampEstimate >= fNoiseThreshold) { // something worth looking at | |
453 | ||
454 | time = timeEstimate; // maxrev in AliCaloRawAnalyzer speak; comes with an offset w.r.t. real timebin | |
455 | Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1); | |
456 | amp = ampEstimate; | |
457 | ||
458 | if ( nsamples > 1 && maxADC<fgkOverflowCut ) { // possibly something to fit | |
459 | FitRaw(first, last, amp, time, chi2, fitDone); | |
460 | time += timebinOffset; | |
461 | timeEstimate += timebinOffset; | |
462 | ndf = nsamples - 2; | |
463 | } | |
464 | ||
465 | } // ampEstimate check | |
466 | } // method selection | |
467 | ||
468 | if ( fitDone ) { // brief sanity check of fit results | |
469 | Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate); | |
470 | Float_t timeDiff = time - timeEstimate; | |
471 | if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) ) { | |
472 | // AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations amp %f time %d", amp, time, ampEstimate, timeEstimate)); | |
473 | ||
474 | // for now just overwrite the fit results with the simple/initial estimate | |
475 | amp = ampEstimate; | |
476 | time = timeEstimate; | |
477 | fitDone = kFALSE; | |
478 | } | |
479 | } // fitDone | |
480 | ||
481 | if (amp >= fNoiseThreshold) { // something to be stored | |
482 | if ( ! fitDone) { // smear ADC with +- 0.5 uniform (avoid discrete effects) | |
483 | amp += (0.5 - gRandom->Rndm()); // Rndm generates a number in ]0,1] | |
484 | } | |
485 | ||
486 | Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; | |
487 | lowGain = in.IsLowGain(); | |
488 | ||
489 | // go from time-bin units to physical time fgtimetrigger | |
490 | time = time * GetRawFormatTimeBinWidth(); // skip subtraction of fgTimeTrigger? | |
491 | // subtract RCU L1 phase (L1Phase is in seconds) w.r.t. L0: | |
492 | time -= in.GetL1Phase(); | |
493 | ||
494 | AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); | |
495 | // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp); | |
496 | ||
497 | AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf); | |
498 | } | |
499 | ||
500 | }//ALTRO | |
501 | else if(fUseFALTRO) | |
502 | {// Fake ALTRO | |
503 | fTriggerRawDigitMaker->Add( bunchlist ); | |
504 | }//Fake ALTRO | |
505 | } // end while over channel | |
32cd4c24 | 506 | } //end while over DDL's, of input stream |
a2339fea | 507 | |
de39a0ff | 508 | fTriggerRawDigitMaker->PostProcess(); |
509 | ||
f4a4dc82 | 510 | TrimDigits(digitsArr); |
511 | ||
ee299369 | 512 | return ; |
513 | } | |
514 | ||
515 | //____________________________________________________________________________ | |
2cd0ffda | 516 | void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf) { |
82cbdfca | 517 | // |
518 | // Add a new digit. | |
519 | // This routine checks whether a digit exists already for this tower | |
520 | // and then decides whether to use the high or low gain info | |
521 | // | |
522 | // Called by Raw2Digits | |
523 | ||
524 | AliEMCALDigit *digit = 0, *tmpdigit = 0; | |
82cbdfca | 525 | TIter nextdigit(digitsArr); |
526 | while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { | |
f4a4dc82 | 527 | if (tmpdigit->GetId() == id) digit = tmpdigit; |
82cbdfca | 528 | } |
529 | ||
530 | if (!digit) { // no digit existed for this tower; create one | |
f4a4dc82 | 531 | Int_t type = AliEMCALDigit::kHG; // use enum in AliEMCALDigit |
532 | if (lowGain) { | |
533 | amp *= fHighLowGainFactor; | |
534 | type = AliEMCALDigit::kLGnoHG; | |
535 | } | |
536 | Int_t idigit = digitsArr->GetEntries(); | |
2cd0ffda | 537 | new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit, chi2, ndf); |
f4a4dc82 | 538 | AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type)); |
539 | }//digit added first time | |
82cbdfca | 540 | else { // a digit already exists, check range |
f4a4dc82 | 541 | // (use high gain if signal < cut value, otherwise low gain) |
542 | if (lowGain) { // new digit is low gain | |
543 | if (digit->GetAmplitude() > fgkOverflowCut) { // use if previously stored (HG) digit is out of range | |
544 | digit->SetAmplitude(fHighLowGainFactor * amp); | |
545 | digit->SetTime(time); | |
546 | digit->SetType(AliEMCALDigit::kLG); | |
547 | AliDebug(2,Form("Add LG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType())); | |
548 | } | |
549 | }//new low gain digit | |
550 | else { // new digit is high gain | |
551 | if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range | |
552 | digit->SetAmplitude(amp); | |
553 | digit->SetTime(time); | |
554 | digit->SetType(AliEMCALDigit::kHG); | |
555 | AliDebug(2,Form("Add HG digit ID %d for the second time, type %d", digit->GetId(), digit->GetType())); | |
556 | } | |
557 | else { // HG out of range, just change flag value to show that HG did exist | |
558 | digit->SetType(AliEMCALDigit::kLG); | |
559 | AliDebug(2,Form("Change LG digit to HG, ID %d, type %d", digit->GetId(), digit->GetType())); | |
560 | } | |
561 | }//new high gain digit | |
562 | }//digit existed replace it | |
563 | ||
564 | } | |
565 | ||
566 | //____________________________________________________________________________ | |
567 | void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr) | |
568 | { | |
569 | // Remove digits with only low gain and large time | |
570 | ||
571 | AliEMCALDigit *digit = 0; | |
572 | Int_t n = 0; | |
573 | Int_t nDigits = digitsArr->GetEntriesFast(); | |
574 | TIter nextdigit(digitsArr); | |
575 | while ((digit = (AliEMCALDigit*) nextdigit())) { | |
a2339fea | 576 | |
f4a4dc82 | 577 | //Check if only LG existed, remove if so |
578 | if (digit->GetType() == AliEMCALDigit::kLGnoHG) { | |
579 | AliDebug(1,Form("Remove digit with id %d, LGnoHG",digit->GetId())); | |
580 | digitsArr->Remove(digit); | |
82cbdfca | 581 | } |
2cd0ffda | 582 | //Check if time is too large or too small, remove if so |
f4a4dc82 | 583 | else if(fTimeMin > digit->GetTime() || fTimeMax < digit->GetTime()) { |
584 | digitsArr->Remove(digit); | |
585 | AliDebug(1,Form("Remove digit with id %d, Bad Time %e",digit->GetId(), digit->GetTime())); | |
82cbdfca | 586 | } |
2cd0ffda | 587 | // Check if Chi2 is undefined |
588 | else if (0 > digit->GetChi2()) { | |
589 | digitsArr->Remove(digit); | |
590 | AliDebug(1,Form("Remove digit with id %d, Bad Chi2 %e",digit->GetId(), digit->GetChi2())); | |
591 | } | |
f4a4dc82 | 592 | //Good digit, just reassign the index of the digit in case there was a previous removal |
593 | else { | |
594 | digit->SetIndexInList(n); | |
595 | n++; | |
596 | } | |
597 | }//while | |
598 | ||
599 | digitsArr->Compress(); | |
600 | AliDebug(1,Form("N Digits before trimming : %d; after array compression %d",nDigits,digitsArr->GetEntriesFast())); | |
601 | ||
82cbdfca | 602 | } |
f4a4dc82 | 603 | |
82cbdfca | 604 | //____________________________________________________________________________ |
2cd0ffda | 605 | void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const |
16605c06 | 606 | { // Fits the raw signal time distribution |
607 | ||
608 | //-------------------------------------------------- | |
609 | //Do the fit, different fitting algorithms available | |
610 | //-------------------------------------------------- | |
611 | int nsamples = lastTimeBin - firstTimeBin + 1; | |
507751ce | 612 | fitDone = kFALSE; |
a2339fea | 613 | |
16605c06 | 614 | switch(fFittingAlgorithm) { |
168c7b3c | 615 | case Algo::kStandard: |
16605c06 | 616 | { |
7683df1d | 617 | if (nsamples < 3) { return; } // nothing much to fit |
16605c06 | 618 | //printf("Standard fitter \n"); |
7683df1d | 619 | |
16605c06 | 620 | // Create Graph to hold data we will fit |
7683df1d | 621 | TGraph *gSig = new TGraph( nsamples); |
622 | for (int i=0; i<nsamples; i++) { | |
a2339fea | 623 | Int_t timebin = firstTimeBin + i; |
624 | gSig->SetPoint(i, timebin, fRawAnalyzer->GetReversed(timebin)); | |
7683df1d | 625 | } |
a2339fea | 626 | |
16605c06 | 627 | TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); |
628 | signalF->SetParameters(10.,5.,fTau,fOrder,0.); //set all defaults once, just to be safe | |
629 | signalF->SetParNames("amp","t0","tau","N","ped"); | |
630 | signalF->FixParameter(2,fTau); // tau in units of time bin | |
631 | signalF->FixParameter(3,fOrder); // order | |
632 | signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here | |
633 | signalF->SetParameter(1, time); | |
634 | signalF->SetParameter(0, amp); | |
507751ce | 635 | // set rather loose parameter limits |
636 | signalF->SetParLimits(0, 0.5*amp, 2*amp ); | |
637 | signalF->SetParLimits(1, time - 4, time + 4); | |
a2339fea | 638 | |
507751ce | 639 | try { |
a2339fea | 640 | gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points |
641 | // assign fit results | |
642 | amp = signalF->GetParameter(0); | |
643 | time = signalF->GetParameter(1); | |
644 | chi2 = signalF->GetChisquare(); | |
645 | fitDone = kTRUE; | |
507751ce | 646 | } |
647 | catch (const std::exception & e) { | |
a2339fea | 648 | AliError( Form("TGraph Fit exception %s", e.what()) ); |
649 | // stay with default amp and time in case of exception, i.e. no special action required | |
650 | fitDone = kFALSE; | |
507751ce | 651 | } |
16605c06 | 652 | delete signalF; |
a2339fea | 653 | |
16605c06 | 654 | //printf("Std : Amp %f, time %g\n",amp, time); |
7683df1d | 655 | delete gSig; // delete TGraph |
a2339fea | 656 | |
16605c06 | 657 | break; |
658 | }//kStandard Fitter | |
a2339fea | 659 | //---------------------------- |
168c7b3c | 660 | case Algo::kLogFit: |
16605c06 | 661 | { |
7683df1d | 662 | if (nsamples < 3) { return; } // nothing much to fit |
663 | //printf("LogFit \n"); | |
a2339fea | 664 | |
7683df1d | 665 | // Create Graph to hold data we will fit |
666 | TGraph *gSigLog = new TGraph( nsamples); | |
667 | for (int i=0; i<nsamples; i++) { | |
a2339fea | 668 | Int_t timebin = firstTimeBin + i; |
669 | gSigLog->SetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) ); | |
7643e728 | 670 | } |
a2339fea | 671 | |
7683df1d | 672 | TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0, GetRawFormatTimeBins(), 5); |
673 | signalFLog->SetParameters(2.3, 5.,fTau,fOrder,0.); //set all defaults once, just to be safe | |
674 | signalFLog->SetParNames("amplog","t0","tau","N","ped"); | |
675 | signalFLog->FixParameter(2,fTau); // tau in units of time bin | |
676 | signalFLog->FixParameter(3,fOrder); // order | |
677 | signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here | |
678 | signalFLog->SetParameter(1, time); | |
679 | if (amp>=1) { | |
a2339fea | 680 | signalFLog->SetParameter(0, TMath::Log(amp)); |
16605c06 | 681 | } |
a2339fea | 682 | |
7683df1d | 683 | gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points |
a2339fea | 684 | |
7683df1d | 685 | // assign fit results |
686 | Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp | |
687 | amp = TMath::Exp(amplog); | |
688 | time = signalFLog->GetParameter(1); | |
507751ce | 689 | fitDone = kTRUE; |
a2339fea | 690 | |
7683df1d | 691 | delete signalFLog; |
692 | //printf("LogFit: Amp %f, time %g\n",amp, time); | |
693 | delete gSigLog; | |
16605c06 | 694 | break; |
7683df1d | 695 | } //kLogFit |
a2339fea | 696 | //---------------------------- |
697 | ||
698 | //---------------------------- | |
16605c06 | 699 | }//switch fitting algorithms |
a2339fea | 700 | |
16605c06 | 701 | return; |
702 | } | |
8cb998bd | 703 | |
16605c06 | 704 | //__________________________________________________________________ |
705 | void AliEMCALRawUtils::FitParabola(const TGraph *gSig, Float_t & amp) const | |
706 | { | |
707 | //BEG YS alternative methods to calculate the amplitude | |
708 | Double_t * ymx = gSig->GetX() ; | |
709 | Double_t * ymy = gSig->GetY() ; | |
710 | const Int_t kN = 3 ; | |
711 | Double_t ymMaxX[kN] = {0., 0., 0.} ; | |
712 | Double_t ymMaxY[kN] = {0., 0., 0.} ; | |
713 | Double_t ymax = 0. ; | |
714 | // find the maximum amplitude | |
715 | Int_t ymiMax = 0 ; | |
716 | for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) { | |
717 | if (ymy[ymi] > ymMaxY[0] ) { | |
718 | ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude | |
719 | ymMaxX[0] = ymx[ymi] ; | |
720 | ymiMax = ymi ; | |
721 | } | |
722 | } | |
723 | // find the maximum by fitting a parabola through the max and the two adjacent samples | |
724 | if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) { | |
725 | ymMaxY[1] = ymy[ymiMax+1] ; | |
726 | ymMaxY[2] = ymy[ymiMax-1] ; | |
727 | ymMaxX[1] = ymx[ymiMax+1] ; | |
728 | ymMaxX[2] = ymx[ymiMax-1] ; | |
729 | if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) { | |
730 | //fit a parabola through the 3 points y= a+bx+x*x*x | |
731 | Double_t sy = 0 ; | |
732 | Double_t sx = 0 ; | |
733 | Double_t sx2 = 0 ; | |
734 | Double_t sx3 = 0 ; | |
735 | Double_t sx4 = 0 ; | |
736 | Double_t sxy = 0 ; | |
737 | Double_t sx2y = 0 ; | |
738 | for (Int_t i = 0; i < kN ; i++) { | |
a2339fea | 739 | sy += ymMaxY[i] ; |
740 | sx += ymMaxX[i] ; | |
741 | sx2 += ymMaxX[i]*ymMaxX[i] ; | |
742 | sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; | |
743 | sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; | |
744 | sxy += ymMaxX[i]*ymMaxY[i] ; | |
745 | sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ; | |
16605c06 | 746 | } |
747 | Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2); | |
748 | Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ; | |
749 | Double_t c = cN / cD ; | |
750 | Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ; | |
751 | Double_t a = (sy-b*sx-c*sx2)/kN ; | |
752 | Double_t xmax = -b/(2*c) ; | |
753 | ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude | |
029fe7a2 | 754 | amp = ymax; |
16605c06 | 755 | } |
756 | } | |
757 | ||
758 | Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ; | |
759 | if (diff > 0.1) | |
760 | amp = ymMaxY[0] ; | |
761 | //printf("Yves : Amp %f, time %g\n",amp, time); | |
762 | //END YS | |
ee299369 | 763 | return; |
764 | } | |
16605c06 | 765 | |
ee299369 | 766 | //__________________________________________________________________ |
767 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) | |
768 | { | |
8cb998bd | 769 | // Matches version used in 2007 beam test |
770 | // | |
ee299369 | 771 | // Shape of the electronics raw reponse: |
772 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
773 | // | |
7643e728 | 774 | // xx = (t - t0 + tau) / tau [xx is just a convenient help variable] |
775 | // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0 | |
776 | // F = 0 for xx < 0 | |
ee299369 | 777 | // |
778 | // parameters: | |
8cb998bd | 779 | // A: par[0] // Amplitude = peak value |
780 | // t0: par[1] | |
781 | // tau: par[2] | |
782 | // N: par[3] | |
783 | // ped: par[4] | |
ee299369 | 784 | // |
de39a0ff | 785 | Double_t signal = 0.; |
786 | Double_t tau = par[2]; | |
787 | Double_t n = par[3]; | |
788 | Double_t ped = par[4]; | |
789 | Double_t xx = ( x[0] - par[1] + tau ) / tau ; | |
ee299369 | 790 | |
5a056daa | 791 | if (xx <= 0) |
8cb998bd | 792 | signal = ped ; |
ee299369 | 793 | else { |
e5bbbc4e | 794 | signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; |
ee299369 | 795 | } |
796 | return signal ; | |
797 | } | |
798 | ||
7683df1d | 799 | //__________________________________________________________________ |
800 | Double_t AliEMCALRawUtils::RawResponseFunctionLog(Double_t *x, Double_t *par) | |
801 | { | |
802 | // Matches version used in 2007 beam test | |
803 | // | |
804 | // Shape of the electronics raw reponse: | |
805 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
806 | // | |
807 | // xx = (t - t0 + tau) / tau [xx is just a convenient help variable] | |
808 | // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0 | |
809 | // F = 0 for xx < 0 | |
810 | // | |
811 | // parameters: | |
812 | // Log[A]: par[0] // Amplitude = peak value | |
813 | // t0: par[1] | |
814 | // tau: par[2] | |
815 | // N: par[3] | |
816 | // ped: par[4] | |
817 | // | |
53e430a3 | 818 | Double_t signal = 0. ; |
de39a0ff | 819 | Double_t tau = par[2]; |
820 | Double_t n = par[3]; | |
7683df1d | 821 | //Double_t ped = par[4]; // not used |
de39a0ff | 822 | Double_t xx = ( x[0] - par[1] + tau ) / tau ; |
7683df1d | 823 | |
824 | if (xx < 0) | |
825 | signal = par[0] - n*TMath::Log(TMath::Abs(xx)) + n * (1 - xx ) ; | |
826 | else { | |
827 | signal = par[0] + n*TMath::Log(xx) + n * (1 - xx ) ; | |
828 | } | |
829 | return signal ; | |
830 | } | |
831 | ||
ee299369 | 832 | //__________________________________________________________________ |
6751f762 | 833 | Bool_t AliEMCALRawUtils::RawSampledResponse(const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL, const Int_t keyErr) const |
ee299369 | 834 | { |
835 | // for a start time dtime and an amplitude damp given by digit, | |
836 | // calculates the raw sampled response AliEMCAL::RawResponseFunction | |
ee299369 | 837 | Bool_t lowGain = kFALSE ; |
a2339fea | 838 | |
48a56166 | 839 | // A: par[0] // Amplitude = peak value |
840 | // t0: par[1] | |
841 | // tau: par[2] | |
842 | // N: par[3] | |
843 | // ped: par[4] | |
a2339fea | 844 | |
56e13066 | 845 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); |
48a56166 | 846 | signalF.SetParameter(0, damp) ; |
56e13066 | 847 | signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ; |
b4133f05 | 848 | signalF.SetParameter(2, fTau) ; |
849 | signalF.SetParameter(3, fOrder); | |
fe93d365 | 850 | signalF.SetParameter(4, fgPedestalValue); |
6751f762 | 851 | |
852 | Double_t signal=0.0, noise=0.0; | |
ee299369 | 853 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { |
916f1e76 | 854 | signal = signalF.Eval(iTime) ; |
7643e728 | 855 | // Next lines commeted for the moment but in principle it is not necessary to add |
ff10f540 | 856 | // extra noise since noise already added at the digits level. |
a2339fea | 857 | |
fe93d365 | 858 | //According to Terry Awes, 13-Apr-2008 |
859 | //add gaussian noise in quadrature to each sample | |
09974781 | 860 | //Double_t noise = gRandom->Gaus(0.,fgFEENoise); |
fe93d365 | 861 | //signal = sqrt(signal*signal + noise*noise); |
a2339fea | 862 | |
e2c2134b | 863 | // March 17,09 for fast fit simulations by Alexei Pavlinov. |
4fe71e02 | 864 | // Get from PHOS analysis. In some sense it is open questions. |
a2339fea | 865 | if(keyErr>0) { |
866 | noise = gRandom->Gaus(0.,fgFEENoise); | |
867 | signal += noise; | |
868 | } | |
6751f762 | 869 | |
ee299369 | 870 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; |
b4133f05 | 871 | if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits |
872 | adcH[iTime] = fgkRawSignalOverflow ; | |
ee299369 | 873 | lowGain = kTRUE ; |
874 | } | |
a2339fea | 875 | |
ee299369 | 876 | signal /= fHighLowGainFactor; |
a2339fea | 877 | |
ee299369 | 878 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; |
b4133f05 | 879 | if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits |
880 | adcL[iTime] = fgkRawSignalOverflow ; | |
a2339fea | 881 | |
ee299369 | 882 | } |
a2339fea | 883 | |
ee299369 | 884 | return lowGain ; |
885 | } | |
4fe71e02 | 886 | |
829ba234 | 887 | //__________________________________________________________________ |
3a187ba0 | 888 | void AliEMCALRawUtils::CalculateChi2(const Double_t* t, const Double_t* y, const Int_t nPoints, |
889 | const Double_t sig, const Double_t tau, const Double_t amp, const Double_t t0, Double_t &chi2) | |
890 | { | |
891 | // Input: | |
892 | // t[] - array of time bins | |
893 | // y[] - array of amplitudes after pedestal subtractions; | |
894 | // nPoints - number of points | |
895 | // sig - error of amplitude measurement (one value for all channels) | |
896 | // if sig<0 that mean sig=1. | |
897 | // tau - filter time response (in timebin units) | |
898 | // amp - amplitude at t0; | |
899 | // t0 - time of max amplitude; | |
900 | // Output: | |
901 | // chi2 - chi2 | |
902 | // ndf = nPoints - 2 when tau fixed | |
903 | // ndf = nPoints - 3 when tau free | |
904 | static Double_t par[5]={0.0, 0.0, 0.0, 2.0, 0.0}; | |
905 | ||
906 | par[0] = amp; | |
907 | par[1] = t0; | |
908 | par[2] = tau; | |
909 | // par[3]=n=2.; par[4]=ped=0.0 | |
910 | ||
911 | Double_t dy = 0.0, x = 0.0, f=0.0; | |
912 | for(Int_t i=0; i<nPoints; i++){ | |
913 | x = t[i]; | |
914 | f = RawResponseFunction(&x, par); | |
915 | dy = y[i] - f; | |
916 | chi2 += dy*dy; | |
a2339fea | 917 | //printf(" AliEMCALRawUtils::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, y[i], f, dy); |
3a187ba0 | 918 | } |
919 | if(sig>0.0) chi2 /= (sig*sig); | |
920 | } | |
921 | ||
4fe71e02 | 922 | //__________________________________________________________________ |
923 | void AliEMCALRawUtils::SetFittingAlgorithm(Int_t fitAlgo) | |
924 | { | |
925 | //Set fitting algorithm and initialize it if this same algorithm was not set before. | |
916f1e76 | 926 | //printf("**** Set Algorithm , number %d ****\n",fitAlgo); |
927 | ||
4fe71e02 | 928 | if(fitAlgo == fFittingAlgorithm && fRawAnalyzer) { |
929 | //Do nothing, this same algorithm already set before. | |
930 | //printf("**** Algorithm already set before, number %d, %s ****\n",fitAlgo, fRawAnalyzer->GetName()); | |
931 | return; | |
932 | } | |
933 | //Initialize the requested algorithm | |
934 | if(fitAlgo != fFittingAlgorithm || !fRawAnalyzer) { | |
935 | //printf("**** Init Algorithm , number %d ****\n",fitAlgo); | |
936 | ||
937 | fFittingAlgorithm = fitAlgo; | |
938 | if (fRawAnalyzer) delete fRawAnalyzer; // delete prev. analyzer if existed. | |
939 | ||
168c7b3c | 940 | if (fitAlgo == Algo::kFastFit) { |
4fe71e02 | 941 | fRawAnalyzer = new AliCaloRawAnalyzerFastFit(); |
942 | } | |
168c7b3c | 943 | else if (fitAlgo == Algo::kNeuralNet) { |
4fe71e02 | 944 | fRawAnalyzer = new AliCaloRawAnalyzerNN(); |
945 | } | |
168c7b3c | 946 | else if (fitAlgo == Algo::kLMS) { |
4fe71e02 | 947 | fRawAnalyzer = new AliCaloRawAnalyzerLMS(); |
948 | } | |
168c7b3c | 949 | else if (fitAlgo == Algo::kPeakFinder) { |
4fe71e02 | 950 | fRawAnalyzer = new AliCaloRawAnalyzerPeakFinder(); |
951 | } | |
168c7b3c | 952 | else if (fitAlgo == Algo::kCrude) { |
4fe71e02 | 953 | fRawAnalyzer = new AliCaloRawAnalyzerCrude(); |
954 | } | |
955 | else { | |
168c7b3c | 956 | // fRawAnalyzer = new AliCaloRawAnalyzer(); |
957 | fRawAnalyzer = 0; | |
4fe71e02 | 958 | } |
959 | } | |
960 | ||
961 | } | |
962 | ||
963 |