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