Coding conventions (F.Carminati)
[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 // SDigits are written to TreeS, branch "EMCAL"
24 // AliEMCALSDigitizer with all current parameters is written 
25 // to TreeS branch "AliEMCALSDigitizer".
26 // Both branches have the same title. If necessary one can produce 
27 // another set of SDigits with different parameters. Two versions
28 // can be distunguished using titles of the branches.
29 // User case:
30 // root [0] AliEMCALSDigitizer * s = new AliEMCALSDigitizer("galice.root")
31 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
32 // root [1] s->ExecuteTask()
33 //             // Makes SDigitis for all events stored in galice.root
34 // root [2] s->SetPedestalParameter(0.001)
35 //             // One can change parameters of digitization
36 // root [3] s->SetSDigitsBranch("Redestal 0.001")
37 //             // and write them into the new branch
38 // root [4] s->ExeciteTask("deb all tim")
39 //             // available parameters:
40 //             deb - print # of produced SDigitis
41 //             deb all  - print # and list of produced SDigits
42 //             tim - print benchmarking information
43 //
44 //*-- Author : Sahal Yacoob (LBL)
45 // based on  : AliPHOSSDigitzer 
46 // Modif: 
47 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
48 //                           of new  IO (à la PHOS)
49 //////////////////////////////////////////////////////////////////////////////
50
51 // --- ROOT system ---
52 #include "TBenchmark.h"
53
54 // --- Standard library ---
55 #include "stdlib.h"
56
57 // --- AliRoot header files ---
58 #include "AliEMCALDigit.h"
59 #include "AliEMCALGetter.h"
60 #include "AliEMCALHit.h"
61 #include "AliEMCALSDigitizer.h"
62 #include "AliEMCALGeometry.h"
63
64 ClassImp(AliEMCALSDigitizer)
65            
66 //____________________________________________________________________________ 
67   AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("","") 
68 {
69   // ctor
70   fFirstEvent = fLastEvent  = 0 ;  
71   fDefaultInit = kTRUE ; 
72 }
73
74 //____________________________________________________________________________ 
75 AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, 
76                                        const char * eventFolderName):
77   TTask("EMCAL"+AliConfig::Instance()->GetSDigitizerTaskName(), alirunFileName),
78   fEventFolderName(eventFolderName)
79 {
80   // ctor
81   fFirstEvent = fLastEvent  = 0 ; // runs one event by defaut  
82   Init();
83   InitParameters() ; 
84   fDefaultInit = kFALSE ; 
85 }
86
87
88 //____________________________________________________________________________ 
89 AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) : TTask(sd) {
90   //cpy ctor 
91
92   fFirstEvent    = sd.fFirstEvent ; 
93   fLastEvent     = sd.fLastEvent ;
94   fA             = sd.fA ;
95   fB             = sd.fB ;
96   fECPrimThreshold  = sd.fECPrimThreshold ;
97   fSDigitsInRun  = sd.fSDigitsInRun ;
98   SetName(sd.GetName()) ; 
99   SetTitle(sd.GetTitle()) ; 
100   fEventFolderName = sd.fEventFolderName;
101 }
102
103
104 //____________________________________________________________________________ 
105 AliEMCALSDigitizer::~AliEMCALSDigitizer() {
106   // dtor
107   AliEMCALGetter * gime = 
108     AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());  
109   gime->EmcalLoader()->CleanSDigitizer();
110 }
111
112 //____________________________________________________________________________ 
113 void AliEMCALSDigitizer::Init(){
114   // Initialization: open root-file, allocate arrays for hits and sdigits,
115   // attach task SDigitizer to the list of EMCAL tasks
116   // 
117   // Initialization can not be done in the default constructor
118   //============================================================= YS
119   //  The initialisation is now done by the getter
120
121   fInit = kTRUE ; 
122    
123   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());  
124   if ( gime == 0 ) {
125     Fatal("Init", "Could not obtain the Getter objectfor file %s and event %s !", GetTitle(), fEventFolderName.Data()) ;  
126     return ;
127   } 
128   
129   TString opt("SDigits") ; 
130   if(gime->VersionExists(opt) ) { 
131     Error( "Init", "Give a version name different from %s", fEventFolderName.Data() ) ;
132     fInit = kFALSE ; 
133   }
134   
135   gime->PostSDigitizer(this);
136   gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
137   
138 }
139
140 //____________________________________________________________________________ 
141 void AliEMCALSDigitizer::InitParameters()
142 {
143   // Initializes parameters
144   fA                      = 0;
145   fB                      = 10000000.;
146
147   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
148   const AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
149   if (geom->GetSampling() == 0.) {
150     Fatal("InitParameters", "Sampling factor not set !") ; 
151   }
152 //   else
153 //     Info("InitParameters", "Sampling factor set to %f", geom->GetSampling()) ; 
154   
155   // this threshold corresponds approximately to 100 MeV
156   fECPrimThreshold     = 100E-3;
157 }
158
159 //____________________________________________________________________________
160 void AliEMCALSDigitizer::Exec(Option_t *option) 
161
162   // Collects all hit of the same tower into digits
163   if (strstr(option, "print") ) {
164     Print() ; 
165     return ; 
166   }
167   
168   if(strstr(option,"tim"))
169     gBenchmark->Start("EMCALSDigitizer");
170
171   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ;
172  
173   //switch off reloading of this task while getting event
174   if (!fInit) { // to prevent overwrite existing file
175     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
176     return ;
177     }
178
179   if (fLastEvent == -1) 
180     fLastEvent = gime->MaxEvent() - 1 ;
181   else 
182     fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // only ine event at the time
183   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
184
185   Int_t ievent ;   
186   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
187     gime->Event(ievent,"H") ;  
188     TTree * treeS = gime->TreeS(); 
189     TClonesArray * hits = gime->Hits() ; 
190     TClonesArray * sdigits = gime->SDigits() ;
191     sdigits->Clear();
192     Int_t nSdigits = 0 ;
193     //Now make SDigits from hits, for EMCAL it is the same, so just copy    
194     Int_t nPrim =  static_cast<Int_t>((gime->TreeH())->GetEntries()) ; 
195     // Attention nPrim is the number of primaries tracked by Geant 
196     // and this number could be different to the number of Primaries in TreeK;
197     Int_t iprim ;
198     for ( iprim = 0 ; iprim < nPrim ; iprim++ ) { 
199       //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
200       gime->Track(iprim) ;
201       Int_t i;
202       for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
203         AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
204         AliEMCALDigit * curSDigit = 0 ;
205         AliEMCALDigit * sdigit = 0 ;
206         Bool_t newsdigit = kTRUE; 
207
208         // Assign primary number only if deposited energy is significant
209         AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
210         if( geom->IsInECA(hit->GetId()) ){
211           if( hit->GetEnergy() >  fECPrimThreshold )
212             curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
213                                             hit->GetIparent(), hit->GetId(), 
214                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
215           else
216             curSDigit =  new AliEMCALDigit( -1               , 
217                                             -1               ,
218                                             hit->GetId(), 
219                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
220         }
221         else{
222           newsdigit = kFALSE;
223         }
224
225         Int_t check = 0 ;
226
227         for(check= 0; check < nSdigits ; check++) {
228           sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
229           if(curSDigit != 0){
230             if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ?              
231             *sdigit = *sdigit + *curSDigit;
232             newsdigit = kFALSE;
233             }
234           }
235         }
236         if (newsdigit) {
237           new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
238           nSdigits++ ;  
239         }
240         delete curSDigit ; 
241       }  // loop over all hits (hit = deposited energy/entering particle)
242     } // loop over iprim
243     sdigits->Sort() ;
244     
245     nSdigits = sdigits->GetEntriesFast() ;
246     fSDigitsInRun += nSdigits ;  
247     sdigits->Expand(nSdigits) ;
248
249     Int_t nPrimarymax = -1 ; 
250     Int_t i ;
251     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
252       AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
253       sdigit->SetIndexInList(i) ;
254     }
255     
256     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {   
257       if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > nPrimarymax)
258         nPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
259     }
260     
261     // Now write SDigits    
262     //First list of sdigits
263
264     Int_t bufferSize = 32000 ;    
265     TBranch * sdigitsBranch = treeS->Branch("EMCAL",&sdigits,bufferSize);
266  
267     sdigitsBranch->Fill() ;
268
269     gime->WriteSDigits("OVERWRITE");
270     
271     //NEXT - SDigitizer
272
273     gime->WriteSDigitizer("OVERWRITE");
274     
275     if(strstr(option,"deb"))
276       PrintSDigits(option) ;  
277   }
278
279   Unload();
280   
281   gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
282   
283   if(strstr(option,"tim")){
284     gBenchmark->Stop("EMCALSDigitizer"); 
285     printf("Exec: took %f seconds for SDigitizing %f seconds per event", 
286          gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ; 
287   }
288 }
289
290
291 //__________________________________________________________________
292 void AliEMCALSDigitizer::Print()const
293
294   // Prints parameters of SDigitizer
295   printf("Print: \n------------------- %s -------------", GetName() ) ; 
296   printf("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()) ;
297   printf("   with digitization parameters  A = %f\n", fA) ; 
298   printf("                                 B = %f\n", fB) ;
299   printf("   Threshold for EC Primary assignment= %f\n", fECPrimThreshold)  ; 
300   printf("---------------------------------------------------\n") ;
301
302 }
303
304 //__________________________________________________________________
305 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
306 {
307   // Equal operator.
308   // SDititizers are equal if their pedestal, slope and threshold are equal
309   if( (fA==sd.fA)&&(fB==sd.fB)&&
310       (fECPrimThreshold==sd.fECPrimThreshold))
311     return kTRUE ;
312   else
313     return kFALSE ;
314 }
315
316 //__________________________________________________________________
317 void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
318 {
319   //Prints list of digits produced at the current pass of AliEMCALDigitizer
320   
321   
322   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
323   const TClonesArray * sdigits = gime->SDigits() ; 
324   
325   printf("\n") ;  
326   printf("event %i", gAlice->GetEvNumber()) ;
327   printf("\n      Number of entries in SDigits list %i", sdigits->GetEntriesFast()); 
328   if(strstr(option,"all")||strstr(option,"EMC")){
329     
330     //loop over digits
331     AliEMCALDigit * digit;
332     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
333     Int_t index ;
334     char * tempo = new char[8192]; 
335     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
336       digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
337       sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d :",
338               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
339       printf(tempo); 
340       
341       Int_t iprimary;
342       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
343         sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
344         printf(tempo); 
345       }          
346     }
347     delete tempo ;
348   }
349 }
350
351 //____________________________________________________________________________ 
352 void AliEMCALSDigitizer::Unload() const
353 {
354   // Unload Hits and SDigits from the folder
355   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
356   AliEMCALLoader * loader = gime->EmcalLoader() ; 
357   loader->UnloadHits() ; 
358   loader->UnloadSDigits() ; 
359 }