classes for the PID construction
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALPIDv1.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 // Implementation version v1 of the EMCAL particle identifier 
20
21 //*-- Author: Yves Schutz (SUBATECH) 
22 //
23 // --- ROOT system ---
24 #include "TROOT.h"
25 #include "TTree.h"
26 #include "TFile.h"
27 #include "TF2.h"
28 #include "TCanvas.h"
29 #include "TFolder.h"
30 #include "TSystem.h"
31 #include "TBenchmark.h"
32
33 #include "TSystem.h"
34
35 // --- Standard library ---
36
37 #include <Riostream.h>
38
39 // --- AliRoot header files ---
40
41 #include "AliRun.h"
42 #include "AliGenerator.h"
43 #include "AliEMCAL.h"
44 #include "AliEMCALPIDv1.h"
45 #include "AliEMCALClusterizerv1.h"
46 #include "AliEMCALTrackSegment.h"
47 #include "AliEMCALTrackSegmentMakerv1.h"
48 #include "AliEMCALRecParticle.h"
49 #include "AliEMCALGeometry.h"
50 #include "AliEMCALGetter.h"
51
52 ClassImp( AliEMCALPIDv1) 
53
54 //____________________________________________________________________________
55 AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID()
56
57   // default ctor
58  
59   InitParameters() ; 
60   fDefaultInit = kTRUE ; 
61
62 }
63
64 //____________________________________________________________________________
65 AliEMCALPIDv1::AliEMCALPIDv1(const char * headerFile,const char * name, const Bool_t toSplit)
66 :AliEMCALPID(headerFile, name,toSplit)
67
68   //ctor with the indication on where to look for the track segments
69  
70   InitParameters() ; 
71
72   Init() ;
73   fDefaultInit = kFALSE ; 
74
75 }
76
77 //____________________________________________________________________________
78 AliEMCALPIDv1::~AliEMCALPIDv1()
79
80   // dtor
81
82   if (!fDefaultInit) {  
83     fSplitFile = 0 ; 
84   }
85 }
86
87 //____________________________________________________________________________
88 const TString AliEMCALPIDv1::BranchName() const 
89 {  
90   TString branchName(GetName() ) ;
91   branchName.Remove(branchName.Index(Version())-1) ;
92   return branchName ;
93 }
94  
95 //____________________________________________________________________________
96 void AliEMCALPIDv1::Init()
97 {
98   // Make all memory allocations that are not possible in default constructor
99   // Add the PID task to the list of EMCAL tasks
100
101   if ( strcmp(GetTitle(), "") == 0 )
102     SetTitle("galice.root") ;
103
104   TString branchname(GetName()) ;
105   branchname.Remove(branchname.Index(Version())-1) ;    
106   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(),branchname.Data(),fToSplit ) ; 
107
108   //  gime->SetRecParticlesTitle(BranchName()) ;
109   if ( gime == 0 ) {
110     Error("Init", "Could not obtain the Getter object !" ) ;  
111     return ;
112   } 
113
114   fSplitFile = 0 ;
115   if(fToSplit){
116     //First - extract full path if necessary
117     TString fileName(GetTitle()) ;
118     Ssiz_t islash = fileName.Last('/') ;
119     if(islash<fileName.Length())
120       fileName.Remove(islash+1,fileName.Length()) ;
121     else
122       fileName="" ;
123     fileName+="EMCAL.RecData." ;
124     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
125       fileName+=branchname.Data() ;
126       fileName+="." ;
127     }
128     fileName+="root" ;
129     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));   
130     if(!fSplitFile)
131       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
132   }
133   
134   gime->PostPID(this) ;
135   gime->PostRecParticles(branchname) ; 
136   
137 }
138
139 //____________________________________________________________________________
140 void AliEMCALPIDv1::InitParameters()
141 {
142
143   fRecParticlesInRun = 0 ; 
144   fNEvent            = 0 ;            
145   fRecParticlesInRun = 0 ;
146   TString pidName( GetName()) ;
147   if (pidName.IsNull() ) 
148     pidName = "Default" ; 
149   pidName.Append(":") ; 
150   pidName.Append(Version()) ; 
151   SetName(pidName) ;
152   fPi0Analysis = kFALSE ;
153 }
154
155 //____________________________________________________________________________
156
157 void  AliEMCALPIDv1::Exec(Option_t * option) 
158 {
159   //Steering method
160   
161   if( strcmp(GetName(), "")== 0 ) 
162     Init() ;
163   
164   if(strstr(option,"tim"))
165     gBenchmark->Start("EMCALPID");
166   
167   if(strstr(option,"print")) {
168     Print("") ; 
169     return ; 
170   }
171   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
172   if(gime->BranchExists("RecParticles") )
173     return ;
174   Int_t nevents = gime->MaxEvent() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
175   Int_t ievent ;
176
177
178   for(ievent = 0; ievent < nevents; ievent++){
179     gime->Event(ievent,"R") ;
180  
181     MakeRecParticles() ;
182     
183     WriteRecParticles(ievent);
184     
185     if(strstr(option,"deb"))
186       PrintRecParticles(option) ;
187
188     //increment the total number of rec particles per run 
189     fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ; 
190
191   }
192   
193   if(strstr(option,"tim")){
194     gBenchmark->Stop("EMCALPID");
195     Info("Exec", "took %f seconds for PID %f seconds per event", 
196          gBenchmark->GetCpuTime("EMCALPID"),  
197          gBenchmark->GetCpuTime("EMCALPID")/nevents) ;
198   } 
199 }
200
201 //____________________________________________________________________________
202 void  AliEMCALPIDv1::MakeRecParticles(){
203
204   // Makes a RecParticle out of a TrackSegment
205   
206   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
207   TObjArray * aECRecPoints = gime->ECALRecPoints() ; 
208   TObjArray * aPRRecPoints = gime->PRERecPoints() ; 
209   TObjArray * aHCRecPoints = gime->HCALRecPoints() ; 
210   TClonesArray * trackSegments = gime->TrackSegments() ; 
211   if ( !aECRecPoints || !aPRRecPoints || !aHCRecPoints || !trackSegments ) {
212     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;  
213   }
214   TClonesArray * recParticles  = gime->RecParticles() ; 
215   recParticles->Clear();
216
217   TIter next(trackSegments) ; 
218   AliEMCALTrackSegment * ts ; 
219   Int_t index = 0 ; 
220   AliEMCALRecParticle * rp ; 
221   while ( (ts = (AliEMCALTrackSegment *)next()) ) {
222     
223     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
224     rp = (AliEMCALRecParticle *)recParticles->At(index) ; 
225     rp->SetTrackSegment(index) ;
226     rp->SetIndexInList(index) ;
227         
228     AliEMCALTowerRecPoint * ecal = 0 ;
229     if(ts->GetECIndex()>=0)
230       ecal = dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(ts->GetECIndex())) ;
231     
232     AliEMCALTowerRecPoint    * pre = 0 ;
233     if(ts->GetPREIndex()>=0)
234       pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRRecPoints->At(ts->GetPREIndex())) ;
235     
236     AliEMCALTowerRecPoint    * hcal = 0 ;
237     if(ts->GetHCIndex()>=0)
238       hcal = dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(ts->GetHCIndex())) ;
239
240     // Now set type (reconstructed) of the particle
241
242     // Choose the cluster energy range
243     
244     if (!ecal) {
245       Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECIndex(), ecal ) ;
246     }
247
248     Float_t    e = ecal->GetEnergy() ;   
249     
250     Float_t  lambda[2] ;
251     ecal->GetElipsAxis(lambda) ;
252     
253     if((lambda[0]>0.01) && (lambda[1]>0.01)){
254       // Looking PCA. Define and calculate the data (X),
255       // introduce in the function X2P that gives the components (P).  
256
257       Float_t  Spher = 0. ;
258       Float_t  Emaxdtotal = 0. ; 
259       
260       if((lambda[0]+lambda[1])!=0) 
261         Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); 
262       
263       Emaxdtotal=ecal->GetMaximalEnergy()/ecal->GetEnergy(); 
264     }
265     
266     //    Float_t time = ecal->GetTime() ;
267       
268     //Set momentum, energy and other parameters 
269     Float_t  encal = e;
270     TVector3 dir(0., 0., 0.) ; 
271     dir.SetMag(encal) ;
272     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
273     rp->SetCalcMass(0);
274     rp->Name(); //If photon sets the particle pdg name to gamma
275     rp->SetProductionVertex(0,0,0,0);
276     rp->SetFirstMother(-1);
277     rp->SetLastMother(-1);
278     rp->SetFirstDaughter(-1);
279     rp->SetLastDaughter(-1);
280     rp->SetPolarisation(0,0,0);
281     index++ ; 
282   }
283   
284 }
285
286 //____________________________________________________________________________
287 void  AliEMCALPIDv1:: Print()
288 {
289   // Print the parameters used for the particle type identification
290
291   TString message ; 
292     message  = "\n=============== AliEMCALPID1 ================\n" ;
293     message += "Making PID\n";
294     message += "    Pricipal analysis file from 0.5 to 100 %s\n" ; 
295     message += "    Name of parameters file     %s\n" ;
296     message += "    Matrix of Parameters: 9x4\n" ;
297     message += "        RCPV 2x3 rows x and z, columns function cut parameters\n" ;
298     message += "        TOF  1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n" ;
299     message += "        PCA  5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n" ;
300     message += "        Energy Calibration  1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n" ;
301     Info("Print", message.Data(), fFileName.Data(), fFileNamePar.Data() ) ; 
302 }
303
304 //____________________________________________________________________________
305 void  AliEMCALPIDv1::WriteRecParticles(Int_t event)
306 {
307  
308   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; 
309
310   TClonesArray * recParticles = gime->RecParticles() ; 
311   recParticles->Expand(recParticles->GetEntriesFast() ) ;
312   TTree * treeR ;
313
314   if(fToSplit){
315     if(!fSplitFile)
316       return ;
317     fSplitFile->cd() ;
318     char name[10] ;
319     sprintf(name,"%s%d", "TreeR",event) ;
320     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); 
321   }
322   else{
323     treeR = gAlice->TreeR();
324   }
325   
326   if(!treeR){
327     gAlice->MakeTree("R", fSplitFile);
328     treeR = gAlice->TreeR() ;
329   }
330   
331   //First rp
332   Int_t bufferSize = 32000 ;    
333   TBranch * rpBranch = treeR->Branch("EMCALRP",&recParticles,bufferSize);
334   rpBranch->SetTitle(BranchName());
335
336   
337   //second, pid
338   Int_t splitlevel = 0 ; 
339   AliEMCALPIDv1 * pid = this ;
340   TBranch * pidBranch = treeR->Branch("AliEMCALPID","AliEMCALPIDv1",&pid,bufferSize,splitlevel);
341   pidBranch->SetTitle(BranchName());
342   
343   rpBranch->Fill() ;
344   pidBranch->Fill() ; 
345   
346   treeR->AutoSave() ; //Write(0,kOverwrite) ;  
347   if(gAlice->TreeR()!=treeR){
348     treeR->Delete();
349   }
350 }
351
352 //____________________________________________________________________________
353 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
354 {
355   // Print table of reconstructed particles
356
357   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; 
358
359   TClonesArray * recParticles = gime->RecParticles(BranchName()) ; 
360
361   TString message ; 
362   message  = "\nevent " ;
363   message += gAlice->GetEvNumber() ; 
364   message += "       found " ; 
365   message += recParticles->GetEntriesFast(); 
366   message += " RecParticles\n" ; 
367
368   if(strstr(option,"all")) {  // printing found TS
369     message += "\n  PARTICLE         Index    \n" ; 
370     
371     Int_t index ;
372     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
373       AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ;       
374       message += "\n" ;
375       message += rp->Name().Data() ;  
376       message += " " ;
377       message += rp->GetIndexInList() ;  
378       message += " " ;
379       message += rp->GetType()  ;
380     }
381   }
382   Info("Print", message.Data() ) ; 
383 }
384
385
386