Line data Source code
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 : #include <TROOT.h>
17 : #include <TTree.h>
18 : #include <TSystem.h>
19 : #include <TBenchmark.h>
20 : #include <TBrowser.h>
21 : #include <TObjectTable.h>
22 : #include <TRandom.h>
23 : #include <TF1.h>
24 : #include <cassert>
25 :
26 : // --- AliRoot header files ---
27 : #include "AliLog.h"
28 : #include "AliRun.h"
29 : #include "AliDigitizationInput.h"
30 : #include "AliRunLoader.h"
31 : #include "AliCDBManager.h"
32 : #include "AliCDBEntry.h"
33 : #include "AliEMCALDigit.h"
34 : #include "AliEMCAL.h"
35 : #include "AliEMCALLoader.h"
36 : #include "AliEMCALDigitizer.h"
37 : #include "AliEMCALSDigitizer.h"
38 : #include "AliEMCALGeometry.h"
39 : #include "AliEMCALCalibData.h"
40 : #include "AliEMCALCalibTime.h"
41 : #include "AliEMCALSimParam.h"
42 : #include "AliEMCALTriggerRawDigit.h"
43 : #include "AliCaloCalibPedestal.h"
44 :
45 : namespace
46 : {
47 : Double_t HeavisideTheta(Double_t x)
48 : {
49 : Double_t signal = 0.;
50 :
51 6256 : if (x > 0.) signal = 1.;
52 :
53 2176 : return signal;
54 : }
55 :
56 : Double_t AnalogFastORFunction(Double_t *x, Double_t *par)
57 : {
58 2176 : Double_t v0 = par[0];
59 1088 : Double_t t0 = par[1];
60 1088 : Double_t tr = par[2];
61 :
62 : Double_t R1 = 1000.;
63 : Double_t C1 = 33e-12;
64 : Double_t R2 = 1800;
65 : Double_t C2 = 22e-12;
66 :
67 1088 : Double_t t = x[0];
68 :
69 5440 : return (((0.8*(-((TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-t + t0)/(C1*R1))*
70 3264 : TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2)) +
71 3264 : C1*C2*R1*R2*(1 - (C2*TMath::Power(TMath::E(),(-t + t0)/(C2*R2))*R2)/(-(C1*R1) + C2*R2)))*v0*
72 2176 : HeavisideTheta(t - t0))/tr
73 3264 : - (0.8*(C1*C2*R1*R2 -
74 2176 : (TMath::Power(C1,2)*C2*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C1*R1))*
75 3264 : TMath::Power(R1,2)*R2)/(C1*R1 - C2*R2) +
76 1088 : (C1*TMath::Power(C2,2)*TMath::Power(TMath::E(),(-1.*t + t0 + 1.25*tr)/(C2*R2))*
77 3264 : R1*TMath::Power(R2,2))/(C1*R1 - C2*R2))*v0*
78 3264 : HeavisideTheta(t - t0 - 1.25*tr))/tr)/(C2*R1));
79 : }
80 : }
81 :
82 42 : ClassImp(AliEMCALDigitizer)
83 :
84 :
85 : //____________________________________________________________________________
86 : AliEMCALDigitizer::AliEMCALDigitizer()
87 0 : : AliDigitizer("",""),
88 0 : fDefaultInit(kTRUE),
89 0 : fDigitsInRun(0),
90 0 : fInit(0),
91 0 : fInput(0),
92 0 : fInputFileNames(0x0),
93 0 : fEventNames(0x0),
94 0 : fDigitThreshold(0),
95 0 : fMeanPhotonElectron(0),
96 0 : fGainFluctuations(0),
97 0 : fPinNoise(0),
98 0 : fTimeNoise(0),
99 0 : fTimeDelay(0),
100 0 : fTimeDelayFromOCDB(0),
101 0 : fTimeResolutionPar0(0),
102 0 : fTimeResolutionPar1(0),
103 0 : fADCchannelEC(0),
104 0 : fADCpedestalEC(0),
105 0 : fADCchannelECDecal(0),
106 0 : fTimeChannel(0),
107 0 : fTimeChannelDecal(0),
108 0 : fNADCEC(0),
109 0 : fEventFolderName(""),
110 0 : fFirstEvent(0),
111 0 : fLastEvent(0),
112 0 : fCalibData(0x0),
113 0 : fCalibTime(0x0),
114 0 : fSDigitizer(0x0)
115 0 : {
116 : // ctor
117 0 : InitParameters() ;
118 0 : fDigInput = 0 ; // We work in the standalone mode
119 0 : }
120 :
121 : //____________________________________________________________________________
122 : AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
123 0 : : AliDigitizer("EMCALDigitizer", alirunFileName),
124 0 : fDefaultInit(kFALSE),
125 0 : fDigitsInRun(0),
126 0 : fInit(0),
127 0 : fInput(0),
128 0 : fInputFileNames(0),
129 0 : fEventNames(0),
130 0 : fDigitThreshold(0),
131 0 : fMeanPhotonElectron(0),
132 0 : fGainFluctuations(0),
133 0 : fPinNoise(0),
134 0 : fTimeNoise(0),
135 0 : fTimeDelay(0),
136 0 : fTimeDelayFromOCDB(0),
137 0 : fTimeResolutionPar0(0),
138 0 : fTimeResolutionPar1(0),
139 0 : fADCchannelEC(0),
140 0 : fADCpedestalEC(0),
141 0 : fADCchannelECDecal(0),
142 0 : fTimeChannel(0),
143 0 : fTimeChannelDecal(0),
144 0 : fNADCEC(0),
145 0 : fEventFolderName(eventFolderName),
146 0 : fFirstEvent(0),
147 0 : fLastEvent(0),
148 0 : fCalibData(0x0),
149 0 : fCalibTime(0x0),
150 0 : fSDigitizer(0x0)
151 0 : {
152 : // ctor
153 0 : InitParameters() ;
154 0 : Init() ;
155 0 : fDigInput = 0 ; // We work in the standalone mode
156 0 : }
157 :
158 : //____________________________________________________________________________
159 : AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d)
160 0 : : AliDigitizer(d.GetName(),d.GetTitle()),
161 0 : fDefaultInit(d.fDefaultInit),
162 0 : fDigitsInRun(d.fDigitsInRun),
163 0 : fInit(d.fInit),
164 0 : fInput(d.fInput),
165 0 : fInputFileNames(d.fInputFileNames),
166 0 : fEventNames(d.fEventNames),
167 0 : fDigitThreshold(d.fDigitThreshold),
168 0 : fMeanPhotonElectron(d.fMeanPhotonElectron),
169 0 : fGainFluctuations(d.fGainFluctuations),
170 0 : fPinNoise(d.fPinNoise),
171 0 : fTimeNoise(d.fTimeNoise),
172 0 : fTimeDelay(d.fTimeDelay),
173 0 : fTimeDelayFromOCDB(d.fTimeDelayFromOCDB),
174 0 : fTimeResolutionPar0(d.fTimeResolutionPar0),
175 0 : fTimeResolutionPar1(d.fTimeResolutionPar1),
176 0 : fADCchannelEC(d.fADCchannelEC),
177 0 : fADCpedestalEC(d.fADCpedestalEC),
178 0 : fADCchannelECDecal(d.fADCchannelECDecal),
179 0 : fTimeChannel(d.fTimeChannel), fTimeChannelDecal(d.fTimeChannelDecal),
180 0 : fNADCEC(d.fNADCEC),
181 0 : fEventFolderName(d.fEventFolderName),
182 0 : fFirstEvent(d.fFirstEvent),
183 0 : fLastEvent(d.fLastEvent),
184 0 : fCalibData(d.fCalibData),
185 0 : fCalibTime(d.fCalibTime),
186 0 : fSDigitizer(d.fSDigitizer ? new AliEMCALSDigitizer(*d.fSDigitizer) : 0)
187 0 : {
188 : // copyy ctor
189 0 : }
190 :
191 : //____________________________________________________________________________
192 : AliEMCALDigitizer::AliEMCALDigitizer(AliDigitizationInput * rd)
193 1 : : AliDigitizer(rd,"EMCALDigitizer"),
194 1 : fDefaultInit(kFALSE),
195 1 : fDigitsInRun(0),
196 1 : fInit(0),
197 1 : fInput(0),
198 1 : fInputFileNames(0),
199 1 : fEventNames(0),
200 1 : fDigitThreshold(0),
201 1 : fMeanPhotonElectron(0),
202 1 : fGainFluctuations(0),
203 1 : fPinNoise(0.),
204 1 : fTimeNoise(0.),
205 1 : fTimeDelay(0.),
206 1 : fTimeDelayFromOCDB(0.),
207 1 : fTimeResolutionPar0(0.),
208 1 : fTimeResolutionPar1(0.),
209 1 : fADCchannelEC(0),
210 1 : fADCpedestalEC(0),
211 1 : fADCchannelECDecal(0),
212 1 : fTimeChannel(0),
213 1 : fTimeChannelDecal(0),
214 1 : fNADCEC(0),
215 1 : fEventFolderName(0),
216 1 : fFirstEvent(0),
217 1 : fLastEvent(0),
218 1 : fCalibData(0x0),
219 1 : fCalibTime(0x0),
220 1 : fSDigitizer(0x0)
221 5 : {
222 : // ctor Init() is called by RunDigitizer
223 1 : fDigInput = rd ;
224 2 : fEventFolderName = fDigInput->GetInputFolderName(0) ;
225 5 : SetTitle(dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
226 1 : InitParameters() ;
227 2 : }
228 :
229 : //____________________________________________________________________________
230 : AliEMCALDigitizer::~AliEMCALDigitizer()
231 6 : {
232 : //dtor
233 4 : delete [] fInputFileNames ;
234 4 : delete [] fEventNames ;
235 3 : if (fSDigitizer) delete fSDigitizer;
236 3 : }
237 :
238 : //____________________________________________________________________________
239 : void AliEMCALDigitizer::Digitize(Int_t event)
240 : {
241 : // Makes the digitization of the collected summable digits
242 : // for this it first creates the array of all EMCAL modules
243 : // filled with noise and after that adds contributions from
244 : // SDigits. This design helps to avoid scanning over the
245 : // list of digits to add contribution of any new SDigit.
246 : //
247 : // JLK 26-Jun-2008
248 : // Note that SDigit energy info is stored as an amplitude, so we
249 : // must call the Calibrate() method of the SDigitizer to convert it
250 : // back to an energy in GeV before adding it to the Digit
251 : //
252 :
253 8 : AliRunLoader *rl = AliRunLoader::Instance();
254 12 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
255 :
256 4 : if(!emcalLoader)
257 : {
258 0 : AliFatal("EMCALLoader is NULL!") ;
259 0 : return; // not needed but in case coverity complains ...
260 : }
261 :
262 : Int_t readEvent = event ;
263 4 : if (fDigInput) // fDigInput is data member from AliDigitizer
264 4 : readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
265 12 : AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
266 : readEvent, GetTitle(), fEventFolderName.Data())) ;
267 4 : rl->GetEvent(readEvent);
268 :
269 4 : TClonesArray * digits = emcalLoader->Digits() ;
270 4 : digits->Delete() ; //correct way to clear array when memory is
271 : //allocated by objects stored in array
272 :
273 : // Load Geometry
274 4 : if (!rl->GetAliRun())
275 : {
276 0 : AliFatal("Could not get AliRun from runLoader");
277 0 : return; // not needed but in case coverity complains ...
278 : }
279 :
280 4 : AliEMCAL * emcal = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
281 4 : AliEMCALGeometry *geom = emcal->GetGeometry();
282 7 : static int nEMC = geom->GetNCells();//max number of digits possible
283 12 : AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
284 :
285 4 : digits->Expand(nEMC) ;
286 :
287 : // RS create a digitizer on the fly
288 9 : if (!fSDigitizer) fSDigitizer = new AliEMCALSDigitizer(rl->GetFileName().Data());
289 4 : fSDigitizer->SetEventRange(0, -1) ;
290 :
291 : //-------------------------------------------------------
292 : //take all the inputs to add together and load the SDigits
293 4 : TObjArray * sdigArray = new TObjArray(fInput) ;
294 4 : sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
295 :
296 : Int_t i ;
297 : Int_t embed = kFALSE;
298 8 : for(i = 1 ; i < fInput ; i++)
299 : {
300 0 : TString tempo(fEventNames[i]) ;
301 0 : tempo += i ;
302 :
303 0 : AliRunLoader * rl2 = AliRunLoader::GetRunLoader(tempo) ;
304 :
305 0 : if (!rl2)
306 0 : rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
307 :
308 0 : if(!rl2)
309 : {
310 0 : AliFatal("Run Loader of second event not available!");
311 0 : return; // not needed but in case coverity complains ...
312 : }
313 :
314 0 : if (fDigInput)
315 0 : readEvent = dynamic_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
316 :
317 0 : AliInfo(Form("Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
318 :
319 0 : rl2->LoadSDigits();
320 : // rl2->LoadDigits();
321 0 : rl2->GetEvent(readEvent);
322 :
323 0 : if(!rl2->GetDetectorLoader("EMCAL"))
324 : {
325 0 : AliFatal("Loader of second input not found");
326 0 : return; // not needed but in case coverity complains ...
327 : }
328 :
329 0 : AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
330 :
331 0 : if(!emcalLoader2)
332 : {
333 0 : AliFatal("EMCAL Loader of second event not available!");
334 0 : return; // not needed but in case coverity complains ...
335 : }
336 :
337 : //Merge 2 simulated sdigits
338 0 : if(!emcalLoader2->SDigits()) continue;
339 :
340 0 : TClonesArray* sdigits2 = emcalLoader2->SDigits();
341 0 : sdigArray->AddAt(sdigits2, i) ;
342 :
343 : // Check if first sdigit is of embedded type, if so, handle the sdigits differently:
344 : // do not smear energy of embedded, do not add noise to any sdigits
345 0 : if( sdigits2->GetEntriesFast() <= 0 ) continue;
346 :
347 : //printf("Merged digit type: %d\n",dynamic_cast<AliEMCALDigit*> (sdigits2->At(0))->GetType());
348 0 : AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*> (sdigits2->At(0));
349 0 : if( digit2 && digit2->GetType()==AliEMCALDigit::kEmbedded ) embed = kTRUE;
350 :
351 0 : }// input loop
352 :
353 : //--------------------------------
354 : //Find the first tower with signal
355 4 : Int_t nextSig = nEMC + 1 ;
356 : TClonesArray * sdigits ;
357 16 : for(i = 0 ; i < fInput ; i++)
358 : {
359 4 : if(i > 0 && embed) continue;
360 :
361 12 : sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
362 4 : if(!sdigits)
363 : {
364 0 : AliDebug(1,"Null sdigit pointer");
365 : continue;
366 : }
367 :
368 4 : if (sdigits->GetEntriesFast() <= 0 )
369 : {
370 0 : AliDebug(1, "No sdigits entries");
371 : continue;
372 : }
373 :
374 12 : AliEMCALDigit *sd = dynamic_cast<AliEMCALDigit *>(sdigits->At(0));
375 4 : if(!sd)
376 : {
377 0 : AliDebug(1, "NULL sdigit pointer");
378 0 : continue;
379 : }
380 :
381 4 : Int_t curNext = sd->GetId() ;
382 4 : if(curNext < nextSig)
383 4 : nextSig = curNext ;
384 12 : AliDebug(1,Form("input %i : #sdigits %i \n",i, sdigits->GetEntriesFast()));
385 :
386 4 : }// input loop
387 :
388 12 : AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
389 :
390 4 : TArrayI index(fInput) ;
391 4 : index.Reset() ; //Set all indexes to zero
392 :
393 : AliEMCALDigit * digit ;
394 : AliEMCALDigit * curSDigit ;
395 :
396 : //---------------------------------------------
397 : //Put Noise contribution, smear time and energy
398 : Float_t timeResolution = 0;
399 : Int_t absID = -1 ;
400 141320 : for(absID = 0; absID < nEMC; absID++)
401 : {
402 : Float_t energy = 0 ;
403 :
404 : // amplitude set to zero, noise will be added later
405 : Float_t noiseTime = 0.;
406 211968 : if(!embed) noiseTime = TimeOfNoise(); //No need for embedded events?
407 211968 : new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0., noiseTime,kFALSE); // absID-1->absID
408 : //look if we have to add signal?
409 282624 : digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
410 :
411 70656 : if (!digit)
412 : {
413 0 : AliDebug(1,"Digit pointer is null");
414 0 : continue;
415 : }
416 :
417 70656 : if(absID==nextSig)
418 : {
419 : // Calculate time as time of the largest digit
420 268 : Float_t time = digit->GetTime() ;
421 268 : Float_t aTime= digit->GetAmplitude() ;
422 :
423 : // loop over input
424 268 : Int_t nInputs = fInput;
425 268 : if(embed) nInputs = 1; // In case of embedding, merge later real digits, do not smear energy and time of data
426 1072 : for(i = 0; i< nInputs ; i++)
427 : {
428 : //loop over (possible) merge sources
429 1072 : TClonesArray* sdtclarr = dynamic_cast<TClonesArray *>(sdigArray->At(i));
430 268 : if(sdtclarr)
431 : {
432 268 : Int_t sDigitEntries = sdtclarr->GetEntriesFast();
433 :
434 1876 : if(sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
435 : else curSDigit = 0 ;
436 :
437 : //May be several digits will contribute from the same input
438 1336 : while(curSDigit && (curSDigit->GetId() == absID))
439 : {
440 : //Shift primary to separate primaries belonging different inputs
441 : Int_t primaryoffset = i ;
442 536 : if(fDigInput) primaryoffset = fDigInput->GetMask(i) ;
443 268 : curSDigit->ShiftPrimary(primaryoffset) ;
444 :
445 268 : if(curSDigit->GetAmplitude()>aTime)
446 : {
447 268 : aTime = curSDigit->GetAmplitude();
448 268 : time = curSDigit->GetTime();
449 268 : }
450 :
451 804 : *digit = *digit + *curSDigit ; //adds amplitudes of each digit
452 :
453 536 : index[i]++ ;
454 :
455 1856 : if( sDigitEntries > index[i] ) curSDigit = dynamic_cast<AliEMCALDigit*>(sdtclarr->At(index[i])) ;
456 : else curSDigit = 0 ;
457 : }// while
458 268 : }// source exists
459 : }// loop over merging sources
460 :
461 : //Here we convert the summed amplitude to an energy in GeV only for simulation or mixing of simulations
462 268 : energy = fSDigitizer->Calibrate(digit->GetAmplitude()) ; // GeV
463 :
464 : // add fluctuations for photo-electron creation
465 : // corrected fluctuations after comparison with beam test, Paraskevi Ganoti (06/11/2011)
466 268 : Float_t fluct = static_cast<Float_t>((energy*fMeanPhotonElectron)/fGainFluctuations);
467 536 : energy *= static_cast<Float_t>(gRandom->Poisson(fluct)) / fluct ;
468 :
469 : //calculate and set time
470 268 : digit->SetTime(time) ;
471 :
472 : //Find next signal module
473 268 : nextSig = nEMC + 1 ;
474 1072 : for(i = 0 ; i < nInputs ; i++)
475 : {
476 1072 : sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
477 :
478 268 : if(sdigits)
479 : {
480 : Int_t curNext = nextSig ;
481 804 : if(sdigits->GetEntriesFast() > index[i])
482 : {
483 1320 : AliEMCALDigit * tmpdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]));
484 528 : if ( tmpdigit ) curNext = tmpdigit->GetId() ;
485 264 : }
486 :
487 532 : if(curNext < nextSig) nextSig = curNext ;
488 268 : }// sdigits exist
489 : } // input loop
490 :
491 268 : }//absID==nextSig
492 :
493 : // add the noise now, no need for embedded events with real data
494 70656 : if(!embed)
495 141312 : energy += gRandom->Gaus(0., fPinNoise) ;
496 :
497 : // JLK 26-June-2008
498 : //Now digitize the energy according to the fSDigitizer method,
499 : //which merely converts the energy to an integer. Later we will
500 : //check that the stored value matches our allowed dynamic ranges
501 141312 : digit->SetAmplitude(fSDigitizer->Digitize(energy)) ;
502 :
503 : //Set time resolution with final energy
504 70656 : timeResolution = GetTimeResolution(energy);
505 141312 : digit->SetTime(gRandom->Gaus(digit->GetTime(),timeResolution) ) ;
506 353280 : AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
507 : absID, energy, nextSig));
508 : //Add delay to time
509 70656 : digit->SetTime(digit->GetTime()) ;
510 :
511 70656 : } // for(absID = 0; absID < nEMC; absID++)
512 :
513 : //---------------------------------------------------------
514 : //Embed simulated amplitude (and time?) to real data digits
515 4 : if ( embed )
516 : {
517 0 : for(Int_t input = 1; input<fInput; input++)
518 : {
519 0 : TClonesArray *realDigits = dynamic_cast<TClonesArray*> (sdigArray->At(input));
520 0 : if(!realDigits)
521 : {
522 0 : AliDebug(1,"No sdigits to merge\n");
523 0 : continue;
524 : }
525 :
526 0 : for(Int_t i2 = 0 ; i2 < realDigits->GetEntriesFast() ; i2++)
527 : {
528 0 : AliEMCALDigit * digit2 = dynamic_cast<AliEMCALDigit*>( realDigits->At(i2) ) ;
529 0 : if ( !digit2 ) continue;
530 :
531 0 : digit = dynamic_cast<AliEMCALDigit*>( digits->At(digit2->GetId()) ) ;
532 0 : if ( !digit ) continue;
533 :
534 : // Put the embedded cell energy in same units as simulated sdigits -> transform to energy units
535 : // and multiply to get a big int.
536 0 : Float_t time2 = digit2->GetTime();
537 0 : Float_t e2 = digit2->GetAmplitude();
538 0 : CalibrateADCTime(e2,time2,digit2->GetId());
539 0 : Float_t amp2 = fSDigitizer->Digitize(e2);
540 :
541 0 : Float_t e0 = digit ->GetAmplitude();
542 0 : if(e0>0.01)
543 0 : AliDebug(1,Form("digit 1: Abs ID %d, amp %f, type %d, time %e; digit2: Abs ID %d, amp %f, type %d, time %e\n",
544 : digit ->GetId(),digit ->GetAmplitude(), digit ->GetType(), digit->GetTime(),
545 : digit2->GetId(),amp2, digit2->GetType(), time2 ));
546 :
547 : // Sum amplitudes, change digit amplitude
548 0 : digit->SetAmplitude( digit->GetAmplitude() + amp2);
549 :
550 : // Rough assumption, assign time of the largest of the 2 digits,
551 : // data or signal to the final digit.
552 0 : if(amp2 > digit->GetAmplitude()) digit->SetTime(time2);
553 :
554 : //Mark the new digit as embedded
555 0 : digit->SetType(AliEMCALDigit::kEmbedded);
556 :
557 0 : if(digit2->GetAmplitude()>0.01 && e0> 0.01 )
558 0 : AliDebug(1,Form("Embedded digit: Abs ID %d, amp %f, type %d\n",
559 : digit->GetId(), digit->GetAmplitude(), digit->GetType()));
560 0 : }//loop digit 2
561 0 : }//input loop
562 0 : }//embed
563 :
564 : //JLK is it better to call Clear() here?
565 8 : delete sdigArray ; //We should not delete its contents
566 :
567 : //------------------------------
568 : // Remove digits below ADC thresholds or
569 : // dead in OCDB, decalibrate them
570 : //
571 :
572 4 : Float_t ampADC = 0;
573 4 : Float_t time = 0;
574 : Int_t idigit = 0;
575 141320 : for(i = 0 ; i < nEMC ; i++)
576 : {
577 282624 : digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
578 70656 : if ( !digit )
579 : {
580 0 : digits->RemoveAt(i) ; // It should not happen, but just in case
581 : continue;
582 : }
583 :
584 : // Get the time and energy before calibration
585 141312 : ampADC = fSDigitizer->Calibrate(digit->GetAmplitude()) ;
586 :
587 70656 : time = digit->GetTime();
588 :
589 : // Then digitize energy (GeV to ADC) and shift time
590 : // using the calibration constants of the OCDB
591 70656 : DigitizeEnergyTime(ampADC, time, digit->GetId()) ;
592 :
593 : // Skip digits with below 3 ADC or found dead in OCDB
594 70772 : if(ampADC < fDigitThreshold || IsDead(digit->GetId()))
595 : {
596 70598 : digits->RemoveAt(i) ;
597 : continue;
598 : }
599 :
600 : // Set digit final values
601 58 : digit->SetIndexInList(idigit++) ;
602 58 : digit->SetAmplitude(ampADC) ;
603 58 : digit->SetTime(time);
604 :
605 58 : } // digit loop
606 :
607 4 : digits->Compress() ;
608 :
609 4 : Int_t ndigits = digits->GetEntriesFast();
610 :
611 4 : if(idigit != ndigits)
612 0 : AliFatal(Form("Total number of digits in array %d different to expected %d",ndigits,idigit));
613 :
614 20 : AliDebug(1,Form("Number of recorded digits is %d",ndigits));
615 :
616 8 : }
617 :
618 : /// JLK 26-June-2008
619 : /// Returns digitized value of the energy and shifted time in a cell absId
620 : /// in the way those parameters are stored in the data digits.
621 : /// This is done using the calibration constants stored in the OCDB
622 : /// or default values if no fCalibData object is found.
623 : /// This effectively converts everything to match the dynamic range
624 : /// of the real data we will collect
625 : ///
626 : /// \param energy: digit energy in GeV
627 : /// \param time: time at generation
628 : /// \param absId: tower ID
629 : ///
630 : //_____________________________________________________________________
631 : void AliEMCALDigitizer::DigitizeEnergyTime(Float_t & energy, Float_t & time, Int_t absId)
632 : {
633 : // Load Geometry and cell indeces
634 141312 : const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
635 :
636 70656 : if (geom==0)
637 : {
638 0 : AliFatal("Did not get geometry from EMCALLoader");
639 0 : return;
640 : }
641 :
642 70656 : Int_t iSupMod = -1;
643 70656 : Int_t nModule = -1;
644 70656 : Int_t nIphi = -1;
645 70656 : Int_t nIeta = -1;
646 70656 : Int_t iphi = -1;
647 70656 : Int_t ieta = -1;
648 70656 : Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
649 70656 : if(!bCell)
650 0 : Error("DigitizeEnergyTime","Wrong cell id number : absId %i ", absId) ;
651 70656 : geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
652 :
653 : // Recover parameters from OCDB for this channel
654 70656 : if(fCalibData)
655 : {
656 : // Energy
657 70656 : fADCpedestalEC = fCalibData->GetADCpedestal (iSupMod,ieta,iphi);
658 70656 : fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
659 70656 : fADCchannelECDecal = fCalibData->GetADCchannelDecal (iSupMod,ieta,iphi);
660 70656 : }
661 :
662 70656 : if(fCalibTime)
663 : {
664 : // Time
665 : // Recover parameters for bunch crossing number equal to 0
666 : // (has simulation different bunch crossings stored? not for the moment)
667 : // Time stored in ns, pass to ns
668 70656 : fTimeChannel = fCalibTime->GetTimeChannel (iSupMod,ieta,iphi,0) * 1e-9;
669 70656 : fTimeChannelDecal = fCalibTime->GetTimeChannelDecal(iSupMod,ieta,iphi) * 1e-9;
670 70656 : }
671 :
672 : // Apply calibration to get ADC counts and partial decalibration as especified in OCDB
673 70656 : energy = (energy + fADCpedestalEC)/fADCchannelEC/fADCchannelECDecal ;
674 :
675 70656 : if ( energy > fNADCEC ) energy = fNADCEC ;
676 :
677 : // Apply shift to time, if requested and calibration parameter is available,
678 : // if not, apply fix shift
679 70656 : if ( fTimeDelayFromOCDB )
680 0 : time += fTimeChannel - fTimeChannelDecal;
681 : else
682 70656 : time += fTimeDelay;
683 141312 : }
684 :
685 : //_____________________________________________________________________
686 : void AliEMCALDigitizer::DecalibrateTrigger(AliEMCALDigit *digit)
687 : {
688 : // Decalibrate, used in Trigger digits
689 :
690 536 : if ( !fCalibData ) return ;
691 :
692 : // Load Geometry
693 268 : const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
694 :
695 268 : if (!geom)
696 : {
697 0 : AliFatal("Did not get geometry from EMCALLoader");
698 0 : return;
699 : }
700 :
701 : // Recover the digit location in row-column-SM
702 268 : Int_t absId = digit->GetId();
703 268 : Int_t iSupMod = -1;
704 268 : Int_t nModule = -1;
705 268 : Int_t nIphi = -1;
706 268 : Int_t nIeta = -1;
707 268 : Int_t iphi = -1;
708 268 : Int_t ieta = -1;
709 :
710 268 : Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
711 268 : if (!bCell) Error("Decalibrate","Wrong cell id number : absId %i ", absId) ;
712 :
713 268 : geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
714 :
715 : // Now decalibrate
716 268 : Float_t adcChannel = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
717 268 : Float_t adcChannelOnline = fCalibData->GetADCchannelOnline(iSupMod,ieta,iphi);
718 : //printf("calib param for (SM%d,col%d,row%d): %1.4f - online param: %1.4f \n",iSupMod,ieta,iphi,adcChannel,adcChannelOnline);
719 268 : Float_t factor = 1./(adcChannel/adcChannelOnline);
720 :
721 536 : *digit = *digit * factor;
722 :
723 536 : }
724 :
725 : //_____________________________________________________________________
726 : void AliEMCALDigitizer::CalibrateADCTime(Float_t & adc, Float_t & time, const Int_t absId)
727 : {
728 : // Returns the energy in a cell absId with a given adc value
729 : // using the calibration constants stored in the OCDB. Time also corrected from parameter in OCDB
730 : // Used in case of embedding, transform ADC counts from real event into energy
731 : // so that we can add the energy of the simulated sdigits which are in energy
732 : // units.
733 : // Same as in AliEMCALClusterizer::Calibrate() but here we do not reject channels being marked as hot
734 : // or with time out of window
735 :
736 : // Load Geometry
737 0 : const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
738 :
739 0 : if (!geom)
740 : {
741 0 : AliFatal("Did not get geometry from EMCALLoader");
742 0 : return;
743 : }
744 :
745 0 : Int_t iSupMod = -1;
746 0 : Int_t nModule = -1;
747 0 : Int_t nIphi = -1;
748 0 : Int_t nIeta = -1;
749 0 : Int_t iphi = -1;
750 0 : Int_t ieta = -1;
751 0 : Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
752 0 : if(!bCell) Error("CalibrateADCTime","Wrong cell id number : absId %i ", absId) ;
753 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
754 :
755 : // Energy calibration
756 0 : if(fCalibData)
757 : {
758 0 : fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
759 0 : fADCchannelEC = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
760 0 : }
761 :
762 0 : adc = adc * fADCchannelEC - fADCpedestalEC;
763 :
764 : // Time calibration
765 : // Assign bunch crossing number equal to 0 (has simulation different bunch crossings? Not for the moment)
766 0 : if(fCalibTime && fTimeDelayFromOCDB)
767 0 : fTimeChannel = fCalibTime->GetTimeChannel(iSupMod,ieta,iphi,0) * 1e-9; // pass from ns to s
768 :
769 0 : time -= fTimeChannel;
770 0 : }
771 :
772 :
773 : //____________________________________________________________________________
774 : void AliEMCALDigitizer::Digitize(Option_t *option)
775 : {
776 : // Steering method to process digitization for events
777 : // in the range from fFirstEvent to fLastEvent.
778 : // This range is optionally set by SetEventRange().
779 : // if fLastEvent=-1, then process events until the end.
780 : // by default fLastEvent = fFirstEvent (process only one event)
781 :
782 8 : if (!fInit)
783 : { // to prevent overwrite existing file
784 0 : Error( "Digitize", "Give a version name different from %s", fEventFolderName.Data() ) ;
785 0 : return ;
786 : }
787 :
788 4 : if (strstr(option,"print"))
789 : {
790 0 : Print();
791 0 : return ;
792 : }
793 :
794 4 : if(strstr(option,"tim"))
795 0 : gBenchmark->Start("EMCALDigitizer");
796 :
797 4 : AliRunLoader *rl = AliRunLoader::Instance();
798 :
799 12 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
800 4 : if(!emcalLoader)
801 : {
802 0 : AliFatal("Did not get the Loader");
803 0 : return; // coverity
804 : }
805 :
806 4 : if (fLastEvent == -1)
807 0 : fLastEvent = rl->GetNumberOfEvents() - 1 ;
808 4 : else if (fDigInput)
809 4 : fLastEvent = fFirstEvent ; // what is this ??
810 :
811 4 : Int_t nEvents = fLastEvent - fFirstEvent + 1;
812 : Int_t ievent = -1;
813 :
814 12 : AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
815 4 : if(!emcal)
816 : {
817 0 : AliFatal("Did not get the AliEMCAL pointer");
818 0 : return; // coverity
819 : }
820 :
821 4 : AliEMCALGeometry *geom = emcal->GetGeometry();
822 4 : if(!geom)
823 : {
824 0 : AliFatal("Geometry pointer null");
825 0 : return; // fix for coverity
826 : }
827 :
828 4 : const Int_t nTRU = geom->GetNTotalTRU();
829 4 : TClonesArray* digitsTMP = new TClonesArray("AliEMCALDigit", nTRU*96);
830 4 : TClonesArray* digitsTRG = new TClonesArray("AliEMCALTriggerRawDigit", nTRU*96);
831 :
832 4 : rl->LoadSDigits("EMCAL");
833 16 : for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++)
834 : {
835 4 : rl->GetEvent(ievent);
836 :
837 4 : Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
838 :
839 4 : WriteDigits() ;
840 :
841 : //Trigger Digits
842 : //-------------------------------------
843 :
844 4 : Digits2FastOR(digitsTMP, digitsTRG);
845 :
846 4 : WriteDigits(digitsTRG);
847 :
848 4 : (emcalLoader->TreeD())->Fill();
849 :
850 4 : emcalLoader->WriteDigits("OVERWRITE");
851 :
852 4 : Unload();
853 :
854 4 : digitsTRG ->Delete();
855 4 : digitsTMP ->Delete();
856 :
857 : //-------------------------------------
858 :
859 4 : if(strstr(option,"deb"))
860 0 : PrintDigits(option);
861 4 : if(strstr(option,"table")) gObjectTable->Print();
862 :
863 : //increment the total number of Digits per run
864 4 : fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;
865 : }//loop
866 :
867 4 : if(strstr(option,"tim"))
868 : {
869 0 : gBenchmark->Stop("EMCALDigitizer");
870 0 : Float_t cputime = gBenchmark->GetCpuTime("EMCALDigitizer");
871 : Float_t avcputime = cputime;
872 0 : if(nEvents==0) avcputime = 0 ;
873 0 : AliInfo(Form("Digitize: took %f seconds for Digitizing %f seconds per event", cputime, avcputime)) ;
874 0 : }
875 8 : }
876 :
877 : //__________________________________________________________________
878 : Float_t AliEMCALDigitizer::GetTimeResolution(Float_t energy) const
879 : {
880 : // Assign a smeared time to the digit, from observed distribution in LED system (?).
881 : // From F. Blanco
882 : Float_t res = -1;
883 141312 : if (energy > 0)
884 : {
885 70934 : res = TMath::Sqrt(fTimeResolutionPar0 +
886 35467 : fTimeResolutionPar1/(energy*energy) );
887 35467 : }
888 : // parametrization above is for ns. Convert to seconds:
889 70656 : return res*1e-9;
890 : }
891 :
892 : //____________________________________________________________________________
893 : void AliEMCALDigitizer::Digits2FastOR(TClonesArray* digitsTMP, TClonesArray* digitsTRG)
894 : {
895 : // FEE digits afterburner to produce TRG digits
896 : // we are only interested in the FEE digit deposited energy
897 : // to be converted later into a voltage value
898 :
899 : // push the FEE digit to its associated FastOR (numbered from 0:95)
900 : // TRU is in charge of summing module digits
901 :
902 8 : AliRunLoader *runLoader = AliRunLoader::Instance();
903 :
904 12 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
905 4 : if (!emcalLoader) AliFatal("Did not get the Loader");
906 :
907 4 : const AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
908 4 : if(!geom)
909 : {
910 0 : AliFatal("Geometry pointer null");
911 0 : return; // fix for coverity
912 : }
913 :
914 : // build FOR from simulated digits
915 : // and xfer to the corresponding TRU input (mapping)
916 :
917 4 : TClonesArray* sdigits = emcalLoader->SDigits();
918 :
919 4 : TClonesArray *digits = (TClonesArray*)sdigits->Clone();
920 :
921 12 : AliDebug(999,Form("=== %d SDigits to trigger digits ===",digits->GetEntriesFast()));
922 :
923 4 : TIter NextDigit(digits);
924 552 : while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
925 : {
926 536 : if (IsDead(digit)) continue;
927 :
928 268 : DecalibrateTrigger(digit);
929 :
930 268 : Int_t id = digit->GetId();
931 :
932 268 : Int_t trgid;
933 536 : if (geom->GetFastORIndexFromCellIndex(id , trgid))
934 : {
935 1340 : AliDebug(1,Form("trigger digit id: %d from cell id: %d\n",trgid,id));
936 :
937 536 : AliEMCALDigit* d = static_cast<AliEMCALDigit*>(digitsTMP->At(trgid));
938 :
939 268 : if (!d)
940 : {
941 474 : new((*digitsTMP)[trgid]) AliEMCALDigit(*digit);
942 316 : d = (AliEMCALDigit*)digitsTMP->At(trgid);
943 158 : d->SetId(trgid);
944 158 : }
945 : else
946 : {
947 330 : *d = *d + *digit;
948 : }
949 268 : }
950 268 : }
951 :
952 16 : if (AliDebugLevel()) printf("Number of TRG digits: %d\n",digitsTMP->GetEntriesFast());
953 :
954 : Int_t nSamples = 32;
955 8 : Int_t *timeSamples = new Int_t[nSamples];
956 :
957 12 : NextDigit = TIter(digitsTMP);
958 332 : while (AliEMCALDigit* digit = (AliEMCALDigit*)NextDigit())
959 : {
960 158 : if (digit)
961 : {
962 158 : Int_t id = digit->GetId();
963 : Float_t time = 50.e-9;
964 :
965 : Double_t depositedEnergy = 0.;
966 418 : for (Int_t j = 1; j <= digit->GetNprimary(); j++) depositedEnergy += digit->GetDEPrimary(j);
967 :
968 632 : if (AliDebugLevel()) printf("Deposited Energy: %f\n", depositedEnergy);
969 :
970 : // FIXME: Check digit time!
971 158 : if (depositedEnergy) {
972 68 : depositedEnergy += gRandom->Gaus(0., .08);
973 34 : DigitalFastOR(time, depositedEnergy, timeSamples, nSamples);
974 :
975 2244 : for (Int_t j=0;j<nSamples;j++) {
976 4352 : if (AliDebugLevel()) printf("timeSamples[%d]: %d\n",j,timeSamples[j]);
977 1088 : timeSamples[j] = ((j << 16) | (timeSamples[j] & 0xFFFF));
978 : }
979 :
980 136 : new((*digitsTRG)[digitsTRG->GetEntriesFast()]) AliEMCALTriggerRawDigit(id, timeSamples, nSamples);
981 :
982 136 : if (AliDebugLevel()) ((AliEMCALTriggerRawDigit*)digitsTRG->At(digitsTRG->GetEntriesFast() - 1))->Print("");
983 : }
984 158 : }
985 158 : }
986 :
987 8 : delete [] timeSamples;
988 :
989 16 : if (digits && digits->GetEntriesFast()) digits->Delete();
990 8 : }
991 :
992 : //____________________________________________________________________________
993 : void AliEMCALDigitizer::DigitalFastOR( Double_t time, Double_t dE, Int_t timeSamples[], Int_t nSamples )
994 : {
995 : // parameters:
996 : // id: 0..95
997 : const Int_t reso = 12; // 11-bit resolution ADC
998 : const Double_t vFSR = 2.; // Full scale input voltage range 2V (p-p)
999 : //const Double_t dNe = 125; // signal of the APD per MeV of energy deposit in a tower: 125 photo-e-/MeV @ M=30
1000 : const Double_t dNe = 125/1.3; // F-ALTRO max V. FEE: factor 4
1001 : const Double_t vA = .136e-6; // CSP output range: 0.136uV/e-
1002 : const Double_t rise = 50e-9; // rise time (10-90%) of the FastOR signal before shaping
1003 :
1004 : const Double_t kTimeBinWidth = 25E-9; // sampling frequency (40MHz)
1005 :
1006 68 : Double_t vV = 1000. * dE * dNe * vA; // GeV 2 MeV
1007 :
1008 34 : TF1 signalF("signal", AnalogFastORFunction, 0, nSamples * kTimeBinWidth, 3);
1009 34 : signalF.SetParameter( 0, vV );
1010 34 : signalF.SetParameter( 1, time ); // FIXME: when does the signal arrive? Might account for cable lengths
1011 34 : signalF.SetParameter( 2, rise );
1012 :
1013 2244 : for (Int_t iTime=0; iTime<nSamples; iTime++)
1014 : {
1015 : // FIXME: add noise (probably not simply Gaussian) according to DA measurements
1016 : // probably plan an access to OCDB
1017 1088 : Double_t sig = signalF.Eval(iTime * kTimeBinWidth);
1018 1088 : if (TMath::Abs(sig) > vFSR/2.) {
1019 0 : AliError("Signal overflow!");
1020 0 : timeSamples[iTime] = (1 << reso) - 1;
1021 0 : } else {
1022 5440 : AliDebug(999,Form("iTime: %d sig: %f\n",iTime,sig));
1023 1088 : timeSamples[iTime] = ((1 << reso) / vFSR) * sig + 0.5;
1024 : }
1025 : }
1026 34 : }
1027 :
1028 : //____________________________________________________________________________
1029 : Bool_t AliEMCALDigitizer::Init()
1030 : {
1031 : // Makes all memory allocations
1032 2 : fInit = kTRUE ;
1033 3 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1034 :
1035 1 : if ( emcalLoader == 0 ) {
1036 0 : AliFatal("Could not obtain the AliEMCALLoader");
1037 0 : return kFALSE;
1038 : }
1039 :
1040 1 : fFirstEvent = 0 ;
1041 1 : fLastEvent = fFirstEvent ;
1042 :
1043 1 : if(fDigInput)
1044 1 : fInput = fDigInput->GetNinputs() ;
1045 : else
1046 0 : fInput = 1 ;
1047 :
1048 5 : fInputFileNames = new TString[fInput] ;
1049 5 : fEventNames = new TString[fInput] ;
1050 1 : fInputFileNames[0] = GetTitle() ;
1051 1 : fEventNames[0] = fEventFolderName.Data() ;
1052 : Int_t index ;
1053 2 : for (index = 1 ; index < fInput ; index++) {
1054 0 : fInputFileNames[index] = dynamic_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
1055 0 : TString tempo = fDigInput->GetInputFolderName(index) ;
1056 0 : fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fDigInput
1057 0 : }
1058 :
1059 : // Calibration instances
1060 1 : fCalibData = emcalLoader->CalibData();
1061 1 : fCalibTime = emcalLoader->CalibTime();
1062 :
1063 1 : return fInit ;
1064 1 : }
1065 :
1066 : //____________________________________________________________________________
1067 : void AliEMCALDigitizer::InitParameters()
1068 : {
1069 : // Parameter initialization for digitizer
1070 :
1071 : // Get the parameters from the OCDB via the loader
1072 2 : AliRunLoader *rl = AliRunLoader::Instance();
1073 3 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1074 : AliEMCALSimParam * simParam = 0x0;
1075 2 : if(emcalLoader) simParam = emcalLoader->SimulationParameters();
1076 :
1077 1 : if(!simParam){
1078 0 : simParam = AliEMCALSimParam::GetInstance();
1079 0 : AliWarning("Simulation Parameters not available in OCDB?");
1080 0 : }
1081 :
1082 1 : fMeanPhotonElectron = simParam->GetMeanPhotonElectron() ; // 4400; // electrons per GeV
1083 1 : fGainFluctuations = simParam->GetGainFluctuations() ; // 15.0;
1084 :
1085 1 : fPinNoise = simParam->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data
1086 1 : if (fPinNoise < 0.0001 )
1087 0 : Warning("InitParameters", "No noise added\n") ;
1088 1 : fTimeNoise = simParam->GetTimeNoise(); // 1.28E-5 s
1089 1 : fDigitThreshold = simParam->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
1090 1 : fTimeResolutionPar0 = simParam->GetTimeResolutionPar0();
1091 1 : fTimeResolutionPar1 = simParam->GetTimeResolutionPar1();
1092 1 : fTimeDelay = simParam->GetTimeDelay(); //600e-9 ; // 600 ns
1093 1 : fTimeDelayFromOCDB = simParam->IsTimeDelayFromOCDB();
1094 :
1095 : // These defaults are normally not used.
1096 : // Values are read from calibration database instead
1097 1 : fADCchannelEC = 0.0162; // Nominal value set online for most of the channels, MIP peak sitting at ~266./16/1024
1098 1 : fADCpedestalEC = 0.0 ; // GeV
1099 1 : fADCchannelECDecal = 1.0; // No decalibration by default
1100 1 : fTimeChannel = 0.0; // No time calibration by default
1101 1 : fTimeChannelDecal = 0.0; // No time decalibration by default
1102 :
1103 1 : fNADCEC = simParam->GetNADCEC();//(Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536
1104 :
1105 3 : AliDebug(2,Form("Mean Photon Electron %d, Gain Fluct. %2.1f; Noise: APD %f, Time %f; Digit Threshold %d,Time Resolution Par0 %g Par1 %g,NADCEC %d",
1106 : fMeanPhotonElectron, fGainFluctuations, fPinNoise,fTimeNoise, fDigitThreshold,fTimeResolutionPar0,fTimeResolutionPar1,fNADCEC));
1107 :
1108 1 : }
1109 :
1110 : //__________________________________________________________________
1111 : void AliEMCALDigitizer::Print1(Option_t * option)
1112 : { // 19-nov-04 - just for convenience
1113 0 : Print();
1114 0 : PrintDigits(option);
1115 0 : }
1116 :
1117 : //__________________________________________________________________
1118 : void AliEMCALDigitizer::Print (Option_t * ) const
1119 : {
1120 : // Print Digitizer's parameters
1121 0 : printf("Print: \n------------------- %s -------------", GetName() ) ;
1122 0 : if( strcmp(fEventFolderName.Data(), "") != 0 ){
1123 0 : printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
1124 :
1125 : Int_t nStreams ;
1126 0 : if (fDigInput)
1127 0 : nStreams = GetNInputStreams() ;
1128 : else
1129 0 : nStreams = fInput ;
1130 :
1131 : AliRunLoader *rl=0;
1132 :
1133 : Int_t index = 0 ;
1134 0 : for (index = 0 ; index < nStreams ; index++) {
1135 0 : TString tempo(fEventNames[index]) ;
1136 0 : tempo += index ;
1137 0 : if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
1138 0 : rl = AliRunLoader::Open(fInputFileNames[index], tempo) ;
1139 0 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
1140 0 : if(emcalLoader){
1141 0 : TString fileName( emcalLoader->GetSDigitsFileName() ) ;
1142 0 : if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1143 0 : fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index] + ".root" ;
1144 0 : printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ;
1145 0 : }//loader
1146 0 : }
1147 :
1148 0 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1149 :
1150 0 : if(emcalLoader) printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
1151 0 : else printf("\nNULL LOADER");
1152 :
1153 0 : printf("\nWith following parameters:\n") ;
1154 0 : printf(" Electronics noise in EMC, APD (fPinNoise) = %f, Time = %f \n", fPinNoise, fTimeNoise) ;
1155 0 : printf(" Threshold in Tower (fDigitThreshold) = %d\n", fDigitThreshold) ;
1156 0 : printf("---------------------------------------------------\n") ;
1157 0 : }
1158 : else
1159 0 : printf("Print: AliEMCALDigitizer not initialized") ;
1160 0 : }
1161 :
1162 : //__________________________________________________________________
1163 : void AliEMCALDigitizer::PrintDigits(Option_t * option)
1164 : {
1165 : //utility method for printing digit information
1166 :
1167 0 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1168 0 : if(emcalLoader){
1169 0 : TClonesArray * digits = emcalLoader->Digits() ;
1170 0 : TClonesArray * sdigits = emcalLoader->SDigits() ;
1171 :
1172 0 : printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ;
1173 0 : printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
1174 :
1175 0 : if(strstr(option,"all")){
1176 : //loop over digits
1177 : AliEMCALDigit * digit;
1178 0 : printf("\nEMC digits (with primaries):\n") ;
1179 0 : printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
1180 : Int_t index ;
1181 0 : for (index = 0 ; index < digits->GetEntries() ; index++) {
1182 0 : digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
1183 0 : if(digit){
1184 0 : printf("\n%6d %8f %6.5e %4d %2d : ",
1185 0 : digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
1186 : Int_t iprimary;
1187 0 : for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
1188 0 : printf("%d ",digit->GetPrimary(iprimary+1) ) ;
1189 : }
1190 0 : }// digit exists
1191 : }// loop
1192 0 : }
1193 0 : printf("\n");
1194 0 : }// loader exists
1195 0 : else printf("NULL LOADER, cannot print\n");
1196 0 : }
1197 :
1198 : //__________________________________________________________________
1199 : Float_t AliEMCALDigitizer::TimeOfNoise(void)
1200 : {
1201 : // Calculates the time signal generated by noise
1202 : //printf("Time noise %e\n",fTimeNoise);
1203 141312 : return gRandom->Rndm() * fTimeNoise;
1204 : }
1205 :
1206 : //__________________________________________________________________
1207 : void AliEMCALDigitizer::Unload()
1208 : {
1209 : // Unloads the SDigits and Digits
1210 : AliRunLoader *rl=0;
1211 :
1212 : Int_t i ;
1213 12 : for(i = 1 ; i < fInput ; i++){
1214 0 : TString tempo(fEventNames[i]) ;
1215 0 : tempo += i;
1216 0 : if ((rl = AliRunLoader::GetRunLoader(tempo)))
1217 0 : rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ;
1218 0 : }
1219 12 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1220 8 : if(emcalLoader)emcalLoader->UnloadDigits() ;
1221 0 : else AliFatal("Did not get the loader");
1222 4 : }
1223 :
1224 : //_________________________________________________________________________________________
1225 : void AliEMCALDigitizer::WriteDigits()
1226 : { // Makes TreeD in the output file.
1227 : // Check if branch already exists:
1228 : // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1229 : // already existing branches.
1230 : // else creates branch with Digits, named "EMCAL", title "...",
1231 : // and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1232 : // and names of files, from which digits are made.
1233 :
1234 16 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1235 :
1236 4 : if(emcalLoader){
1237 4 : const TClonesArray * digits = emcalLoader->Digits() ;
1238 4 : TTree * treeD = emcalLoader->TreeD();
1239 4 : if ( !treeD ) {
1240 4 : emcalLoader->MakeDigitsContainer();
1241 4 : treeD = emcalLoader->TreeD();
1242 4 : }
1243 :
1244 : // -- create Digits branch
1245 : Int_t bufferSize = 32000 ;
1246 : TBranch * digitsBranch = 0;
1247 4 : if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
1248 0 : digitsBranch->SetAddress(&digits);
1249 0 : AliWarning("Digits Branch already exists. Not all digits will be visible");
1250 0 : }
1251 : else
1252 4 : treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
1253 : //digitsBranch->SetTitle(fEventFolderName);
1254 :
1255 : // treeD->Fill() ;
1256 : /*
1257 : emcalLoader->WriteDigits("OVERWRITE");
1258 : emcalLoader->WriteDigitizer("OVERWRITE");
1259 :
1260 : Unload() ;
1261 : */
1262 :
1263 4 : }// loader exists
1264 0 : else AliFatal("Loader not available");
1265 4 : }
1266 :
1267 : //__________________________________________________________________
1268 : void AliEMCALDigitizer::WriteDigits(TClonesArray* digits, const char* branchName)
1269 : {
1270 : // overloaded method
1271 12 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
1272 4 : if(emcalLoader){
1273 :
1274 4 : TTree* treeD = emcalLoader->TreeD();
1275 4 : if (!treeD)
1276 : {
1277 0 : emcalLoader->MakeDigitsContainer();
1278 0 : treeD = emcalLoader->TreeD();
1279 0 : }
1280 :
1281 : // -- create Digits branch
1282 : Int_t bufferSize = 32000;
1283 :
1284 4 : if (TBranch* triggerBranch = treeD->GetBranch(branchName))
1285 : {
1286 0 : triggerBranch->SetAddress(&digits);
1287 0 : }
1288 : else
1289 : {
1290 4 : treeD->Branch(branchName,"TClonesArray",&digits,bufferSize);
1291 : }
1292 :
1293 : // treeD->Fill();
1294 4 : }// loader exists
1295 0 : else AliFatal("Loader not available");
1296 4 : }
1297 :
1298 : //__________________________________________________________________
1299 : Bool_t AliEMCALDigitizer::IsDead(AliEMCALDigit *digit)
1300 : {
1301 : // Check if cell is defined as dead, so that it is not included
1302 : // input is digit
1303 536 : Int_t absId = digit->GetId();
1304 :
1305 268 : return IsDead(absId);
1306 :
1307 : }
1308 :
1309 :
1310 : //__________________________________________________________________
1311 : Bool_t AliEMCALDigitizer::IsDead(Int_t absId)
1312 : {
1313 : // Check if cell absID is defined as dead, so that it is not included
1314 :
1315 652 : AliRunLoader *runLoader = AliRunLoader::Instance();
1316 978 : AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
1317 326 : if (!emcalLoader)
1318 : {
1319 0 : AliFatal("Did not get the Loader");
1320 0 : return kTRUE;
1321 : }
1322 :
1323 326 : AliCaloCalibPedestal *caloPed = emcalLoader->PedestalData();
1324 326 : if (!caloPed)
1325 : {
1326 0 : AliWarning("Could not access pedestal data! No dead channel removal applied");
1327 0 : return kFALSE;
1328 : }
1329 :
1330 : // Load Geometry
1331 326 : const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
1332 326 : if (!geom)
1333 : {
1334 0 : AliFatal("Did not get geometry from EMCALLoader");
1335 0 : return kTRUE; //coverity
1336 : }
1337 :
1338 326 : Int_t iSupMod = -1;
1339 326 : Int_t nModule = -1;
1340 326 : Int_t nIphi = -1;
1341 326 : Int_t nIeta = -1;
1342 326 : Int_t iphi = -1;
1343 326 : Int_t ieta = -1;
1344 :
1345 326 : Bool_t bCell = geom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
1346 :
1347 326 : if (!bCell) Error("IsDead","Wrong cell id number : absId %i ", absId) ;
1348 326 : geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1349 :
1350 326 : Int_t channelStatus = (Int_t)(caloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
1351 :
1352 326 : if (channelStatus == AliCaloCalibPedestal::kDead)
1353 0 : return kTRUE;
1354 : else
1355 326 : return kFALSE;
1356 652 : }
1357 :
|
