6cd2cea2a933bf5513f454c2239ea851b2d6fd49
[u/mrichter/AliRoot.git] / PHOS / AliPHOSSDigitizer.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 primary in one active cell
22 // A threshold for assignment of the primary to SDigit is applied 
23 // SDigits are written to TreeS, branch "PHOS"
24 // AliPHOSSDigitizer with all current parameters is written 
25 // to TreeS branch "AliPHOSSDigitizer".
26 // Both branches, "PHOS" and "AliPHOSSDigitizer", are written to the same
27 // file, and therefore, changing branch file name one can produce several
28 // versions of SDigitization from the same hits.
29 // 
30 //
31 //*-- Author :  Dmitri Peressounko (SUBATECH & KI) 
32 //////////////////////////////////////////////////////////////////////////////
33
34 // --- ROOT system ---
35 #include "TFile.h"
36 #include "TTask.h"
37 #include "TTree.h"
38 #include "TSystem.h"
39 #include "TROOT.h"
40 #include "TFolder.h"
41 #include "TBenchmark.h"
42 // --- Standard library ---
43 #include <iomanip.h>
44
45 // --- AliRoot header files ---
46 #include "AliRun.h"
47 #include "AliPHOSDigit.h"
48 #include "AliPHOSHit.h"
49 #include "AliPHOSSDigitizer.h"
50
51
52 ClassImp(AliPHOSSDigitizer)
53
54            
55 //____________________________________________________________________________ 
56   AliPHOSSDigitizer::AliPHOSSDigitizer():TTask("AliPHOSSDigitizer","") 
57 {
58   // ctor
59   fA = 0;
60   fB = 10000000. ;
61   fPrimThreshold = 0.01 ;
62   fNevents = 0 ;     
63   fSDigits = 0 ;
64   fHits = 0 ;
65   fIsInitialized = kFALSE ;
66
67 }
68
69 //____________________________________________________________________________ 
70 AliPHOSSDigitizer::AliPHOSSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliPHOSSDigitizer","")
71 {
72   // ctor
73   fA = 0;
74   fB = 10000000.;
75   fPrimThreshold = 0.01 ;
76   fNevents = 0 ;      
77   fSDigitsTitle = sDigitsTitle ;
78   fHeadersFile = headerFile ;
79   fSDigits = new TClonesArray("AliPHOSDigit",1000);
80   fHits    = new TClonesArray("AliPHOSHit",1000);
81
82   TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
83   
84   //File was not opened yet
85   if(file == 0){
86     file = new TFile(fHeadersFile.Data(),"update") ;
87     gAlice = (AliRun *) file->Get("gAlice") ;
88   }
89   
90   //add Task to //root/Tasks folder
91   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
92   roottasks->Add(this) ; 
93   
94   fIsInitialized = kTRUE ;
95 }
96
97 //____________________________________________________________________________ 
98 AliPHOSSDigitizer::~AliPHOSSDigitizer()
99 {
100   // dtor
101   if(fSDigits)
102     delete fSDigits ;
103   if(fHits)
104     delete fHits ;
105 }
106 //____________________________________________________________________________ 
107 void AliPHOSSDigitizer::Init(){
108   //Initialization can not be done in the default constructor
109
110   if(!fIsInitialized){
111
112     if(fHeadersFile.IsNull())
113       fHeadersFile="galice.root" ;
114
115     TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ;
116     
117     //if file was not opened yet, read gAlice
118     if(file == 0){
119       file = new TFile(fHeadersFile.Data(),"update") ;
120       gAlice = (AliRun *) file->Get("gAlice") ;
121     }
122     
123     fHits    = new TClonesArray("AliPHOSHit",1000);
124     fSDigits = new TClonesArray("AliPHOSDigit",1000);
125     
126     // add Task to //root/Tasks folder
127     TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
128     roottasks->Add(this) ; 
129     
130     fIsInitialized = kTRUE ;
131   }
132 }
133 //____________________________________________________________________________
134 void AliPHOSSDigitizer::Exec(Option_t *option) { 
135   //Collects all hits in the same active volume into digit
136   
137   if(!fIsInitialized)
138     Init() ;
139
140   if(strstr(option,"tim"))
141     gBenchmark->Start("PHOSSDigitizer");
142   
143   fNevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
144   
145   Int_t ievent ;
146   for(ievent = 0; ievent < fNevents; ievent++){
147     gAlice->GetEvent(ievent) ;
148     gAlice->SetEvent(ievent) ;
149     
150     if(gAlice->TreeH()==0){
151       cout << "AliPHOSSDigitizer: There is no Hit Tree" << endl;
152       return ;
153     }
154     
155     //set address of the hits 
156     TBranch * branch = gAlice->TreeH()->GetBranch("PHOS");
157     if (branch) 
158       branch->SetAddress(&fHits);
159     else{
160       cout << "ERROR in AliPHOSSDigitizer: "<< endl ;
161       cout << "      no branch PHOS in TreeH"<< endl ;
162       cout << "      do nothing " << endl ;
163       return ;
164     }
165     
166     fSDigits->Clear();
167     Int_t nSdigits = 0 ;
168     
169     
170     //Now made SDigits from hits, for PHOS it is the same, so just copy    
171     Int_t itrack ;
172     for (itrack=0; itrack<gAlice->GetNtrack(); itrack++){
173       
174       //=========== Get the Hits Tree for the Primary track itrack
175       gAlice->ResetHits();    
176       gAlice->TreeH()->GetEvent(itrack);
177       
178       Int_t i;
179       for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
180         AliPHOSHit * hit = (AliPHOSHit*)fHits->At(i) ;
181         AliPHOSDigit * newdigit ;
182
183         // Assign primary number only if contribution is significant
184         if( hit->GetEnergy() > fPrimThreshold)
185           newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
186         else
187           newdigit = new AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
188         
189         new((*fSDigits)[nSdigits]) AliPHOSDigit(* newdigit) ;
190         nSdigits++ ;  
191         
192         delete newdigit ;    
193       } 
194       
195     } // loop over tracks
196     
197     fSDigits->Sort() ;
198     
199     nSdigits = fSDigits->GetEntriesFast() ;
200     fSDigits->Expand(nSdigits) ;
201     
202     Int_t i ;
203     for (i = 0 ; i < nSdigits ; i++) { 
204       AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(i) ; 
205       digit->SetIndexInList(i) ;     
206     }
207
208     if(gAlice->TreeS() == 0)
209       gAlice->MakeTree("S") ;
210     
211     //check, if this branch already exits?
212     TBranch * sdigitsBranch = 0;
213     TBranch * sdigitizerBranch = 0;
214     
215     TObjArray * branches = gAlice->TreeS()->GetListOfBranches() ;
216     Int_t ibranch;
217     Bool_t phosNotFound = kTRUE ;
218     Bool_t sdigitizerNotFound = kTRUE ;
219     
220     for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
221       
222       if(phosNotFound){
223         sdigitsBranch=(TBranch *) branches->At(ibranch) ;
224         if( (strcmp("PHOS",sdigitsBranch->GetName())==0 ) &&
225             (fSDigitsTitle.CompareTo(sdigitsBranch->GetTitle()) == 0) )
226           phosNotFound = kFALSE ;
227       }
228       if(sdigitizerNotFound){
229         sdigitizerBranch = (TBranch *) branches->At(ibranch) ;
230         if( (strcmp(sdigitizerBranch->GetName(),"AliPHOSSDigitizer") == 0)&&
231             (fSDigitsTitle.CompareTo(sdigitizerBranch->GetTitle()) == 0) )
232           sdigitizerNotFound = kFALSE ;
233       }
234     }
235
236     if(!(sdigitizerNotFound && phosNotFound)){
237       cout << "AliPHOSSdigitizer error:" << endl ;
238       cout << "Can not overwrite existing branches: do not write" << endl ;
239       return ;
240     }
241     
242     //Make (if necessary) branches    
243     char * file =0;
244     if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
245       file = new char[strlen(gAlice->GetBaseFile())+20] ;
246       sprintf(file,"%s/PHOS.SDigits.root",gAlice->GetBaseFile()) ;
247     }
248     
249     TDirectory *cwd = gDirectory;
250     
251     //First list of sdigits
252     Int_t bufferSize = 32000 ;    
253     sdigitsBranch = gAlice->TreeS()->Branch("PHOS",&fSDigits,bufferSize);
254     sdigitsBranch->SetTitle(fSDigitsTitle.Data());
255     if (file) {
256       sdigitsBranch->SetFile(file);
257       TIter next( sdigitsBranch->GetListOfBranches());
258       while ((sdigitsBranch=(TBranch*)next())) {
259         sdigitsBranch->SetFile(file);
260       }   
261       cwd->cd();
262     } 
263       
264     //second - SDigitizer
265     Int_t splitlevel = 0 ;
266     AliPHOSSDigitizer * sd = this ;
267     sdigitizerBranch = gAlice->TreeS()->Branch("AliPHOSSDigitizer","AliPHOSSDigitizer",
268                                                &sd,bufferSize,splitlevel); 
269     sdigitizerBranch->SetTitle(fSDigitsTitle.Data());
270     if (file) {
271       sdigitizerBranch->SetFile(file);
272       TIter next( sdigitizerBranch->GetListOfBranches());
273       while ((sdigitizerBranch=(TBranch*)next())) {
274         sdigitizerBranch->SetFile(file);
275       }   
276       cwd->cd();
277       delete file;
278     }
279
280     gAlice->TreeS()->Fill() ;
281     gAlice->TreeS()->Write(0,TObject::kOverwrite) ;
282     
283     if(strstr(option,"deb"))
284       PrintSDigits(option) ;
285     
286   }
287   
288   if(strstr(option,"tim")){
289     gBenchmark->Stop("PHOSSDigitizer");
290     cout << "AliPHOSSDigitizer:" << endl ;
291     cout << "   took " << gBenchmark->GetCpuTime("PHOSSDigitizer") << " seconds for SDigitizing " 
292          <<  gBenchmark->GetCpuTime("PHOSSDigitizer")/fNevents << " seconds per event " << endl ;
293     cout << endl ;
294   }
295   
296   
297 }
298 //__________________________________________________________________
299 void AliPHOSSDigitizer::SetSDigitsBranch(const char * title ){
300   //Seting title to branch SDigits 
301   if(!fSDigitsTitle.IsNull())
302     cout << "AliPHOSSdigitizer: changing SDigits file from " <<fSDigitsTitle.Data() << " to " << title << endl ;
303   fSDigitsTitle=title ;
304 }
305 //__________________________________________________________________
306 void AliPHOSSDigitizer::Print(Option_t* option)const{
307   cout << "------------------- "<< GetName() << " -------------" << endl ;
308   cout << "   Writing SDigitis to branch with title  " << fSDigitsTitle.Data() << endl ;
309   cout << "   with digitization parameters  A = " << fA << endl ;
310   cout << "                                 B = " << fB << endl ;
311   cout << "   Threshold for Primary assignment= " << fPrimThreshold << endl ; 
312   cout << "---------------------------------------------------"<<endl ;
313   
314 }
315 //__________________________________________________________________
316 Bool_t AliPHOSSDigitizer::operator==( AliPHOSSDigitizer const &sd )const{
317   if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
318     return kTRUE ;
319   else
320     return kFALSE ;
321 }
322 //__________________________________________________________________
323 void AliPHOSSDigitizer::PrintSDigits(Option_t * option){
324   //Prints list of digits produced at the current pass of AliPHOSDigitizer
325   
326   cout << "AliPHOSSDigitizer: " << endl ;
327   cout << "       Number of entries in SDigits list  " << fSDigits->GetEntriesFast() << endl ;
328   cout << endl ;
329   
330   if(strstr(option,"all")){// print all digits
331     
332     //loop over digits
333     AliPHOSDigit * digit;
334     cout << "SDigit Id " << " Amplitude " <<  " Index "  <<  " Nprim " << " Primaries list " <<  endl;    
335     Int_t index ;
336     for (index = 0 ; index < fSDigits->GetEntries() ; index++) {
337       digit = (AliPHOSDigit * )  fSDigits->At(index) ;
338       cout << setw(8)  <<  digit->GetId() << " "  <<    setw(3)  <<  digit->GetAmp() <<   "  "  
339            << setw(6)  <<  digit->GetIndexInList() << "  "   
340            << setw(5)  <<  digit->GetNprimary() <<"  ";
341       
342       Int_t iprimary;
343       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
344         cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << "  ";
345       cout << endl;      
346     }
347     
348   }
349 }