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