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