Transition to NewIO
[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 #include "TSystem.h"
33   
34   // --- Standard library ---
35   
36   // --- AliRoot header files ---
37   
38 #include "AliGenerator.h"
39 #include "AliEMCAL.h"
40 #include "AliEMCALPIDv1.h"
41 #include "AliEMCALClusterizerv1.h"
42 #include "AliEMCALTrackSegment.h"
43 #include "AliEMCALTrackSegmentMakerv1.h"
44 #include "AliEMCALRecParticle.h"
45 #include "AliEMCALGeometry.h"
46 #include "AliEMCALGetter.h"
47   
48   ClassImp( AliEMCALPIDv1) 
49
50 //____________________________________________________________________________
51 AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID()
52
53   // default ctor
54   
55   InitParameters() ; 
56   fDefaultInit = kTRUE ; 
57   
58 }
59
60 //____________________________________________________________________________
61 AliEMCALPIDv1::AliEMCALPIDv1(const AliEMCALPIDv1 & pid ):AliEMCALPID(pid)
62
63   // ctor
64   InitParameters() ; 
65   Init() ;
66   
67 }
68
69 //____________________________________________________________________________
70 AliEMCALPIDv1::AliEMCALPIDv1(const TString alirunFileName, const TString eventFolderName):AliEMCALPID(alirunFileName, eventFolderName)
71
72   //ctor with the indication on where to look for the track segments
73  
74   InitParameters() ; 
75   Init() ;
76   fDefaultInit = kFALSE ; 
77
78 }
79
80 //____________________________________________________________________________
81 AliEMCALPIDv1::~AliEMCALPIDv1()
82
83   // dtor
84 }
85
86 //____________________________________________________________________________
87 const TString AliEMCALPIDv1::BranchName() const 
88 {  
89
90   return GetName() ;
91 }
92  
93 //____________________________________________________________________________
94 void AliEMCALPIDv1::Init()
95 {
96   // Make all memory allocations that are not possible in default constructor
97   // Add the PID task to the list of EMCAL tasks
98
99
100   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data()) ; 
101
102   if ( !gime->PID() ) 
103     gime->PostPID(this) ;
104 }
105
106 //____________________________________________________________________________
107 void AliEMCALPIDv1::InitParameters()
108 {
109
110   fRecParticlesInRun = 0 ; 
111   fNEvent            = 0 ;            
112   fRecParticlesInRun = 0 ;
113 }
114
115 //____________________________________________________________________________
116
117 void  AliEMCALPIDv1::Exec(Option_t * option) 
118 {
119   //Steering method
120
121   if(strstr(option,"tim"))
122     gBenchmark->Start("EMCALPID");
123   
124   if(strstr(option,"print")) {
125     Print("") ; 
126     return ; 
127   }
128   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
129
130   Int_t nevents = gime->MaxEvent() ;      
131   Int_t ievent ;
132
133
134   for(ievent = 0; ievent < nevents; ievent++){
135     gime->Event(ievent,"TR") ;
136     if(gime->TrackSegments() && //Skip events, where no track segments made
137        gime->TrackSegments()->GetEntriesFast()) {   
138       MakeRecParticles() ;
139       WriteRecParticles(ievent);
140       if(strstr(option,"deb"))
141         PrintRecParticles(option) ;
142       //increment the total number of rec particles per run 
143       fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; 
144     }
145   }
146   if(strstr(option,"tim")){
147     gBenchmark->Stop("EMCALPID");
148     Info("Exec", "took %f seconds for PID %f seconds per event", 
149          gBenchmark->GetCpuTime("EMCALPID"),  
150          gBenchmark->GetCpuTime("EMCALPID")/nevents) ;
151   } 
152
153   Unload();
154 }
155
156 //____________________________________________________________________________
157 void  AliEMCALPIDv1::MakeRecParticles(){
158
159   // Makes a RecParticle out of a TrackSegment
160   
161   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
162   TObjArray * aECARecPoints = gime->ECARecPoints() ; 
163   TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
164   TObjArray * aHCARecPoints = gime->HCARecPoints() ; 
165   TClonesArray * trackSegments = gime->TrackSegments() ; 
166   if ( !aECARecPoints || !aPRERecPoints || !aHCARecPoints || !trackSegments ) {
167     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;  
168   }
169   TClonesArray * recParticles  = gime->RecParticles() ; 
170   recParticles->Clear();
171
172   TIter next(trackSegments) ; 
173   AliEMCALTrackSegment * ts ; 
174   Int_t index = 0 ; 
175   AliEMCALRecParticle * rp ; 
176   while ( (ts = (AliEMCALTrackSegment *)next()) ) {
177     
178     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
179     rp = (AliEMCALRecParticle *)recParticles->At(index) ; 
180     rp->SetTrackSegment(index) ;
181     rp->SetIndexInList(index) ;
182         
183     AliEMCALTowerRecPoint * eca = 0 ;
184     if(ts->GetECAIndex()>=0)
185       eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(ts->GetECAIndex())) ;
186     
187     AliEMCALTowerRecPoint * pre = 0 ;
188     if(ts->GetPREIndex()>=0)
189       pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(ts->GetPREIndex())) ;
190     
191     AliEMCALTowerRecPoint * hca = 0 ;
192     if(ts->GetHCAIndex()>=0)
193       hca = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(ts->GetHCAIndex())) ;
194
195     // Now set type (reconstructed) of the particle
196
197     // Choose the cluster energy range
198     
199     if (!eca) {
200       Fatal("MakeRecParticles", "-> emcal(%d) = %d", ts->GetECAIndex(), eca ) ;
201     }
202
203     Float_t    e = eca->GetEnergy() ;   
204     
205     Float_t  lambda[2] ;
206     eca->GetElipsAxis(lambda) ;
207     
208     if((lambda[0]>0.01) && (lambda[1]>0.01)){
209       // Looking PCA. Define and calculate the data (X),
210       // introduce in the function X2P that gives the components (P).  
211
212       Float_t  Spher = 0. ;
213       Float_t  Emaxdtotal = 0. ; 
214       
215       if((lambda[0]+lambda[1])!=0) 
216         Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); 
217       
218       Emaxdtotal=eca->GetMaximalEnergy()/eca->GetEnergy(); 
219     }
220     
221     //    Float_t time = ecal->GetTime() ;
222       
223     //Set momentum, energy and other parameters 
224     Float_t  enca = e;
225     TVector3 dir(0., 0., 0.) ; 
226     dir.SetMag(enca) ;
227     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),enca) ;
228     rp->SetCalcMass(0);
229     rp->Name(); //If photon sets the particle pdg name to gamma
230     rp->SetProductionVertex(0,0,0,0);
231     rp->SetFirstMother(-1);
232     rp->SetLastMother(-1);
233     rp->SetFirstDaughter(-1);
234     rp->SetLastDaughter(-1);
235     rp->SetPolarisation(0,0,0);
236     index++ ; 
237   }
238   
239 }
240
241 //____________________________________________________________________________
242 void  AliEMCALPIDv1:: Print()
243 {
244   // Print the parameters used for the particle type identification
245
246     Info("Print", "=============== AliEMCALPID1 ================") ;
247     printf("Making PID\n");
248     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileName.Data() ) ; 
249     printf("    Name of parameters file     %s\n", fFileNamePar.Data() )  ;
250 }
251
252 //____________________________________________________________________________
253 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
254 {
255   // Print table of reconstructed particles
256
257   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
258
259   TClonesArray * recParticles = gime->RecParticles() ; 
260
261   TString message ; 
262   message  = "\nevent " ;
263   message += gAlice->GetEvNumber() ; 
264   message += "       found " ; 
265   message += recParticles->GetEntriesFast(); 
266   message += " RecParticles\n" ; 
267
268   if(strstr(option,"all")) {  // printing found TS
269     message += "\n  PARTICLE         Index    \n" ; 
270     
271     Int_t index ;
272     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
273       AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ;       
274       message += "\n" ;
275       message += rp->Name().Data() ;  
276       message += " " ;
277       message += rp->GetIndexInList() ;  
278       message += " " ;
279       message += rp->GetType()  ;
280     }
281   }
282   Info("Print", message.Data() ) ; 
283 }
284
285 //____________________________________________________________________________
286 void AliEMCALPIDv1::Unload() 
287 {
288   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;  
289   gime->EmcalLoader()->UnloadRecPoints() ;
290   gime->EmcalLoader()->UnloadTracks() ;
291   gime->EmcalLoader()->UnloadRecParticles() ;
292 }
293
294 //____________________________________________________________________________
295 void  AliEMCALPIDv1::WriteRecParticles(Int_t event)
296 {
297  
298   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
299
300   TClonesArray * recParticles = gime->RecParticles() ; 
301   recParticles->Expand(recParticles->GetEntriesFast() ) ;
302   TTree * treeP =  gime->TreeP() ;
303
304   
305   
306   //First rp
307   Int_t bufferSize = 32000 ;    
308   TBranch * rpBranch = treeP->Branch("EMCALRP",&recParticles,bufferSize);
309   rpBranch->SetTitle(BranchName());
310
311   rpBranch->Fill() ;
312  
313   gime->WriteRecParticles("OVERWRITE");
314   gime->WritePID("OVERWRITE");
315
316 }
317