]>
Commit | Line | Data |
---|---|---|
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$ */ | |
17 | ||
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. | |
27 | //*-- Author: Marco van Leeuwen (LBL) | |
28 | ||
29 | #include "AliEMCALRawUtils.h" | |
30 | ||
31 | #include "TF1.h" | |
32 | #include "TGraph.h" | |
33 | #include <TRandom.h> | |
34 | class TSystem; | |
35 | ||
36 | class AliLog; | |
37 | #include "AliRun.h" | |
38 | #include "AliRunLoader.h" | |
39 | class AliCaloAltroMapping; | |
40 | #include "AliAltroBuffer.h" | |
41 | #include "AliRawReader.h" | |
42 | #include "AliCaloRawStreamV3.h" | |
43 | #include "AliDAQ.h" | |
44 | ||
45 | #include "AliEMCALRecParam.h" | |
46 | #include "AliEMCALLoader.h" | |
47 | #include "AliEMCALGeometry.h" | |
48 | class AliEMCALDigitizer; | |
49 | #include "AliEMCALDigit.h" | |
50 | #include "AliEMCAL.h" | |
51 | ||
52 | ClassImp(AliEMCALRawUtils) | |
53 | ||
54 | // Signal shape parameters | |
55 | Int_t AliEMCALRawUtils::fgTimeBins = 256; // number of time bins for EMCAL | |
56 | Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns | |
57 | Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec | |
58 | ||
59 | // some digitization constants | |
60 | Int_t AliEMCALRawUtils::fgThreshold = 1; | |
61 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule | |
62 | Int_t AliEMCALRawUtils::fgPedestalValue = 32; // pedestal value for digits2raw | |
63 | Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled) | |
64 | ||
65 | AliEMCALRawUtils::AliEMCALRawUtils() | |
66 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), | |
67 | fNPedSamples(0), fGeom(0), fOption("") | |
68 | { | |
69 | ||
70 | //These are default parameters. | |
71 | //Can be re-set from without with setter functions | |
72 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) | |
73 | fOrder = 2; // order of gamma fn | |
74 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
75 | fNoiseThreshold = 3; | |
76 | fNPedSamples = 5; | |
77 | ||
78 | //Get Mapping RCU files from the AliEMCALRecParam | |
79 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
80 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
81 | ||
82 | for(Int_t i = 0; i < 4; i++) { | |
83 | fMapping[i] = (AliAltroMapping*)maps->At(i); | |
84 | } | |
85 | ||
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 | |
88 | AliRunLoader *rl = AliRunLoader::Instance(); | |
89 | if(!rl) AliError("Cannot find RunLoader!"); | |
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")); | |
94 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); | |
95 | } | |
96 | ||
97 | if(!fGeom) AliFatal(Form("Could not get geometry!")); | |
98 | ||
99 | } | |
100 | ||
101 | //____________________________________________________________________________ | |
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 | |
118 | fNoiseThreshold = 3; | |
119 | fNPedSamples = 5; | |
120 | ||
121 | //Get Mapping RCU files from the AliEMCALRecParam | |
122 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
123 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
124 | ||
125 | for(Int_t i = 0; i < 4; i++) { | |
126 | fMapping[i] = (AliAltroMapping*)maps->At(i); | |
127 | } | |
128 | ||
129 | if(!fGeom) AliFatal(Form("Could not get geometry!")); | |
130 | ||
131 | } | |
132 | ||
133 | //____________________________________________________________________________ | |
134 | AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) | |
135 | : TObject(), | |
136 | fHighLowGainFactor(rawU.fHighLowGainFactor), | |
137 | fOrder(rawU.fOrder), | |
138 | fTau(rawU.fTau), | |
139 | fNoiseThreshold(rawU.fNoiseThreshold), | |
140 | fNPedSamples(rawU.fNPedSamples), | |
141 | fGeom(rawU.fGeom), | |
142 | fOption(rawU.fOption) | |
143 | { | |
144 | //copy ctor | |
145 | fMapping[0] = rawU.fMapping[0]; | |
146 | fMapping[1] = rawU.fMapping[1]; | |
147 | fMapping[2] = rawU.fMapping[2]; | |
148 | fMapping[3] = rawU.fMapping[3]; | |
149 | } | |
150 | ||
151 | //____________________________________________________________________________ | |
152 | AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) | |
153 | { | |
154 | //assignment operator | |
155 | ||
156 | if(this != &rawU) { | |
157 | fHighLowGainFactor = rawU.fHighLowGainFactor; | |
158 | fOrder = rawU.fOrder; | |
159 | fTau = rawU.fTau; | |
160 | fNoiseThreshold = rawU.fNoiseThreshold; | |
161 | fNPedSamples = rawU.fNPedSamples; | |
162 | fGeom = rawU.fGeom; | |
163 | fOption = rawU.fOption; | |
164 | fMapping[0] = rawU.fMapping[0]; | |
165 | fMapping[1] = rawU.fMapping[1]; | |
166 | fMapping[2] = rawU.fMapping[2]; | |
167 | fMapping[3] = rawU.fMapping[3]; | |
168 | } | |
169 | ||
170 | return *this; | |
171 | ||
172 | } | |
173 | ||
174 | //____________________________________________________________________________ | |
175 | AliEMCALRawUtils::~AliEMCALRawUtils() { | |
176 | //dtor | |
177 | ||
178 | } | |
179 | ||
180 | //____________________________________________________________________________ | |
181 | void AliEMCALRawUtils::Digits2Raw() | |
182 | { | |
183 | // convert digits of the current event to raw data | |
184 | ||
185 | AliRunLoader *rl = AliRunLoader::Instance(); | |
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 | } | |
197 | ||
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 | ||
203 | TArrayI adcValuesLow(fgTimeBins); | |
204 | TArrayI adcValuesHigh(fgTimeBins); | |
205 | ||
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; | |
219 | fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); | |
220 | fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; | |
221 | ||
222 | //Check which is the RCU, 0 or 1, of the cell. | |
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 | |
232 | ||
233 | if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same | |
234 | ||
235 | if (iRCU<0) | |
236 | Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); | |
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); | |
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]); | |
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 { | |
266 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh.GetArray(), adcValuesLow.GetArray()) ; | |
267 | if (lowgain) | |
268 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow.GetArray(), fgThreshold); | |
269 | else | |
270 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh.GetArray(), fgThreshold); | |
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 | } | |
282 | ||
283 | loader->UnloadDigits(); | |
284 | } | |
285 | ||
286 | //____________________________________________________________________________ | |
287 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr) | |
288 | { | |
289 | // convert raw data of the current event to digits | |
290 | ||
291 | digitsArr->Clear(); | |
292 | ||
293 | if (!digitsArr) { | |
294 | Error("Raw2Digits", "no digits found !"); | |
295 | return; | |
296 | } | |
297 | if (!reader) { | |
298 | Error("Raw2Digits", "no raw reader found !"); | |
299 | return; | |
300 | } | |
301 | ||
302 | AliCaloRawStreamV3 in(reader,"EMCAL",fMapping); | |
303 | // Select EMCAL DDL's; | |
304 | reader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL's | |
305 | ||
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); | |
312 | signalF->SetParameters(10.,0.,fTau,fOrder,5.); //set all defaults once, just to be safe | |
313 | signalF->SetParNames("amp","t0","tau","N","ped"); | |
314 | signalF->SetParameter(2,fTau); // tau in units of time bin | |
315 | signalF->SetParLimits(2,2,-1); | |
316 | signalF->SetParameter(3,fOrder); // order | |
317 | signalF->SetParLimits(3,2,-1); | |
318 | ||
319 | Int_t id = -1; | |
320 | Float_t time = 0. ; | |
321 | Float_t amp = 0. ; | |
322 | Int_t i = 0; | |
323 | Int_t startBin = 0; | |
324 | ||
325 | //Graph to hold data we will fit (should be converted to an array | |
326 | //later to speed up processing | |
327 | TGraph * gSig = new TGraph(GetRawFormatTimeBins()); | |
328 | ||
329 | Int_t lowGain = 0; | |
330 | Int_t caloFlag = 0; // low, high gain, or TRU, or LED ref. | |
331 | ||
332 | // start loop over input stream | |
333 | while (in.NextDDL()) { | |
334 | while (in.NextChannel()) { | |
335 | ||
336 | //Check if the signal is high or low gain and then do the fit, | |
337 | //if it is from TRU do not fit | |
338 | caloFlag = in.GetCaloFlag(); | |
339 | if (caloFlag != 0 && caloFlag != 1) continue; | |
340 | ||
341 | // There can be zero-suppression in the raw data, | |
342 | // so set up the TGraph in advance | |
343 | for (i=0; i < GetRawFormatTimeBins(); i++) { | |
344 | gSig->SetPoint(i, i , 0); | |
345 | } | |
346 | ||
347 | Int_t maxTime = 0; | |
348 | Int_t min = 0x3ff; // init to 10-bit max | |
349 | Int_t max = 0; // init to 10-bit min | |
350 | int nsamples = 0; | |
351 | while (in.NextBunch()) { | |
352 | const UShort_t *sig = in.GetSignals(); | |
353 | startBin = in.GetStartTimeBin(); | |
354 | ||
355 | if (((UInt_t) maxTime) < in.GetStartTimeBin()) { | |
356 | maxTime = in.GetStartTimeBin(); // timebins come in reverse order | |
357 | } | |
358 | ||
359 | if (maxTime < 0 || maxTime >= GetRawFormatTimeBins()) { | |
360 | AliWarning(Form("Invalid time bin %d",maxTime)); | |
361 | maxTime = GetRawFormatTimeBins(); | |
362 | } | |
363 | nsamples += in.GetBunchLength(); | |
364 | for (i = 0; i < in.GetBunchLength(); i++) { | |
365 | time = startBin--; | |
366 | gSig->SetPoint(time, time, (Double_t) sig[i]) ; | |
367 | if (max < sig[i]) max= sig[i]; | |
368 | if (min > sig[i]) min = sig[i]; | |
369 | } | |
370 | } // loop over bunches | |
371 | ||
372 | if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout | |
373 | ||
374 | id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; | |
375 | lowGain = in.IsLowGain(); | |
376 | ||
377 | gSig->Set(maxTime+1); | |
378 | if ( (max - min) > fNoiseThreshold) FitRaw(gSig, signalF, amp, time) ; | |
379 | ||
380 | //if (caloFlag == 0 || caloFlag == 1) { // low gain or high gain | |
381 | if (amp > 0 && amp < 2000) { //check both high and low end of | |
382 | //result, 2000 is somewhat arbitrary - not nice with magic numbers in the code.. | |
383 | AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); | |
384 | ||
385 | AddDigit(digitsArr, id, lowGain, (Int_t)amp, time); | |
386 | } | |
387 | ||
388 | //} | |
389 | ||
390 | // Reset graph | |
391 | for (Int_t index = 0; index < gSig->GetN(); index++) { | |
392 | gSig->SetPoint(index, index, 0) ; | |
393 | } | |
394 | // Reset starting parameters for fit function | |
395 | signalF->SetParameters(10.,0.,fTau,fOrder,5.); //reset all defaults just to be safe | |
396 | ||
397 | } // nsamples>0 check, some data found for this channel; not only trailer/header | |
398 | } // end while over channel | |
399 | } //end while over DDL's, of input stream | |
400 | ||
401 | delete signalF ; | |
402 | delete gSig; | |
403 | ||
404 | return ; | |
405 | } | |
406 | ||
407 | //____________________________________________________________________________ | |
408 | void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) { | |
409 | // | |
410 | // Add a new digit. | |
411 | // This routine checks whether a digit exists already for this tower | |
412 | // and then decides whether to use the high or low gain info | |
413 | // | |
414 | // Called by Raw2Digits | |
415 | ||
416 | AliEMCALDigit *digit = 0, *tmpdigit = 0; | |
417 | ||
418 | TIter nextdigit(digitsArr); | |
419 | while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { | |
420 | if (tmpdigit->GetId() == id) | |
421 | digit = tmpdigit; | |
422 | } | |
423 | ||
424 | if (!digit) { // no digit existed for this tower; create one | |
425 | if (lowGain) | |
426 | amp = Int_t(fHighLowGainFactor * amp); | |
427 | Int_t idigit = digitsArr->GetEntries(); | |
428 | new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ; | |
429 | } | |
430 | else { // a digit already exists, check range | |
431 | // (use high gain if signal < cut value, otherwise low gain) | |
432 | if (lowGain) { // new digit is low gain | |
433 | if (digit->GetAmp() > fgkOverflowCut) { // use if stored digit is out of range | |
434 | digit->SetAmp(Int_t(fHighLowGainFactor * amp)); | |
435 | digit->SetTime(time); | |
436 | } | |
437 | } | |
438 | else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range | |
439 | digit->SetAmp(amp); | |
440 | digit->SetTime(time); | |
441 | } | |
442 | } | |
443 | } | |
444 | ||
445 | //____________________________________________________________________________ | |
446 | void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time) const | |
447 | { | |
448 | // Fits the raw signal time distribution; from AliEMCALGetter | |
449 | amp = time = 0. ; | |
450 | Double_t ped = 0; | |
451 | Int_t nPed = 0; | |
452 | ||
453 | for (Int_t index = 0; index < fNPedSamples; index++) { | |
454 | Double_t ttime, signal; | |
455 | gSig->GetPoint(index, ttime, signal) ; | |
456 | if (signal > 0) { | |
457 | ped += signal; | |
458 | nPed++; | |
459 | } | |
460 | } | |
461 | ||
462 | if (nPed > 0) | |
463 | ped /= nPed; | |
464 | else { | |
465 | AliWarning("Could not determine pedestal"); | |
466 | ped = 10; // put some small value as first guess | |
467 | } | |
468 | ||
469 | ||
470 | Int_t maxFound = 0; | |
471 | Int_t iMax = 0; | |
472 | Float_t max = -1; | |
473 | Float_t maxFit = gSig->GetN(); | |
474 | Float_t minAfterSig = 9999; | |
475 | Int_t tminAfterSig = gSig->GetN(); | |
476 | Int_t nPedAfterSig = 0; | |
477 | Int_t plateauWidth = 0; | |
478 | Int_t plateauStart = 9999; | |
479 | Float_t cut = 0.3; | |
480 | ||
481 | for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) { | |
482 | Double_t ttime, signal; | |
483 | gSig->GetPoint(i, ttime, signal) ; | |
484 | if (!maxFound && signal > max) { | |
485 | iMax = i; | |
486 | max = signal; | |
487 | } | |
488 | else if ( max > ped + fNoiseThreshold ) { | |
489 | maxFound = 1; | |
490 | minAfterSig = signal; | |
491 | tminAfterSig = i; | |
492 | } | |
493 | if (maxFound) { | |
494 | if ( signal < minAfterSig) { | |
495 | minAfterSig = signal; | |
496 | tminAfterSig = i; | |
497 | } | |
498 | if (i > tminAfterSig + 5) { // Two close peaks; end fit at minimum | |
499 | maxFit = tminAfterSig; | |
500 | break; | |
501 | } | |
502 | if ( signal < cut*max){ //stop fit at 30% amplitude(avoid the pulse shape falling edge) | |
503 | maxFit = i; | |
504 | break; | |
505 | } | |
506 | if ( signal < ped + fNoiseThreshold) | |
507 | nPedAfterSig++; | |
508 | if (nPedAfterSig >= 5) { // include 5 pedestal bins after peak | |
509 | maxFit = i; | |
510 | break; | |
511 | } | |
512 | } | |
513 | //Add check on plateau | |
514 | if (signal >= fgkRawSignalOverflow - fNoiseThreshold) { | |
515 | if(plateauWidth == 0) plateauStart = i; | |
516 | plateauWidth++; | |
517 | } | |
518 | } | |
519 | ||
520 | if(plateauWidth > 0) { | |
521 | for(int j = 0; j < plateauWidth; j++) { | |
522 | //Note, have to remove the same point N times because after each | |
523 | //remove, the positions of all subsequent points have shifted down | |
524 | gSig->RemovePoint(plateauStart); | |
525 | } | |
526 | } | |
527 | ||
528 | if ( max - ped > fNoiseThreshold ) { // else its noise | |
529 | AliDebug(2,Form("Fitting max %d ped %d", max, ped)); | |
530 | signalF->SetRange(0,maxFit); | |
531 | ||
532 | if(max-ped > 50) | |
533 | signalF->SetParLimits(2,1,3); | |
534 | ||
535 | signalF->SetParameter(4, ped) ; | |
536 | signalF->SetParameter(1, iMax); | |
537 | signalF->SetParameter(0, max); | |
538 | ||
539 | gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points | |
540 | amp = signalF->GetParameter(0); | |
541 | time = signalF->GetParameter(1)*GetRawFormatTimeBinWidth() - fgTimeTrigger; | |
542 | } | |
543 | return; | |
544 | } | |
545 | //__________________________________________________________________ | |
546 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) | |
547 | { | |
548 | // Matches version used in 2007 beam test | |
549 | // | |
550 | // Shape of the electronics raw reponse: | |
551 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
552 | // | |
553 | // t' = (t - t0 + tau) / tau | |
554 | // F = A * t**N * exp( N * ( 1 - t) ) for t >= 0 | |
555 | // F = 0 for t < 0 | |
556 | // | |
557 | // parameters: | |
558 | // A: par[0] // Amplitude = peak value | |
559 | // t0: par[1] | |
560 | // tau: par[2] | |
561 | // N: par[3] | |
562 | // ped: par[4] | |
563 | // | |
564 | Double_t signal ; | |
565 | Double_t tau =par[2]; | |
566 | Double_t n =par[3]; | |
567 | Double_t ped = par[4]; | |
568 | Double_t xx = ( x[0] - par[1] + tau ) / tau ; | |
569 | ||
570 | if (xx <= 0) | |
571 | signal = ped ; | |
572 | else { | |
573 | signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; | |
574 | } | |
575 | return signal ; | |
576 | } | |
577 | ||
578 | //__________________________________________________________________ | |
579 | Bool_t AliEMCALRawUtils::RawSampledResponse( | |
580 | const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const | |
581 | { | |
582 | // for a start time dtime and an amplitude damp given by digit, | |
583 | // calculates the raw sampled response AliEMCAL::RawResponseFunction | |
584 | ||
585 | Bool_t lowGain = kFALSE ; | |
586 | ||
587 | // A: par[0] // Amplitude = peak value | |
588 | // t0: par[1] | |
589 | // tau: par[2] | |
590 | // N: par[3] | |
591 | // ped: par[4] | |
592 | ||
593 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); | |
594 | signalF.SetParameter(0, damp) ; | |
595 | signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ; | |
596 | signalF.SetParameter(2, fTau) ; | |
597 | signalF.SetParameter(3, fOrder); | |
598 | signalF.SetParameter(4, fgPedestalValue); | |
599 | ||
600 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { | |
601 | Double_t signal = signalF.Eval(iTime) ; | |
602 | ||
603 | //According to Terry Awes, 13-Apr-2008 | |
604 | //add gaussian noise in quadrature to each sample | |
605 | //Double_t noise = gRandom->Gaus(0.,fgFEENoise); | |
606 | //signal = sqrt(signal*signal + noise*noise); | |
607 | ||
608 | // March 17,09 for fast fit simulations by Alexei Pavlinov. | |
609 | // Get from PHOS analysis. In some sense it is open questions. | |
610 | Double_t noise = gRandom->Gaus(0.,fgFEENoise); | |
611 | signal += noise; | |
612 | ||
613 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
614 | if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits | |
615 | adcH[iTime] = fgkRawSignalOverflow ; | |
616 | lowGain = kTRUE ; | |
617 | } | |
618 | ||
619 | signal /= fHighLowGainFactor; | |
620 | ||
621 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
622 | if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits | |
623 | adcL[iTime] = fgkRawSignalOverflow ; | |
624 | } | |
625 | return lowGain ; | |
626 | } |