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