Change digitization threshold to >=3 ACD units, now cut on ADC value of digit after...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //_________________________________________________________________________
19 // 
20 //////////////////////////////////////////////////////////////////////////////
21 // Class performs digitization of Summable digits from simulated data
22 //  
23 // In addition it performs mixing of summable digits from different events.
24 //
25 // For each event two branches are created in TreeD:
26 //   "EMCAL" - list of digits
27 //   "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
28 //
29 // Note, that one cset title for new digits branch, and repeat digitization with
30 // another set of parameters.
31 //
32 // Examples of use:
33 // root[0] AliEMCALDigitizer * d = new AliEMCALDigitizer() ;
34 // root[1] d->ExecuteTask()             
35 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
36 //                       //Digitizes SDigitis in all events found in file galice.root 
37 //
38 // root[2] AliEMCALDigitizer * d1 = new AliEMCALDigitizer("galice1.root") ;  
39 //                       // Will read sdigits from galice1.root
40 // root[3] d1->MixWith("galice2.root")       
41 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
42 //                       // Reads another portion of sdigits from galice2.root
43 // root[3] d1->MixWith("galice3.root")       
44 //                       // Reads another portion of sdigits from galice3.root
45 // root[4] d->ExecuteTask("deb timing")    
46 //                       // Reads SDigits from files galice1.root, galice2.root ....
47 //                       // mixes them and stores produced Digits in file galice1.root          
48 //                       // deb - prints number of produced digits
49 //                       // deb all - prints list of produced digits
50 //                       // timing  - prints time used for digitization
51 ////////////////////////////////////////////////////////////////////////////////////
52 //
53 //*-- Author: Sahal Yacoob (LBL)
54 // based on : AliEMCALDigitizer
55 // Modif: 
56 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
57 //                           of new  IO (à la PHOS)
58 //  November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry 
59 //_________________________________________________________________________________
60
61 // --- ROOT system ---
62 #include <TROOT.h>
63 #include <TTree.h>
64 #include <TSystem.h>
65 #include <TBenchmark.h>
66 #include <TBrowser.h>
67 #include <TObjectTable.h>
68 #include <TRandom.h>
69 #include <cassert>
70
71 // --- AliRoot header files ---
72 #include "AliLog.h"
73 #include "AliRun.h"
74 #include "AliRunDigitizer.h"
75 #include "AliRunLoader.h"
76 #include "AliCDBManager.h"
77 #include "AliCDBEntry.h"
78 #include "AliEMCALDigit.h"
79 #include "AliEMCAL.h"
80 #include "AliEMCALLoader.h"
81 #include "AliEMCALDigitizer.h"
82 #include "AliEMCALSDigitizer.h"
83 #include "AliEMCALGeometry.h"
84 #include "AliEMCALTick.h"
85 #include "AliEMCALCalibData.h"
86 #include "AliEMCALSimParam.h"
87
88 ClassImp(AliEMCALDigitizer)
89
90
91 //____________________________________________________________________________ 
92 AliEMCALDigitizer::AliEMCALDigitizer()
93   : AliDigitizer("",""),
94     fDefaultInit(kTRUE),
95     fDigitsInRun(0),
96     fInit(0),
97     fInput(0),
98     fInputFileNames(0x0),
99     fEventNames(0x0),
100     fDigitThreshold(0),
101     fMeanPhotonElectron(0),
102 //    fPedestal(0), //Not used, remove?
103 //    fSlope(0),    //Not used, remove?
104     fPinNoise(0),
105     fTimeResolution(0),
106 //    fTimeThreshold(0),    //Not used, remove?
107 //    fTimeSignalLength(0), //Not used, remove?
108     fADCchannelEC(0),
109     fADCpedestalEC(0),
110     fNADCEC(0),
111     fEventFolderName(""),
112     fFirstEvent(0),
113     fLastEvent(0),
114     fCalibData(0x0)
115 {
116   // ctor
117   InitParameters() ; 
118   fManager = 0 ;                     // We work in the standalone mode
119 }
120
121 //____________________________________________________________________________ 
122 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
123   : AliDigitizer("EMCAL"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
124     fDefaultInit(kFALSE),
125     fDigitsInRun(0),
126     fInit(0),
127     fInput(0),
128     fInputFileNames(0), 
129     fEventNames(0), 
130     fDigitThreshold(0),
131     fMeanPhotonElectron(0),
132 //    fPedestal(0),//Not used, remove?
133 //    fSlope(0),   //Not used, remove?
134     fPinNoise(0),
135     fTimeResolution(0),
136 //    fTimeThreshold(0),    //Not used, remove?
137 //    fTimeSignalLength(0), //Not used, remove?
138     fADCchannelEC(0),
139     fADCpedestalEC(0),
140     fNADCEC(0),
141     fEventFolderName(eventFolderName),
142     fFirstEvent(0),
143     fLastEvent(0),
144     fCalibData(0x0)
145 {
146   // ctor
147   InitParameters() ; 
148   Init() ;
149   fManager = 0 ;                     // We work in the standalone mode
150 }
151
152 //____________________________________________________________________________ 
153 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d) 
154   : AliDigitizer(d.GetName(),d.GetTitle()),
155     fDefaultInit(d.fDefaultInit),
156     fDigitsInRun(d.fDigitsInRun),
157     fInit(d.fInit),
158     fInput(d.fInput),
159     fInputFileNames(d.fInputFileNames),
160     fEventNames(d.fEventNames),
161     fDigitThreshold(d.fDigitThreshold),
162     fMeanPhotonElectron(d.fMeanPhotonElectron),
163 //    fPedestal(d.fPedestal), //Not used, remove?
164 //    fSlope(d.fSlope),       //Not used, remove?
165     fPinNoise(d.fPinNoise),
166     fTimeResolution(d.fTimeResolution),
167 //    fTimeThreshold(d.fTimeThreshold),       //Not used, remove?
168 //    fTimeSignalLength(d.fTimeSignalLength), //Not used, remove?
169     fADCchannelEC(d.fADCchannelEC),
170     fADCpedestalEC(d.fADCpedestalEC),
171     fNADCEC(d.fNADCEC),
172     fEventFolderName(d.fEventFolderName),
173     fFirstEvent(d.fFirstEvent),
174     fLastEvent(d.fLastEvent),
175     fCalibData(d.fCalibData)
176 {
177   // copyy ctor 
178  }
179
180 //____________________________________________________________________________ 
181 AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * rd)
182   : AliDigitizer(rd,"EMCAL"+AliConfig::Instance()->GetDigitizerTaskName()),
183     fDefaultInit(kFALSE),
184     fDigitsInRun(0),
185     fInit(0),
186     fInput(0),
187     fInputFileNames(0),
188     fEventNames(0),
189     fDigitThreshold(0.),
190     fMeanPhotonElectron(0),
191 //    fPedestal(0), //Not used, remove?
192 //    fSlope(0.),   //Not used, remove?
193     fPinNoise(0),
194     fTimeResolution(0.),
195 //    fTimeThreshold(0),    //Not used, remove?
196 //    fTimeSignalLength(0), //Not used, remove?
197     fADCchannelEC(0),
198     fADCpedestalEC(0),
199     fNADCEC(0),
200     fEventFolderName(0),
201     fFirstEvent(0),
202     fLastEvent(0),
203     fCalibData(0x0)
204 {
205   // ctor Init() is called by RunDigitizer
206   fManager = rd ; 
207   fEventFolderName = fManager->GetInputFolderName(0) ;
208   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
209   InitParameters() ; 
210 }
211
212 //____________________________________________________________________________ 
213   AliEMCALDigitizer::~AliEMCALDigitizer()
214 {
215   //dtor
216   if (AliRunLoader::Instance()) {
217     AliLoader *emcalLoader=0;
218     if ((emcalLoader = AliRunLoader::Instance()->GetDetectorLoader("EMCAL")))
219       emcalLoader->CleanDigitizer();
220   }
221   else
222     AliDebug(1," no runloader present");
223   delete [] fInputFileNames ; 
224   delete [] fEventNames ; 
225
226 }
227
228 //____________________________________________________________________________
229 void AliEMCALDigitizer::Digitize(Int_t event) 
230
231
232   // Makes the digitization of the collected summable digits
233   // for this it first creates the array of all EMCAL modules
234   // filled with noise and after that adds contributions from 
235   // SDigits. This design helps to avoid scanning over the 
236   // list of digits to add  contribution of any new SDigit.
237   //
238   // JLK 26-Jun-2008
239   // Note that SDigit energy info is stored as an amplitude, so we
240   // must call the Calibrate() method of the SDigitizer to convert it
241   // back to an energy in GeV before adding it to the Digit
242   //
243   static int nEMC=0; //max number of digits possible
244
245   AliRunLoader *rl = AliRunLoader::Instance();
246   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
247   Int_t readEvent = event ; 
248   // fManager is data member from AliDigitizer
249   if (fManager) 
250     readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ; 
251   AliDebug(1,Form("Adding event %d from input stream 0 %s %s", 
252                   readEvent, GetTitle(), fEventFolderName.Data())) ; 
253   rl->GetEvent(readEvent);
254
255   TClonesArray * digits = emcalLoader->Digits() ; 
256   digits->Delete() ;  //correct way to clear array when memory is
257                       //allocated by objects stored in array
258
259   // Load Geometry
260   AliEMCALGeometry *geom = 0;
261   if (rl->GetAliRun()) {
262     AliEMCAL * emcal  = (AliEMCAL*)rl->GetAliRun()->GetDetector("EMCAL");
263     geom = emcal->GetGeometry();
264   }
265   else 
266     AliFatal("Could not get AliRun from runLoader");
267
268   nEMC = geom->GetNCells();
269   AliDebug(1,Form("nEMC %i (number cells in EMCAL) | %s \n", nEMC, geom->GetName()));
270   
271   Int_t absID ;
272
273   digits->Expand(nEMC) ;
274
275   // get first the sdigitizer from the tasks list (must have same name as the digitizer)
276   if ( !emcalLoader->SDigitizer() ) 
277     emcalLoader->LoadSDigitizer();
278   AliEMCALSDigitizer * sDigitizer = dynamic_cast<AliEMCALSDigitizer*>(emcalLoader->SDigitizer()); 
279   
280   if ( !sDigitizer )
281     Fatal("Digitize", "SDigitizer with name %s %s not found", fEventFolderName.Data(), GetTitle() ) ; 
282
283   //take all the inputs to add together and load the SDigits
284   TObjArray * sdigArray = new TObjArray(fInput) ;
285   sdigArray->AddAt(emcalLoader->SDigits(), 0) ;
286   Int_t i ;
287
288   for(i = 1 ; i < fInput ; i++){
289     TString tempo(fEventNames[i]) ; 
290     tempo += i ;
291
292     AliRunLoader *  rl2 = AliRunLoader::GetRunLoader(tempo) ; 
293
294     if (rl2==0) 
295       rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ; 
296
297     if (fManager) 
298       readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ; 
299     Info("Digitize", "Adding event %d from input stream %d %s %s", readEvent, i, fInputFileNames[i].Data(), tempo.Data()) ; 
300     rl2->LoadSDigits();
301     rl2->GetEvent(readEvent);
302     AliEMCALLoader *emcalLoader2 = dynamic_cast<AliEMCALLoader*>(rl2->GetDetectorLoader("EMCAL"));
303     sdigArray->AddAt(emcalLoader2->SDigits(), i) ;
304   }
305   
306   //Find the first tower with signal
307   Int_t nextSig = nEMC + 1 ; 
308   TClonesArray * sdigits ;  
309   for(i = 0 ; i < fInput ; i++){
310     sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
311     if ( !sdigits->GetEntriesFast() )
312       continue ; 
313     Int_t curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(0))->GetId() ;
314      if(curNext < nextSig) 
315        nextSig = curNext ;
316      AliDebug(1,Form("input %i : #sdigits %i \n",
317                      i, sdigits->GetEntriesFast()));
318   }
319   AliDebug(1,Form("FIRST tower with signal %i \n", nextSig));
320
321   TArrayI index(fInput) ;
322   index.Reset() ;  //Set all indexes to zero
323
324   AliEMCALDigit * digit ;
325   AliEMCALDigit * curSDigit ;
326
327   //  TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
328
329   //Put Noise contribution
330   for(absID = 0; absID < nEMC; absID++){ // Nov 30, 2006 by PAI; was from 1 to nEMC
331     Float_t energy = 0 ;
332     // amplitude set to zero, noise will be added later
333     new((*digits)[absID]) AliEMCALDigit( -1, -1, absID, 0, TimeOfNoise() ); // absID-1->absID
334     //look if we have to add signal?
335     digit = dynamic_cast<AliEMCALDigit *>(digits->At(absID)); // absID-1->absID
336     
337     if(absID==nextSig){
338       //Add SDigits from all inputs    
339       //      ticks->Clear() ;
340       //Int_t contrib = 0 ;
341
342       //Follow PHOS and comment out this timing model til a better one
343       //can be developed - JLK 28-Apr-2008
344
345       //Float_t a = digit->GetAmp() ;
346       //Float_t b = TMath::Abs( a /fTimeSignalLength) ;
347       //Mark the beginning of the signal
348       //new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);  
349       //Mark the end of the signal     
350       //new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
351
352       // Calculate time as time of the largest digit
353       Float_t time = digit->GetTime() ;
354       Float_t aTime= digit->GetAmp() ;
355       
356       // loop over input
357       for(i = 0; i< fInput ; i++){  //loop over (possible) merge sources
358         if(dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
359           curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;      
360         else
361           curSDigit = 0 ;
362         //May be several digits will contribute from the same input
363         while(curSDigit && (curSDigit->GetId() == absID)){         
364           //Shift primary to separate primaries belonging different inputs
365           Int_t primaryoffset ;
366           if(fManager)
367             primaryoffset = fManager->GetMask(i) ; 
368           else
369             primaryoffset = i ;
370           curSDigit->ShiftPrimary(primaryoffset) ;
371
372           //Remove old timing model - JLK 28-April-2008
373           //a = curSDigit->GetAmp() ;
374           //b = a /fTimeSignalLength ;
375           //new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);  
376           //new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
377           if(curSDigit->GetAmp()>aTime) {
378             aTime = curSDigit->GetAmp();
379             time = curSDigit->GetTime();
380           }
381
382           *digit = *digit + *curSDigit ;  //adds amplitudes of each digit
383
384           index[i]++ ;
385           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
386             curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
387           else
388             curSDigit = 0 ;
389         }
390       }
391       //Here we convert the summed amplitude to an energy in GeV
392       energy = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
393       // add fluctuations for photo-electron creation
394       energy *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
395   
396       //calculate and set time
397       //New timing model needed - JLK 28-April-2008
398       //Float_t time = FrontEdgeTime(ticks) ;
399       digit->SetTime(time) ;
400
401       //Find next signal module
402       nextSig = nEMC + 1 ;
403       for(i = 0 ; i < fInput ; i++){
404         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
405         Int_t curNext = nextSig ;
406         if(sdigits->GetEntriesFast() > index[i] ){
407           curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]))->GetId() ;       
408         }
409         if(curNext < nextSig) nextSig = curNext ;
410       }
411     }
412     // add the noise now
413     energy += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
414     // JLK 26-June-2008
415     //Now digitize the energy according to the sDigitizer method,
416     //which merely converts the energy to an integer.  Later we will
417     //check that the stored value matches our allowed dynamic ranges
418     digit->SetAmp(sDigitizer->Digitize(energy)) ;  
419     AliDebug(10,Form(" absID %5i energy %f nextSig %5i\n",
420                      absID, energy, nextSig));
421   } // for(absID = 0; absID < nEMC; absID++)
422   
423   //ticks->Delete() ;
424   //delete ticks ;
425
426   //JLK is it better to call Clear() here?
427   delete sdigArray ; //We should not delete its contents
428
429   //remove digits below thresholds
430   // until 10-02-2010 remove digits with energy smaller than fDigitThreshold 3*fPinNoise
431   // now, remove digits with Digitized ADC smaller than fDigitThreshold = 3
432   Float_t energy=0;
433   for(i = 0 ; i < nEMC ; i++){
434     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
435     //First get the energy in GeV units.
436     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
437     //Then digitize using the calibration constants of the ocdb
438     Int_t ampADC = DigitizeEnergy(energy, digit->GetId())  ;      
439     //if(ampADC>2)printf("Digit energy %f, amp %d, amp cal %d, threshold %d\n",energy,digit->GetAmp(),ampADC,fDigitThreshold);
440     if(ampADC < fDigitThreshold)
441       digits->RemoveAt(i) ;
442     else 
443       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
444   }
445   
446   digits->Compress() ;  
447   
448   Int_t ndigits = digits->GetEntriesFast() ; 
449   
450   //JLK 26-June-2008
451   //After we have done the summing and digitizing to create the
452   //digits, now we want to calibrate the resulting amplitude to match
453   //the dynamic range of our real data.  
454   for (i = 0 ; i < ndigits ; i++) { 
455     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ; 
456     digit->SetIndexInList(i) ; 
457     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
458     digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ;
459   }
460
461 }
462
463 // //_____________________________________________________________________
464 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
465
466   // JLK 26-June-2008
467   // Returns digitized value of the energy in a cell absId
468   // using the calibration constants stored in the OCDB
469   // or default values if no CalibData object is found.
470   // This effectively converts everything to match the dynamic range
471   // of the real data we will collect
472   //
473   // Load Geometry
474   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
475
476   if (geom==0)
477     AliFatal("Did not get geometry from EMCALLoader");
478
479   Int_t iSupMod = -1;
480   Int_t nModule  = -1;
481   Int_t nIphi   = -1;
482   Int_t nIeta   = -1;
483   Int_t iphi    = -1;
484   Int_t ieta    = -1;
485   Int_t channel = -999; 
486
487   Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
488   if(!bCell)
489     Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
490   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
491   
492   if(fCalibData) {
493     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
494     fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
495   }
496   
497   channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
498   
499   if(channel > fNADCEC ) 
500     channel =  fNADCEC ; 
501   return channel ;
502   
503 }
504
505 //____________________________________________________________________________
506 void AliEMCALDigitizer::Exec(Option_t *option) 
507
508   // Steering method to process digitization for events
509   // in the range from fFirstEvent to fLastEvent.
510   // This range is optionally set by SetEventRange().
511   // if fLastEvent=-1, then process events until the end.
512   // by default fLastEvent = fFirstEvent (process only one event)
513
514   if (!fInit) { // to prevent overwrite existing file
515     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
516     return ;
517   } 
518
519   if (strstr(option,"print")) {
520
521     Print();
522     return ; 
523   }
524   
525   if(strstr(option,"tim"))
526     gBenchmark->Start("EMCALDigitizer");
527
528   AliRunLoader *rl = AliRunLoader::Instance();
529   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
530
531   // Post Digitizer to the white board
532   emcalLoader->PostDigitizer(this) ;
533   
534   if (fLastEvent == -1)  {
535     fLastEvent = rl->GetNumberOfEvents() - 1 ;
536   }
537   else if (fManager) 
538     fLastEvent = fFirstEvent ; // what is this ??
539
540   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
541   Int_t ievent;
542
543   rl->LoadSDigits("EMCAL");
544   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
545     
546     rl->GetEvent(ievent);
547
548     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
549
550     WriteDigits() ;
551
552     if(strstr(option,"deb"))
553       PrintDigits(option);
554     if(strstr(option,"table")) gObjectTable->Print();
555
556     //increment the total number of Digits per run 
557     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;  
558   }
559   
560   emcalLoader->CleanDigitizer() ;
561
562   if(strstr(option,"tim")){
563     gBenchmark->Stop("EMCALDigitizer");
564     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event", 
565          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
566   } 
567 }
568
569 //____________________________________________________________________________ 
570 //Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks) 
571 //{ 
572 //  //  Returns the shortest time among all time ticks
573 //
574 //  ticks->Sort() ; //Sort in accordance with times of ticks
575 //  TIter it(ticks) ;
576 //  AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
577 //  Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
578 //  
579 //  AliEMCALTick * t ;  
580 //  while((t=(AliEMCALTick*) it.Next())){
581 //    if(t->GetTime() < time)  //This tick starts before crossing
582 //      *ctick+=*t ;
583 //    else
584 //      return time ;
585 //    
586 //    time = ctick->CrossingTime(fTimeThreshold) ;    
587 //  }
588 //  return time ;
589 //}
590 //
591
592 //____________________________________________________________________________ 
593 Bool_t AliEMCALDigitizer::Init()
594 {
595   // Makes all memory allocations
596   fInit = kTRUE ; 
597   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
598
599   if ( emcalLoader == 0 ) {
600     Fatal("Init", "Could not obtain the AliEMCALLoader");  
601     return kFALSE;
602   } 
603
604   fFirstEvent = 0 ; 
605   fLastEvent = fFirstEvent ; 
606   
607   if(fManager)
608     fInput = fManager->GetNinputs() ; 
609   else 
610     fInput           = 1 ;
611
612   fInputFileNames  = new TString[fInput] ; 
613   fEventNames      = new TString[fInput] ; 
614   fInputFileNames[0] = GetTitle() ; 
615   fEventNames[0]     = fEventFolderName.Data() ; 
616   Int_t index ; 
617   for (index = 1 ; index < fInput ; index++) {
618     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0); 
619     TString tempo = fManager->GetInputFolderName(index) ;
620     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager 
621   }
622   
623   //to prevent cleaning of this object while GetEvent is called
624   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
625
626   //Calibration instance
627   fCalibData = emcalLoader->CalibData();
628   return fInit ;    
629 }
630
631 //____________________________________________________________________________ 
632 void AliEMCALDigitizer::InitParameters()
633
634   // Parameter initialization for digitizer
635  
636   fMeanPhotonElectron = AliEMCALSimParam::GetInstance()->GetMeanPhotonElectron();//4400;  // electrons per GeV 
637   fPinNoise           = AliEMCALSimParam::GetInstance()->GetPinNoise();//0.012; // pin noise in GeV from analysis test beam data 
638   if (fPinNoise < 0.0001 ) 
639     Warning("InitParameters", "No noise added\n") ; 
640   fDigitThreshold     = AliEMCALSimParam::GetInstance()->GetDigitThreshold(); //fPinNoise * 3; // 3 * sigma
641   fTimeResolution     = AliEMCALSimParam::GetInstance()->GetTimeResolution(); //0.6e-9 ; // 600 psc
642
643   // These defaults are normally not used. 
644   // Values are read from calibration database instead
645   fADCchannelEC       = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
646   fADCpedestalEC      = 0.0 ;  // GeV
647
648   fNADCEC          = AliEMCALSimParam::GetInstance()->GetNADCEC();//(Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
649
650   AliDebug(2,Form("Mean Photon Electron %d, Noise %f, E Threshold %f,Time Resolution %g,NADCEC %d",
651                 fMeanPhotonElectron,fPinNoise,fDigitThreshold,fTimeResolution,fNADCEC));
652
653   // Not used anymore, remove?
654   // fTimeSignalLength   = 1.0e-9 ;
655   // fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
656
657 }
658
659 //__________________________________________________________________
660 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
661 {
662   // Allows to produce digits by superimposing background and signal event.
663   // It is assumed, that headers file with SIGNAL events is opened in 
664   // the constructor. 
665   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed 
666   // Thus we avoid writing (changing) huge and expensive 
667   // backgound files: all output will be writen into SIGNAL, i.e. 
668   // opened in constructor file. 
669   //
670   // One can open as many files to mix with as one needs.
671   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
672   // can be mixed.
673
674   if( strcmp(GetName(), "") == 0 )
675     Init() ;
676   
677   if(fManager){
678     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
679     return ;
680   } 
681   // looking for file which contains AliRun
682   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
683     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
684     return ; 
685   }
686   // looking for the file which contains SDigits
687   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
688   TString fileName( emcalLoader->GetSDigitsFileName() ) ; 
689     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
690       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
691     if ( (gSystem->AccessPathName(fileName)) ) { 
692       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
693       return ;
694     }
695     // need to increase the arrays
696     // MvL: This code only works when fInput == 1, I think.
697     TString tempo = fInputFileNames[fInput-1] ; 
698     delete [] fInputFileNames ; 
699     fInputFileNames = new TString[fInput+1] ; 
700     fInputFileNames[fInput-1] = tempo ; 
701  
702     tempo = fEventNames[fInput-1] ; 
703     delete [] fEventNames ; 
704     fEventNames = new TString[fInput+1] ; 
705     fEventNames[fInput-1] = tempo ; 
706
707     fInputFileNames[fInput] = alirunFileName ; 
708     fEventNames[fInput]     = eventFolderName ;
709     fInput++ ;
710 }  
711
712 //__________________________________________________________________
713 void AliEMCALDigitizer::Print1(Option_t * option)
714 { // 19-nov-04 - just for convinience
715   Print(); 
716   PrintDigits(option);
717 }
718
719 //__________________________________________________________________
720 void AliEMCALDigitizer::Print(Option_t*)const 
721 {
722   // Print Digitizer's parameters
723   printf("Print: \n------------------- %s -------------", GetName() ) ; 
724   if( strcmp(fEventFolderName.Data(), "") != 0 ){
725     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
726     
727     Int_t nStreams ; 
728     if (fManager) 
729       nStreams =  GetNInputStreams() ;
730     else 
731       nStreams = fInput ; 
732     
733     Int_t index = 0 ;  
734
735     AliRunLoader *rl=0;
736
737     for (index = 0 ; index < nStreams ; index++) {  
738       TString tempo(fEventNames[index]) ; 
739       tempo += index ;
740       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
741         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ; 
742       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
743       TString fileName( emcalLoader->GetSDigitsFileName() ) ; 
744       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
745         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
746       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ; 
747     }
748
749     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
750
751     printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
752     
753     printf("\nWith following parameters:\n") ;
754     
755     printf("    Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
756     printf("    Threshold  in Tower  (fDigitThreshold) = %d\n", fDigitThreshold)  ;
757     printf("---------------------------------------------------\n")  ;
758   }
759   else
760     printf("Print: AliEMCALDigitizer not initialized") ; 
761 }
762
763 //__________________________________________________________________
764 void AliEMCALDigitizer::PrintDigits(Option_t * option)
765 {
766   //utility method for printing digit information
767
768   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
769   TClonesArray * digits  = emcalLoader->Digits() ;
770   TClonesArray * sdigits = emcalLoader->SDigits() ;
771   
772   printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ; 
773   printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
774   
775   if(strstr(option,"all")){  
776     //loop over digits
777     AliEMCALDigit * digit;
778     printf("\nEMC digits (with primaries):\n")  ;
779     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
780     Int_t index ;
781     for (index = 0 ; index < digits->GetEntries() ; index++) {
782       digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
783       printf("\n%6d  %8d    %6.5e %4d      %2d : ",
784               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
785       Int_t iprimary;
786       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
787         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
788       }
789     }   
790   }
791   printf("\n");
792 }
793
794 //__________________________________________________________________
795 Float_t AliEMCALDigitizer::TimeOfNoise(void)
796 {  
797   // Calculates the time signal generated by noise
798   //PH  Info("TimeOfNoise", "Change me") ; 
799   return gRandom->Rndm() * 1.28E-5;
800 }
801
802 //__________________________________________________________________
803 void AliEMCALDigitizer::Unload() 
804 {  
805   // Unloads the SDigits and Digits
806   AliRunLoader *rl=0;
807     
808   Int_t i ; 
809   for(i = 1 ; i < fInput ; i++){
810     TString tempo(fEventNames[i]) ; 
811     tempo += i;
812     if ((rl = AliRunLoader::GetRunLoader(tempo))) 
813       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ; 
814   }
815   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
816   emcalLoader->UnloadDigits() ; 
817 }
818
819 //_________________________________________________________________________________________
820 void AliEMCALDigitizer::WriteDigits()
821 {
822
823   // Makes TreeD in the output file. 
824   // Check if branch already exists: 
825   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
826   //      already existing branches. 
827   //   else creates branch with Digits, named "EMCAL", title "...",
828   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
829   //      and names of files, from which digits are made.
830
831   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
832
833   const TClonesArray * digits = emcalLoader->Digits() ; 
834   TTree * treeD = emcalLoader->TreeD(); 
835   if ( !treeD ) {
836     emcalLoader->MakeDigitsContainer();
837     treeD = emcalLoader->TreeD(); 
838   }
839
840   // -- create Digits branch
841   Int_t bufferSize = 32000 ;    
842   TBranch * digitsBranch = 0;
843   if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
844     digitsBranch->SetAddress(&digits);
845     AliWarning("Digits Branch already exists. Not all digits will be visible");
846   }
847   else
848     treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
849   //digitsBranch->SetTitle(fEventFolderName);
850   treeD->Fill() ;
851   
852   emcalLoader->WriteDigits("OVERWRITE");
853   emcalLoader->WriteDigitizer("OVERWRITE");
854
855   Unload() ; 
856
857 }
858
859 //__________________________________________________________________
860 void AliEMCALDigitizer::Browse(TBrowser* b)
861 {
862   TTask::Browse(b);
863 }