First version of standalone Makefile for test of PDC04
[u/mrichter/AliRoot.git] / ANALYSIS / AliReaderESD.cxx
CommitLineData
a5556ea5 1#include "AliReaderESD.h"
2//____________________________________________________________________
3//////////////////////////////////////////////////////////////////////
4// //
5// class AliReaderESD //
6// //
7// reader for ALICE Event Summary Data (ESD). //
8// //
9// Piotr.Skowronski@cern.ch //
10// //
11//////////////////////////////////////////////////////////////////////
12
13#include <TPDGCode.h>
14#include <TString.h>
15#include <TObjString.h>
16#include <TTree.h>
17#include <TFile.h>
18#include <TKey.h>
19#include <TParticle.h>
20#include <TH1.h>
21
22#include <AliRun.h>
23#include <AliRunLoader.h>
24#include <AliStack.h>
25#include <AliESDtrack.h>
7814b246 26#include <AliKalmanTrack.h>
a5556ea5 27#include <AliESD.h>
28
29#include "AliAnalysis.h"
30#include "AliAODRun.h"
31#include "AliAOD.h"
afa8b37b 32#include "AliAODParticle.h"
a5556ea5 33#include "AliAODParticleCut.h"
34#include "AliTrackPoints.h"
35#include "AliClusterMap.h"
36
37ClassImp(AliReaderESD)
38
39AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
40 fESDFileName(esdfilename),
41 fGAlFileName(galfilename),
42 fFile(0x0),
43 fRunLoader(0x0),
44 fKeyIterator(0x0),
45 fReadSim(kFALSE),
46 fCheckParticlePID(kFALSE),
47 fReadMostProbableOnly(kFALSE),
48 fNTrackPoints(0),
49 fdR(0.0),
50 fClusterMap(kFALSE),
7814b246 51 fITSTrackPoints(kFALSE),
52 fMustTPC(kFALSE),
beb1c41d 53 fReadCentralBarrel(kFALSE),
54 fReadMuon(kFALSE),
55 fReadPHOS(kFALSE),
a5556ea5 56 fNTPCClustMin(0),
57 fNTPCClustMax(150),
58 fTPCChi2PerClustMin(0.0),
59 fTPCChi2PerClustMax(10e5),
60 fChi2Min(0.0),
61 fChi2Max(10e5),
62 fC00Min(0.0),
63 fC00Max(10e5),
64 fC11Min(0.0),
65 fC11Max(10e5),
66 fC22Min(0.0),
67 fC22Max(10e5),
68 fC33Min(0.0),
69 fC33Max(10e5),
70 fC44Min(0.0),
71 fC44Max(10e5),
72 fTPCC00Min(0.0),
73 fTPCC00Max(10e5),
74 fTPCC11Min(0.0),
75 fTPCC11Max(10e5),
76 fTPCC22Min(0.0),
77 fTPCC22Max(10e5),
78 fTPCC33Min(0.0),
79 fTPCC33Max(10e5),
80 fTPCC44Min(0.0),
81 fTPCC44Max(10e5)
82
83{
84 //cosntructor
85 if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
86 Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
87}
88/********************************************************************/
89
90AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
91 AliReader(dirs),
92 fESDFileName(esdfilename),
93 fGAlFileName(galfilename),
94 fFile(0x0),
95 fRunLoader(0x0),
96 fKeyIterator(0x0),
97 fReadSim(kFALSE),
98 fCheckParticlePID(kFALSE),
99 fReadMostProbableOnly(kFALSE),
100 fNTrackPoints(0),
101 fdR(0.0),
102 fClusterMap(kFALSE),
7814b246 103 fITSTrackPoints(kFALSE),
104 fMustTPC(kFALSE),
beb1c41d 105 fReadCentralBarrel(kFALSE),
106 fReadMuon(kFALSE),
107 fReadPHOS(kFALSE),
a5556ea5 108 fNTPCClustMin(0),
109 fNTPCClustMax(150),
110 fTPCChi2PerClustMin(0.0),
111 fTPCChi2PerClustMax(10e5),
112 fChi2Min(0.0),
113 fChi2Max(10e5),
114 fC00Min(0.0),
115 fC00Max(10e5),
116 fC11Min(0.0),
117 fC11Max(10e5),
118 fC22Min(0.0),
119 fC22Max(10e5),
120 fC33Min(0.0),
121 fC33Max(10e5),
122 fC44Min(0.0),
123 fC44Max(10e5),
124 fTPCC00Min(0.0),
125 fTPCC00Max(10e5),
126 fTPCC11Min(0.0),
127 fTPCC11Max(10e5),
128 fTPCC22Min(0.0),
129 fTPCC22Max(10e5),
130 fTPCC33Min(0.0),
131 fTPCC33Max(10e5),
132 fTPCC44Min(0.0),
133 fTPCC44Max(10e5)
134{
135 //cosntructor
136 if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
137 Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
138}
139/********************************************************************/
140
141AliReaderESD::~AliReaderESD()
142{
143 //desctructor
beb1c41d 144 delete fRunLoader;
145 delete fKeyIterator;
146 delete fFile;
a5556ea5 147}
beb1c41d 148
a5556ea5 149/**********************************************************/
150Int_t AliReaderESD::ReadNext()
151{
152//reads next event from fFile
153 //fRunLoader is for reading Kine
154
afa8b37b 155 if (AliVAODParticle::GetDebug())
a5556ea5 156 Info("ReadNext","Entered");
157
158 if (fEventSim == 0x0) fEventSim = new AliAOD();
159 if (fEventRec == 0x0) fEventRec = new AliAOD();
160
161 fEventSim->Reset();
162 fEventRec->Reset();
163
164 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
beb1c41d 165 {
166 if (fFile == 0x0)
167 {
168 fFile = OpenFile(fCurrentDir);//rl is opened here
169 if (fFile == 0x0)
170 {
171 Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
172 fCurrentDir++;
173 continue;
174 }
175 fCurrentEvent = 0;
176 }
a5556ea5 177 TString esdname = "ESD";
178 esdname+=fCurrentEvent;
179 AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
180 if (esd == 0x0)
181 {
afa8b37b 182 if (AliVAODParticle::GetDebug() > 2 )
a5556ea5 183 {
184 Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
185 }
186 fCurrentDir++;
a5556ea5 187 delete fFile;//we have to assume there is no more ESD objects in the fFile
188 fFile = 0x0;
189 delete fRunLoader;
190 fRunLoader = 0x0;
191 continue;
192 }
beb1c41d 193 ReadESD(esd);
a5556ea5 194
beb1c41d 195 fCurrentEvent++;
196 fNEventsRead++;
197 delete esd;
198 return 0;//success -> read one event
199 }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
a5556ea5 200
201 return 1; //no more directories to read
202}
a5556ea5 203
beb1c41d 204/**********************************************************/
a5556ea5 205Int_t AliReaderESD::ReadESD(AliESD* esd)
206{
beb1c41d 207//Reads esd data
208 if (esd == 0x0)
209 {
210 Error("ReadESD","ESD is NULL");
211 return 1;
212 }
213
214 // seperate each method
215 if (fReadCentralBarrel) ReadESDCentral(esd);
216
217 if (fReadMuon) ReadESDMuon(esd);
218
219 if (fReadPHOS) ReadESDPHOS(esd);
220
221 return 1;
222}
223
224/**********************************************************/
225Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
226{
a5556ea5 227 //****** Tentative particle type "concentrations"
228 static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
229
230 Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
231 Double_t w[kNSpecies];
232 Double_t mom[3];//momentum
233 Double_t pos[3];//position
234 Double_t vertexpos[3];//vertex position
235 //Reads one ESD
beb1c41d 236
a5556ea5 237 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
238 if (pdgdb == 0x0)
239 {
240 Error("ReadESD","Can not get PDG Database Instance.");
241 return 1;
242 }
243
244 Float_t mf = esd->GetMagneticField();
245
7814b246 246 if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
a5556ea5 247 {
248 Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
249 return 1;
250 }
251
7814b246 252 if (fITSTrackPoints)
253 {
254 Info("ReadESD","Magnetic Field is %f",mf/10.);
255 AliKalmanTrack::SetMagneticField(mf/10.);
256 }
257
a5556ea5 258 AliStack* stack = 0x0;
259 if (fReadSim && fRunLoader)
260 {
261 fRunLoader->GetEvent(fCurrentEvent);
262 stack = fRunLoader->Stack();
263 }
264
265 const AliESDVertex* vertex = esd->GetVertex();
266 if (vertex == 0x0)
267 {
268 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
269 vertexpos[0] = 0.0;
270 vertexpos[1] = 0.0;
271 vertexpos[2] = 0.0;
272 }
273 else
274 {
275 vertex->GetXYZ(vertexpos);
276 }
277
afa8b37b 278 if (AliVAODParticle::GetDebug() > 0)
a5556ea5 279 {
280 Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]);
281 }
282
283 Info("ReadESD","Reading Event %d",fCurrentEvent);
284
285 Int_t ntr = esd->GetNumberOfTracks();
286 Info("ReadESD","Found %d tracks.",ntr);
287 for (Int_t i = 0;i<ntr; i++)
288 {
289 AliESDtrack *esdtrack = esd->GetTrack(i);
290 if (esdtrack == 0x0)
291 {
292 Error("Next","Can not get track %d", i);
293 continue;
294 }
295
296 //if (esdtrack->HasVertexParameters() == kFALSE)
297 if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
298 {
afa8b37b 299 if (AliVAODParticle::GetDebug() > 2)
a5556ea5 300 Info("ReadNext","Particle skipped: Data at vertex not available.");
301 continue;
302 }
303
7814b246 304 if (fMustTPC)
305 {
306 if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
307 {
308 if (AliVAODParticle::GetDebug() > 2)
309 Info("ReadNext","Particle skipped: Was not reconstructed in TPC.");
310 continue;
311 }
312 }
313
a5556ea5 314 if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
315 {
afa8b37b 316 if (AliVAODParticle::GetDebug() > 2)
a5556ea5 317 Info("ReadNext","Particle skipped: PID BIT is not set.");
318 continue;
319 }
320
321
322 Double_t extx;
323 Double_t extp[5];
324 esdtrack->GetConstrainedExternalParameters(extx,extp);
325 if (extp[4] == 0.0)
326 {
afa8b37b 327 if (AliVAODParticle::GetDebug() > 2)
a5556ea5 328 Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
329 continue;
330 }
331 esdtrack->GetESDpid(pidtable);
332 esdtrack->GetConstrainedPxPyPz(mom);
333 esdtrack->GetConstrainedXYZ(pos);
334 pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
335 pos[1] -= vertexpos[1];
336 pos[2] -= vertexpos[2];
337
338 Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
339
340 //Particle from kinematics
afa8b37b 341 AliAODParticle* particle = 0;
a5556ea5 342 Bool_t keeppart = kFALSE;
343 if ( fReadSim && stack )
344 {
345 if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track
346 TParticle *p = stack->Particle(esdtrack->GetLabel());
347 if (p==0x0)
348 {
349 Error("ReadNext","Can not find track with such label.");
350 continue;
351 }
7814b246 352
353 if (fCheckParticlePID)
a5556ea5 354 {
cea0a066 355 if(Rejected(p->GetPdgCode()))
7814b246 356 {
357 if ( AliVAODParticle::GetDebug() > 5 )
358 Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode());
359 continue; //check if we are intersted with particles of this type
360 }
361 }
a5556ea5 362// if(p->GetPdgCode()<0) charge = -1;
afa8b37b 363 particle = new AliAODParticle(*p,i);
a5556ea5 364
365 }
366
367 if(CheckTrack(esdtrack)) continue;
368
369 //Here we apply Bayes' formula
370 Double_t rc=0.;
371 for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
372 if (rc==0.0)
373 {
afa8b37b 374 if (AliVAODParticle::GetDebug() > 2)
a5556ea5 375 Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
376 continue;
377 }
378
379 for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
380
afa8b37b 381 if (AliVAODParticle::GetDebug() > 4)
a5556ea5 382 {
383 Info("ReadNext","###########################################################################");
384 Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
385 Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
386 TString msg("Pid list got from track:");
387 for (Int_t s = 0;s<kNSpecies;s++)
388 {
389 msg+="\n ";
390 msg+=s;
391 msg+="(";
392 msg+=charge*GetSpeciesPdgCode((ESpecies)s);
393 msg+="): ";
394 msg+=w[s];
395 msg+=" (";
396 msg+=pidtable[s];
397 msg+=")";
398 }
399 Info("ReadNext","%s",msg.Data());
afa8b37b 400 }//if (AliVAODParticle::GetDebug()>4)
a5556ea5 401
402 AliTrackPoints* tpts = 0x0;
403 if (fNTrackPoints > 0)
404 {
405 tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
0d8a4589 406// tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
7814b246 407 }
408
409 AliTrackPoints* itstpts = 0x0;
410 if (fITSTrackPoints)
411 {
412 itstpts = new AliTrackPoints(AliTrackPoints::kITS,esdtrack);
0d8a4589 413// itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
a5556ea5 414 }
415
416 AliClusterMap* cmap = 0x0;
417 if ( fClusterMap )
418 {
419 cmap = new AliClusterMap(esdtrack);
420 }
421
422 //If this flag fReadMostProbableOnly is false the
423 //loop over species (see "LOOPPIDS") is over all possible PIDs
424 //in other case the most probablle PID is searched
425 //and the loop is limited to that PID only
426
427 Int_t firstspecie = 0;
428 Int_t lastspecie = kNSpecies;
429
430 if (fReadMostProbableOnly)
431 {
432 //find the most probable PID
433 Int_t spec = 0;
434 Float_t maxprob = w[0];
435 for (Int_t s=1; s<AliESDtrack::kSPECIES; s++)
436 {
437 if (w[s]>maxprob)
438 {
439 maxprob = w[s];
440 spec = s;
441 }
442 }
443 firstspecie = spec;
444 lastspecie = spec + 1;
445 }
446
447 for (Int_t s = firstspecie; s<lastspecie; s++)//LOOPPIDS
448 {
449 Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
450 Float_t pp = w[s];
451 if (pp == 0.0)
452 {
afa8b37b 453 if ( AliVAODParticle::GetDebug() > 5 )
a5556ea5 454 Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
455 continue;
456 }
457
cea0a066 458 if(Rejected(pdgcode))
a5556ea5 459 {
afa8b37b 460 if ( AliVAODParticle::GetDebug() > 5 )
a5556ea5 461 Info("ReadNext","PID (%d) did not pass the cut.",pdgcode);
462 continue; //check if we are intersted with particles of this type
463 }
464
465 Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
466 Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
467
afa8b37b 468 AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i,
a5556ea5 469 mom[0], mom[1], mom[2], tEtot,
470 pos[0], pos[1], pos[2], 0.);
471 //copy probabilitis of other species (if not zero)
472 for (Int_t k = 0; k<kNSpecies; k++)
473 {
474 if (k == s) continue;
475 if (w[k] == 0.0) continue;
476 track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
477 }
478
cea0a066 479 if(Rejected(track))//check if meets all criteria of any of our cuts
a5556ea5 480 //if it does not delete it and take next good track
481 {
afa8b37b 482 if ( AliVAODParticle::GetDebug() > 4 )
a5556ea5 483 Info("ReadNext","Track did not pass the cut");
484 delete track;
485 continue;
486 }
487
488 //Single Particle cuts on cluster map and track points rather do not have sense
489 if (tpts)
490 {
78d7c6d3 491 track->SetTPCTrackPoints(tpts);
a5556ea5 492 }
7814b246 493
494 if (itstpts)
495 {
496 track->SetITSTrackPoints(itstpts);
497 }
a5556ea5 498
499 if (cmap)
500 {
501 track->SetClusterMap(cmap);
502 }
503
504 fEventRec->AddParticle(track);
505 if (particle) fEventSim->AddParticle(particle);
506 keeppart = kTRUE;
507
afa8b37b 508 if (AliVAODParticle::GetDebug() > 4 )
a5556ea5 509 {
510 Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
511 track->Print();
512 if (particle) particle->Print();
513 Info("ReadNext","\n----------------------------------------------\n");
514 }
515 }//for (Int_t s = 0; s<kNSpecies; s++)
516
517 if (keeppart == kFALSE)
518 {
519 delete particle;//particle was not stored in event
520 delete tpts;
7814b246 521 delete itstpts;
a5556ea5 522 delete cmap;
523 }
7814b246 524 else
525 {
f4b459f7 526 if ( fReadSim && stack )
7814b246 527 {
f4b459f7 528 if (particle->P() < 0.00001)
529 {
530 Info("ReadNext","###################################");
531 Info("ReadNext","###################################");
532 Info("ReadNext","Track Label %d",esdtrack->GetLabel());
533 TParticle *p = stack->Particle(esdtrack->GetLabel());
534 Info("ReadNext","");
535 p->Print();
536 Info("ReadNext","");
537 particle->Print();
538 }
539 }
7814b246 540 }
a5556ea5 541
542 }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
543
544 Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
545 fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(),
546 fNEventsRead,fCurrentEvent,fCurrentDir);
547 fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
0d8a4589 548
549 /******************************************************/
550 /****** Setting glevet properties *************/
551 /******************************************************/
552
553 if (fEventRec->GetNumberOfParticles() > 0)
554 {
555 fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
556 }
a5556ea5 557 return 0;
558}
559
560/**********************************************************/
beb1c41d 561Int_t AliReaderESD::ReadESDMuon(AliESD* esd)
562{
563
564 Double_t vertexpos[3];//vertex position, assuming no secondary decay
565
566 const AliESDVertex* vertex = esd->GetVertex();
567
568 if (vertex == 0x0) {
569 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
570 vertexpos[0] = 0.0;
571 vertexpos[1] = 0.0;
572 vertexpos[2] = 0.0;
573 } else {
574 vertex->GetXYZ(vertexpos);
575 }
576
577 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
578
579 if (AliVAODParticle::GetDebug() > 0) {
580 Info("ReadESD","Reading Event %d",fCurrentEvent);
581 Info("ReadESD","Found %d tracks.",nTracks);
582 }
583 // settings
584 Float_t Chi2Cut = 100.;
585 Float_t PtCutMin = 1.;
586 Float_t PtCutMax = 10000.;
587 Float_t muonMass = 0.105658389;
588 Int_t pdgcode = -13;
589 Double_t thetaX, thetaY, pYZ;
590 Double_t pxRec1, pyRec1, pzRec1, E1;
591 Int_t charge;
592
593 Int_t ntrackhits;
594 Double_t fitfmin;
595
596 TLorentzVector fV1;
597 fEventRec->Reset();
598 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
599
600 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
601
602 thetaX = muonTrack->GetThetaX();
603 thetaY = muonTrack->GetThetaY();
604
605 pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
606 pzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
607 pxRec1 = pzRec1 * TMath::Tan(thetaX);
608 pyRec1 = pzRec1 * TMath::Tan(thetaY);
609 charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
610 E1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1);
611 fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, E1);
612
613 ntrackhits = muonTrack->GetNHit();
614 fitfmin = muonTrack->GetChi2();
615
616 // transverse momentum
617 Float_t pt1 = fV1.Pt();
618
619 // chi2 per d.o.f.
620 Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5);
621
622 if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
623 AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack,
624 pxRec1, pyRec1,pzRec1, E1,
625 vertexpos[0], vertexpos[1], vertexpos[2], 0.);
626 fEventRec->AddParticle(track);
627 }
628
629 }
630 fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
631 return 0;
632}
633
634/**********************************************************/
a5556ea5 635
636void AliReaderESD::Rewind()
637{
638 //rewinds reading
beb1c41d 639 // delete fKeyIterator;
a5556ea5 640 delete fFile;
641 fFile = 0x0;
beb1c41d 642 // fKeyIterator = 0x0;
a5556ea5 643 delete fRunLoader;
644 fRunLoader = 0x0;
645 fCurrentDir = 0;
646 fNEventsRead = 0;
647 if (fTrackCounter) fTrackCounter->Reset();
648}
649/**********************************************************/
650
651TFile* AliReaderESD::OpenFile(Int_t n)
652{
beb1c41d 653//opens fFile with tree
a5556ea5 654
655 const TString& dirname = GetDirName(n);
656 if (dirname == "")
657 {
658 Error("OpenFiles","Can not get directory name");
659 return 0x0;
660 }
661 TString filename = dirname +"/"+ fESDFileName;
662 TFile *ret = TFile::Open(filename.Data());
663
664 if ( ret == 0x0)
665 {
666 Error("OpenFiles","Can't open fFile %s",filename.Data());
667 return 0x0;
668 }
669 if (!ret->IsOpen())
670 {
671 Error("OpenFiles","Can't open fFile %s",filename.Data());
672 return 0x0;
673 }
674
675 if (fReadSim )
676 {
677 fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
678 if (fRunLoader == 0x0)
679 {
680 Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
681 delete ret;
682 return 0x0;
683 }
684
685 fRunLoader->LoadHeader();
686 if (fRunLoader->LoadKinematics())
687 {
688 Error("Next","Error occured while loading kinematics.");
689 delete fRunLoader;
690 delete ret;
691 return 0x0;
692 }
693 }
694
695 return ret;
696}
697/**********************************************************/
698
699Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
700{
701 //returns pdg code from the PID index
702 //ask jura about charge
703 switch (spec)
704 {
705 case kESDElectron:
706 return kPositron;
707 break;
708 case kESDMuon:
709 return kMuonPlus;
710 break;
711 case kESDPion:
712 return kPiPlus;
713 break;
714 case kESDKaon:
715 return kKPlus;
716 break;
717 case kESDProton:
718 return kProton;
719 break;
720 default:
721 ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
722 break;
723 }
724 return 0;
725}
726/********************************************************************/
727Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const
728{
729 //Performs check of the track
730
731 if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE;
732
733 if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE;
734
735 if (t->GetTPCclusters(0x0) > 0)
736 {
737 Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0));
738 if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE;
739 }
740
741 Double_t cc[15];
742 t->GetConstrainedExternalCovariance(cc);
743
744 if ( (cc[0] < fC00Min) || (cc[0] > fC00Max) ) return kTRUE;
745 if ( (cc[2] < fC11Min) || (cc[2] > fC11Max) ) return kTRUE;
746 if ( (cc[5] < fC22Min) || (cc[5] > fC22Max) ) return kTRUE;
747 if ( (cc[9] < fC33Min) || (cc[9] > fC33Max) ) return kTRUE;
748 if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;
749
750
751 t->GetInnerExternalCovariance(cc);
752
753 if ( (cc[0] < fTPCC00Min) || (cc[0] > fTPCC00Max) ) return kTRUE;
754 if ( (cc[2] < fTPCC11Min) || (cc[2] > fTPCC11Max) ) return kTRUE;
755 if ( (cc[5] < fTPCC22Min) || (cc[5] > fTPCC22Max) ) return kTRUE;
756 if ( (cc[9] < fTPCC33Min) || (cc[9] > fTPCC33Max) ) return kTRUE;
757 if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE;
758
759 return kFALSE;
760
761}
762/********************************************************************/
763
764void AliReaderESD::SetChi2Range(Float_t min, Float_t max)
765{
766 //sets range of Chi2 per Cluster
767 fChi2Min = min;
768 fChi2Max = max;
769}
770/********************************************************************/
771
772void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
773{
774 //sets range of Number Of Clusters that tracks have to have
775 fNTPCClustMin = min;
776 fNTPCClustMax = max;
777}
778/********************************************************************/
779
780void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
781{
782 //sets range of Chi2 per Cluster
783 fTPCChi2PerClustMin = min;
784 fTPCChi2PerClustMax = max;
785}
786/********************************************************************/
787
788void AliReaderESD::SetC00Range(Float_t min, Float_t max)
789{
790 //Sets range of C00 parameter of covariance matrix of the track
791 //it defines uncertainty of the momentum
792 fC00Min = min;
793 fC00Max = max;
794}
795/********************************************************************/
796
797void AliReaderESD::SetC11Range(Float_t min, Float_t max)
798{
799 //Sets range of C11 parameter of covariance matrix of the track
800 //it defines uncertainty of the momentum
801 fC11Min = min;
802 fC11Max = max;
803}
804/********************************************************************/
805
806void AliReaderESD::SetC22Range(Float_t min, Float_t max)
807{
808 //Sets range of C22 parameter of covariance matrix of the track
809 //it defines uncertainty of the momentum
810 fC22Min = min;
811 fC22Max = max;
812}
813/********************************************************************/
814
815void AliReaderESD::SetC33Range(Float_t min, Float_t max)
816{
817 //Sets range of C33 parameter of covariance matrix of the track
818 //it defines uncertainty of the momentum
819 fC33Min = min;
820 fC33Max = max;
821}
822/********************************************************************/
823
824void AliReaderESD::SetC44Range(Float_t min, Float_t max)
825{
826 //Sets range of C44 parameter of covariance matrix of the track
827 //it defines uncertainty of the momentum
828 fC44Min = min;
829 fC44Max = max;
830}