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