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