Initialize arrays in constructor (Ivana)
[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 "TBenchmark.h"
27 #include "TSystem.h"
28   
29   // --- Standard library ---
30   
31   // --- AliRoot header files ---
32 #include "AliGenerator.h"
33 #include "AliEMCALPIDv1.h"
34 #include "AliEMCALRecParticle.h"
35 #include "AliEMCALGetter.h"
36   
37   ClassImp( AliEMCALPIDv1) 
38
39 //____________________________________________________________________________
40 AliEMCALPIDv1::AliEMCALPIDv1():AliEMCALPID()
41
42   // default ctor
43   
44   InitParameters() ; 
45   fDefaultInit = kTRUE ; 
46   
47 }
48
49 //____________________________________________________________________________
50 AliEMCALPIDv1::AliEMCALPIDv1(const AliEMCALPIDv1 & pid ):AliEMCALPID(pid)
51
52   // ctor
53   InitParameters() ; 
54   Init() ;
55   
56 }
57
58 //____________________________________________________________________________
59 AliEMCALPIDv1::AliEMCALPIDv1(const TString alirunFileName, const TString eventFolderName):AliEMCALPID(alirunFileName, eventFolderName)
60
61   //ctor with the indication on where to look for the track segments
62  
63   InitParameters() ; 
64   Init() ;
65   fDefaultInit = kFALSE ; 
66
67 }
68
69 //____________________________________________________________________________
70 AliEMCALPIDv1::~AliEMCALPIDv1()
71
72   // dtor
73 }
74
75 //____________________________________________________________________________
76 const TString AliEMCALPIDv1::BranchName() const 
77 {  
78
79   return GetName() ;
80 }
81  
82 //____________________________________________________________________________
83 Float_t  AliEMCALPIDv1::GetCalibratedEnergy(Float_t e) const
84 {
85 //      It calibrates Energy depending on the recpoint energy.
86 //      The energy of the reconstructed cluster is corrected with 
87 //      the formula A + B* E  + C* E^2, whose parameters where obtained 
88 //      through the study of the reconstructed energy distribution of 
89 //      monoenergetic photons.
90  
91   //Float_t p[]={0.,0.,0.};
92   //for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i);
93   Float_t enerec = e ; // p[0] +  p[1]*e + p[2]*e*e;
94   return enerec ;
95
96 }
97 //____________________________________________________________________________
98 TVector3 AliEMCALPIDv1::GetMomentumDirection(AliEMCALRecPoint * emc)const 
99
100   // Calculates the momentum direction:
101   // direction is given by IP and this RecPoint
102   
103
104   TVector3 dir(0,0,0) ; 
105   TVector3 emcglobalpos ;
106   // TMatrix  dummy ;
107   
108   emc->GetGlobalPosition(emcglobalpos) ;
109   
110
111   dir = emcglobalpos ;  
112   // dir.SetMag(1.) ; Removed to avoid warings !!!!!!!!!!!!!! TO BE REVISED
113
114   //account correction to the position of IP
115   Float_t xo,yo,zo ; //Coordinates of the origin
116   gAlice->Generator()->GetOrigin(xo,yo,zo) ;
117   TVector3 origin(xo,yo,zo);
118   dir = dir - origin ;
119
120   return dir ;  
121 }
122
123 //____________________________________________________________________________
124 void AliEMCALPIDv1::Init()
125 {
126   // Make all memory allocations that are not possible in default constructor
127   // Add the PID task to the list of EMCAL tasks
128
129
130   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data()) ; 
131
132   if ( !gime->PID() ) 
133     gime->PostPID(this) ;
134 }
135
136 //____________________________________________________________________________
137 void AliEMCALPIDv1::InitParameters()
138 {
139   // Initialize the parameters
140   fRecParticlesInRun = 0 ; 
141   fNEvent            = 0 ;            
142   fRecParticlesInRun = 0 ;
143 }
144
145 //____________________________________________________________________________
146
147 void  AliEMCALPIDv1::Exec(Option_t * option) 
148 {
149   //Steering method
150
151   if(strstr(option,"tim"))
152     gBenchmark->Start("EMCALPID");
153   
154   if(strstr(option,"print")) {
155     Print("") ; 
156     return ; 
157   }
158   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
159   
160   if (fLastEvent == -1) 
161     fLastEvent = gime->MaxEvent() - 1 ;
162   else 
163     fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent());
164   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
165
166   Int_t ievent ;
167
168   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
169     gime->Event(ievent,"R") ;
170     MakeRecParticles() ;
171     WriteRecParticles();
172     if(strstr(option,"deb"))
173       PrintRecParticles(option) ;
174     //increment the total number of rec particles per run 
175     fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ;  
176   }
177   if(strstr(option,"tim")){
178     gBenchmark->Stop("EMCALPID");
179     printf("Exec: took %f seconds for PID %f seconds per event", 
180          gBenchmark->GetCpuTime("EMCALPID"),  
181          gBenchmark->GetCpuTime("EMCALPID")/nEvents) ;
182   } 
183
184   Unload();
185 }
186
187 //____________________________________________________________________________
188 void  AliEMCALPIDv1::MakeRecParticles(){
189
190   // Makes a RecParticle out of a TrackSegment
191   
192   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
193   TObjArray * aECARecPoints = gime->ECARecPoints() ; 
194   if ( !aECARecPoints ) {
195     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;  
196   }
197   TClonesArray * recParticles  = gime->RecParticles() ; 
198   recParticles->Clear();
199
200   TIter next(aECARecPoints) ; 
201   AliEMCALRecPoint * eca ; 
202   Int_t index = 0 ; 
203   AliEMCALRecParticle * rp ; 
204   while ( (eca = (AliEMCALRecPoint *)next()) ) {
205     
206     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
207     rp = (AliEMCALRecParticle *)recParticles->At(index) ; 
208     rp->SetRecPoint(index) ;
209     rp->SetIndexInList(index) ;
210         
211     // Now set type (reconstructed) of the particle
212
213     // Choose the cluster energy range
214     
215     Float_t  lambda[2] ;
216     eca->GetElipsAxis(lambda) ;
217     
218     if((lambda[0]>0.01) && (lambda[1]>0.01)){
219       // Looking PCA. Define and calculate the data (X),
220       // introduce in the function X2P that gives the components (P).  
221
222       Float_t  spher = 0. ;
223       Float_t  emaxdtotal = 0. ; 
224       
225       if((lambda[0]+lambda[1])!=0) 
226         spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); 
227       
228       emaxdtotal=eca->GetMaximalEnergy()/eca->GetEnergy(); 
229     }
230     
231     //    Float_t time = eca->GetTime() ;
232       
233     //Set momentum, energy and other parameters 
234     Float_t  encal = GetCalibratedEnergy(eca->GetEnergy());
235     TVector3 dir   = GetMomentumDirection(eca) ; 
236     // dir.SetMag(encal) ;Removed to avoid warings !!!!!!!!!!!!!! TO BE REVISED
237     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
238     rp->SetCalcMass(0);
239     rp->Name(); //If photon sets the particle pdg name to gamma
240     rp->SetProductionVertex(0,0,0,0);
241     rp->SetFirstMother(-1);
242     rp->SetLastMother(-1);
243     rp->SetFirstDaughter(-1);
244     rp->SetLastDaughter(-1);
245     rp->SetPolarisation(0,0,0);
246     //Set the position in global coordinate system from the RecPoint
247     //AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
248     //AliEMCALTowerRecPoint  * erp = gime->ECARecPoint(rp->GetEMCALRPIndex()) ; 
249     TVector3 pos ; 
250     //geom->GetGlobal(erp, pos) ; !!!!!!!!!! to check 
251     rp->SetPos(pos);
252
253
254     index++ ; 
255   }
256   
257 }
258
259 //____________________________________________________________________________
260 void  AliEMCALPIDv1:: Print(Option_t * /*option*/) const
261 {
262   // Print the parameters used for the particle type identification
263
264     printf("Print: =============== AliEMCALPID1 ================") ;
265     printf("Making PID\n");
266     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileName.Data() ) ; 
267     printf("    Name of parameters file     %s\n", fFileNamePar.Data() )  ;
268 }
269
270 //____________________________________________________________________________
271 void  AliEMCALPIDv1::Print() const
272 {
273   // Print the parameters used for the particle type identification
274
275     Info("Print", "=============== AliEMCALPIDv1 ================") ;
276 }
277
278 //____________________________________________________________________________
279 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
280 {
281   // Print table of reconstructed particles
282
283   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
284
285   TClonesArray * recParticles = gime->RecParticles() ; 
286
287   printf("\nevent %i", gAlice->GetEvNumber()); 
288   printf("       found %i", recParticles->GetEntriesFast()); 
289   printf(" RecParticles\n"); 
290
291   if(strstr(option,"all")) {  // printing found TS
292     printf("\n  PARTICLE         Index    \n"); 
293     
294     Int_t index ;
295     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
296       AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ;       
297       printf("\n");
298       printf(rp->Name().Data());  
299       printf(" %i", rp->GetIndexInList());  
300       printf(" %i", rp->GetType());
301     }
302   }
303 }
304
305 //____________________________________________________________________________
306 void AliEMCALPIDv1::Unload() 
307 {
308   // Unloads RecPoints, TrackSegments and RecParticles from the folder 
309   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;  
310   gime->EmcalLoader()->UnloadRecPoints() ;
311   gime->EmcalLoader()->UnloadTracks() ;
312   gime->EmcalLoader()->UnloadRecParticles() ;
313 }
314
315 //____________________________________________________________________________
316 void  AliEMCALPIDv1::WriteRecParticles()
317 {
318   // Write RecParticles array to file
319   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
320
321   TClonesArray * recParticles = gime->RecParticles() ; 
322   recParticles->Expand(recParticles->GetEntriesFast() ) ;
323   TTree * treeP =  gime->TreeP() ;
324
325   
326   
327   //First rp
328   Int_t bufferSize = 32000 ;    
329   TBranch * rpBranch = treeP->Branch("EMCALRP",&recParticles,bufferSize);
330   rpBranch->SetTitle(BranchName());
331
332   rpBranch->Fill() ;
333  
334   gime->WriteRecParticles("OVERWRITE");
335   gime->WritePID("OVERWRITE");
336
337 }
338