Enlarging array
[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(GetTitle()) ; 
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
169   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
170     gime->Event(ievent,"R") ;
171     MakeRecParticles() ;
172     WriteRecParticles();
173     if(strstr(option,"deb"))
174       PrintRecParticles(option) ;
175     //increment the total number of rec particles per run 
176     fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; 
177     
178   }
179   if(strstr(option,"tim")){
180     gBenchmark->Stop("EMCALPID");
181     printf("Exec: took %f seconds for PID %f seconds per event", 
182          gBenchmark->GetCpuTime("EMCALPID"),  
183          gBenchmark->GetCpuTime("EMCALPID")/nEvents) ;
184   } 
185
186   Unload();
187 }
188
189 //____________________________________________________________________________
190 void  AliEMCALPIDv1::MakeRecParticles(){
191
192   // Makes a RecParticle out of a TrackSegment
193   
194   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
195   TObjArray * aECARecPoints = gime->ECARecPoints() ; 
196   if ( !aECARecPoints ) {
197     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;  
198   }
199   TClonesArray * recParticles  = gime->RecParticles() ; 
200   recParticles->Clear();
201
202   TIter next(aECARecPoints) ; 
203   AliEMCALRecPoint * eca ; 
204   Int_t index = 0 ; 
205   AliEMCALRecParticle * rp ; 
206   while ( (eca = (AliEMCALRecPoint *)next()) ) {
207     
208     new( (*recParticles)[index] ) AliEMCALRecParticle() ;
209     rp = (AliEMCALRecParticle *)recParticles->At(index) ; 
210     rp->SetRecPoint(index) ;
211     rp->SetIndexInList(index) ;
212         
213     // Now set type (reconstructed) of the particle
214
215     // Choose the cluster energy range
216     
217     Float_t  lambda[2] ;
218     eca->GetElipsAxis(lambda) ;
219     
220     if((lambda[0]>0.01) && (lambda[1]>0.01)){
221       // Looking PCA. Define and calculate the data (X),
222       // introduce in the function X2P that gives the components (P).  
223
224       Float_t  spher = 0. ;
225       Float_t  emaxdtotal = 0. ; 
226       
227       if((lambda[0]+lambda[1])!=0) 
228         spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); 
229       
230       emaxdtotal=eca->GetMaximalEnergy()/eca->GetEnergy(); 
231     }
232     
233     //    Float_t time = eca->GetTime() ;
234       
235     //Set momentum, energy and other parameters 
236     Float_t  encal = GetCalibratedEnergy(eca->GetEnergy());
237     TVector3 dir   = GetMomentumDirection(eca) ; 
238     // dir.SetMag(encal) ;Removed to avoid warings !!!!!!!!!!!!!! TO BE REVISED
239     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
240     rp->SetCalcMass(0);
241     rp->Name(); //If photon sets the particle pdg name to gamma
242     rp->SetProductionVertex(0,0,0,0);
243     rp->SetFirstMother(-1);
244     rp->SetLastMother(-1);
245     rp->SetFirstDaughter(-1);
246     rp->SetLastDaughter(-1);
247     rp->SetPolarisation(0,0,0);
248     //Set the position in global coordinate system from the RecPoint
249     //AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
250     //AliEMCALTowerRecPoint  * erp = gime->ECARecPoint(rp->GetEMCALRPIndex()) ; 
251     TVector3 pos ; 
252     //geom->GetGlobal(erp, pos) ; !!!!!!!!!! to check 
253     rp->SetPos(pos);
254
255
256     index++ ; 
257   }
258   
259 }
260
261 //____________________________________________________________________________
262 void  AliEMCALPIDv1:: Print(Option_t * /*option*/) const
263 {
264   // Print the parameters used for the particle type identification
265
266     printf("Print: =============== AliEMCALPID1 ================") ;
267     printf("Making PID\n");
268     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileName.Data() ) ; 
269     printf("    Name of parameters file     %s\n", fFileNamePar.Data() )  ;
270 }
271
272 //____________________________________________________________________________
273 void  AliEMCALPIDv1::Print() const
274 {
275   // Print the parameters used for the particle type identification
276
277     Info("Print", "=============== AliEMCALPIDv1 ================") ;
278 }
279
280 //____________________________________________________________________________
281 void AliEMCALPIDv1::PrintRecParticles(Option_t * option)
282 {
283   // Print table of reconstructed particles
284
285   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
286
287   TClonesArray * recParticles = gime->RecParticles() ; 
288
289   printf("\nevent %i", gAlice->GetEvNumber()); 
290   printf("       found %i", recParticles->GetEntriesFast()); 
291   printf(" RecParticles\n"); 
292
293   if(strstr(option,"all")) {  // printing found TS
294     printf("\n  PARTICLE         Index    \n"); 
295     
296     Int_t index ;
297     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
298       AliEMCALRecParticle * rp = (AliEMCALRecParticle * ) recParticles->At(index) ;       
299       printf("\n");
300       printf(rp->Name().Data());  
301       printf(" %i", rp->GetIndexInList());  
302       printf(" %i", rp->GetType());
303     }
304   }
305 }
306
307 //____________________________________________________________________________
308 void AliEMCALPIDv1::Unload() 
309 {
310   // Unloads RecPoints, TrackSegments and RecParticles from the folder 
311   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;  
312   gime->EmcalLoader()->UnloadRecPoints() ;
313   gime->EmcalLoader()->UnloadTracks() ;
314   gime->EmcalLoader()->UnloadRecParticles() ;
315 }
316
317 //____________________________________________________________________________
318 void  AliEMCALPIDv1::WriteRecParticles()
319 {
320   // Write RecParticles array to file
321   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
322
323   TClonesArray * recParticles = gime->RecParticles() ; 
324   recParticles->Expand(recParticles->GetEntriesFast() ) ;
325   TTree * treeP =  gime->TreeP() ;
326
327   
328   
329   //First rp
330   Int_t bufferSize = 32000 ;    
331   TBranch * rpBranch = treeP->Branch("EMCALRP",&recParticles,bufferSize);
332   rpBranch->SetTitle(BranchName());
333
334   rpBranch->Fill() ;
335  
336   gime->WriteRecParticles("OVERWRITE");
337   gime->WritePID("OVERWRITE");
338
339 }
340