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