]>
Commit | Line | Data |
---|---|---|
ee299369 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | /* History of cvs commits: | |
18 | * | |
c47157cd | 19 | * $Log$ |
6d0b6861 | 20 | * Revision 1.6 2007/11/01 01:23:51 mvl |
21 | * Removed call to SetOldRCUFormat, which is only needed for testbeam data | |
22 | * | |
f31c5945 | 23 | * Revision 1.5 2007/11/01 01:20:33 mvl |
24 | * Further improvement of peak finding; more robust fit | |
25 | * | |
e9dbb64a | 26 | * Revision 1.4 2007/10/31 17:15:24 mvl |
27 | * Fixed bug in raw data unpacking; Added pedestal to signal fit; Added logic to deal with high/low gain | |
28 | * | |
82cbdfca | 29 | * Revision 1.3 2007/09/27 08:36:46 mvl |
30 | * More robust setting of fit range in FitRawSignal (P. Hristov) | |
31 | * | |
5a056daa | 32 | * Revision 1.2 2007/09/03 20:55:35 jklay |
33 | * EMCAL e-by-e reconstruction methods from Cvetan | |
34 | * | |
c47157cd | 35 | * Revision 1.1 2007/03/17 19:56:38 mvl |
36 | * Moved signal shape routines from AliEMCAL to separate class AliEMCALRawUtils to streamline raw data reconstruction code. | |
37 | * */ | |
ee299369 | 38 | |
39 | //*-- Author: Marco van Leeuwen (LBL) | |
40 | #include "AliEMCALRawUtils.h" | |
41 | ||
42 | #include "TF1.h" | |
43 | #include "TGraph.h" | |
44 | #include "TSystem.h" | |
45 | ||
46 | #include "AliLog.h" | |
47 | #include "AliRunLoader.h" | |
48 | #include "AliCaloAltroMapping.h" | |
49 | #include "AliAltroBuffer.h" | |
50 | #include "AliRawReader.h" | |
51 | #include "AliCaloRawStream.h" | |
52 | #include "AliDAQ.h" | |
53 | ||
54 | #include "AliEMCALLoader.h" | |
55 | #include "AliEMCALGeometry.h" | |
56 | #include "AliEMCALDigitizer.h" | |
57 | #include "AliEMCALDigit.h" | |
58 | ||
59 | ||
60 | ClassImp(AliEMCALRawUtils) | |
61 | ||
62 | // Signal shape parameters | |
82cbdfca | 63 | Int_t AliEMCALRawUtils::fgOrder = 2 ; // Order of gamma function |
64 | Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns | |
65 | Double_t AliEMCALRawUtils::fgTau = 235E-9 ; // 235 ns (from CERN testbeam; not very accurate) | |
ee299369 | 66 | Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec |
67 | ||
68 | // some digitization constants | |
69 | Int_t AliEMCALRawUtils::fgThreshold = 1; | |
70 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule | |
71 | ||
72 | AliEMCALRawUtils::AliEMCALRawUtils(): fHighLowGainFactor(0.) { | |
73 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) | |
74 | } | |
75 | //____________________________________________________________________________ | |
76 | AliEMCALRawUtils::~AliEMCALRawUtils() { | |
77 | } | |
78 | //____________________________________________________________________________ | |
79 | void AliEMCALRawUtils::Digits2Raw() | |
80 | { | |
81 | // convert digits of the current event to raw data | |
82 | ||
83 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); | |
84 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); | |
85 | ||
86 | // get the digits | |
87 | loader->LoadDigits("EMCAL"); | |
88 | loader->GetEvent(); | |
89 | TClonesArray* digits = loader->Digits() ; | |
90 | ||
91 | if (!digits) { | |
92 | Warning("Digits2Raw", "no digits found !"); | |
93 | return; | |
94 | } | |
95 | ||
96 | // get the geometry | |
97 | AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(); | |
98 | if (!geom) { | |
99 | AliError(Form("No geometry found !")); | |
100 | return; | |
101 | } | |
102 | ||
103 | static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here | |
104 | AliAltroBuffer* buffers[nDDL]; | |
105 | for (Int_t i=0; i < nDDL; i++) | |
106 | buffers[i] = 0; | |
107 | ||
108 | Int_t adcValuesLow[fgkTimeBins]; | |
109 | Int_t adcValuesHigh[fgkTimeBins]; | |
110 | ||
111 | //Load Mapping RCU files once | |
112 | TString path = gSystem->Getenv("ALICE_ROOT"); | |
113 | path += "/EMCAL/mapping/RCU"; | |
114 | TString path0 = path+"0.data";//This file will change in future | |
115 | TString path1 = path+"1.data";//This file will change in future | |
116 | AliAltroMapping * mapping[2] ; // For the moment only 2 | |
117 | mapping[0] = new AliCaloAltroMapping(path0.Data()); | |
118 | mapping[1] = new AliCaloAltroMapping(path1.Data()); | |
119 | ||
120 | // loop over digits (assume ordered digits) | |
121 | for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { | |
122 | AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ; | |
123 | if (digit->GetAmp() < fgThreshold) | |
124 | continue; | |
125 | ||
126 | //get cell indices | |
127 | Int_t nSM = 0; | |
128 | Int_t nIphi = 0; | |
129 | Int_t nIeta = 0; | |
130 | Int_t iphi = 0; | |
131 | Int_t ieta = 0; | |
132 | Int_t nModule = 0; | |
133 | geom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); | |
134 | geom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; | |
135 | ||
136 | //Check which is the RCU of the cell. | |
137 | Int_t iRCU = -111; | |
138 | //RCU0 | |
139 | if (0<=iphi&&iphi<8) iRCU=0; // first cable row | |
140 | else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; | |
141 | //second cable row | |
142 | //RCU1 | |
143 | else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; | |
144 | //second cable row | |
145 | else if(16<=iphi&&iphi<24) iRCU=1; // third cable row | |
146 | ||
147 | //Which DDL? | |
148 | Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; | |
149 | if (iDDL >= nDDL) | |
150 | Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); | |
151 | ||
152 | if (buffers[iDDL] == 0) { | |
153 | // open new file and write dummy header | |
154 | TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); | |
155 | buffers[iDDL] = new AliAltroBuffer(fileName.Data(),mapping[iRCU]); | |
156 | buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; | |
157 | } | |
158 | ||
159 | // out of time range signal (?) | |
160 | if (digit->GetTimeR() > GetRawFormatTimeMax() ) { | |
161 | AliInfo("Signal is out of time range.\n"); | |
162 | buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp()); | |
163 | buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin | |
164 | buffers[iDDL]->FillBuffer(3); // bunch length | |
165 | buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer | |
166 | // calculate the time response function | |
167 | } else { | |
168 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; | |
169 | if (lowgain) | |
170 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold); | |
171 | else | |
172 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold); | |
173 | } | |
174 | } | |
175 | ||
176 | // write headers and close files | |
177 | for (Int_t i=0; i < nDDL; i++) { | |
178 | if (buffers[i]) { | |
179 | buffers[i]->Flush(); | |
180 | buffers[i]->WriteDataHeader(kFALSE, kFALSE); | |
181 | delete buffers[i]; | |
182 | } | |
183 | } | |
184 | mapping[0]->Delete(); | |
185 | mapping[1]->Delete(); | |
186 | loader->UnloadDigits(); | |
187 | } | |
188 | ||
189 | //____________________________________________________________________________ | |
c47157cd | 190 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr) |
ee299369 | 191 | { |
192 | // convert raw data of the current event to digits | |
193 | AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); | |
194 | if (!geom) { | |
195 | AliError(Form("No geometry found !")); | |
196 | return; | |
197 | } | |
198 | ||
c47157cd | 199 | digitsArr->Clear(); |
ee299369 | 200 | |
c47157cd | 201 | if (!digitsArr) { |
ee299369 | 202 | Error("Raw2Digits", "no digits found !"); |
203 | return; | |
204 | } | |
205 | if (!reader) { | |
206 | Error("Raw2Digits", "no raw reader found !"); | |
207 | return; | |
208 | } | |
209 | ||
210 | // Use AliAltroRawStream to read the ALTRO format. No need to | |
211 | // reinvent the wheel :-) | |
212 | AliCaloRawStream in(reader,"EMCAL"); | |
213 | // Select EMCAL DDL's; | |
214 | reader->Select("EMCAL"); | |
f31c5945 | 215 | //in.SetOldRCUFormat(kTRUE); // Needed for testbeam data |
e9dbb64a | 216 | |
ee299369 | 217 | // reading is from previously existing AliEMCALGetter.cxx |
218 | // ReadRaw method | |
219 | TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4); | |
ee299369 | 220 | |
221 | Int_t id = -1; | |
82cbdfca | 222 | Float_t time = 0. ; |
223 | Float_t amp = 0. ; | |
ee299369 | 224 | |
225 | TGraph * gSig = new TGraph(GetRawFormatTimeBins()) ; | |
226 | ||
82cbdfca | 227 | Int_t readOk = 1; |
ee299369 | 228 | Int_t lowGain = 0; |
229 | ||
82cbdfca | 230 | while (readOk && in.GetModule() < 0) |
231 | readOk = in.Next(); // Go to first digit | |
e9dbb64a | 232 | |
82cbdfca | 233 | while (readOk) { |
ee299369 | 234 | id = geom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; |
235 | lowGain = in.IsLowGain(); | |
82cbdfca | 236 | Int_t maxTime = in.GetTime(); // timebins come in reverse order |
e9dbb64a | 237 | if (maxTime < 0 || maxTime >= GetRawFormatTimeBins()) { |
238 | AliWarning(Form("Invalid time bin %d",maxTime)); | |
239 | maxTime = GetRawFormatTimeBins(); | |
240 | } | |
82cbdfca | 241 | gSig->Set(maxTime+1); |
242 | // There is some kind of zero-suppression in the raw data, | |
243 | // so set up the TGraph in advance | |
244 | for (Int_t i=0; i < maxTime; i++) { | |
245 | gSig->SetPoint(i, i * GetRawFormatTimeBinWidth(), 0); | |
246 | } | |
ee299369 | 247 | |
82cbdfca | 248 | Int_t iTime = 0; |
ee299369 | 249 | do { |
82cbdfca | 250 | if (in.GetTime() >= gSig->GetN()) { |
251 | AliWarning("Too many time bins"); | |
252 | gSig->Set(in.GetTime()); | |
ee299369 | 253 | } |
82cbdfca | 254 | gSig->SetPoint(in.GetTime(), |
255 | in.GetTime() * GetRawFormatTimeBinWidth(), | |
256 | in.GetSignal()) ; | |
257 | if (in.GetTime() > maxTime) | |
258 | maxTime = in.GetTime(); | |
ee299369 | 259 | iTime++; |
82cbdfca | 260 | } while ((readOk = in.Next()) && !in.IsNewHWAddress()); |
261 | signalF->SetRange(0,(Float_t)maxTime*GetRawFormatTimeBinWidth()); | |
ee299369 | 262 | |
263 | FitRaw(gSig, signalF, amp, time) ; | |
ee299369 | 264 | |
265 | if (amp > 0) { | |
82cbdfca | 266 | AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); |
267 | AddDigit(digitsArr, id, lowGain, (Int_t)amp, time); | |
ee299369 | 268 | } |
ee299369 | 269 | |
270 | // Reset graph | |
82cbdfca | 271 | for (Int_t index = 0; index < gSig->GetN(); index++) { |
272 | gSig->SetPoint(index, index * GetRawFormatTimeBinWidth(), 0) ; | |
ee299369 | 273 | } |
82cbdfca | 274 | }; // EMCAL entries loop |
ee299369 | 275 | |
276 | delete signalF ; | |
277 | delete gSig; | |
278 | ||
279 | return ; | |
280 | } | |
281 | ||
282 | //____________________________________________________________________________ | |
82cbdfca | 283 | void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) { |
284 | // | |
285 | // Add a new digit. | |
286 | // This routine checks whether a digit exists already for this tower | |
287 | // and then decides whether to use the high or low gain info | |
288 | // | |
289 | // Called by Raw2Digits | |
290 | ||
291 | AliEMCALDigit *digit = 0, *tmpdigit = 0; | |
292 | ||
293 | TIter nextdigit(digitsArr); | |
294 | while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { | |
295 | if (tmpdigit->GetId() == id) | |
296 | digit = tmpdigit; | |
297 | } | |
298 | ||
299 | if (!digit) { // no digit existed for this tower; create one | |
300 | if (lowGain) | |
301 | amp = Int_t(fHighLowGainFactor * amp); | |
302 | Int_t idigit = digitsArr->GetEntries(); | |
303 | new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ; | |
304 | } | |
305 | else { // a digit already exists, check range | |
306 | // (use high gain if signal < 800, otherwise low gain) | |
307 | if (lowGain) { // new digit is low gain | |
308 | if (digit->GetAmp() > 800) { // use if stored digit is out of range | |
309 | digit->SetAmp(Int_t(fHighLowGainFactor * amp)); | |
310 | digit->SetTime(time); | |
311 | } | |
312 | } | |
313 | else if (amp < 800) { // new digit is high gain; use if not out of range | |
314 | digit->SetAmp(amp); | |
315 | digit->SetTime(time); | |
316 | } | |
317 | } | |
318 | } | |
319 | ||
320 | //____________________________________________________________________________ | |
321 | void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time) | |
ee299369 | 322 | { |
323 | // Fits the raw signal time distribution; from AliEMCALGetter | |
324 | ||
e9dbb64a | 325 | const Int_t kNoiseThreshold = 5; |
326 | const Int_t kNPedSamples = 10; | |
ee299369 | 327 | amp = time = 0. ; |
82cbdfca | 328 | Double_t ped = 0; |
329 | Int_t nPed = 0; | |
ee299369 | 330 | |
e9dbb64a | 331 | for (Int_t index = 0; index < kNPedSamples; index++) { |
332 | Double_t ttime, signal; | |
82cbdfca | 333 | gSig->GetPoint(index, ttime, signal) ; |
334 | if (signal > 0) { | |
335 | ped += signal; | |
336 | nPed++; | |
ee299369 | 337 | } |
338 | } | |
e9dbb64a | 339 | |
82cbdfca | 340 | if (nPed > 0) |
341 | ped /= nPed; | |
e9dbb64a | 342 | else { |
343 | AliWarning("Could determine pedestal"); | |
82cbdfca | 344 | ped = 10; // put some small value as first guess |
e9dbb64a | 345 | } |
82cbdfca | 346 | |
e9dbb64a | 347 | Int_t max_found = 0; |
348 | Int_t i_max = 0; | |
349 | Float_t max = -1; | |
350 | Float_t tmax = 0; | |
351 | Float_t max_fit = gSig->GetN()*GetRawFormatTimeBinWidth(); | |
352 | Float_t min_after_sig = 9999; | |
353 | Int_t imin_after_sig = gSig->GetN(); | |
354 | Float_t tmin_after_sig = gSig->GetN()*GetRawFormatTimeBinWidth(); | |
355 | Int_t n_ped_after_sig = 0; | |
356 | ||
357 | for (Int_t i=kNPedSamples; i < gSig->GetN(); i++) { | |
358 | Double_t ttime, signal; | |
359 | gSig->GetPoint(i, ttime, signal) ; | |
360 | if (!max_found && signal > max) { | |
361 | i_max = i; | |
362 | tmax = ttime; | |
363 | max = signal; | |
364 | } | |
365 | else if ( max > ped + kNoiseThreshold ) { | |
366 | max_found = 1; | |
367 | min_after_sig = signal; | |
368 | imin_after_sig = i; | |
369 | tmin_after_sig = ttime; | |
370 | } | |
371 | if (max_found) { | |
372 | if ( signal < min_after_sig) { | |
373 | min_after_sig = signal; | |
374 | imin_after_sig = i; | |
375 | tmin_after_sig = ttime; | |
376 | } | |
377 | if (i > tmin_after_sig + 5) { // Two close peaks; end fit at minimum | |
378 | max_fit = tmin_after_sig; | |
379 | break; | |
380 | } | |
381 | if ( signal < ped + kNoiseThreshold) | |
382 | n_ped_after_sig++; | |
383 | if (n_ped_after_sig >= 5) { // include 5 pedestal bins after peak | |
384 | max_fit = ttime; | |
385 | break; | |
386 | } | |
5a056daa | 387 | } |
82cbdfca | 388 | } |
5a056daa | 389 | |
e9dbb64a | 390 | if ( max - ped > kNoiseThreshold ) { // else its noise |
391 | AliDebug(2,Form("Fitting max %d ped %d", max, ped)); | |
82cbdfca | 392 | signalF->SetParameter(0, ped) ; |
e9dbb64a | 393 | signalF->SetParameter(1, max - ped) ; |
394 | signalF->SetParameter(2, tmax) ; | |
395 | signalF->SetParLimits(2, 0, max_fit) ; | |
396 | gSig->Fit(signalF, "QRON", "", 0., max_fit); //, "QRON") ; | |
82cbdfca | 397 | amp = signalF->GetParameter(1); |
398 | time = signalF->GetParameter(2) - fgTimeTrigger; | |
ee299369 | 399 | } |
400 | return; | |
401 | } | |
402 | //__________________________________________________________________ | |
403 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) | |
404 | { | |
405 | // Shape of the electronics raw reponse: | |
406 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
407 | // | |
408 | // t' = (t - t0 + tau) / tau | |
409 | // F = A * t**N * exp( N * ( 1 - t) ) for t >= 0 | |
410 | // F = 0 for t < 0 | |
411 | // | |
412 | // parameters: | |
82cbdfca | 413 | // ped: par[0] |
414 | // A: par[1] // Amplitude = peak value | |
415 | // t0: par[2] | |
ee299369 | 416 | // tau: fgTau |
417 | // N: fgOrder | |
418 | // | |
419 | Double_t signal ; | |
82cbdfca | 420 | Double_t xx = ( x[0] - par[2] + fgTau ) / fgTau ; |
ee299369 | 421 | |
5a056daa | 422 | if (xx <= 0) |
82cbdfca | 423 | signal = par[0] ; |
ee299369 | 424 | else { |
82cbdfca | 425 | signal = par[0] + par[1] * TMath::Power(xx , fgOrder) * TMath::Exp(fgOrder * (1 - xx )) ; |
426 | ||
ee299369 | 427 | } |
428 | return signal ; | |
429 | } | |
430 | ||
431 | //__________________________________________________________________ | |
432 | Bool_t AliEMCALRawUtils::RawSampledResponse( | |
433 | const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const | |
434 | { | |
435 | // for a start time dtime and an amplitude damp given by digit, | |
436 | // calculates the raw sampled response AliEMCAL::RawResponseFunction | |
437 | ||
438 | const Int_t kRawSignalOverflow = 0x3FF ; | |
82cbdfca | 439 | const Int_t pedVal = 32; |
ee299369 | 440 | Bool_t lowGain = kFALSE ; |
441 | ||
442 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4); | |
82cbdfca | 443 | signalF.SetParameter(0, pedVal) ; |
444 | signalF.SetParameter(1, damp) ; | |
445 | signalF.SetParameter(2, dtime + fgTimeTrigger) ; | |
ee299369 | 446 | |
447 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { | |
82cbdfca | 448 | Double_t time = iTime * GetRawFormatTimeBinWidth() ; |
ee299369 | 449 | Double_t signal = signalF.Eval(time) ; |
450 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
451 | if ( adcH[iTime] > kRawSignalOverflow ){ // larger than 10 bits | |
452 | adcH[iTime] = kRawSignalOverflow ; | |
453 | lowGain = kTRUE ; | |
454 | } | |
455 | ||
456 | signal /= fHighLowGainFactor; | |
457 | ||
458 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
459 | if ( adcL[iTime] > kRawSignalOverflow) // larger than 10 bits | |
460 | adcL[iTime] = kRawSignalOverflow ; | |
461 | } | |
462 | return lowGain ; | |
463 | } |