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