Access simulation parameters in digitizers via file in OCDB, simulation parameter...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSDigitizer.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 // This is a TTask that makes SDigits out of Hits
20 // A Summable Digits is the sum of all hits originating 
21 // from one in one tower of the EMCAL 
22 // A threshold for assignment of the primary to SDigit is applied 
23 //
24 // JLK 26-Jun-2008 Added explanation:
25 // SDigits need to hold the energy sum of the hits, but AliEMCALDigit
26 // can (should) only store amplitude.  Therefore, the SDigit energy is
27 // "digitized" before being stored and must be "calibrated" back to an
28 // energy before SDigits are summed to form true Digits
29 //
30 // SDigits are written to TreeS, branch "EMCAL"
31 // AliEMCALSDigitizer with all current parameters is written 
32 // to TreeS branch "AliEMCALSDigitizer".
33 // Both branches have the same title. If necessary one can produce 
34 // another set of SDigits with different parameters. Two versions
35 // can be distunguished using titles of the branches.
36 // User case:
37 // root [0] AliEMCALSDigitizer * s = new AliEMCALSDigitizer("galice.root")
38 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
39 // root [1] s->ExecuteTask()
40 //             // Makes SDigitis for all events stored in galice.root
41 // root [2] s->SetPedestalParameter(0.001)
42 //             // One can change parameters of digitization
43 // root [3] s->SetSDigitsBranch("Redestal 0.001")
44 //             // and write them into the new branch
45 // root [4] s->ExeciteTask("deb all tim")
46 //             // available parameters:
47 //             deb - print # of produced SDigitis
48 //             deb all  - print # and list of produced SDigits
49 //             tim - print benchmarking information
50 //
51 //*-- Author : Sahal Yacoob (LBL)
52 // based on  : AliPHOSSDigitzer 
53 //////////////////////////////////////////////////////////////////////////////
54
55 // --- ROOT system ---
56 #include <TBenchmark.h>
57 #include <TBrowser.h>
58 //#include <Riostream.h>
59 #include <TMath.h>
60 #include <TROOT.h>
61
62 // --- Standard library ---
63 #include "stdlib.h"
64
65 // --- AliRoot header files ---
66 #include "AliLog.h"
67 #include "AliRunLoader.h"
68 #include "AliStack.h"
69 #include "AliEMCALDigit.h"
70 #include "AliEMCALLoader.h"
71 #include "AliEMCALHit.h"
72 #include "AliEMCALSDigitizer.h"
73 #include "AliEMCALGeometry.h"
74 #include "AliEMCALSimParam.h"
75
76 ClassImp(AliEMCALSDigitizer)
77            
78 //____________________________________________________________________________ 
79 AliEMCALSDigitizer::AliEMCALSDigitizer()
80   : TTask("",""),
81     fA(0.),fB(0.),fECPrimThreshold(0.),
82     fDefaultInit(kTRUE),
83     fEventFolderName(0),
84     fInit(0),
85     fSDigitsInRun(0),
86     fFirstEvent(0),
87     fLastEvent(0),
88     fSampling(0.)
89 {
90   // ctor
91   InitParameters();
92 }
93
94 //____________________________________________________________________________ 
95 AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, 
96                                        const char * eventFolderName)
97   : TTask("EMCAL"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
98     fA(0.),fB(0.),fECPrimThreshold(0.),
99     fDefaultInit(kFALSE),
100     fEventFolderName(eventFolderName),
101     fInit(0),
102     fSDigitsInRun(0),
103     fFirstEvent(0),
104     fLastEvent(0),
105     fSampling(0.)
106 {
107   // ctor
108   Init();
109   InitParameters() ; 
110
111 }
112
113
114 //____________________________________________________________________________ 
115 AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) 
116   : TTask(sd.GetName(),sd.GetTitle()),
117     fA(sd.fA),
118     fB(sd.fB),
119     fECPrimThreshold(sd.fECPrimThreshold),
120     fDefaultInit(sd.fDefaultInit),
121     fEventFolderName(sd.fEventFolderName),
122     fInit(sd.fInit),
123     fSDigitsInRun(sd.fSDigitsInRun),
124     fFirstEvent(sd.fFirstEvent),
125     fLastEvent(sd.fLastEvent),
126     fSampling(sd.fSampling)
127 {
128   //cpy ctor 
129 }
130
131
132 //____________________________________________________________________________ 
133 AliEMCALSDigitizer::~AliEMCALSDigitizer() {
134   //dtor
135   AliLoader *emcalLoader=0;
136   if ((emcalLoader = AliRunLoader::Instance()->GetDetectorLoader("EMCAL")))
137       emcalLoader->CleanSDigitizer();
138 }
139
140 //____________________________________________________________________________ 
141 void AliEMCALSDigitizer::Init(){
142   // Initialization: open root-file, allocate arrays for hits and sdigits,
143   // attach task SDigitizer to the list of EMCAL tasks
144   // 
145   // Initialization can not be done in the default constructor
146   //============================================================= YS
147   //  The initialisation is now done by the getter
148
149   fInit = kTRUE; 
150    
151   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
152
153   if ( emcalLoader == 0 ) {
154     Fatal("Init", "Could not obtain the AliEMCALLoader");
155     return ;
156   } 
157   
158   emcalLoader->PostSDigitizer(this);
159   emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
160   
161 }
162
163 //____________________________________________________________________________ 
164 void AliEMCALSDigitizer::InitParameters()
165 {
166   //initialize parameters for sdigitization
167
168   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
169   if (geom->GetSampling() == 0.) {
170     Fatal("InitParameters", "Sampling factor not set !") ; 
171   }
172
173   // Get the parameters from the OCDB via the loader
174   AliRunLoader *rl = AliRunLoader::Instance();
175   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
176   AliEMCALSimParam * simParam = 0x0;
177   if(emcalLoader) simParam = emcalLoader->SimulationParameters();
178         
179   if(!simParam){
180         simParam = AliEMCALSimParam::GetInstance();
181         AliWarning("Simulation Parameters not available in OCDB?");
182   }
183         
184   //
185   //JLK 26-Jun-2008 THIS SHOULD HAVE BEEN EXPLAINED AGES AGO:
186   //
187   //In order to be able to store SDigit Energy info into
188   //AliEMCALDigit, we need to convert it temporarily to an ADC amplitude
189   //and later when summing SDigits to form digits, convert it back to
190   //energy.  These fA and fB parameters accomplish this through the
191   //Digitize() and Calibrate() methods
192   //
193   // Initializes parameters
194   fA         = simParam->GetA(); //0;
195   fB         = simParam->GetB(); //1.e+6;  // Changed 24 Apr 2007. Dynamic range now 2 TeV
196   fSampling  = geom->GetSampling();
197
198   // threshold for deposit energy of hit
199   fECPrimThreshold  = simParam->GetECPrimaryThreshold();//0.05;// GeV // 22-may-07 was 0// 24-nov-04 - was 1.e-6;
200   
201   AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
202   AliDebug(2,Form("   fInit                                 %i\n", int(fInit)));
203   AliDebug(2,Form("   fFirstEvent                           %i\n", fFirstEvent));
204   AliDebug(2,Form("   fLastEvent                            %i\n", fLastEvent));
205   AliDebug(2,Form("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()));
206   AliDebug(2,Form("   with digitization parameters       A = %f\n", fA));
207   AliDebug(2,Form("                                      B = %f\n", fB));
208   AliDebug(2,Form("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold));
209   AliDebug(2,Form("   Sampling                             = %f\n", fSampling));
210   AliDebug(2,Form("---------------------------------------------------\n"));
211
212 }
213
214 //____________________________________________________________________________
215 void AliEMCALSDigitizer::Exec(Option_t *option) 
216
217   // Collects all hit of the same tower into digits
218   TString o(option); o.ToUpper();
219   if (strstr(option, "print") ) {
220
221     AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName()));
222     AliDebug(2,Form("   fInit                                 %i\n", int(fInit)));
223     AliDebug(2,Form("   fFirstEvent                           %i\n", fFirstEvent));
224     AliDebug(2,Form("   fLastEvent                            %i\n", fLastEvent));
225     AliDebug(2,Form("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()));
226     AliDebug(2,Form("   with digitization parameters       A = %f\n", fA));
227     AliDebug(2,Form("                                      B = %f\n", fB));
228     AliDebug(2,Form("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold));
229     AliDebug(2,Form("   Sampling                             = %f\n", fSampling));
230     AliDebug(2,Form("---------------------------------------------------\n"));
231
232     return ; 
233   }
234  
235
236   if(strstr(option,"tim"))
237     gBenchmark->Start("EMCALSDigitizer");
238
239   AliRunLoader *rl = AliRunLoader::Instance();
240   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
241
242   //switch off reloading of this task while getting event
243   if (!fInit) { // to prevent overwrite existing file
244     AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ;
245     return ;
246     }
247
248   if (fLastEvent == -1) 
249     fLastEvent = rl->GetNumberOfEvents() - 1 ;
250   else {
251     fLastEvent = TMath::Min(fLastEvent, rl->GetNumberOfEvents()-1);
252   }
253   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
254
255   Int_t ievent;
256   Float_t energy=0.; // de * fSampling - 23-nov-04
257   rl->LoadKinematics();
258   rl->LoadHits("EMCAL");
259
260   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
261     rl->GetEvent(ievent);
262     TTree * treeS = emcalLoader->TreeS();
263     if ( !treeS ) { 
264       emcalLoader->MakeSDigitsContainer();
265       treeS = emcalLoader->TreeS();
266     }
267     TClonesArray * hits = emcalLoader->Hits() ; 
268     TClonesArray * sdigits = emcalLoader->SDigits() ;
269     sdigits->Clear();
270
271     Int_t nSdigits = 0 ;
272     Int_t i;
273
274         AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); 
275     for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
276       AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
277       AliEMCALDigit * curSDigit = 0 ;
278       AliEMCALDigit * sdigit = 0 ;
279       Bool_t newsdigit = kTRUE; 
280       
281       // hit->GetId() - Absolute Id number EMCAL segment
282       if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId())
283         energy = hit->GetEnergy() * fSampling; // 23-nov-04
284         if(energy >  fECPrimThreshold )
285           // Assign primary number only if deposited energy is significant
286           curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
287                                           hit->GetIparent(), hit->GetId(), 
288                                           Digitize(energy), hit->GetTime(), 
289                                           -1, energy ) ;
290           else
291             curSDigit =  new AliEMCALDigit( -1               , 
292                                             -1               ,
293                                             hit->GetId(), 
294                                             Digitize(energy), hit->GetTime(), 
295                                             -1, energy ) ;
296       } else {
297         Warning("Exec"," abs id %i is bad \n", hit->GetId());
298         newsdigit = kFALSE;
299         curSDigit = 0;
300       }
301       
302       if(curSDigit != 0){
303         for(Int_t check= 0; check < nSdigits ; check++) {
304           sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
305           
306           if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ?              
307             *sdigit = *sdigit + *curSDigit;
308             newsdigit = kFALSE;
309           }
310         }
311       }
312       if (newsdigit) {
313         new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
314         nSdigits++ ;  
315       }
316       delete curSDigit ; 
317     }  // loop over all hits (hit = deposited energy/entering particle)
318     sdigits->Sort() ;
319     
320     nSdigits = sdigits->GetEntriesFast() ;
321     fSDigitsInRun += nSdigits ;  
322
323     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
324       AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
325       sdigit->SetIndexInList(i) ;
326     }
327     
328     // Now write SDigits    
329     
330     Int_t bufferSize = 32000 ;    
331     TBranch * sdigitsBranch = treeS->GetBranch("EMCAL");
332     if (sdigitsBranch)
333       sdigitsBranch->SetAddress(&sdigits);
334     else
335       treeS->Branch("EMCAL",&sdigits,bufferSize);
336     
337     treeS->Fill();
338     
339     emcalLoader->WriteSDigits("OVERWRITE");
340     
341     //NEXT - SDigitizer
342     emcalLoader->WriteSDigitizer("OVERWRITE");  // why in event cycle ?
343     
344     if(strstr(option,"deb"))
345       PrintSDigits(option) ;  
346   }
347   
348   Unload();
349   
350   emcalLoader->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
351   
352   if(strstr(option,"tim")){
353     gBenchmark->Stop("EMCALSDigitizer"); 
354     printf("\n Exec: took %f seconds for SDigitizing %f seconds per event\n", 
355            gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ; 
356   }
357 }
358
359 //__________________________________________________________________
360 Int_t AliEMCALSDigitizer::Digitize(Float_t energy)const {
361   // Digitize the energy
362   //
363   //JLK 26-Jun-2008 EXPLANATION LONG OVERDUE:
364   //
365   //We have to digitize the SDigit energy so that it can be stored in
366   //AliEMCALDigit, which has only an ADC amplitude member and
367   //(rightly) no energy member.  This method converts the energy to an
368   //integer which can be re-converted back to an energy with the
369   //Calibrate(energy) method when it is time to create Digits from SDigits 
370   //
371   Double_t aSignal = fA + energy*fB;
372   if (TMath::Abs(aSignal)>2147483647.0) { 
373     //PH 2147483647 is the max. integer
374     //PH This apparently is a problem which needs investigation
375     AliWarning(Form("Too big or too small energy %f",aSignal));
376     aSignal = TMath::Sign((Double_t)2147483647,aSignal);
377   }
378   return (Int_t ) aSignal;
379 }
380
381 //__________________________________________________________________
382 Float_t AliEMCALSDigitizer::Calibrate(Int_t amp)const {
383   //
384   // Convert the amplitude back to energy in GeV
385   //
386   //JLK 26-Jun-2008 EXPLANATION LONG OVERDUE:
387   //
388   //We have to digitize the SDigit energy with the method Digitize() 
389   //so that it can be stored in AliEMCALDigit, which has only an ADC 
390   //amplitude member and (rightly) no energy member.  This method is
391   //just the reverse of Digitize(): it converts the stored amplitude 
392   //back to an energy value in GeV so that the SDigit energies can be 
393   //summed before adding noise and creating digits out of them
394   //
395   return (Float_t)(amp - fA)/fB;
396
397 }
398  
399
400 //__________________________________________________________________
401 void AliEMCALSDigitizer::Print1(Option_t * option)
402 {
403   Print(); 
404   PrintSDigits(option);
405 }
406
407 //__________________________________________________________________
408 void AliEMCALSDigitizer::Print(Option_t *option) const
409
410   // Prints parameters of SDigitizer
411   printf("Print: \n------------------- %s ------------- option %s\n", GetName() , option) ; 
412   printf("   fInit                                 %i\n", int(fInit));
413   printf("   fFirstEvent                           %i\n", fFirstEvent);
414   printf("   fLastEvent                            %i\n", fLastEvent);
415   printf("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()) ;
416   printf("   with digitization parameters       A = %f\n", fA) ; 
417   printf("                                      B = %f\n", fB) ;
418   printf("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold)  ;
419   printf("   Sampling                             = %f\n", fSampling);
420   printf("---------------------------------------------------\n") ;
421 }
422
423 //__________________________________________________________________
424 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
425 {
426   // Equal operator.
427   // SDititizers are equal if their pedestal, slope and threshold are equal
428   if( (fA==sd.fA)&&(fB==sd.fB)&&
429       (fECPrimThreshold==sd.fECPrimThreshold))
430     return kTRUE ;
431   else
432     return kFALSE ;
433 }
434
435 //__________________________________________________________________
436 void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
437 {
438   //Prints list of digits produced at the current pass of AliEMCALDigitizer
439     
440   AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
441   const TClonesArray * sdigits = rl->SDigits() ; 
442   
443   printf("\n") ;  
444   printf("event %i", rl->GetRunLoader()->GetEventNumber());
445   printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast()); 
446   if(strstr(option,"all")||strstr(option,"EMC")){
447     
448     //loop over digits
449     AliEMCALDigit * digit;
450     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
451     Int_t index, isum=0;
452     char * tempo = new char[8192]; 
453     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
454       digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
455       sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d :",
456               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
457       printf("%s",tempo);
458       isum += digit->GetAmp();
459       
460       Int_t iprimary;
461       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
462                   sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
463         printf("%s",tempo); 
464       }          
465     }
466     delete tempo ;
467     printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, Calibrate(isum));
468   } else printf("\n");
469 }
470
471 //____________________________________________________________________________ 
472 void AliEMCALSDigitizer::Unload() const
473 {
474   // Unload Hits and SDigits from the folder
475   AliEMCALLoader *rl = dynamic_cast<AliEMCALLoader*>(AliRunLoader::Instance()->GetDetectorLoader("EMCAL"));
476   rl->UnloadHits() ; 
477   rl->UnloadSDigits() ; 
478 }
479
480 //____________________________________________________________________________ 
481 void AliEMCALSDigitizer::Browse(TBrowser* b)
482 {
483   TTask::Browse(b);
484 }