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