updating apd and fee parameters
[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() ;  
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       Float_t a = digit->GetAmp() ;
342       Float_t b = TMath::Abs( a /fTimeSignalLength) ;
343       //Mark the beginning of the signal
344       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);  
345       //Mark the end of the signal     
346       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
347       
348       // loop over input
349       for(i = 0; i< fInput ; i++){  //loop over (possible) merge sources
350         if(dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
351           curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;      
352         else
353           curSDigit = 0 ;
354         //May be several digits will contribute from the same input
355         while(curSDigit && (curSDigit->GetId() == absID)){         
356           //Shift primary to separate primaries belonging different inputs
357           Int_t primaryoffset ;
358           if(fManager)
359             primaryoffset = fManager->GetMask(i) ; 
360           else
361             primaryoffset = i ;
362           curSDigit->ShiftPrimary(primaryoffset) ;
363           
364           a = curSDigit->GetAmp() ;
365           b = a /fTimeSignalLength ;
366           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);  
367           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
368
369           *digit = *digit + *curSDigit ;  //add energies
370
371           index[i]++ ;
372           if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
373             curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
374           else
375             curSDigit = 0 ;
376         }
377       }
378       // add fluctuations for photo-electron creation
379       amp = sDigitizer->Calibrate(digit->GetAmp()) ; // GeV
380       amp *= static_cast<Float_t>(gRandom->Poisson(fMeanPhotonElectron)) / static_cast<Float_t>(fMeanPhotonElectron) ;
381   
382       //calculate and set time
383       Float_t time = FrontEdgeTime(ticks) ;
384       digit->SetTime(time) ;
385
386       //Find next signal module
387       nextSig = nEMC + 1 ;
388       for(i = 0 ; i < fInput ; i++){
389         sdigits = dynamic_cast<TClonesArray *>(sdigArray->At(i)) ;
390         Int_t curNext = nextSig ;
391         if(sdigits->GetEntriesFast() > index[i] ){
392           curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(index[i]))->GetId() ;       
393         }
394         if(curNext < nextSig) nextSig = curNext ;
395       }
396     }
397     // add the noise now
398     amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
399     digit->SetAmp(sDigitizer->Digitize(amp)) ;  
400     AliDebug(10,Form(" absID %5i amp %f nextSig %5i\n",
401                      absID, amp, nextSig));
402   } // for(absID = 1; absID <= nEMC; absID++)
403   
404   ticks->Delete() ;
405   delete ticks ;
406
407   delete sdigArray ; //We should not delete its contents
408
409   //remove digits below thresholds
410   for(i = 0 ; i < nEMC ; i++){
411     digit = dynamic_cast<AliEMCALDigit*>( digits->At(i) ) ;
412     Float_t threshold = fDigitThreshold ; 
413     if(sDigitizer->Calibrate( digit->GetAmp() ) < threshold)
414       digits->RemoveAt(i) ;
415     else 
416       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
417   }
418   
419   digits->Compress() ;  
420   
421   Int_t ndigits = digits->GetEntriesFast() ; 
422   
423   //Set indexes in list of digits and fill hists.
424   AliEMCALHistoUtilities::FillH1(fHists, 0, Double_t(ndigits));
425   Float_t energy=0., esum=0.;
426   for (i = 0 ; i < ndigits ; i++) { 
427     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ; 
428     digit->SetIndexInList(i) ; 
429     energy = sDigitizer->Calibrate(digit->GetAmp()) ;
430     esum += energy;
431     digit->SetAmp(DigitizeEnergy(energy, digit->GetId()) ) ; // for what ??
432     AliEMCALHistoUtilities::FillH1(fHists, 2, double(digit->GetAmp()));
433     AliEMCALHistoUtilities::FillH1(fHists, 3, double(energy));
434     AliEMCALHistoUtilities::FillH1(fHists, 4, double(digit->GetId()));
435     //    if(digit->GetId() == nEMC) {
436     //  printf(" i %i \n", i );
437     //  digit->Dump();
438     //  assert(0);
439     //}
440   }
441   AliEMCALHistoUtilities::FillH1(fHists, 1, esum);
442 }
443
444 // //_____________________________________________________________________
445 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t AbsId)
446
447   // Returns digitized value of the energy in a cell absId
448
449   // Load Geometry
450   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
451
452   if (geom==0)
453     AliFatal("Did not get geometry from EMCALLoader");
454
455   Int_t iSupMod = -1;
456   Int_t nModule  = -1;
457   Int_t nIphi   = -1;
458   Int_t nIeta   = -1;
459   Int_t iphi    = -1;
460   Int_t ieta    = -1;
461   Int_t channel = -999; 
462
463   Bool_t bCell = geom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
464   if(!bCell)
465     Error("DigitizeEnergy","Wrong cell id number : AbsId %i ", AbsId) ;
466   geom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
467   
468   if(fCalibData) {
469     fADCpedestalEC = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
470     fADCchannelEC = fCalibData->GetADCchannel(iSupMod,ieta,iphi);
471   }
472   
473   channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalEC)/fADCchannelEC ))  ;
474   
475   if(channel > fNADCEC ) 
476     channel =  fNADCEC ; 
477   return channel ;
478   
479 }
480
481 //____________________________________________________________________________
482 void AliEMCALDigitizer::Exec(Option_t *option) 
483
484   // Steering method to process digitization for events
485   // in the range from fFirstEvent to fLastEvent.
486   // This range is optionally set by SetEventRange().
487   // if fLastEvent=-1, then process events until the end.
488   // by default fLastEvent = fFirstEvent (process only one event)
489
490   if (!fInit) { // to prevent overwrite existing file
491     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
492     return ;
493   } 
494
495   if (strstr(option,"print")) {
496
497     Print();
498     return ; 
499   }
500   
501   if(strstr(option,"tim"))
502     gBenchmark->Start("EMCALDigitizer");
503
504   AliRunLoader *rl = AliRunLoader::GetRunLoader();
505   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
506
507   // Post Digitizer to the white board
508   emcalLoader->PostDigitizer(this) ;
509   
510   if (fLastEvent == -1)  {
511     fLastEvent = rl->GetNumberOfEvents() - 1 ;
512   }
513   else if (fManager) 
514     fLastEvent = fFirstEvent ; // what is this ??
515
516   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
517   Int_t ievent;
518
519   rl->LoadSDigits("EMCAL");
520   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
521     
522     rl->GetEvent(ievent);
523
524     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
525
526     WriteDigits() ;
527
528     if(strstr(option,"deb"))
529       PrintDigits(option);
530     if(strstr(option,"table")) gObjectTable->Print();
531
532     //increment the total number of Digits per run 
533     fDigitsInRun += emcalLoader->Digits()->GetEntriesFast() ;  
534   }
535   
536   emcalLoader->CleanDigitizer() ;
537
538   if(strstr(option,"tim")){
539     gBenchmark->Stop("EMCALDigitizer");
540     AliInfo(Form("Exec: took %f seconds for Digitizing %f seconds per event", 
541          gBenchmark->GetCpuTime("EMCALDigitizer"), gBenchmark->GetCpuTime("EMCALDigitizer")/nEvents )) ;
542   } 
543 }
544
545 //____________________________________________________________________________ 
546 Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks) 
547
548   //  Returns the shortest time among all time ticks
549
550   ticks->Sort() ; //Sort in accordance with times of ticks
551   TIter it(ticks) ;
552   AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
553   Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
554   
555   AliEMCALTick * t ;  
556   while((t=(AliEMCALTick*) it.Next())){
557     if(t->GetTime() < time)  //This tick starts before crossing
558       *ctick+=*t ;
559     else
560       return time ;
561     
562     time = ctick->CrossingTime(fTimeThreshold) ;    
563   }
564   return time ;
565 }
566
567 //____________________________________________________________________________ 
568 Bool_t AliEMCALDigitizer::Init()
569 {
570   // Makes all memory allocations
571   fInit = kTRUE ; 
572   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
573
574   if ( emcalLoader == 0 ) {
575     Fatal("Init", "Could not obtain the AliEMCALLoader");  
576     return kFALSE;
577   } 
578
579   fFirstEvent = 0 ; 
580   fLastEvent = fFirstEvent ; 
581   
582   if(fManager)
583     fInput = fManager->GetNinputs() ; 
584   else 
585     fInput           = 1 ;
586
587   fInputFileNames  = new TString[fInput] ; 
588   fEventNames      = new TString[fInput] ; 
589   fInputFileNames[0] = GetTitle() ; 
590   fEventNames[0]     = fEventFolderName.Data() ; 
591   Int_t index ; 
592   for (index = 1 ; index < fInput ; index++) {
593     fInputFileNames[index] = dynamic_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0); 
594     TString tempo = fManager->GetInputFolderName(index) ;
595     fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added bt fManager 
596   }
597   
598   //to prevent cleaning of this object while GetEvent is called
599   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
600
601   //PH  Print();
602   //Calibration instance
603   fCalibData = emcalLoader->CalibData();
604   return fInit ;    
605 }
606
607 //____________________________________________________________________________ 
608 void AliEMCALDigitizer::InitParameters()
609
610   // Parameter initialization for digitizer
611   // Tune parameters - 24-nov-04; Apr 29, 2007
612   // New parameters JLK 14-Apr-2008
613
614   fMeanPhotonElectron = 4400;  // electrons per GeV 
615   fPinNoise           = 0.037; // pin noise in GEV from analysis test beam data 
616   if (fPinNoise == 0. ) 
617     Warning("InitParameters", "No noise added\n") ; 
618   fDigitThreshold     = fPinNoise * 3; // 3 * sigma
619   fTimeResolution     = 0.3e-9 ; // 300 psc
620   fTimeSignalLength   = 1.0e-9 ;
621
622   // These defaults are normally not used. 
623   // Values are read from calibration database instead
624   fADCchannelEC    = 0.0153; // Update 24 Apr 2007: 250./16/1024 - width of one ADC channel in GeV
625   fADCpedestalEC   = 0.0 ;  // GeV
626   fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
627
628   fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
629   // hists. for control; no hists on default
630   fControlHists = 0;
631   fHists        = 0;
632 }
633
634 //__________________________________________________________________
635 void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
636 {
637   // Allows to produce digits by superimposing background and signal event.
638   // It is assumed, that headers file with SIGNAL events is opened in 
639   // the constructor. 
640   // Sets the BACKGROUND event, with which the SIGNAL event is to be mixed 
641   // Thus we avoid writing (changing) huge and expensive 
642   // backgound files: all output will be writen into SIGNAL, i.e. 
643   // opened in constructor file. 
644   //
645   // One can open as many files to mix with as one needs.
646   // However only Sdigits with the same name (i.e. constructed with the same SDigitizer)
647   // can be mixed.
648
649   if( strcmp(GetName(), "") == 0 )
650     Init() ;
651   
652   if(fManager){
653     Error("MixWith", "Cannot use this method under AliRunDigitizer") ;
654     return ;
655   } 
656   // looking for file which contains AliRun
657   if (gSystem->AccessPathName(alirunFileName)) {// file does not exist
658     Error("MixWith", "File %s does not exist!", alirunFileName.Data()) ;
659     return ; 
660   }
661   // looking for the file which contains SDigits
662   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
663   TString fileName( emcalLoader->GetSDigitsFileName() ) ; 
664     if ( eventFolderName != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
665       fileName = fileName.ReplaceAll(".root", "") + "_" + eventFolderName + ".root" ;
666     if ( (gSystem->AccessPathName(fileName)) ) { 
667       Error("MixWith", "The file %s does not exist!", fileName.Data()) ;
668       return ;
669     }
670     // need to increase the arrays
671     // MvL: This code only works when fInput == 1, I think.
672     TString tempo = fInputFileNames[fInput-1] ; 
673     delete [] fInputFileNames ; 
674     fInputFileNames = new TString[fInput+1] ; 
675     fInputFileNames[fInput-1] = tempo ; 
676  
677     tempo = fEventNames[fInput-1] ; 
678     delete [] fEventNames ; 
679     fEventNames = new TString[fInput+1] ; 
680     fEventNames[fInput-1] = tempo ; 
681
682     fInputFileNames[fInput] = alirunFileName ; 
683     fEventNames[fInput]     = eventFolderName ;
684     fInput++ ;
685 }  
686
687 //__________________________________________________________________
688 void AliEMCALDigitizer::Print1(Option_t * option)
689 { // 19-nov-04 - just for convinience
690   Print(); 
691   PrintDigits(option);
692 }
693
694 //__________________________________________________________________
695 void AliEMCALDigitizer::Print(Option_t*)const 
696 {
697   // Print Digitizer's parameters
698   printf("Print: \n------------------- %s -------------", GetName() ) ; 
699   if( strcmp(fEventFolderName.Data(), "") != 0 ){
700     printf(" Writing Digits to branch with title  %s\n", fEventFolderName.Data()) ;
701     
702     Int_t nStreams ; 
703     if (fManager) 
704       nStreams =  GetNInputStreams() ;
705     else 
706       nStreams = fInput ; 
707     
708     Int_t index = 0 ;  
709
710     AliRunLoader *rl=0;
711
712     for (index = 0 ; index < nStreams ; index++) {  
713       TString tempo(fEventNames[index]) ; 
714       tempo += index ;
715       if ((rl = AliRunLoader::GetRunLoader(tempo)) == 0)
716         rl = AliRunLoader::Open(fInputFileNames[index], tempo) ; 
717       AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
718       TString fileName( emcalLoader->GetSDigitsFileName() ) ; 
719       if ( fEventNames[index] != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
720         fileName = fileName.ReplaceAll(".root", "") + "_" + fEventNames[index]  + ".root" ;
721       printf ("Adding SDigits from %s %s\n", fInputFileNames[index].Data(), fileName.Data()) ; 
722     }
723
724     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
725
726     printf("\nWriting digits to %s", emcalLoader->GetDigitsFileName().Data()) ;
727     
728     printf("\nWith following parameters:\n") ;
729     
730     printf("    Electronics noise in EMC (fPinNoise) = %f\n", fPinNoise) ;
731     printf("    Threshold  in EMC  (fDigitThreshold) = %f\n", fDigitThreshold)  ;
732     printf("---------------------------------------------------\n")  ;
733   }
734   else
735     printf("Print: AliEMCALDigitizer not initialized") ; 
736 }
737
738 //__________________________________________________________________
739 void AliEMCALDigitizer::PrintDigits(Option_t * option)
740 {
741   //utility method for printing digit information
742
743   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
744   TClonesArray * digits  = emcalLoader->Digits() ;
745   TClonesArray * sdigits = emcalLoader->SDigits() ;
746   
747   printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ; 
748   printf("\n event %d", emcalLoader->GetRunLoader()->GetEventNumber());
749   
750   if(strstr(option,"all")){  
751     //loop over digits
752     AliEMCALDigit * digit;
753     printf("\nEMC digits (with primaries):\n")  ;
754     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
755     Int_t index ;
756     for (index = 0 ; index < digits->GetEntries() ; index++) {
757       digit = dynamic_cast<AliEMCALDigit *>(digits->At(index)) ;
758       printf("\n%6d  %8d    %6.5e %4d      %2d : ",
759               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
760       Int_t iprimary;
761       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
762         printf("%d ",digit->GetPrimary(iprimary+1) ) ; 
763       }
764     }   
765   }
766   printf("\n");
767 }
768
769 //__________________________________________________________________
770 Float_t AliEMCALDigitizer::TimeOfNoise(void)
771 {  
772   // Calculates the time signal generated by noise
773   //PH  Info("TimeOfNoise", "Change me") ; 
774   return gRandom->Rndm() * 1.28E-5;
775 }
776
777 //__________________________________________________________________
778 void AliEMCALDigitizer::Unload() 
779 {  
780   // Unloads the SDigits and Digits
781   AliRunLoader *rl=0;
782     
783   Int_t i ; 
784   for(i = 1 ; i < fInput ; i++){
785     TString tempo(fEventNames[i]) ; 
786     tempo += i;
787     if ((rl = AliRunLoader::GetRunLoader(tempo))) 
788       rl->GetDetectorLoader("EMCAL")->UnloadSDigits() ; 
789   }
790   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
791   emcalLoader->UnloadDigits() ; 
792 }
793
794 //_________________________________________________________________________________________
795 void AliEMCALDigitizer::WriteDigits()
796 {
797
798   // Makes TreeD in the output file. 
799   // Check if branch already exists: 
800   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
801   //      already existing branches. 
802   //   else creates branch with Digits, named "EMCAL", title "...",
803   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
804   //      and names of files, from which digits are made.
805
806   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
807
808   const TClonesArray * digits = emcalLoader->Digits() ; 
809   TTree * treeD = emcalLoader->TreeD(); 
810   if ( !treeD ) {
811     emcalLoader->MakeDigitsContainer();
812     treeD = emcalLoader->TreeD(); 
813   }
814
815   // -- create Digits branch
816   Int_t bufferSize = 32000 ;    
817   TBranch * digitsBranch = 0;
818   if ((digitsBranch = treeD->GetBranch("EMCAL"))) {
819     digitsBranch->SetAddress(&digits);
820     AliWarning("Digits Branch already exists. Not all digits will be visible");
821   }
822   else
823     treeD->Branch("EMCAL","TClonesArray",&digits,bufferSize);
824   //digitsBranch->SetTitle(fEventFolderName);
825   treeD->Fill() ;
826   
827   emcalLoader->WriteDigits("OVERWRITE");
828   emcalLoader->WriteDigitizer("OVERWRITE");
829
830   Unload() ; 
831
832 }
833
834 //__________________________________________________________________
835 void AliEMCALDigitizer::Browse(TBrowser* b)
836 {
837   if(fHists) b->Add(fHists);
838   TTask::Browse(b);
839 }
840
841 //__________________________________________________________________
842 TList *AliEMCALDigitizer::BookControlHists(int var)
843
844   // 22-nov-04
845   // histograms for monitoring digitizer performance
846
847   Info("BookControlHists"," started ");
848   gROOT->cd();
849   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
850   if(var>=1){
851     new TH1F("hDigiN",  "#EMCAL digits with fAmp > fDigitThreshold", 
852     fNADCEC+1, -0.5, Double_t(fNADCEC));
853     new TH1F("HDigiSumEnergy","Sum.EMCAL energy from digi", 1000, 0.0, 200.);
854     new TH1F("hDigiAmp",  "EMCAL digital amplitude", fNADCEC+1, -0.5, Double_t(fNADCEC));
855     new TH1F("hDigiEnergy","EMCAL cell energy", 2000, 0.0, 200.);
856     new TH1F("hDigiAbsId","EMCAL absId cells with fAmp > fDigitThreshold ",
857     geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
858   }
859
860   fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalDigiControlHists", kFALSE);
861   fHists = 0; //huh? JLK 03-Mar-2006
862   return fHists;
863 }
864
865 //__________________________________________________________________
866 void AliEMCALDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
867 {
868   AliEMCALHistoUtilities::SaveListOfHists(fHists, name, kSingleKey, opt); 
869 }