code updates suggested by Federico to improve memory management
[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 
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 <TList.h>
67 #include <TH1.h>
68 #include <TBrowser.h>
69 #include <TObjectTable.h>
70 #include <TRandom.h>
71 #include <cassert>
72
73 // --- AliRoot header files ---
74 #include "AliLog.h"
75 #include "AliRun.h"
76 #include "AliRunDigitizer.h"
77 #include "AliRunLoader.h"
78 #include "AliCDBManager.h"
79 #include "AliCDBEntry.h"
80 #include "AliEMCALDigit.h"
81 #include "AliEMCAL.h"
82 #include "AliEMCALLoader.h"
83 #include "AliEMCALDigitizer.h"
84 #include "AliEMCALSDigitizer.h"
85 #include "AliEMCALGeometry.h"
86 #include "AliEMCALTick.h"
87 #include "AliEMCALHistoUtilities.h"
88
89 ClassImp(AliEMCALDigitizer)
90
91
92 //____________________________________________________________________________ 
93 AliEMCALDigitizer::AliEMCALDigitizer()
94   : AliDigitizer("",""),
95     fDefaultInit(kTRUE),
96     fDigitsInRun(0),
97     fInit(0),
98     fInput(0),
99     fInputFileNames(0x0),
100     fEventNames(0x0),
101     fDigitThreshold(0),
102     fMeanPhotonElectron(0),
103     fPedestal(0),
104     fSlope(0),
105     fPinNoise(0),
106     fTimeResolution(0),
107     fTimeThreshold(0),    
108     fTimeSignalLength(0),
109     fADCchannelEC(0),
110     fADCpedestalEC(0),
111     fNADCEC(0),
112     fEventFolderName(""),
113     fFirstEvent(0),
114     fLastEvent(0),
115     fControlHists(0),
116     fHists(0),fCalibData(0x0)
117 {
118   // ctor
119   InitParameters() ; 
120   fManager = 0 ;                     // We work in the standalong mode
121 }
122
123 //____________________________________________________________________________ 
124 AliEMCALDigitizer::AliEMCALDigitizer(TString alirunFileName, TString eventFolderName)
125   : AliDigitizer("EMCAL"+AliConfig::Instance()->GetDigitizerTaskName(), alirunFileName),
126     fDefaultInit(kFALSE),
127     fDigitsInRun(0),
128     fInit(0),
129     fInput(0),
130     fInputFileNames(0), 
131     fEventNames(0), 
132     fDigitThreshold(0),
133     fMeanPhotonElectron(0),
134     fPedestal(0),
135     fSlope(0),
136     fPinNoise(0),
137     fTimeResolution(0),
138     fTimeThreshold(0),
139     fTimeSignalLength(0),
140     fADCchannelEC(0),
141     fADCpedestalEC(0),
142     fNADCEC(0),
143     fEventFolderName(eventFolderName),
144     fFirstEvent(0),
145     fLastEvent(0),
146     fControlHists(0),
147     fHists(0),fCalibData(0x0)
148 {
149   // ctor
150   InitParameters() ; 
151   Init() ;
152   fManager = 0 ;                     // We work in the standalong mode
153 }
154
155 //____________________________________________________________________________ 
156 AliEMCALDigitizer::AliEMCALDigitizer(const AliEMCALDigitizer & d) 
157   : AliDigitizer(d.GetName(),d.GetTitle()),
158     fDefaultInit(d.fDefaultInit),
159     fDigitsInRun(d.fDigitsInRun),
160     fInit(d.fInit),
161     fInput(d.fInput),
162     fInputFileNames(d.fInputFileNames),
163     fEventNames(d.fEventNames),
164     fDigitThreshold(d.fDigitThreshold),
165     fMeanPhotonElectron(d.fMeanPhotonElectron),
166     fPedestal(d.fPedestal),
167     fSlope(d.fSlope),
168     fPinNoise(d.fPinNoise),
169     fTimeResolution(d.fTimeResolution),
170     fTimeThreshold(d.fTimeThreshold),
171     fTimeSignalLength(d.fTimeSignalLength),
172     fADCchannelEC(d.fADCchannelEC),
173     fADCpedestalEC(d.fADCpedestalEC),
174     fNADCEC(d.fNADCEC),
175     fEventFolderName(d.fEventFolderName),
176     fFirstEvent(d.fFirstEvent),
177     fLastEvent(d.fLastEvent),
178     fControlHists(d.fControlHists),
179     fHists(d.fHists),fCalibData(d.fCalibData)
180 {
181   // copyy ctor 
182  }
183
184 //____________________________________________________________________________ 
185 AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * rd)
186   : AliDigitizer(rd,"EMCAL"+AliConfig::Instance()->GetDigitizerTaskName()),
187     fDefaultInit(kFALSE),
188     fDigitsInRun(0),
189     fInit(0),
190     fInput(0),
191     fInputFileNames(0),
192     fEventNames(0),
193     fDigitThreshold(0.),
194     fMeanPhotonElectron(0),
195     fPedestal(0),
196     fSlope(0.),
197     fPinNoise(0),
198     fTimeResolution(0.),
199     fTimeThreshold(0),
200     fTimeSignalLength(0),
201     fADCchannelEC(0),
202     fADCpedestalEC(0),
203     fNADCEC(0),
204     fEventFolderName(0),
205     fFirstEvent(0),
206     fLastEvent(0),
207     fControlHists(0),
208     fHists(0),fCalibData(0x0)
209 {
210   // ctor Init() is called by RunDigitizer
211   fManager = rd ; 
212   fEventFolderName = fManager->GetInputFolderName(0) ;
213   SetTitle(dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
214   InitParameters() ; 
215 }
216
217 //____________________________________________________________________________ 
218   AliEMCALDigitizer::~AliEMCALDigitizer()
219 {
220   //dtor
221   if (AliRunLoader::GetRunLoader()) {
222     AliLoader *emcalLoader=0;
223     if ((emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")))
224       emcalLoader->CleanDigitizer();
225   }
226   else
227     AliDebug(1," no runloader present");
228   delete [] fInputFileNames ; 
229   delete [] fEventNames ; 
230
231   if(fHists) delete fHists;
232 }
233
234 //____________________________________________________________________________
235 void AliEMCALDigitizer::Digitize(Int_t event) 
236
237
238   // Makes the digitization of the collected summable digits
239   // for this it first creates the array of all EMCAL modules
240   // filled with noise (different for EMC, CPV and PPSD) and
241   // after that adds contributions from SDigits. This design 
242   // helps to avoid scanning over the list of digits to add 
243   // contribution of any new SDigit.
244   static int nEMC=0; //max number of digits possible
245
246   AliRunLoader *rl = AliRunLoader::GetRunLoader();
247   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
248   Int_t readEvent = event ; 
249   // fManager is data member from AliDigitizer
250   if (fManager) 
251     readEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ; 
252   AliDebug(1,Form("Adding event %d from input stream 0 %s %s", 
253                   readEvent, GetTitle(), fEventFolderName.Data())) ; 
254   rl->GetEvent(readEvent);
255
256   TClonesArray * digits = emcalLoader->Digits() ; 
257   digits->Delete() ;  //JLK why is this created then deleted?
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 amp = 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 eTime= 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()>eTime) {
378             eTime = curSDigit->GetAmp();
379             time = curSDigit->GetTime();
380           }
381
382           *digit = *digit + *curSDigit ;  //add energies
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       // add fluctuations for photo-electron creation
392       amp = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
393       amp *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
394   
395       //calculate and set time
396       //New timing model needed - JLK 28-April-2008
397       //Float_t time = FrontEdgeTime(ticks) ;
398       digit->SetTime(time) ;
399
400       //Find next signal module
401       nextSig = nEMC + 1 ;
402       for(i = 0 ; i < fInput ; i++){
403         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
404         Int_t curNext = nextSig ;
405         if(sdigits->GetEntriesFast() > index[i] ){
406           curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]))->GetId() ;       
407         }
408         if(curNext < nextSig) nextSig = curNext ;
409       }
410     }
411     // add the noise now
412     amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
413     digit->SetAmp(sDigitizer->Digitize(amp)) ;  
414     AliDebug(10,Form(" absID %5i amp %f nextSig %5i\n",
415                      absID, amp, nextSig));
416   } // for(absID = 1; absID <= nEMC; absID++)
417   
418   //ticks->Delete() ;
419   //delete ticks ;
420
421   delete sdigArray ; //We should not delete its contents
422
423   //remove digits below thresholds
424   for(i = 0 ; i < nEMC ; i++){
425     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
426     Float_t threshold = fDigitThreshold ; 
427     if(sDigitizer->Calibrate( digit->GetAmp() ) < threshold)
428       digits->RemoveAt(i) ;
429     else 
430       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
431   }
432   
433   digits->Compress() ;  
434   
435   Int_t ndigits = digits->GetEntriesFast() ; 
436   
437   //Set indexes in list of digits and fill hists.
438   AliEMCALHistoUtilities::FillH1(fHists, 0, Double_t(ndigits));
439   Float_t energy=0., esum=0.;
440   for (i = 0 ; i < ndigits ; i++) { 
441     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ; 
442     digit->SetIndexInList(i) ; 
443     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
444     esum += energy;
445     digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ; // for what ??
446     AliEMCALHistoUtilities::FillH1(fHists, 2, double(digit->GetAmp()));
447     AliEMCALHistoUtilities::FillH1(fHists, 3, double(energy));
448     AliEMCALHistoUtilities::FillH1(fHists, 4, double(digit->GetId()));
449     //    if(digit->GetId() == nEMC) {
450     //  printf(" i %i \n", i );
451     //  digit->Dump();
452     //  assert(0);
453     //}
454   }
455   AliEMCALHistoUtilities::FillH1(fHists, 1, esum);
456 }
457
458 // //_____________________________________________________________________
459 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
460
461   // Returns digitized value of the energy in a cell absId
462
463   // Load Geometry
464   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
465
466   if (geom==0)
467     AliFatal("Did not get geometry from EMCALLoader");
468
469   Int_t iSupMod = -1;
470   Int_t nModule  = -1;
471   Int_t nIphi   = -1;
472   Int_t nIeta   = -1;
473   Int_t iphi    = -1;
474   Int_t ieta    = -1;
475   Int_t channel = -999; 
476
477   Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
478   if(!bCell)
479     Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
480   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
481   
482   if(fCalibData) {
483     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
484     fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
485   }
486   
487   channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
488   
489   if(channel > fNADCEC ) 
490     channel =  fNADCEC ; 
491   return channel ;
492   
493 }
494
495 //____________________________________________________________________________
496 void AliEMCALDigitizer::Exec(Option_t *option) 
497
498   // Steering method to process digitization for events
499   // in the range from fFirstEvent to fLastEvent.
500   // This range is optionally set by SetEventRange().
501   // if fLastEvent=-1, then process events until the end.
502   // by default fLastEvent = fFirstEvent (process only one event)
503
504   if (!fInit) { // to prevent overwrite existing file
505     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
506     return ;
507   } 
508
509   if (strstr(option,"print")) {
510
511     Print();
512     return ; 
513   }
514   
515   if(strstr(option,"tim"))
516     gBenchmark->Start("EMCALDigitizer");
517
518   AliRunLoader *rl = AliRunLoader::GetRunLoader();
519   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
520
521   // Post Digitizer to the white board
522   emcalLoader->PostDigitizer(this) ;
523   
524   if (fLastEvent == -1)  {
525     fLastEvent = rl->GetNumberOfEvents() - 1 ;
526   }
527   else if (fManager) 
528     fLastEvent = fFirstEvent ; // what is this ??
529
530   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
531   Int_t ievent;
532
533   rl->LoadSDigits("EMCAL");
534   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
535     
536     rl->GetEvent(ievent);
537
538     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
539
540     WriteDigits() ;
541
542     if(strstr(option,"deb"))
543       PrintDigits(option);
544     if(strstr(option,"table")) gObjectTable->Print();
545
546     //increment the total number of Digits per run 
547     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;  
548   }
549   
550   emcalLoader->CleanDigitizer() ;
551
552   if(strstr(option,"tim")){
553     gBenchmark->Stop("EMCALDigitizer");
554     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event", 
555          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
556   } 
557 }
558
559 //____________________________________________________________________________ 
560 Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks) 
561
562   //  Returns the shortest time among all time ticks
563
564   ticks->Sort() ; //Sort in accordance with times of ticks
565   TIter it(ticks) ;
566   AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
567   Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
568   
569   AliEMCALTick * t ;  
570   while((t=(AliEMCALTick*) it.Next())){
571     if(t->GetTime() < time)  //This tick starts before crossing
572       *ctick+=*t ;
573     else
574       return time ;
575     
576     time = ctick->CrossingTime(fTimeThreshold) ;    
577   }
578   return time ;
579 }
580
581 //____________________________________________________________________________ 
582 Bool_t AliEMCALDigitizer::Init()
583 {
584   // Makes all memory allocations
585   fInit = kTRUE ; 
586   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
587
588   if ( emcalLoader == 0 ) {
589     Fatal("Init", "Could not obtain the AliEMCALLoader");  
590     return kFALSE;
591   } 
592
593   fFirstEvent = 0 ; 
594   fLastEvent = fFirstEvent ; 
595   
596   if(fManager)
597     fInput = fManager->GetNinputs() ; 
598   else 
599     fInput           = 1 ;
600
601   fInputFileNames  = new TString[fInput] ; 
602   fEventNames      = new TString[fInput] ; 
603   fInputFileNames[0] = GetTitle() ; 
604   fEventNames[0]     = fEventFolderName.Data() ; 
605   Int_t index ; 
606   for (index = 1 ; index < fInput ; index++) {
607     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0); 
608     TString tempo = fManager->GetInputFolderName(index) ;
609     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager 
610   }
611   
612   //to prevent cleaning of this object while GetEvent is called
613   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
614
615   //PH  Print();
616   //Calibration instance
617   fCalibData = emcalLoader->CalibData();
618   return fInit ;    
619 }
620
621 //____________________________________________________________________________ 
622 void AliEMCALDigitizer::InitParameters()
623
624   // Parameter initialization for digitizer
625   // Tune parameters - 24-nov-04; Apr 29, 2007
626   // New parameters JLK 14-Apr-2008
627
628   fMeanPhotonElectron = 4400;  // electrons per GeV 
629   fPinNoise           = 0.037; // pin noise in GEV from analysis test beam data 
630   if (fPinNoise == 0. ) 
631     Warning("InitParameters", "No noise added\n") ; 
632   fDigitThreshold     = fPinNoise * 3; // 3 * sigma
633   fTimeResolution     = 0.3e-9 ; // 300 psc
634   fTimeSignalLength   = 1.0e-9 ;
635
636   // These defaults are normally not used. 
637   // Values are read from calibration database instead
638   fADCchannelEC    = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
639   fADCpedestalEC   = 0.0 ;  // GeV
640   fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
641
642   fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
643   // hists. for control; no hists on default
644   fControlHists = 0;
645   fHists        = 0;
646 }
647
648 //__________________________________________________________________
649 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
650 {
651   // Allows to produce digits by superimposing background and signal event.
652   // It is assumed, that headers file with SIGNAL events is opened in 
653   // the constructor. 
654   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed 
655   // Thus we avoid writing (changing) huge and expensive 
656   // backgound files: all output will be writen into SIGNAL, i.e. 
657   // opened in constructor file. 
658   //
659   // One can open as many files to mix with as one needs.
660   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
661   // can be mixed.
662
663   if( strcmp(GetName(), "") == 0 )
664     Init() ;
665   
666   if(fManager){
667     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
668     return ;
669   } 
670   // looking for file which contains AliRun
671   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
672     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
673     return ; 
674   }
675   // looking for the file which contains SDigits
676   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
677   TString fileName( emcalLoader->GetSDigitsFileName() ) ; 
678     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
679       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
680     if ( (gSystem->AccessPathName(fileName)) ) { 
681       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
682       return ;
683     }
684     // need to increase the arrays
685     // MvL: This code only works when fInput == 1, I think.
686     TString tempo = fInputFileNames[fInput-1] ; 
687     delete [] fInputFileNames ; 
688     fInputFileNames = new TString[fInput+1] ; 
689     fInputFileNames[fInput-1] = tempo ; 
690  
691     tempo = fEventNames[fInput-1] ; 
692     delete [] fEventNames ; 
693     fEventNames = new TString[fInput+1] ; 
694     fEventNames[fInput-1] = tempo ; 
695
696     fInputFileNames[fInput] = alirunFileName ; 
697     fEventNames[fInput]     = eventFolderName ;
698     fInput++ ;
699 }  
700
701 //__________________________________________________________________
702 void AliEMCALDigitizer::Print1(Option_t * option)
703 { // 19-nov-04 - just for convinience
704   Print(); 
705   PrintDigits(option);
706 }
707
708 //__________________________________________________________________
709 void AliEMCALDigitizer::Print(Option_t*)const 
710 {
711   // Print Digitizer's parameters
712   printf("Print: \n------------------- %s -------------", GetName() ) ; 
713   if( strcmp(fEventFolderName.Data(), "") != 0 ){
714     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
715     
716     Int_t nStreams ; 
717     if (fManager) 
718       nStreams =  GetNInputStreams() ;
719     else 
720       nStreams = fInput ; 
721     
722     Int_t index = 0 ;  
723
724     AliRunLoader *rl=0;
725
726     for (index = 0 ; index < nStreams ; index++) {  
727       TString tempo(fEventNames[index]) ; 
728       tempo += index ;
729       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
730         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ; 
731       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
732       TString fileName( emcalLoader->GetSDigitsFileName() ) ; 
733       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
734         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
735       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ; 
736     }
737
738     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
739
740     printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
741     
742     printf("\nWith following parameters:\n") ;
743     
744     printf("    Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
745     printf("    Threshold  in EMC  (fDigitThreshold) = %f\n", fDigitThreshold)  ;
746     printf("---------------------------------------------------\n")  ;
747   }
748   else
749     printf("Print: AliEMCALDigitizer not initialized") ; 
750 }
751
752 //__________________________________________________________________
753 void AliEMCALDigitizer::PrintDigits(Option_t * option)
754 {
755   //utility method for printing digit information
756
757   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
758   TClonesArray * digits  = emcalLoader->Digits() ;
759   TClonesArray * sdigits = emcalLoader->SDigits() ;
760   
761   printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ; 
762   printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
763   
764   if(strstr(option,"all")){  
765     //loop over digits
766     AliEMCALDigit * digit;
767     printf("\nEMC digits (with primaries):\n")  ;
768     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
769     Int_t index ;
770     for (index = 0 ; index < digits->GetEntries() ; index++) {
771       digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
772       printf("\n%6d  %8d    %6.5e %4d      %2d : ",
773               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
774       Int_t iprimary;
775       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
776         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
777       }
778     }   
779   }
780   printf("\n");
781 }
782
783 //__________________________________________________________________
784 Float_t AliEMCALDigitizer::TimeOfNoise(void)
785 {  
786   // Calculates the time signal generated by noise
787   //PH  Info("TimeOfNoise", "Change me") ; 
788   return gRandom->Rndm() * 1.28E-5;
789 }
790
791 //__________________________________________________________________
792 void AliEMCALDigitizer::Unload() 
793 {  
794   // Unloads the SDigits and Digits
795   AliRunLoader *rl=0;
796     
797   Int_t i ; 
798   for(i = 1 ; i < fInput ; i++){
799     TString tempo(fEventNames[i]) ; 
800     tempo += i;
801     if ((rl = AliRunLoader::GetRunLoader(tempo))) 
802       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ; 
803   }
804   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
805   emcalLoader->UnloadDigits() ; 
806 }
807
808 //_________________________________________________________________________________________
809 void AliEMCALDigitizer::WriteDigits()
810 {
811
812   // Makes TreeD in the output file. 
813   // Check if branch already exists: 
814   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
815   //      already existing branches. 
816   //   else creates branch with Digits, named "EMCAL", title "...",
817   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
818   //      and names of files, from which digits are made.
819
820   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
821
822   const TClonesArray * digits = emcalLoader->Digits() ; 
823   TTree * treeD = emcalLoader->TreeD(); 
824   if ( !treeD ) {
825     emcalLoader->MakeDigitsContainer();
826     treeD = emcalLoader->TreeD(); 
827   }
828
829   // -- create Digits branch
830   Int_t bufferSize = 32000 ;    
831   TBranch * digitsBranch = 0;
832   if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
833     digitsBranch->SetAddress(&digits);
834     AliWarning("Digits Branch already exists. Not all digits will be visible");
835   }
836   else
837     treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
838   //digitsBranch->SetTitle(fEventFolderName);
839   treeD->Fill() ;
840   
841   emcalLoader->WriteDigits("OVERWRITE");
842   emcalLoader->WriteDigitizer("OVERWRITE");
843
844   Unload() ; 
845
846 }
847
848 //__________________________________________________________________
849 void AliEMCALDigitizer::Browse(TBrowser* b)
850 {
851   if(fHists) b->Add(fHists);
852   TTask::Browse(b);
853 }
854
855 //__________________________________________________________________
856 TList *AliEMCALDigitizer::BookControlHists(int var)
857
858   // 22-nov-04
859   // histograms for monitoring digitizer performance
860
861   Info("BookControlHists"," started ");
862   gROOT->cd();
863   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
864   if(var>=1){
865     new TH1F("hDigiN",  "#EMCAL digits with fAmp > fDigitThreshold", 
866     fNADCEC+1, -0.5, Double_t(fNADCEC));
867     new TH1F("HDigiSumEnergy","Sum.EMCAL energy from digi", 1000, 0.0, 200.);
868     new TH1F("hDigiAmp",  "EMCAL digital amplitude", fNADCEC+1, -0.5, Double_t(fNADCEC));
869     new TH1F("hDigiEnergy","EMCAL cell energy", 2000, 0.0, 200.);
870     new TH1F("hDigiAbsId","EMCAL absId cells with fAmp > fDigitThreshold ",
871     geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
872   }
873
874   fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalDigiControlHists", kFALSE);
875   fHists = 0; //huh? JLK 03-Mar-2006
876   return fHists;
877 }
878
879 //__________________________________________________________________
880 void AliEMCALDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
881 {
882   AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt); 
883 }