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