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