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