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