Coding Convention
[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   InitParameters() ; 
71   fDefaultInit = kTRUE ; 
72 }
73
74 //____________________________________________________________________________ 
75 AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, const char * eventFolderName):
76   TTask("EMCAL"+AliConfig::fgkSDigitizerTaskName, alirunFileName),
77   fEventFolderName(eventFolderName)
78 {
79   // ctor
80   Init();
81   InitParameters() ; 
82   fDefaultInit = kFALSE ; 
83 }
84
85
86 //____________________________________________________________________________ 
87 AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) : TTask(sd) {
88   //cpy ctor 
89
90   fA             = sd.fA ;
91   fB             = sd.fB ;
92   fPREPrimThreshold = sd.fPREPrimThreshold ;
93   fECPrimThreshold  = sd.fECPrimThreshold ;
94   fHCPrimThreshold  = sd.fHCPrimThreshold ;
95   fSDigitsInRun  = sd.fSDigitsInRun ;
96   SetName(sd.GetName()) ; 
97   SetTitle(sd.GetTitle()) ; 
98   fEventFolderName = sd.fEventFolderName;
99 }
100
101
102 //____________________________________________________________________________ 
103 void AliEMCALSDigitizer::Init(){
104   // Initialization: open root-file, allocate arrays for hits and sdigits,
105   // attach task SDigitizer to the list of EMCAL tasks
106   // 
107   // Initialization can not be done in the default constructor
108   //============================================================= YS
109   //  The initialisation is now done by the getter
110
111   fInit = kTRUE ; 
112    
113   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());  
114   if ( gime == 0 ) {
115     Fatal("Init", "Could not obtain the Getter objectfor file %s and event %s !", GetTitle(), fEventFolderName.Data()) ;  
116     return ;
117   } 
118   
119   TString opt("SDigits") ; 
120   if(gime->VersionExists(opt) ) { 
121     Error( "Init", "Give a version name different from %s", fEventFolderName.Data() ) ;
122     fInit = kFALSE ; 
123   }
124   
125   gime->PostSDigitizer(this);
126   gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
127   
128 }
129
130 //____________________________________________________________________________ 
131 void AliEMCALSDigitizer::InitParameters()
132 {
133   // Initializes parameters
134   fA                      = 0;
135   fB                      = 10000000.;
136
137   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
138   const AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
139   if (geom->GetSampling() == 0.) {
140     Error("InitParameters", "Sampling factor not set !") ; 
141     abort() ;
142   }
143   else
144     Info("InitParameters", "Sampling factor set to %f\n", geom->GetSampling()) ; 
145   
146   // this threshold corresponds approximately to 100 MeV
147   fECPrimThreshold     = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNECLayers() ;
148   fPREPrimThreshold    = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNPRLayers() ; 
149   fHCPrimThreshold     = fECPrimThreshold/5. ; // 5 is totally arbitrary
150
151 }
152
153 //____________________________________________________________________________
154 void AliEMCALSDigitizer::Exec(Option_t *option) 
155
156   // Collects all hits in the section (PRE/ECAL/HCAL) of the same tower into digit
157   
158   if (strstr(option, "print") ) {
159     Print() ; 
160     return ; 
161   }
162   
163   if(strstr(option,"tim"))
164     gBenchmark->Start("EMCALSDigitizer");
165
166   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
167  
168   //switch off reloading of this task while getting event
169   if (!fInit) { // to prevent overwrite existing file
170     Error( "Exec", "Give a version name different from %s", fEventFolderName.Data() ) ;
171     return ;
172   }
173   
174   Int_t nevents = gime->MaxEvent() ; 
175   Int_t ievent ;   
176   for(ievent = 0; ievent < nevents; ievent++){     
177   
178     gime->Event(ievent,"H") ;  
179
180     TTree * treeS = gime->TreeS(); 
181     TClonesArray * hits = gime->Hits() ; 
182     TClonesArray * sdigits = gime->SDigits() ;
183     sdigits->Clear();
184     Int_t nSdigits = 0 ;
185     //Now make SDigits from hits, for EMCAL it is the same, so just copy    
186     Int_t nPrim =  static_cast<Int_t>((gime->TreeH())->GetEntries()) ; 
187     // Attention nPrim is the number of primaries tracked by Geant 
188     // and this number could be different to the number of Primaries in TreeK;
189     Int_t iprim ;
190
191     for ( iprim = 0 ; iprim < nPrim ; iprim++ ) { 
192       //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
193       gime->Track(iprim) ;
194       Int_t i;
195       for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
196         AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
197         AliEMCALDigit * curSDigit = 0 ;
198         AliEMCALDigit * sdigit = 0 ;
199         Bool_t newsdigit = kTRUE; 
200
201         // Assign primary number only if deposited energy is significant
202
203         AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
204
205         if( geom->IsInPRE(hit->GetId()) )  
206           if( hit->GetEnergy() > fPREPrimThreshold )
207             curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
208                                             hit->GetIparent(), hit->GetId(), 
209                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
210           else
211             curSDigit =  new AliEMCALDigit( -1               , 
212                                             -1               ,
213                                             hit->GetId(), 
214                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
215         else if( geom->IsInECA(hit->GetId()) )
216           if( hit->GetEnergy() >  fECPrimThreshold )
217             curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
218                                             hit->GetIparent(), hit->GetId(), 
219                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
220           else
221             curSDigit =  new AliEMCALDigit( -1               , 
222                                             -1               ,
223                                             hit->GetId(), 
224                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
225         else if( geom->IsInHCA(hit->GetId()) )
226           if( hit->GetEnergy() >  fHCPrimThreshold )
227             
228             curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
229                                             hit->GetIparent(), hit->GetId(), 
230                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
231           else
232             curSDigit =  new AliEMCALDigit( -1               , 
233                                             -1               ,
234                                             hit->GetId(), 
235                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
236         
237         Int_t check = 0 ;
238         for(check= 0; check < nSdigits ; check++) {
239           sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
240           if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL/HCAL/preshower tower ?              
241             *sdigit = *sdigit + *curSDigit;
242             newsdigit = kFALSE;
243           }
244         }
245         if (newsdigit) { 
246           new((*sdigits)[nSdigits])  AliEMCALDigit(*curSDigit);
247           nSdigits++ ;  
248         }
249         delete curSDigit ; 
250       }  // loop over all hits (hit = deposited energy/layer/entering particle)
251     } // loop over iprim
252     
253     sdigits->Sort() ;
254     
255     nSdigits = sdigits->GetEntriesFast() ;
256     fSDigitsInRun += nSdigits ;  
257     sdigits->Expand(nSdigits) ;
258
259     Int_t nPrimarymax = -1 ; 
260     Int_t i ;
261     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
262       AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
263       sdigit->SetIndexInList(i) ;
264     }
265     
266     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {   
267       if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > nPrimarymax)
268         nPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
269     }
270     
271     // Now write SDigits    
272     
273     //First list of sdigits
274
275     Int_t bufferSize = 32000 ;    
276     TBranch * sdigitsBranch = treeS->Branch("EMCAL",&sdigits,bufferSize);
277  
278     sdigitsBranch->Fill() ;
279
280     gime->WriteSDigits("OVERWRITE");
281     
282     //NEXT - SDigitizer
283
284     gime->WriteSDigitizer("OVERWRITE");
285     
286     if(strstr(option,"deb"))
287       PrintSDigits(option) ;  
288   }
289
290   Unload();
291   
292   gime->EmcalLoader()->GetSDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
293   
294   if(strstr(option,"tim")){
295     gBenchmark->Stop("EMCALSDigitizer"); 
296     Info("Exec", "took %f seconds for SDigitizing %f seconds per event", 
297          gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ; 
298   }
299 }
300
301
302 //__________________________________________________________________
303 void AliEMCALSDigitizer::Print()const
304
305   // Prints parameters of SDigitizer
306   Info("Print", "\n------------------- %s -------------", GetName() ) ; 
307   printf("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()) ;
308   printf("   with digitization parameters  A = %f\n", fA) ; 
309   printf("                                 B = %f\n", fB) ;
310   printf("   Threshold for PRE Primary assignment= %f\n", fPREPrimThreshold)  ; 
311   printf("   Threshold for EC Primary assignment= %f\n", fECPrimThreshold)  ; 
312   printf("   Threshold for HC Primary assignment= %f\n", fHCPrimThreshold)  ; 
313   printf("---------------------------------------------------\n") ;
314
315 }
316
317 //__________________________________________________________________
318 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
319 {
320   // Equal operator.
321   // SDititizers are equal if their pedestal, slope and threshold are equal
322
323   if( (fA==sd.fA)&&(fB==sd.fB)&&
324       (fECPrimThreshold==sd.fECPrimThreshold) &&
325       (fHCPrimThreshold==sd.fHCPrimThreshold) &&
326       (fPREPrimThreshold==sd.fPREPrimThreshold))
327     return kTRUE ;
328   else
329     return kFALSE ;
330 }
331
332 //__________________________________________________________________
333 void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
334 {
335   //Prints list of digits produced at the current pass of AliEMCALDigitizer
336   
337   
338   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
339   const TClonesArray * sdigits = gime->SDigits() ; 
340   
341   TString message("\n") ;  
342   message += "event " ; 
343   message += gAlice->GetEvNumber() ;
344   message += "\n      Number of entries in SDigits list " ;
345   message += sdigits->GetEntriesFast() ; 
346   if(strstr(option,"all")||strstr(option,"EMC")){
347     
348     //loop over digits
349     AliEMCALDigit * digit;
350     message += "\n   Id  Amplitude    Time          Index Nprim: Primaries list \n" ;    
351     Int_t index ;
352     char * tempo = new char[8192]; 
353     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
354       digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
355       sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d :",
356               digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
357       message += tempo ; 
358       
359       Int_t iprimary;
360       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
361         sprintf(tempo, "%d ",digit->GetPrimary(iprimary+1) ) ; 
362         message += tempo ; 
363       }          
364     }
365     delete tempo ;
366   }
367   Info("PrintSDigits", message.Data() ) ; 
368 }
369
370 //____________________________________________________________________________ 
371 void AliEMCALSDigitizer::Unload() const
372 {
373   // Unload Hits and SDigits from the folder
374   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
375   AliEMCALLoader * loader = gime->EmcalLoader() ; 
376   loader->UnloadHits() ; 
377   loader->UnloadSDigits() ; 
378 }