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