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