ANALYSIS/AliReaderESD.cxx

   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
  18 //____________________________________________________________________
  19 //////////////////////////////////////////////////////////////////////
  20 //                                                                  //
  21 // class AliReaderESD                                               //
  22 //                                                                  //
  23 // reader for ALICE Event Summary Data (ESD).                       //
  24 //                                                                  //
  25 // Piotr.Skowronski@cern.ch                                         //
  26 //                                                                  //
  27 //////////////////////////////////////////////////////////////////////
  28
  29 #include <TFile.h>
  30 #include <TGliteXmlEventlist.h>
  31 #include <TH1.h>
  32 #include <TPDGCode.h>
  33 #include <TParticle.h>
  34 #include <TString.h>
  35
  36 #include "AliAOD.h"
  37 #include "AliAODParticle.h"
  38 #include "AliClusterMap.h"
  39 #include "AliESD.h"
  40 #include "AliESDtrack.h"
  41 #include "AliLog.h"
  42 #include "AliReaderESD.h"
  43 #include "AliRunLoader.h"
  44 #include "AliStack.h"
  45
  46 ClassImp(AliReaderESD)
  47
  48 AliReaderESD::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),
  60  fITSTrackPoints(kFALSE),
  61  fITSTrackPointsType(AliTrackPoints::kITS),
  62  fMustTPC(kFALSE),
  63  fReadCentralBarrel(kTRUE),
  64  fReadMuon(kFALSE),
  65  fReadPHOS(kFALSE),
  66  fNTPCClustMin(0),
  67  fNTPCClustMax(1500),
  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
  95   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
  96     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
  97 }
  98 /********************************************************************/
  99
 100 AliReaderESD::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),
 113  fITSTrackPoints(kFALSE),
 114  fITSTrackPointsType(AliTrackPoints::kITS),
 115  fMustTPC(kFALSE),
 116  fReadCentralBarrel(kTRUE),
 117  fReadMuon(kFALSE),
 118  fReadPHOS(kFALSE),
 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
 147   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
 148     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
 149 }
 150 /********************************************************************/
 151
 152 AliReaderESD::~AliReaderESD()
 153 {
 154  //desctructor
 155     delete fRunLoader;
 156     delete fKeyIterator;
 157     delete fFile;
 158 }
 159
 160 /**********************************************************/
 161 Int_t AliReaderESD::ReadNext()
 162 {
 163 //reads next event from fFile
 164   //fRunLoader is for reading Kine
 165
 166   AliDebug(1,"Entered");
 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 "./"
 175     {
 176       if (fFile == 0x0)
 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        }
 187      TString esdname = "ESD";
 188      esdname+=fCurrentEvent;
 189      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
 190      if (esd == 0x0)
 191       {
 192         AliDebug(3,Form("Can not find AliESD object named %s",esdname.Data()));
 193         fCurrentDir++;
 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       }
 200       ReadESD(esd);
 201
 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
 207
 208   return 1; //no more directories to read
 209 }
 210
 211 /**********************************************************/
 212 Int_t AliReaderESD::ReadESD(AliESD* esd)
 213 {
 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 /**********************************************************/
 232 Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
 233 {
 234   //****** Tentative particle type "concentrations"
 235   static const Double_t kConcentr[5]={0.05, 0., 0.85, 0.10, 0.05};
 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
 243
 244   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
 245   if (pdgdb == 0x0)
 246    {
 247      Error("ReadESD","Can not get PDG Database Instance.");
 248      return 1;
 249    }
 250
 251   Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG
 252
 253   if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
 254    {
 255       Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
 256
 257    }
 258
 259   if (fITSTrackPoints)
 260    {
 261      Info("ReadESD","Magnetic Field is %f",mf);
 262      //AliKalmanTrack::SetMagneticField(mf);
 263    }
 264
 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
 285   AliDebug(1,Form("Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]));
 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       {
 303         AliDebug(3,Form("Particle skipped: Data at vertex not available."));
 304         continue;
 305       }
 306
 307      if (fMustTPC)
 308       {
 309        if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
 310         {
 311           AliDebug(3,"Particle skipped: Was not reconstructed in TPC.");
 312           continue;
 313         }
 314       }
 315
 316      if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
 317       {
 318         AliDebug(3,"Particle skipped: PID BIT is not set.");
 319         continue;
 320       }
 321
 322
 323      Double_t alpha,extx;
 324      Double_t extp[5];
 325      esdtrack->GetConstrainedExternalParameters(alpha,extx,extp);
 326      if (extp[4] == 0.0)
 327       {
 328         AliDebug(3,"Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
 329         continue;
 330       }
 331      esdtrack->GetESDpid(pidtable);
 332      esdtrack->GetConstrainedPxPyPz(mom);
 333      esdtrack->GetConstrainedXYZ(pos);
 334      pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
 335      pos[1] -= vertexpos[1];
 336      pos[2] -= vertexpos[2];
 337
 338      Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
 339
 340      //Particle from kinematics
 341      AliAODParticle* particle = 0;
 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          }
 352
 353         if (fCheckParticlePID)
 354          {
 355            if(Rejected(p->GetPdgCode()))
 356             {
 357               AliDebug(6,Form("Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode()));
 358               continue; //check if we are intersted with particles of this type
 359             }
 360          }
 361 //           if(p->GetPdgCode()<0) charge = -1;
 362         particle = new AliAODParticle(*p,i);
 363
 364       }
 365
 366      if(CheckTrack(esdtrack)) continue;
 367
 368      //Here we apply Bayes' formula
 369      Double_t rc=0.;
 370      for (Int_t s=0; s<AliPID::kSPECIES; s++) rc+=kConcentr[s]*pidtable[s];
 371      if (rc==0.0)
 372       {
 373         AliDebug(3,"Particle rejected since total bayessian PID probab. is zero.");
 374         continue;
 375       }
 376
 377      for (Int_t s=0; s<AliPID::kSPECIES; s++) w[s]=kConcentr[s]*pidtable[s]/rc;
 378
 379      if (AliDebugLevel() > 4)
 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());
 398       }//if (AliDebugLevel()>4)
 399
 400       AliTrackPoints* tpts = 0x0;
 401       if (fNTrackPoints > 0)
 402        {
 403          tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR);
 404 //         tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
 405        }
 406
 407       AliTrackPoints* itstpts = 0x0;
 408       if (fITSTrackPoints)
 409        {
 410          itstpts = new AliTrackPoints(fITSTrackPointsType,esdtrack,mf*10.0);
 411 //         itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
 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];
 433         for (Int_t s=1; s<AliPID::kSPECIES; s++)
 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          {
 451            AliDebug(6,Form("Probability of being PID %d is zero. Continuing.",pdgcode));
 452            continue;
 453          }
 454
 455         if(Rejected(pdgcode))
 456          {
 457            AliDebug(6,Form("PID (%d) did not pass the cut.",pdgcode));
 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
 464         AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i,
 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
 475         if(Rejected(track))//check if meets all criteria of any of our cuts
 476                        //if it does not delete it and take next good track
 477          {
 478            AliDebug(5,"Track did not pass the cut");
 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          {
 486            track->SetTPCTrackPoints(tpts);
 487          }
 488
 489         if (itstpts)
 490          {
 491            track->SetITSTrackPoints(itstpts);
 492          }
 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
 503         if (AliDebugLevel() > 4 )
 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;
 516         delete itstpts;
 517         delete cmap;
 518       }
 519      else
 520       {
 521         if ( fReadSim && stack  )
 522          {
 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          }
 535       }
 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());
 543
 544   /******************************************************/
 545   /******     Setting glevet properties     *************/
 546   /******************************************************/
 547
 548   if (fEventRec->GetNumberOfParticles() > 0)
 549    {
 550      fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
 551    }
 552   return 0;
 553 }
 554
 555 /**********************************************************/
 556 Int_t AliReaderESD::ReadESDMuon(AliESD* esd)
 557 {
 558   // Reads the muon tracks from the ESD
 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
 574  AliDebug(1,Form("Reading Event %d \nFound %d tracks.",fCurrentEvent,nTracks));
 575
 576  // settings
 577   Float_t chi2Cut = 100.;
 578   Float_t ptCutMin = 1.;
 579   Float_t ptCutMax = 10000.;
 580   Float_t muonMass = 0.105658389;
 581   Int_t pdgcode = -13;
 582   Double_t thetaX, thetaY, pYZ;
 583   Double_t pxRec1, pyRec1, pzRec1, e1;
 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()));
 603       e1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1);
 604       fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, e1);
 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
 615       if ((ch1 < chi2Cut) && (pt1 > ptCutMin) && (pt1 < ptCutMax)) {
 616         AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack,
 617                                                    pxRec1, pyRec1,pzRec1, e1,
 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 /**********************************************************/
 628
 629 void AliReaderESD::Rewind()
 630 {
 631   //rewinds reading
 632   //  delete fKeyIterator;
 633   delete fFile;
 634   fFile = 0x0;
 635   // fKeyIterator = 0x0;
 636   delete fRunLoader;
 637   fRunLoader = 0x0;
 638   fCurrentDir = 0;
 639   fNEventsRead = 0;
 640   if (fEventList) fEventList->Reset();
 641   if (fTrackCounter) fTrackCounter->Reset();
 642 }
 643 /**********************************************************/
 644
 645 TFile* AliReaderESD::OpenFile(Int_t n)
 646 {
 647 //opens fFile with  tree
 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
 664  const TString& dirname = GetDirName(n);
 665  if (dirname == "")
 666   {
 667    Error("OpenFiles","Can not get directory name");
 668    return 0x0;
 669   }
 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());
 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   {
 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
 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();
 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
 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
 738 Int_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 /********************************************************************/
 766 Bool_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
 803 void 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
 811 void 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
 819 void 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
 827 void 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
 836 void 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
 845 void 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
 854 void 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
 863 void 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 }