removed iostream
[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 "TFile.h"
53 #include "TTask.h"
54 #include "TTree.h"
55 #include "TSystem.h"
56 #include "TROOT.h"
57 #include "TFolder.h"
58 #include "TBenchmark.h"
59 #include "TGeometry.h"
60
61 // --- Standard library ---
62
63 // --- AliRoot header files ---
64 #include "AliRun.h"
65 #include "AliHeader.h"
66 #include "AliEMCALDigit.h"
67 #include "AliEMCALGeometry.h"
68 #include "AliEMCALGetter.h"
69 #include "AliEMCALHit.h"
70 #include "AliEMCALSDigitizer.h"
71
72 ClassImp(AliEMCALSDigitizer)
73            
74 //____________________________________________________________________________ 
75   AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","") 
76 {
77   // ctor
78   InitParameters() ; 
79   fDefaultInit = kTRUE ; 
80 }
81
82 //____________________________________________________________________________ 
83 AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle, const Bool_t toSplit):
84 TTask(sDigitsTitle, headerFile)
85 {
86   // ctor
87   InitParameters() ; 
88   fToSplit = toSplit ;
89   Init();
90   fDefaultInit = kFALSE ; 
91 }
92
93 //____________________________________________________________________________ 
94 AliEMCALSDigitizer::~AliEMCALSDigitizer()
95 {
96   // dtor
97   fSplitFile = 0 ; 
98 }
99
100 //____________________________________________________________________________ 
101 void AliEMCALSDigitizer::Init(){
102   // Initialization: open root-file, allocate arrays for hits and sdigits,
103   // attach task SDigitizer to the list of EMCAL tasks
104   // 
105   // Initialization can not be done in the default constructor
106   //============================================================= YS
107   //  The initialisation is now done by the getter
108
109   if( strcmp(GetTitle(), "") == 0 )
110     SetTitle("galice.root") ;
111    
112   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), fToSplit) ; 
113   if ( gime == 0 ) {
114     Error("Init", "Could not obtain the Getter object !" ) ;  
115     return ;
116   } 
117   
118   gime->PostSDigits( GetName(), GetTitle() ) ; 
119   
120   fSplitFile = 0 ;
121   if(fToSplit){
122     // construct the name of the file as /path/EMCAL.SDigits.root
123     // First - extract full path if necessary
124     TString sDigitsFileName(GetTitle()) ;
125     Ssiz_t islash = sDigitsFileName.Last('/') ;
126     if(islash<sDigitsFileName.Length())
127       sDigitsFileName.Remove(islash+1,sDigitsFileName.Length()) ;
128     else
129       sDigitsFileName="" ;
130
131     // Next - append the file name 
132     sDigitsFileName+="EMCAL.SDigits." ;
133     if((strcmp(GetName(),"Default")!=0)&&(strcmp(GetName(),"")!=0)){
134       sDigitsFileName+=GetName() ;
135       sDigitsFileName+="." ;
136     }
137     sDigitsFileName+="root" ;
138
139     // Finally - check if the file already opened or open the file
140     fSplitFile = static_cast<TFile*>(gROOT->GetFile(sDigitsFileName.Data()));   
141     if(!fSplitFile)
142       fSplitFile =  TFile::Open(sDigitsFileName.Data(),"update") ;
143   }
144
145   TString sdname(GetName() );
146   sdname.Append(":") ;
147   sdname.Append(GetTitle() ) ;
148   SetName(sdname) ;
149   gime->PostSDigitizer(this) ;
150 }
151
152 //____________________________________________________________________________ 
153 void AliEMCALSDigitizer::InitParameters()
154 {
155   fA                      = 0;
156   fB                      = 10000000.;
157   fTowerPrimThreshold     = 0.01 ;
158   fPreShowerPrimThreshold = 0.0001 ; 
159   fPhotonElectronFactor   = 5000. ; // photoelectrons per GeV 
160   fSplitFile              = 0 ; 
161   fToSplit                = kFALSE ;
162 }
163
164 //____________________________________________________________________________
165 void AliEMCALSDigitizer::Exec(Option_t *option) 
166
167   // Collects all hits in the same active volume into digit
168
169   if( strcmp(GetName(), "") == 0 )
170     Init() ;
171   
172   if (strstr(option, "print") ) {
173     Print("") ; 
174     return ; 
175   }
176   
177   if(strstr(option,"tim"))
178     gBenchmark->Start("EMCALSDigitizer");
179
180   //Check, if this branch already exits
181   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
182   if(gime->BranchExists("SDigits") ) 
183     return;   
184
185   TString sdname(GetName()) ;
186   sdname.Remove(sdname.Index(GetTitle())-1) ;
187
188   Int_t nevents = gime->MaxEvent() ; 
189   Int_t ievent ;
190    
191     for(ievent = 0; ievent < nevents; ievent++){     
192       gime->Event(ievent,"H") ;
193      
194       const TClonesArray * hits = gime->Hits() ;
195    
196     TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
197     sdigits->Clear();
198     Int_t nSdigits = 0 ;
199
200     //Collects all hits in the same active volume into digit
201     
202     //Now make SDigits from hits, for EMCAL it is the same, so just copy    
203     Int_t nPrim =  static_cast<Int_t>((gAlice->TreeH())->GetEntries()) ; 
204     // Attention nPrim is the number of primaries tracked by Geant 
205     // and this number could be different to the number of Primaries in TreeK;
206     Int_t iprim ;
207       for ( iprim = 0 ; iprim < nPrim ; iprim++ ) { 
208         //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
209         gime->Track(iprim) ;
210         Int_t i;
211         for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
212           AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
213           AliEMCALDigit * curSDigit = 0 ;
214           AliEMCALDigit * sdigit = 0 ;
215           Bool_t newsdigit = kTRUE; 
216           // Assign primary number only if deposited energy is significant
217    
218           if( (!hit->IsInPreShower() && hit->GetEnergy() > fTowerPrimThreshold) || 
219               (hit->IsInPreShower() && hit->GetEnergy() > fPreShowerPrimThreshold)) 
220             curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
221                                             hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()), 
222                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;
223           else 
224             curSDigit =  new AliEMCALDigit( -1               , 
225                                             -1               ,
226                                             Layer2TowerID(hit->GetId(),hit->IsInPreShower()), 
227                                             Digitize(hit->GetEnergy()), hit->GetTime() ) ;      
228           Int_t check = 0 ;
229           for(check= 0; check < nSdigits ; check++) {
230             sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(check)) ;
231             if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same tower or the same preshower ?              
232               *sdigit = *sdigit + *curSDigit;
233               newsdigit = kFALSE;
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/layer/entering particle)
242       } // loop over iprim
243       
244       sdigits->Sort() ;
245       
246       nSdigits = sdigits->GetEntriesFast() ;
247       fSDigitsInRun += nSdigits ;  
248       sdigits->Expand(nSdigits) ;
249         
250       const AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
251       
252       if (nSdigits != 0 ) {
253         Int_t lastPreShowerIndex = nSdigits - 1 ;
254        
255         
256         if (!(dynamic_cast<AliEMCALDigit *>(sdigits->At(lastPreShowerIndex))->IsInPreShower()))
257           
258           lastPreShowerIndex = -2; 
259         
260         Int_t firstPreShowerIndex = 100000 ; 
261         Int_t index ; 
262         AliEMCALDigit * sdigit = 0 ;
263         
264         for ( index = 0; index < nSdigits ; index++) {    
265           sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ) ;
266           if (sdigit->IsInPreShower() ){ 
267             firstPreShowerIndex = index ;
268             break ;
269           }
270         }
271         
272         AliEMCALDigit * preshower ;
273         AliEMCALDigit * tower ;
274         Int_t lastIndex = lastPreShowerIndex +1 ; 
275         
276         for (index = firstPreShowerIndex ; index <= lastPreShowerIndex; index++) {  
277           preshower = dynamic_cast<AliEMCALDigit *>(sdigits->At(index) ); 
278           Bool_t towerFound = kFALSE ;
279           Int_t jndex ;
280           for (jndex = 0; jndex < firstPreShowerIndex; jndex++) {
281             tower  = dynamic_cast<AliEMCALDigit *>(sdigits->At(jndex) ); 
282             if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) {          
283               Float_t towerEnergy  = static_cast<Float_t>(tower->GetAmp()) ; 
284               Float_t preshoEnergy = static_cast<Float_t>(preshower->GetAmp()) ; 
285               towerEnergy +=preshoEnergy ; 
286               *tower = *tower + *preshower    ; // and add preshower multiplied by layer ratio to tower
287               tower->SetAmp(static_cast<Int_t>(TMath::Ceil(towerEnergy))) ; 
288               towerFound = kTRUE ;
289             }
290           }
291           if ( !towerFound ) {  
292             new((*sdigits)[lastIndex])  AliEMCALDigit(*preshower);
293             AliEMCALDigit * temp = dynamic_cast<AliEMCALDigit *>(sdigits->At(lastIndex)) ;
294             temp->SetId(temp->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) ;
295             lastIndex++ ; 
296           }
297         }
298         sdigits->Sort() ;
299         Int_t NPrimarymax = -1 ; 
300         Int_t i ;
301         for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
302           sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
303           sdigit->SetIndexInList(i) ;
304         }
305         
306         for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {   
307           if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > NPrimarymax)
308             NPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
309         }
310       }
311       
312       //Now write SDigits
313       
314       if(gAlice->TreeS() == 0 || (fSplitFile))  //<--- To be checked: we should not create TreeS if it is already here
315         gAlice->MakeTree("S",fSplitFile);
316      
317       if(fSplitFile)
318         fSplitFile->cd() ;
319       
320       //First list of sdigits
321       Int_t bufferSize = 32000 ;    
322       TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize);
323       sdigitsBranch->SetTitle(sdname);
324       
325       //NEXT - SDigitizer
326       Int_t splitlevel = 0 ;
327       AliEMCALSDigitizer * sd = this ;
328       TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer",
329                                                            &sd,bufferSize,splitlevel); 
330       sdigitizerBranch->SetTitle(sdname);
331       
332       sdigitsBranch->Fill() ; 
333       sdigitizerBranch->Fill() ; 
334       gAlice->TreeS()->AutoSave() ;
335
336       if(strstr(option,"deb"))
337         PrintSDigits(option) ;
338       
339   }
340
341   if(strstr(option,"tim")){
342     gBenchmark->Stop("EMCALSDigitizer"); 
343     Info("Exec", "took %f seconds for SDigitizing %f seconds per event", 
344          gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ;
345   }   
346 }
347
348 //__________________________________________________________________
349 void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){
350  
351   // Setting title to branch SDigits 
352
353   TString stitle(title) ;
354
355   // check if branch with title already exists
356   TBranch * sdigitsBranch    = 
357     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("EMCAL")) ; 
358   TBranch * sdigitizerBranch =  
359     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliEMCALSDigitizer")) ;
360   const char * sdigitsTitle    = sdigitsBranch ->GetTitle() ;  
361   const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
362   if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
363     Error("SetSDigitsBranch", "Cannot overwrite existing branch with title %s", title) ;
364     return ;
365   }
366   
367   Info("SetSDigitsBranch", "Changing SDigits file from %s to %s", GetName(), title) ;
368
369   SetName(title) ; 
370     
371   // Post to the WhiteBoard
372   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
373   gime->PostSDigits( title, GetTitle()) ; 
374 }
375
376
377 //__________________________________________________________________
378 void AliEMCALSDigitizer::Print(Option_t* option)const
379
380   // Prints parameters of SDigitizer
381
382   TString message("\n") ; 
383   message += "------------------- "; 
384   message += GetName() ; 
385   message += " -------------\n" ;
386   message += "   Writing SDigitis to branch with title  " ; 
387   message += GetName() ;
388   message += "\n   with digitization parameters  A               = " ; 
389   message += fA ;
390   message += "\n                                 B               = " ; 
391   message += fB ; 
392   message += "\n   Threshold for Primary assignment in Tower     = " ; 
393   message += fTowerPrimThreshold ; 
394   message += "\n   Threshold for Primary assignment in PreShower = " ; 
395   message += fPreShowerPrimThreshold ; 
396   message += "\n---------------------------------------------------" ;
397   
398   Info("Print", message.Data() ) ; 
399 }
400
401 //__________________________________________________________________
402 Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const
403 {
404   // Equal operator.
405   // SDititizers are equal if their pedestal, slope and threshold are equal
406
407   if( (fA==sd.fA)&&(fB==sd.fB)&&
408       (fTowerPrimThreshold==sd.fTowerPrimThreshold) &&
409       (fPreShowerPrimThreshold==sd.fPreShowerPrimThreshold))
410     return kTRUE ;
411   else
412     return kFALSE ;
413 }
414
415 //__________________________________________________________________
416 void AliEMCALSDigitizer::PrintSDigits(Option_t * option){
417   //Prints list of digits produced at the current pass of AliEMCALDigitizer
418   
419   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
420   TString sdname(GetName()) ;
421   sdname.Remove(sdname.Index(GetTitle())-1) ;
422   const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ; 
423   
424   TString message("\n") ;  
425   message += "event " ; 
426   message += gAlice->GetEvNumber() ;
427   message += "\n      Number of entries in SDigits list " ;
428   message += sdigits->GetEntriesFast() ; 
429   if(strstr(option,"all")||strstr(option,"EMC")){
430     
431     //loop over digits
432     AliEMCALDigit * digit;
433     message += "\nSDigit Id  Amplitude      Time      Index  Nprim  Primaries list \n" ;    
434     Int_t index ;
435     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
436       digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
437        message += digit->GetId() ; 
438        message += "   " ; 
439        message += digit->GetAmp() ; 
440        message += "    "  ; 
441        message += digit->GetTime() ; 
442        message += "    ";
443        message += digit->GetIndexInList() ; 
444        message += "    " ;  
445        message += digit->GetNprimary() ; 
446        message += " : ";
447       
448       Int_t iprimary;
449       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
450         message += digit->GetPrimary(iprimary+1) ; 
451         message += "  ";
452       }          
453     }
454   }
455   Info("PrintSDigits", message.Data() ) ; 
456 }
457
458 //________________________________________________________________________
459 const Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t ihit, Bool_t preshower)
460 {
461   // Method to Transform from Hit Id to Digit Id
462   // This function should be one to one
463   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
464   const AliEMCALGeometry * geom = gime->EMCALGeometry();
465   Int_t ieta  = ((ihit-1)/geom->GetNPhi())%geom->GetNZ(); // eta Tower Index
466   Int_t iphi = (ihit-1)%(geom->GetNPhi())+1; //phi Tower Index
467   Int_t it = -10;
468   Int_t ipre = 0;
469
470   if (preshower)ipre = 1;
471   if (iphi > 0 && ieta >= 0){
472     it = iphi+ieta*geom->GetNPhi() + ipre*geom->GetNPhi()*geom->GetNZ();
473     return it;
474   }else{
475     Error("Layer2TowerID", "there is an error: Eta number = %f Phi number = %f", ieta, iphi) ;
476     return it;
477   } // end if iphi>0 && ieta>0
478 }
479 //_______________________________________________________________________________________
480 // void AliEMCALSDigitizer::TestTowerID(void)
481 // {
482 //   Int_t j;
483
484 //   Bool_t preshower = kFALSE;
485 //   for (j = 0 ; j < 10 ; j++){  // loop over hit id
486 //     Int_t i;
487 //    for (i = 0 ; i <= 2 ; i++){  // loop over 
488 //      Int_t k = i*96*144+j*144+1;
489 //       Info("TestTowerID", " Hit Index = %d  %d   TOWERID = %d", k, j*10, Layer2TowerID(k, preshower) ) ;
490 //     }
491 //   }
492 // }
493
494 //____________________________________________________________________________ 
495 void AliEMCALSDigitizer::UseHitsFrom(const char * filename)
496 {
497   SetTitle(filename) ; 
498   Init() ; 
499 }