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