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$ */ |
17 | /* History of cvs commits: |
18 | * |
c47157cd |
19 | * $Log$ |
48a56166 |
20 | * Revision 1.10 2007/12/06 13:58:11 hristov |
21 | * Additional pritection. Do not delete the mapping, it is owned by another class |
22 | * |
e36e3bcf |
23 | * Revision 1.9 2007/12/06 02:19:51 jklay |
24 | * incorporated fitting procedure from testbeam analysis into AliRoot |
25 | * |
8cb998bd |
26 | * Revision 1.8 2007/12/05 02:30:51 jklay |
27 | * modification to read Altro mappings into AliEMCALRecParam and pass to AliEMCALRawUtils from AliEMCALReconstructor; add option to AliEMCALRawUtils to set old RCU format (for testbeam) or not |
28 | * |
feedcab9 |
29 | * Revision 1.7 2007/11/14 15:51:46 gustavo |
30 | * Take out few unnecessary prints |
31 | * |
6d0b6861 |
32 | * Revision 1.6 2007/11/01 01:23:51 mvl |
33 | * Removed call to SetOldRCUFormat, which is only needed for testbeam data |
34 | * |
f31c5945 |
35 | * Revision 1.5 2007/11/01 01:20:33 mvl |
36 | * Further improvement of peak finding; more robust fit |
37 | * |
e9dbb64a |
38 | * Revision 1.4 2007/10/31 17:15:24 mvl |
39 | * Fixed bug in raw data unpacking; Added pedestal to signal fit; Added logic to deal with high/low gain |
40 | * |
82cbdfca |
41 | * Revision 1.3 2007/09/27 08:36:46 mvl |
42 | * More robust setting of fit range in FitRawSignal (P. Hristov) |
43 | * |
5a056daa |
44 | * Revision 1.2 2007/09/03 20:55:35 jklay |
45 | * EMCAL e-by-e reconstruction methods from Cvetan |
46 | * |
c47157cd |
47 | * Revision 1.1 2007/03/17 19:56:38 mvl |
48 | * Moved signal shape routines from AliEMCAL to separate class AliEMCALRawUtils to streamline raw data reconstruction code. |
49 | * */ |
ee299369 |
50 | |
51 | //*-- Author: Marco van Leeuwen (LBL) |
52 | #include "AliEMCALRawUtils.h" |
53 | |
54 | #include "TF1.h" |
55 | #include "TGraph.h" |
56 | #include "TSystem.h" |
57 | |
58 | #include "AliLog.h" |
72c58de0 |
59 | #include "AliRun.h" |
ee299369 |
60 | #include "AliRunLoader.h" |
61 | #include "AliCaloAltroMapping.h" |
62 | #include "AliAltroBuffer.h" |
63 | #include "AliRawReader.h" |
64 | #include "AliCaloRawStream.h" |
65 | #include "AliDAQ.h" |
66 | |
feedcab9 |
67 | #include "AliEMCALRecParam.h" |
ee299369 |
68 | #include "AliEMCALLoader.h" |
69 | #include "AliEMCALGeometry.h" |
70 | #include "AliEMCALDigitizer.h" |
71 | #include "AliEMCALDigit.h" |
72 | |
73 | |
74 | ClassImp(AliEMCALRawUtils) |
75 | |
76 | // Signal shape parameters |
82cbdfca |
77 | Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns |
ee299369 |
78 | Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec |
79 | |
80 | // some digitization constants |
81 | Int_t AliEMCALRawUtils::fgThreshold = 1; |
82 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule |
83 | |
8cb998bd |
84 | AliEMCALRawUtils::AliEMCALRawUtils() |
b4133f05 |
85 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), |
86 | fNPedSamples(0), fGeom(0), fOption("") |
8cb998bd |
87 | { |
b4133f05 |
88 | |
89 | //These are default parameters. |
90 | //Can be re-set from without with setter functions |
ee299369 |
91 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) |
b4133f05 |
92 | fOrder = 2; // order of gamma fn |
93 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam |
94 | fNoiseThreshold = 3; |
95 | fNPedSamples = 5; |
65bdc82f |
96 | |
97 | //Get Mapping RCU files from the AliEMCALRecParam |
98 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); |
99 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); |
100 | |
101 | for(Int_t i = 0; i < 2; i++) { |
102 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
103 | } |
104 | |
72c58de0 |
105 | //To make sure we match with the geometry in a simulation file, |
106 | //let's try to get it first. If not, take the default geometry |
107 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
108 | if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) { |
109 | fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); |
110 | } else { |
111 | AliInfo(Form("Using default geometry in raw reco")); |
937d0661 |
112 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); |
65bdc82f |
113 | } |
114 | |
72c58de0 |
115 | if(!fGeom) AliFatal(Form("Could not get geometry!")); |
116 | |
65bdc82f |
117 | } |
118 | |
119 | //____________________________________________________________________________ |
120 | AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) |
121 | : TObject(), |
122 | fHighLowGainFactor(rawU.fHighLowGainFactor), |
b4133f05 |
123 | fOrder(rawU.fOrder), |
124 | fTau(rawU.fTau), |
125 | fNoiseThreshold(rawU.fNoiseThreshold), |
126 | fNPedSamples(rawU.fNPedSamples), |
65bdc82f |
127 | fGeom(rawU.fGeom), |
128 | fOption(rawU.fOption) |
129 | { |
130 | //copy ctor |
131 | fMapping[0] = rawU.fMapping[0]; |
132 | fMapping[1] = rawU.fMapping[1]; |
133 | } |
134 | |
135 | //____________________________________________________________________________ |
136 | AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) |
137 | { |
138 | //assignment operator |
139 | |
140 | if(this != &rawU) { |
141 | fHighLowGainFactor = rawU.fHighLowGainFactor; |
b4133f05 |
142 | fOrder = rawU.fOrder; |
143 | fTau = rawU.fTau; |
144 | fNoiseThreshold = rawU.fNoiseThreshold; |
145 | fNPedSamples = rawU.fNPedSamples; |
65bdc82f |
146 | fGeom = rawU.fGeom; |
147 | fOption = rawU.fOption; |
148 | fMapping[0] = rawU.fMapping[0]; |
149 | fMapping[1] = rawU.fMapping[1]; |
150 | } |
151 | |
152 | return *this; |
153 | |
ee299369 |
154 | } |
65bdc82f |
155 | |
ee299369 |
156 | //____________________________________________________________________________ |
157 | AliEMCALRawUtils::~AliEMCALRawUtils() { |
65bdc82f |
158 | |
ee299369 |
159 | } |
65bdc82f |
160 | |
ee299369 |
161 | //____________________________________________________________________________ |
65bdc82f |
162 | void AliEMCALRawUtils::Digits2Raw() |
ee299369 |
163 | { |
164 | // convert digits of the current event to raw data |
165 | |
166 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
167 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); |
168 | |
169 | // get the digits |
170 | loader->LoadDigits("EMCAL"); |
171 | loader->GetEvent(); |
172 | TClonesArray* digits = loader->Digits() ; |
173 | |
174 | if (!digits) { |
175 | Warning("Digits2Raw", "no digits found !"); |
176 | return; |
177 | } |
65bdc82f |
178 | |
ee299369 |
179 | static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here |
180 | AliAltroBuffer* buffers[nDDL]; |
181 | for (Int_t i=0; i < nDDL; i++) |
182 | buffers[i] = 0; |
183 | |
184 | Int_t adcValuesLow[fgkTimeBins]; |
185 | Int_t adcValuesHigh[fgkTimeBins]; |
186 | |
ee299369 |
187 | // loop over digits (assume ordered digits) |
188 | for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { |
189 | AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ; |
190 | if (digit->GetAmp() < fgThreshold) |
191 | continue; |
192 | |
193 | //get cell indices |
194 | Int_t nSM = 0; |
195 | Int_t nIphi = 0; |
196 | Int_t nIeta = 0; |
197 | Int_t iphi = 0; |
198 | Int_t ieta = 0; |
199 | Int_t nModule = 0; |
65bdc82f |
200 | fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); |
201 | fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; |
ee299369 |
202 | |
203 | //Check which is the RCU of the cell. |
204 | Int_t iRCU = -111; |
205 | //RCU0 |
206 | if (0<=iphi&&iphi<8) iRCU=0; // first cable row |
207 | else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; |
208 | //second cable row |
209 | //RCU1 |
210 | else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; |
211 | //second cable row |
212 | else if(16<=iphi&&iphi<24) iRCU=1; // third cable row |
e36e3bcf |
213 | if (iRCU<0) |
214 | Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); |
ee299369 |
215 | |
216 | //Which DDL? |
217 | Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; |
218 | if (iDDL >= nDDL) |
219 | Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); |
220 | |
221 | if (buffers[iDDL] == 0) { |
222 | // open new file and write dummy header |
223 | TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); |
65bdc82f |
224 | buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCU]); |
ee299369 |
225 | buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; |
226 | } |
227 | |
228 | // out of time range signal (?) |
229 | if (digit->GetTimeR() > GetRawFormatTimeMax() ) { |
230 | AliInfo("Signal is out of time range.\n"); |
231 | buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp()); |
232 | buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin |
233 | buffers[iDDL]->FillBuffer(3); // bunch length |
234 | buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer |
235 | // calculate the time response function |
236 | } else { |
237 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; |
238 | if (lowgain) |
239 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold); |
240 | else |
241 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold); |
242 | } |
243 | } |
244 | |
245 | // write headers and close files |
246 | for (Int_t i=0; i < nDDL; i++) { |
247 | if (buffers[i]) { |
248 | buffers[i]->Flush(); |
249 | buffers[i]->WriteDataHeader(kFALSE, kFALSE); |
250 | delete buffers[i]; |
251 | } |
252 | } |
65bdc82f |
253 | |
ee299369 |
254 | loader->UnloadDigits(); |
255 | } |
256 | |
257 | //____________________________________________________________________________ |
65bdc82f |
258 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr) |
ee299369 |
259 | { |
65bdc82f |
260 | // convert raw data of the current event to digits |
ee299369 |
261 | |
c47157cd |
262 | digitsArr->Clear(); |
ee299369 |
263 | |
c47157cd |
264 | if (!digitsArr) { |
ee299369 |
265 | Error("Raw2Digits", "no digits found !"); |
266 | return; |
267 | } |
268 | if (!reader) { |
269 | Error("Raw2Digits", "no raw reader found !"); |
270 | return; |
271 | } |
272 | |
65bdc82f |
273 | AliCaloRawStream in(reader,"EMCAL",fMapping); |
ee299369 |
274 | // Select EMCAL DDL's; |
275 | reader->Select("EMCAL"); |
feedcab9 |
276 | |
277 | TString option = GetOption(); |
278 | if (option.Contains("OldRCUFormat")) |
279 | in.SetOldRCUFormat(kTRUE); // Needed for testbeam data |
280 | else |
281 | in.SetOldRCUFormat(kFALSE); |
282 | |
283 | |
8cb998bd |
284 | //Updated fitting routine from 2007 beam test takes into account |
285 | //possibility of two peaks in data and selects first one for fitting |
286 | //Also sets some of the starting parameters based on the shape of the |
287 | //given raw signal being fit |
288 | |
289 | TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); |
b4133f05 |
290 | signalF->SetParameters(10.,0.,fTau,fOrder,5.); //set all defaults once, just to be safe |
8cb998bd |
291 | signalF->SetParNames("amp","t0","tau","N","ped"); |
b4133f05 |
292 | signalF->SetParameter(2,fTau); // tau in units of time bin |
8cb998bd |
293 | signalF->SetParLimits(2,2,-1); |
b4133f05 |
294 | signalF->SetParameter(3,fOrder); // order |
8cb998bd |
295 | signalF->SetParLimits(3,2,-1); |
48a56166 |
296 | |
ee299369 |
297 | Int_t id = -1; |
82cbdfca |
298 | Float_t time = 0. ; |
299 | Float_t amp = 0. ; |
ee299369 |
300 | |
8cb998bd |
301 | //Graph to hold data we will fit (should be converted to an array |
302 | //later to speed up processing |
303 | TGraph * gSig = new TGraph(GetRawFormatTimeBins()); |
ee299369 |
304 | |
82cbdfca |
305 | Int_t readOk = 1; |
ee299369 |
306 | Int_t lowGain = 0; |
307 | |
82cbdfca |
308 | while (readOk && in.GetModule() < 0) |
309 | readOk = in.Next(); // Go to first digit |
e9dbb64a |
310 | |
8cb998bd |
311 | Int_t col = 0; |
312 | Int_t row = 0; |
313 | |
82cbdfca |
314 | while (readOk) { |
65bdc82f |
315 | id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; |
ee299369 |
316 | lowGain = in.IsLowGain(); |
82cbdfca |
317 | Int_t maxTime = in.GetTime(); // timebins come in reverse order |
e9dbb64a |
318 | if (maxTime < 0 || maxTime >= GetRawFormatTimeBins()) { |
319 | AliWarning(Form("Invalid time bin %d",maxTime)); |
320 | maxTime = GetRawFormatTimeBins(); |
321 | } |
82cbdfca |
322 | gSig->Set(maxTime+1); |
323 | // There is some kind of zero-suppression in the raw data, |
324 | // so set up the TGraph in advance |
325 | for (Int_t i=0; i < maxTime; i++) { |
8cb998bd |
326 | gSig->SetPoint(i, i , 0); |
82cbdfca |
327 | } |
ee299369 |
328 | |
82cbdfca |
329 | Int_t iTime = 0; |
ee299369 |
330 | do { |
82cbdfca |
331 | if (in.GetTime() >= gSig->GetN()) { |
332 | AliWarning("Too many time bins"); |
333 | gSig->Set(in.GetTime()); |
ee299369 |
334 | } |
8cb998bd |
335 | col = in.GetColumn(); |
336 | row = in.GetRow(); |
337 | |
338 | gSig->SetPoint(in.GetTime(), in.GetTime(), in.GetSignal()) ; |
82cbdfca |
339 | if (in.GetTime() > maxTime) |
340 | maxTime = in.GetTime(); |
ee299369 |
341 | iTime++; |
82cbdfca |
342 | } while ((readOk = in.Next()) && !in.IsNewHWAddress()); |
ee299369 |
343 | |
344 | FitRaw(gSig, signalF, amp, time) ; |
ee299369 |
345 | |
3ced962c |
346 | if (amp > 0 && amp < 10000) { //check both high and low end of |
347 | //result, 10000 is somewhat arbitrary |
82cbdfca |
348 | AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); |
8cb998bd |
349 | //cout << "col " << col-40 << " row " << row-8 << " lowGain " << lowGain << " amp " << amp << endl; |
82cbdfca |
350 | AddDigit(digitsArr, id, lowGain, (Int_t)amp, time); |
ee299369 |
351 | } |
ee299369 |
352 | |
353 | // Reset graph |
82cbdfca |
354 | for (Int_t index = 0; index < gSig->GetN(); index++) { |
8cb998bd |
355 | gSig->SetPoint(index, index, 0) ; |
ee299369 |
356 | } |
48a56166 |
357 | // Reset starting parameters for fit function |
b4133f05 |
358 | signalF->SetParameters(10.,0.,fTau,fOrder,5.); //reset all defaults just to be safe |
48a56166 |
359 | |
82cbdfca |
360 | }; // EMCAL entries loop |
ee299369 |
361 | |
362 | delete signalF ; |
363 | delete gSig; |
364 | |
365 | return ; |
366 | } |
367 | |
368 | //____________________________________________________________________________ |
82cbdfca |
369 | void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) { |
370 | // |
371 | // Add a new digit. |
372 | // This routine checks whether a digit exists already for this tower |
373 | // and then decides whether to use the high or low gain info |
374 | // |
375 | // Called by Raw2Digits |
376 | |
377 | AliEMCALDigit *digit = 0, *tmpdigit = 0; |
378 | |
379 | TIter nextdigit(digitsArr); |
380 | while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { |
381 | if (tmpdigit->GetId() == id) |
382 | digit = tmpdigit; |
383 | } |
384 | |
385 | if (!digit) { // no digit existed for this tower; create one |
386 | if (lowGain) |
387 | amp = Int_t(fHighLowGainFactor * amp); |
388 | Int_t idigit = digitsArr->GetEntries(); |
389 | new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ; |
390 | } |
391 | else { // a digit already exists, check range |
b4133f05 |
392 | // (use high gain if signal < cut value, otherwise low gain) |
82cbdfca |
393 | if (lowGain) { // new digit is low gain |
b4133f05 |
394 | if (digit->GetAmp() > fgkOverflowCut) { // use if stored digit is out of range |
82cbdfca |
395 | digit->SetAmp(Int_t(fHighLowGainFactor * amp)); |
396 | digit->SetTime(time); |
397 | } |
398 | } |
b4133f05 |
399 | else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range |
82cbdfca |
400 | digit->SetAmp(amp); |
401 | digit->SetTime(time); |
402 | } |
403 | } |
404 | } |
405 | |
406 | //____________________________________________________________________________ |
407 | void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time) |
ee299369 |
408 | { |
409 | // Fits the raw signal time distribution; from AliEMCALGetter |
410 | |
ee299369 |
411 | amp = time = 0. ; |
82cbdfca |
412 | Double_t ped = 0; |
413 | Int_t nPed = 0; |
ee299369 |
414 | |
b4133f05 |
415 | for (Int_t index = 0; index < fNPedSamples; index++) { |
e9dbb64a |
416 | Double_t ttime, signal; |
82cbdfca |
417 | gSig->GetPoint(index, ttime, signal) ; |
418 | if (signal > 0) { |
419 | ped += signal; |
420 | nPed++; |
ee299369 |
421 | } |
422 | } |
e9dbb64a |
423 | |
82cbdfca |
424 | if (nPed > 0) |
425 | ped /= nPed; |
e9dbb64a |
426 | else { |
8cb998bd |
427 | AliWarning("Could not determine pedestal"); |
82cbdfca |
428 | ped = 10; // put some small value as first guess |
e9dbb64a |
429 | } |
82cbdfca |
430 | |
e9dbb64a |
431 | Int_t max_found = 0; |
432 | Int_t i_max = 0; |
433 | Float_t max = -1; |
8cb998bd |
434 | Float_t max_fit = gSig->GetN(); |
e9dbb64a |
435 | Float_t min_after_sig = 9999; |
8cb998bd |
436 | Int_t tmin_after_sig = gSig->GetN(); |
e9dbb64a |
437 | Int_t n_ped_after_sig = 0; |
438 | |
b4133f05 |
439 | for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) { |
e9dbb64a |
440 | Double_t ttime, signal; |
441 | gSig->GetPoint(i, ttime, signal) ; |
442 | if (!max_found && signal > max) { |
443 | i_max = i; |
e9dbb64a |
444 | max = signal; |
445 | } |
b4133f05 |
446 | else if ( max > ped + fNoiseThreshold ) { |
e9dbb64a |
447 | max_found = 1; |
448 | min_after_sig = signal; |
8cb998bd |
449 | tmin_after_sig = i; |
e9dbb64a |
450 | } |
451 | if (max_found) { |
452 | if ( signal < min_after_sig) { |
453 | min_after_sig = signal; |
8cb998bd |
454 | tmin_after_sig = i; |
e9dbb64a |
455 | } |
456 | if (i > tmin_after_sig + 5) { // Two close peaks; end fit at minimum |
457 | max_fit = tmin_after_sig; |
458 | break; |
459 | } |
b4133f05 |
460 | if ( signal < ped + fNoiseThreshold) |
e9dbb64a |
461 | n_ped_after_sig++; |
462 | if (n_ped_after_sig >= 5) { // include 5 pedestal bins after peak |
8cb998bd |
463 | max_fit = i; |
e9dbb64a |
464 | break; |
465 | } |
5a056daa |
466 | } |
82cbdfca |
467 | } |
5a056daa |
468 | |
b4133f05 |
469 | if ( max - ped > fNoiseThreshold ) { // else its noise |
e9dbb64a |
470 | AliDebug(2,Form("Fitting max %d ped %d", max, ped)); |
8cb998bd |
471 | signalF->SetRange(0,max_fit); |
472 | |
473 | if(max-ped > 50) |
474 | signalF->SetParLimits(2,1,3); |
475 | |
476 | signalF->SetParameter(4, ped) ; |
477 | signalF->SetParameter(1, i_max); |
478 | signalF->SetParameter(0, max); |
479 | |
480 | gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points |
481 | amp = signalF->GetParameter(0); |
482 | time = signalF->GetParameter(1)*GetRawFormatTimeBinWidth() - fgTimeTrigger; |
ee299369 |
483 | } |
484 | return; |
485 | } |
486 | //__________________________________________________________________ |
487 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) |
488 | { |
8cb998bd |
489 | // Matches version used in 2007 beam test |
490 | // |
ee299369 |
491 | // Shape of the electronics raw reponse: |
492 | // It is a semi-gaussian, 2nd order Gamma function of the general form |
493 | // |
494 | // t' = (t - t0 + tau) / tau |
495 | // F = A * t**N * exp( N * ( 1 - t) ) for t >= 0 |
496 | // F = 0 for t < 0 |
497 | // |
498 | // parameters: |
8cb998bd |
499 | // A: par[0] // Amplitude = peak value |
500 | // t0: par[1] |
501 | // tau: par[2] |
502 | // N: par[3] |
503 | // ped: par[4] |
ee299369 |
504 | // |
505 | Double_t signal ; |
8cb998bd |
506 | Double_t tau =par[2]; |
507 | Double_t N =par[3]; |
508 | Double_t ped = par[4]; |
509 | Double_t xx = ( x[0] - par[1] + tau ) / tau ; |
ee299369 |
510 | |
5a056daa |
511 | if (xx <= 0) |
8cb998bd |
512 | signal = ped ; |
ee299369 |
513 | else { |
8cb998bd |
514 | signal = ped + par[0] * TMath::Power(xx , N) * TMath::Exp(N * (1 - xx )) ; |
ee299369 |
515 | } |
516 | return signal ; |
517 | } |
518 | |
519 | //__________________________________________________________________ |
520 | Bool_t AliEMCALRawUtils::RawSampledResponse( |
521 | const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const |
522 | { |
523 | // for a start time dtime and an amplitude damp given by digit, |
524 | // calculates the raw sampled response AliEMCAL::RawResponseFunction |
525 | |
82cbdfca |
526 | const Int_t pedVal = 32; |
ee299369 |
527 | Bool_t lowGain = kFALSE ; |
528 | |
48a56166 |
529 | // A: par[0] // Amplitude = peak value |
530 | // t0: par[1] |
531 | // tau: par[2] |
532 | // N: par[3] |
533 | // ped: par[4] |
534 | |
535 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 5); |
536 | signalF.SetParameter(0, damp) ; |
537 | signalF.SetParameter(1, dtime + fgTimeTrigger) ; |
b4133f05 |
538 | signalF.SetParameter(2, fTau) ; |
539 | signalF.SetParameter(3, fOrder); |
48a56166 |
540 | signalF.SetParameter(4, pedVal); |
ee299369 |
541 | |
542 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { |
82cbdfca |
543 | Double_t time = iTime * GetRawFormatTimeBinWidth() ; |
ee299369 |
544 | Double_t signal = signalF.Eval(time) ; |
545 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; |
b4133f05 |
546 | if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits |
547 | adcH[iTime] = fgkRawSignalOverflow ; |
ee299369 |
548 | lowGain = kTRUE ; |
549 | } |
550 | |
551 | signal /= fHighLowGainFactor; |
552 | |
553 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; |
b4133f05 |
554 | if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits |
555 | adcL[iTime] = fgkRawSignalOverflow ; |
ee299369 |
556 | } |
557 | return lowGain ; |
558 | } |