ANALYSIS/AliReaderESD.cxx

   1 #include "AliReaderESD.h"
   2 //____________________________________________________________________
   3 //////////////////////////////////////////////////////////////////////
   4 //                                                                  //
   5 // class AliReaderESD                                               //
   6 //                                                                  //
   7 // reader for ALICE Event Summary Data (ESD).                       //
   8 //                                                                  //
   9 // Piotr.Skowronski@cern.ch                                         //
  10 //                                                                  //
  11 //////////////////////////////////////////////////////////////////////
  12
  13 #include <TPDGCode.h>
  14 #include <TString.h>
  15 #include <TObjString.h>
  16 #include <TTree.h>
  17 #include <TFile.h>
  18 #include <TKey.h>
  19 #include <TParticle.h>
  20 #include <TH1.h>
  21
  22 #include <TGliteXmlEventlist.h>
  23
  24 #include <AliRun.h>
  25 #include <AliRunLoader.h>
  26 #include <AliStack.h>
  27 #include <AliESDtrack.h>
  28 #include <AliKalmanTrack.h>
  29 #include <AliESD.h>
  30
  31 #include "AliAnalysis.h"
  32 #include "AliAODRun.h"
  33 #include "AliAOD.h"
  34 #include "AliAODParticle.h"
  35 #include "AliAODParticleCut.h"
  36 #include "AliClusterMap.h"
  37
  38 ClassImp(AliReaderESD)
  39
  40 AliReaderESD::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),
  52  fITSTrackPoints(kFALSE),
  53  fITSTrackPointsType(AliTrackPoints::kITS),
  54  fMustTPC(kFALSE),
  55  fReadCentralBarrel(kTRUE),
  56  fReadMuon(kFALSE),
  57  fReadPHOS(kFALSE),
  58  fNTPCClustMin(0),
  59  fNTPCClustMax(150),
  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
  87   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
  88     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
  89 }
  90 /********************************************************************/
  91
  92 AliReaderESD::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),
 105  fITSTrackPoints(kFALSE),
 106  fITSTrackPointsType(AliTrackPoints::kITS),
 107  fMustTPC(kFALSE),
 108  fReadCentralBarrel(kTRUE),
 109  fReadMuon(kFALSE),
 110  fReadPHOS(kFALSE),
 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
 139   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
 140     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
 141 }
 142 /********************************************************************/
 143
 144 AliReaderESD::~AliReaderESD()
 145 {
 146  //desctructor
 147     delete fRunLoader;
 148     delete fKeyIterator;
 149     delete fFile;
 150 }
 151
 152 /**********************************************************/
 153 Int_t AliReaderESD::ReadNext()
 154 {
 155 //reads next event from fFile
 156   //fRunLoader is for reading Kine
 157
 158   if (AliVAODParticle::GetDebug())
 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 "./"
 168     {
 169       if (fFile == 0x0)
 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        }
 180      TString esdname = "ESD";
 181      esdname+=fCurrentEvent;
 182      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
 183      if (esd == 0x0)
 184       {
 185         if (AliVAODParticle::GetDebug() > 2 )
 186           {
 187             Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
 188           }
 189         fCurrentDir++;
 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       }
 196       ReadESD(esd);
 197
 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
 203
 204   return 1; //no more directories to read
 205 }
 206
 207 /**********************************************************/
 208 Int_t AliReaderESD::ReadESD(AliESD* esd)
 209 {
 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 /**********************************************************/
 228 Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
 229 {
 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
 239
 240   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
 241   if (pdgdb == 0x0)
 242    {
 243      Error("ReadESD","Can not get PDG Database Instance.");
 244      return 1;
 245    }
 246
 247   Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG
 248
 249   if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
 250    {
 251       Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
 252
 253    }
 254
 255   if (fITSTrackPoints)
 256    {
 257      Info("ReadESD","Magnetic Field is %f",mf);
 258      AliKalmanTrack::SetMagneticField(mf);
 259    }
 260
 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
 281   if (AliVAODParticle::GetDebug() > 0)
 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       {
 302         if (AliVAODParticle::GetDebug() > 2)
 303           Info("ReadNext","Particle skipped: Data at vertex not available.");
 304         continue;
 305       }
 306
 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
 317      if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
 318       {
 319         if (AliVAODParticle::GetDebug() > 2)
 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       {
 330         if (AliVAODParticle::GetDebug() > 2)
 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
 344      AliAODParticle* particle = 0;
 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          }
 355
 356         if (fCheckParticlePID)
 357          {
 358            if(Rejected(p->GetPdgCode()))
 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          }
 365 //           if(p->GetPdgCode()<0) charge = -1;
 366         particle = new AliAODParticle(*p,i);
 367
 368       }
 369
 370      if(CheckTrack(esdtrack)) continue;
 371
 372      //Here we apply Bayes' formula
 373      Double_t rc=0.;
 374      for (Int_t s=0; s<AliPID::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
 375      if (rc==0.0)
 376       {
 377         if (AliVAODParticle::GetDebug() > 2)
 378           Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
 379         continue;
 380       }
 381
 382      for (Int_t s=0; s<AliPID::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
 383
 384      if (AliVAODParticle::GetDebug() > 4)
 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());
 403       }//if (AliVAODParticle::GetDebug()>4)
 404
 405       AliTrackPoints* tpts = 0x0;
 406       if (fNTrackPoints > 0)
 407        {
 408          tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR);
 409 //         tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
 410        }
 411
 412       AliTrackPoints* itstpts = 0x0;
 413       if (fITSTrackPoints)
 414        {
 415          itstpts = new AliTrackPoints(fITSTrackPointsType,esdtrack,mf*10.0);
 416 //         itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
 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];
 438         for (Int_t s=1; s<AliPID::kSPECIES; s++)
 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          {
 456            if ( AliVAODParticle::GetDebug() > 5 )
 457              Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
 458            continue;
 459          }
 460
 461         if(Rejected(pdgcode))
 462          {
 463            if ( AliVAODParticle::GetDebug() > 5 )
 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
 471         AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i,
 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
 482         if(Rejected(track))//check if meets all criteria of any of our cuts
 483                        //if it does not delete it and take next good track
 484          {
 485            if ( AliVAODParticle::GetDebug() > 4 )
 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          {
 494            track->SetTPCTrackPoints(tpts);
 495          }
 496
 497         if (itstpts)
 498          {
 499            track->SetITSTrackPoints(itstpts);
 500          }
 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
 511         if (AliVAODParticle::GetDebug() > 4 )
 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;
 524         delete itstpts;
 525         delete cmap;
 526       }
 527      else
 528       {
 529         if ( fReadSim && stack  )
 530          {
 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          }
 543       }
 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());
 551
 552   /******************************************************/
 553   /******     Setting glevet properties     *************/
 554   /******************************************************/
 555
 556   if (fEventRec->GetNumberOfParticles() > 0)
 557    {
 558      fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
 559    }
 560   return 0;
 561 }
 562
 563 /**********************************************************/
 564 Int_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 /**********************************************************/
 638
 639 void AliReaderESD::Rewind()
 640 {
 641   //rewinds reading
 642   //  delete fKeyIterator;
 643   delete fFile;
 644   fFile = 0x0;
 645   // fKeyIterator = 0x0;
 646   delete fRunLoader;
 647   fRunLoader = 0x0;
 648   fCurrentDir = 0;
 649   fNEventsRead = 0;
 650   if (fEventList) fEventList->Reset();
 651   if (fTrackCounter) fTrackCounter->Reset();
 652 }
 653 /**********************************************************/
 654
 655 TFile* AliReaderESD::OpenFile(Int_t n)
 656 {
 657 //opens fFile with  tree
 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
 674  const TString& dirname = GetDirName(n);
 675  if (dirname == "")
 676   {
 677    Error("OpenFiles","Can not get directory name");
 678    return 0x0;
 679   }
 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());
 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   {
 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
 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();
 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
 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
 748 Int_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 /********************************************************************/
 776 Bool_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
 813 void 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
 821 void 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
 829 void 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
 837 void 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
 846 void 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
 855 void 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
 864 void 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
 873 void 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 }