Bug fix in the constructor (thanks to A.Marin)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
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
108 ClassImp(AliEMCALDigitizer)
109   
110   
111 //____________________________________________________________________________ 
112 AliEMCALDigitizer::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 //____________________________________________________________________________ 
144 AliEMCALDigitizer::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 //____________________________________________________________________________ 
177 AliEMCALDigitizer::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 //____________________________________________________________________________ 
206 AliEMCALDigitizer::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 //____________________________________________________________________________
257 void 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 //_____________________________________________________________________
619 void 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 //_____________________________________________________________________
664 void 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 //____________________________________________________________________________
707 void 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 //__________________________________________________________________
800 Float_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 //____________________________________________________________________________ 
813 void 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 //____________________________________________________________________________ 
906 void 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 //____________________________________________________________________________ 
935 Bool_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 //____________________________________________________________________________ 
974 void 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 //__________________________________________________________________
1015 void AliEMCALDigitizer::Print1(Option_t * option)
1016 { // 19-nov-04 - just for convenience
1017   Print(); 
1018   PrintDigits(option);
1019 }
1020
1021 //__________________________________________________________________
1022 void 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 //__________________________________________________________________
1067 void 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 //__________________________________________________________________
1103 Float_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 //__________________________________________________________________
1111 void 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 //_________________________________________________________________________________________
1129 void 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 //__________________________________________________________________
1172 void 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