Added a protection in the dtor. When the tasks is created by default ctor (to access...
[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
17 /* $Id$ */
18
19 //_________________________________________________________________________
20 // This is a TTask that makes SDigits out of Hits
21 // The name of the TTask is also the title of the branch that will contain 
22 // the created SDigits
23 // The title of the TTAsk is the name of the file that contains the hits from
24 // which the SDigits are created
25 // A Summable Digits is the sum of all hits originating 
26 // from one primary in one active cell
27 // A threshold for assignment of the primary to SDigit is applied 
28 // SDigits are written to TreeS, branch "PHOS"
29 // AliPHOSSDigitizer with all current parameters is written 
30 // to TreeS branch "AliPHOSSDigitizer".
31 // Both branches have the same title. If necessary one can produce 
32 // another set of SDigits with different parameters. Two versions
33 // can be distunguished using titles of the branches.
34 // User case:
35 //  root [0] AliPHOSSDigitizer * s = new AliPHOSSDigitizer("galice.root")
36 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
37 //  root [1] s->ExecuteTask()
38 //             // Makes SDigitis for all events stored in galice.root
39 //  root [2] s->SetPedestalParameter(0.001)
40 //             // One can change parameters of digitization
41 // root [3] s->SetSDigitsBranch("Pedestal 0.001")
42 //             // and write them into the new branch
43 // root [4] s->ExecuteTask("deb all tim")
44 //             // available parameters:
45 //             deb - print # of produced SDigitis
46 //             deb all  - print # and list of produced SDigits
47 //             tim - print benchmarking information
48 //
49 //*-- Author :  Dmitri Peressounko (SUBATECH & KI) 
50 //////////////////////////////////////////////////////////////////////////////
51
52
53 // --- ROOT system ---
54 #include "TFile.h"
55 #include "TTask.h"
56 #include "TTree.h"
57 #include "TSystem.h"
58 #include "TROOT.h"
59 #include "TFolder.h"
60 #include "TBenchmark.h"
61 #include "TGeometry.h"
62
63 // --- Standard library ---
64 #include <iomanip.h>
65
66 // --- AliRoot header files ---
67 #include "AliRun.h"
68 #include "AliHeader.h"
69 #include "AliPHOSDigit.h"
70 #include "AliPHOSGeometry.h"
71 #include "AliPHOSGetter.h"
72 #include "AliPHOSHit.h"
73 #include "AliPHOSSDigitizer.h"
74
75
76 ClassImp(AliPHOSSDigitizer)
77
78            
79 //____________________________________________________________________________ 
80   AliPHOSSDigitizer::AliPHOSSDigitizer():TTask("","") {
81   // ctor
82   InitParameters() ; 
83 }
84
85 //____________________________________________________________________________ 
86 AliPHOSSDigitizer::AliPHOSSDigitizer(const char * headerFile, const char * sDigitsTitle):TTask(sDigitsTitle, headerFile)
87 {
88   // ctor
89   InitParameters() ; 
90   Init();
91 }
92
93 //____________________________________________________________________________ 
94 AliPHOSSDigitizer::~AliPHOSSDigitizer()
95 {
96   // dtor
97   // gime=0 if Digitizer created by default ctor (to get just the parameters)
98
99   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
100   
101   if (gime) {
102     // remove the task from the folder list
103     gime->RemoveTask("S",GetName()) ;
104     
105     TString name(GetName()) ; 
106     if (! name.IsNull() ) 
107       name.Remove(name.Index(":")) ; 
108     
109     // remove the Hits from the folder list
110     gime->RemoveObjects("H",name) ;
111     
112     // remove the SDigits from the folder list
113     gime->RemoveObjects("S", name) ;
114     
115     // Delete gAlice
116     gime->CloseFile() ; 
117     
118   }
119   fSplitFile = 0 ; 
120 }
121
122 //____________________________________________________________________________ 
123 void AliPHOSSDigitizer::Init()
124 {
125   // Initialization: open root-file, allocate arrays for hits and sdigits,
126   // attach task SDigitizer to the list of PHOS tasks
127   // 
128   // Initialization can not be done in the default constructor
129   //============================================================= YS
130   //  The initialisation is now done by AliPHOSGetter
131   
132   if( strcmp(GetTitle(), "") == 0 )
133     SetTitle("galice.root") ;
134
135   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), GetName()) ;     
136   if ( gime == 0 ) {
137     cerr << "ERROR: AliPHOSSDigitizer::Init -> Could not obtain the Getter object !" << endl ; 
138     return ;
139   } 
140   
141   gime->PostSDigits( GetName(), GetTitle() ) ; 
142
143   TString sdname(GetName() );
144   sdname.Append(":") ;
145   sdname.Append(GetTitle() ) ;
146   SetName(sdname) ;
147   gime->PostSDigitizer(this) ;
148 }
149
150 //____________________________________________________________________________ 
151 void AliPHOSSDigitizer::InitParameters()
152
153   fA             = 0;
154   fB             = 10000000.;
155   fPrimThreshold = 0.01 ;
156   fSDigitsInRun  = 0 ;
157   fSplitFile     = 0 ; 
158 }
159
160 //____________________________________________________________________________
161 void AliPHOSSDigitizer::Exec(Option_t *option) 
162
163   // Collects all hits in the same active volume into digit
164   if( strcmp(GetName(), "") == 0 )
165     Init() ;
166   
167   if (strstr(option, "print") ) {
168     Print("") ; 
169     return ; 
170   }
171
172   if(strstr(option,"tim"))
173     gBenchmark->Start("PHOSSDigitizer");
174
175   //Check, if this branch already exits
176   gAlice->GetEvent(0) ;
177
178   TString sdname(GetName()) ;
179   sdname.Remove(sdname.Index(GetTitle())-1) ;
180   
181   if(gAlice->TreeS() ) {
182     TObjArray * lob = static_cast<TObjArray*>(gAlice->TreeS()->GetListOfBranches()) ;
183     TIter next(lob) ; 
184     TBranch * branch = 0 ;  
185     Bool_t phosfound = kFALSE, sdigitizerfound = kFALSE ; 
186     
187     while ( (branch = (static_cast<TBranch*>(next()))) && (!phosfound || !sdigitizerfound) ) {
188       TString thisName( GetName() ) ; 
189       TString branchName( branch->GetTitle() ) ; 
190       branchName.Append(":") ; 
191       if ( (strcmp(branch->GetName(), "PHOS")==0) && thisName.BeginsWith(branchName) )  
192         phosfound = kTRUE ;
193       else if ( (strcmp(branch->GetName(), "AliPHOSSDigitizer")==0) &&  thisName.BeginsWith(branchName) )
194         sdigitizerfound = kTRUE ; 
195     }
196     
197     if ( phosfound || sdigitizerfound ) {
198       cerr << "WARNING: AliPHOSSDigitizer::Exec -> SDigits and/or SDigitizer branch with name " << GetName() 
199            << " already exits" << endl ;
200       return ; 
201     }   
202   }  
203
204   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
205
206   Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; 
207   
208   Int_t ievent ;
209   for(ievent = 0; ievent < nevents; ievent++){
210     gime->Event(ievent,"H") ;
211     const TClonesArray * hits = gime->Hits() ;
212     TClonesArray * sdigits = gime->SDigits(sdname.Data()) ;
213     sdigits->Clear();
214     Int_t nSdigits = 0 ;
215     
216     
217     //Now make SDigits from hits, for PHOS it is the same, so just copy    
218
219     Int_t Nprim =  (Int_t)  (gAlice->TreeH())->GetEntries(); 
220     // Attention Nprim is the number of primaries tracked by Geant and this number could be different to the number of Primaries in TreeK;
221     Int_t iprim;
222     for (iprim=0; iprim<Nprim; iprim++) { 
223       //=========== Get the PHOS branch from Hits Tree for the Primary iprim
224       gime->Track(iprim) ;
225       Int_t i;
226       for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
227         AliPHOSHit * hit = dynamic_cast<AliPHOSHit *>(hits->At(i)) ;
228         // Assign primary number only if contribution is significant
229         
230         if( hit->GetEnergy() > fPrimThreshold)
231           new((*sdigits)[nSdigits]) AliPHOSDigit(hit->GetPrimary(),hit->GetId(),
232                                                   Digitize(hit->GetEnergy()), hit->GetTime()) ;
233         else
234           new((*sdigits)[nSdigits]) AliPHOSDigit( -1              , hit->GetId(), 
235                                                    Digitize(hit->GetEnergy()), hit->GetTime()) ;
236         nSdigits++ ;    
237         
238       }
239
240     } // loop over iprim
241     
242     sdigits->Sort() ;
243     
244     nSdigits = sdigits->GetEntriesFast() ;
245     fSDigitsInRun += nSdigits ;  
246     sdigits->Expand(nSdigits) ;
247
248     Int_t i ;
249     for (i = 0 ; i < nSdigits ; i++) { 
250       AliPHOSDigit * digit = dynamic_cast<AliPHOSDigit *>(sdigits->At(i)) ; 
251       digit->SetIndexInList(i) ;     
252     }
253
254     if(gAlice->TreeS() == 0)
255       gAlice->MakeTree("S", fSplitFile);
256     
257     //First list of sdigits
258     Int_t bufferSize = 32000 ;    
259     TBranch * sdigitsBranch = gAlice->TreeS()->Branch("PHOS",&sdigits,bufferSize);
260     sdigitsBranch->SetTitle(sdname);
261       
262     //Next - SDigitizer
263     Int_t splitlevel = 0 ;
264     AliPHOSSDigitizer * sd = this ;
265     TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliPHOSSDigitizer","AliPHOSSDigitizer",
266                                                &sd,bufferSize,splitlevel); 
267     sdigitizerBranch->SetTitle(sdname);
268   
269     sdigitsBranch->Fill() ; 
270     sdigitizerBranch->Fill() ;
271     gAlice->TreeS()->AutoSave() ;
272     
273     if(strstr(option,"deb"))
274       PrintSDigits(option) ;
275   
276   }
277
278   if(strstr(option,"tim")){
279     gBenchmark->Stop("PHOSSDigitizer");
280     cout << "AliPHOSSDigitizer:" << endl ;
281     cout << "   took " << gBenchmark->GetCpuTime("PHOSSDigitizer") << " seconds for SDigitizing " 
282          <<  gBenchmark->GetCpuTime("PHOSSDigitizer")/nevents << " seconds per event " << endl ;
283     cout << endl ;
284   }
285   
286   
287 }
288
289 //__________________________________________________________________
290 void AliPHOSSDigitizer::SetSDigitsBranch(const char * title )
291 {
292   // Setting title to branch SDigits 
293
294   TString stitle(title) ;
295
296   // check if branch with title already exists
297   TBranch * sdigitsBranch    = 
298     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("PHOS")) ; 
299   TBranch * sdigitizerBranch =  
300     static_cast<TBranch*>(gAlice->TreeS()->GetListOfBranches()->FindObject("AliPHOSSDigitizer")) ;
301   const char * sdigitsTitle    = sdigitsBranch ->GetTitle() ;  
302   const char * sdigitizerTitle = sdigitizerBranch ->GetTitle() ;
303   if ( stitle.CompareTo(sdigitsTitle)==0 || stitle.CompareTo(sdigitizerTitle)==0 ){
304     cerr << "ERROR: AliPHOSSdigitizer::SetSDigitsBranch -> Cannot overwrite existing branch with title " << title << endl ;
305     return ;
306   }
307   
308   cout << "AliPHOSSdigitizer::SetSDigitsBranch -> Changing SDigits file from " << GetName() << " to " << title << endl ;
309
310   SetName(title) ; 
311     
312   // Post to the WhiteBoard
313   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
314   gime->PostSDigits( title, GetTitle()) ; 
315 }
316
317 //__________________________________________________________________
318 void AliPHOSSDigitizer::SetSplitFile(const TString splitFileName) 
319 {
320   // Diverts the SDigits in a file separate from the hits file
321   
322   TDirectory * cwd = gDirectory ;
323   
324   if ( !(gAlice->GetTreeSFileName() == splitFileName) ) {
325     if (gAlice->GetTreeSFile() ) 
326       gAlice->GetTreeSFile()->Close() ; 
327   }
328   
329   fSplitFile = gAlice->InitTreeFile("S",splitFileName.Data());
330   fSplitFile->cd() ; 
331   gAlice->Write(0, TObject::kOverwrite);
332   
333   TTree *treeE  = gAlice->TreeE();
334   if (!treeE) {
335     cerr << "ERROR: AliPHOSSDigitizer::SetSPlitFile -> No TreeE found "<<endl;
336     abort() ;
337   }      
338   
339   // copy TreeE
340     AliHeader *header = new AliHeader();
341     treeE->SetBranchAddress("Header", &header);
342     treeE->SetBranchStatus("*",1);
343     TTree *treeENew =  treeE->CloneTree();
344     treeENew->Write(0, TObject::kOverwrite);
345   
346
347   // copy AliceGeom
348     TGeometry *AliceGeom = static_cast<TGeometry*>(cwd->Get("AliceGeom"));
349     if (!AliceGeom) {
350       cerr << "ERROR: AliPHOSSDigitizer::SetSPlitFile -> AliceGeom was not found in the input file "<<endl;
351       abort() ;
352     }
353     AliceGeom->Write(0, TObject::kOverwrite);
354
355   gAlice->MakeTree("S",fSplitFile);
356   cwd->cd() ; 
357   cout << "INFO: AliPHOSSDigitizer::SetSPlitMode -> SDigits will be stored in " << splitFileName.Data() << endl ; 
358
359 }
360
361 //__________________________________________________________________
362 void AliPHOSSDigitizer::Print(Option_t* option)const
363 {
364   // Prints parameters of SDigitizer
365   cout << "------------------- "<< GetName() << " -------------" << endl ;
366   cout << "   Writing SDigits to branch with title  " << GetName() << endl ;
367   cout << "   with digitization parameters  A = " << fA << endl ;
368   cout << "                                 B = " << fB << endl ;
369   cout << "   Threshold for Primary assignment= " << fPrimThreshold << endl ; 
370   cout << "---------------------------------------------------"<<endl ;
371   
372 }
373 //__________________________________________________________________
374 Bool_t AliPHOSSDigitizer::operator==( AliPHOSSDigitizer const &sd )const
375 {
376   // Equal operator.
377   // SDititizers are equal if their pedestal, slope and threshold are equal
378
379   if( (fA==sd.fA)&&(fB==sd.fB)&&(fPrimThreshold==sd.fPrimThreshold))
380     return kTRUE ;
381   else
382     return kFALSE ;
383 }
384 //__________________________________________________________________
385 //__________________________________________________________________
386 void AliPHOSSDigitizer::PrintSDigits(Option_t * option)
387 {
388   // Prints list of digits produced in the current pass of AliPHOSDigitizer
389
390
391   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
392   TString sdname(GetName()) ;
393   sdname.Remove(sdname.Index(GetTitle())-1) ;
394   const TClonesArray * sdigits = gime->SDigits(sdname.Data()) ; 
395
396   cout << "AliPHOSSDigitiser: event " << gAlice->GetEvNumber() << endl ;
397   cout << "      Number of entries in SDigits list " << sdigits->GetEntriesFast() << endl ;
398   cout << endl ;
399   if(strstr(option,"all")||strstr(option,"EMC")){
400     
401     //loop over digits
402     AliPHOSDigit * digit;
403     cout << "EMC sdigits " << endl ;
404     cout << "Digit Id    Amplitude     Index     Nprim  Primaries list " <<  endl;      
405     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
406     Int_t index ;
407     for (index = 0 ; (index < sdigits->GetEntriesFast()) && 
408          (((AliPHOSDigit * )  sdigits->At(index))->GetId() <= maxEmc) ; index++) {
409       digit = (AliPHOSDigit * )  sdigits->At(index) ;
410       if(digit->GetNprimary() == 0) continue;
411       cout << setw(6)  <<  digit->GetId() << "   "  <<  setw(10)  <<  digit->GetAmp() <<   "    "  
412            << setw(6)  <<  digit->GetIndexInList() << "    "   
413            << setw(5)  <<  digit->GetNprimary() <<"    ";
414       
415       Int_t iprimary;
416       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
417         cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << "    ";
418       cout << endl;      
419     }    
420     cout << endl;
421   }
422
423   if(strstr(option,"all")||strstr(option,"CPV")){
424     
425     //loop over CPV digits
426     AliPHOSDigit * digit;
427     cout << "CPV sdigits " << endl ;
428     cout << "Digit Id    Amplitude     Index     Nprim  Primaries list " <<  endl;      
429     Int_t maxEmc = gime->PHOSGeometry()->GetNModules()*gime->PHOSGeometry()->GetNCristalsInModule() ;
430     Int_t index ;
431     for (index = 0 ; index < sdigits->GetEntriesFast(); index++) {
432       digit = (AliPHOSDigit * )  sdigits->At(index) ;
433       if(digit->GetId() > maxEmc){
434         cout << setw(6)  <<  digit->GetId() << "   "  <<        setw(10)  <<  digit->GetAmp() <<   "    "  
435              << setw(6)  <<  digit->GetIndexInList() << "    "   
436              << setw(5)  <<  digit->GetNprimary() <<"    ";
437         
438         Int_t iprimary;
439         for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
440           cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << "    ";
441         cout << endl;    
442       }    
443     }
444   }
445
446 }
447
448 //____________________________________________________________________________ 
449 void AliPHOSSDigitizer::UseHitsFrom(const char * filename)
450 {
451   SetTitle(filename) ; 
452   Init() ; 
453 }