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