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