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 <TKey.h>
  33 #include <TObjString.h>
  34 #include <TPDGCode.h>
  35 #include <TParticle.h>
  36 #include <TString.h>
  37 #include <TTree.h>
  38
  39 #include "AliAnalysis.h"
  40 #include "AliAODRun.h"
  41
  42 #include "AliAOD.h"
  43 #include "AliAODParticle.h"
  44 #include "AliAODParticleCut.h"
  45 #include "AliAODRun.h"
  46 #include "AliAnalysis.h"
  47 #include "AliClusterMap.h"
  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"
  56
  57 ClassImp(AliReaderESD)
  58
  59 AliReaderESD::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),
  71  fITSTrackPoints(kFALSE),
  72  fITSTrackPointsType(AliTrackPoints::kITS),
  73  fMustTPC(kFALSE),
  74  fReadCentralBarrel(kTRUE),
  75  fReadMuon(kFALSE),
  76  fReadPHOS(kFALSE),
  77  fNTPCClustMin(0),
  78  fNTPCClustMax(1500),
  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
 106   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
 107     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
 108 }
 109 /********************************************************************/
 110
 111 AliReaderESD::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),
 124  fITSTrackPoints(kFALSE),
 125  fITSTrackPointsType(AliTrackPoints::kITS),
 126  fMustTPC(kFALSE),
 127  fReadCentralBarrel(kTRUE),
 128  fReadMuon(kFALSE),
 129  fReadPHOS(kFALSE),
 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
 158   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
 159     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
 160 }
 161 /********************************************************************/
 162
 163 AliReaderESD::~AliReaderESD()
 164 {
 165  //desctructor
 166     delete fRunLoader;
 167     delete fKeyIterator;
 168     delete fFile;
 169 }
 170
 171 /**********************************************************/
 172 Int_t AliReaderESD::ReadNext()
 173 {
 174 //reads next event from fFile
 175   //fRunLoader is for reading Kine
 176
 177   AliDebug(1,"Entered");
 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 "./"
 186     {
 187       if (fFile == 0x0)
 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        }
 198      TString esdname = "ESD";
 199      esdname+=fCurrentEvent;
 200      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
 201      if (esd == 0x0)
 202       {
 203         AliDebug(3,Form("Can not find AliESD object named %s",esdname.Data()));
 204         fCurrentDir++;
 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       }
 211       ReadESD(esd);
 212
 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
 218
 219   return 1; //no more directories to read
 220 }
 221
 222 /**********************************************************/
 223 Int_t AliReaderESD::ReadESD(AliESD* esd)
 224 {
 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 /**********************************************************/
 243 Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
 244 {
 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
 254
 255   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
 256   if (pdgdb == 0x0)
 257    {
 258      Error("ReadESD","Can not get PDG Database Instance.");
 259      return 1;
 260    }
 261
 262   Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG
 263
 264   if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
 265    {
 266       Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
 267
 268    }
 269
 270   if (fITSTrackPoints)
 271    {
 272      Info("ReadESD","Magnetic Field is %f",mf);
 273      //AliKalmanTrack::SetMagneticField(mf);
 274    }
 275
 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
 296   AliDebug(1,Form("Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]));
 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       {
 314         AliDebug(3,Form("Particle skipped: Data at vertex not available."));
 315         continue;
 316       }
 317
 318      if (fMustTPC)
 319       {
 320        if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
 321         {
 322           AliDebug(3,"Particle skipped: Was not reconstructed in TPC.");
 323           continue;
 324         }
 325       }
 326
 327      if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
 328       {
 329         AliDebug(3,"Particle skipped: PID BIT is not set.");
 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       {
 339         AliDebug(3,"Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
 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
 352      AliAODParticle* particle = 0;
 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          }
 363
 364         if (fCheckParticlePID)
 365          {
 366            if(Rejected(p->GetPdgCode()))
 367             {
 368               AliDebug(6,Form("Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode()));
 369               continue; //check if we are intersted with particles of this type
 370             }
 371          }
 372 //           if(p->GetPdgCode()<0) charge = -1;
 373         particle = new AliAODParticle(*p,i);
 374
 375       }
 376
 377      if(CheckTrack(esdtrack)) continue;
 378
 379      //Here we apply Bayes' formula
 380      Double_t rc=0.;
 381      for (Int_t s=0; s<AliPID::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
 382      if (rc==0.0)
 383       {
 384         AliDebug(3,"Particle rejected since total bayessian PID probab. is zero.");
 385         continue;
 386       }
 387
 388      for (Int_t s=0; s<AliPID::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
 389
 390      if (AliDebugLevel() > 4)
 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());
 409       }//if (AliDebugLevel()>4)
 410
 411       AliTrackPoints* tpts = 0x0;
 412       if (fNTrackPoints > 0)
 413        {
 414          tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR);
 415 //         tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
 416        }
 417
 418       AliTrackPoints* itstpts = 0x0;
 419       if (fITSTrackPoints)
 420        {
 421          itstpts = new AliTrackPoints(fITSTrackPointsType,esdtrack,mf*10.0);
 422 //         itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
 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];
 444         for (Int_t s=1; s<AliPID::kSPECIES; s++)
 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          {
 462            AliDebug(6,Form("Probability of being PID %d is zero. Continuing.",pdgcode));
 463            continue;
 464          }
 465
 466         if(Rejected(pdgcode))
 467          {
 468            AliDebug(6,Form("PID (%d) did not pass the cut.",pdgcode));
 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
 475         AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i,
 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
 486         if(Rejected(track))//check if meets all criteria of any of our cuts
 487                        //if it does not delete it and take next good track
 488          {
 489            AliDebug(5,"Track did not pass the cut");
 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          {
 497            track->SetTPCTrackPoints(tpts);
 498          }
 499
 500         if (itstpts)
 501          {
 502            track->SetITSTrackPoints(itstpts);
 503          }
 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
 514         if (AliDebugLevel() > 4 )
 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;
 527         delete itstpts;
 528         delete cmap;
 529       }
 530      else
 531       {
 532         if ( fReadSim && stack  )
 533          {
 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          }
 546       }
 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());
 554
 555   /******************************************************/
 556   /******     Setting glevet properties     *************/
 557   /******************************************************/
 558
 559   if (fEventRec->GetNumberOfParticles() > 0)
 560    {
 561      fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
 562    }
 563   return 0;
 564 }
 565
 566 /**********************************************************/
 567 Int_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
 585  AliDebug(1,Form("Reading Event %d \nFound %d tracks.",fCurrentEvent,nTracks));
 586
 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 /**********************************************************/
 639
 640 void AliReaderESD::Rewind()
 641 {
 642   //rewinds reading
 643   //  delete fKeyIterator;
 644   delete fFile;
 645   fFile = 0x0;
 646   // fKeyIterator = 0x0;
 647   delete fRunLoader;
 648   fRunLoader = 0x0;
 649   fCurrentDir = 0;
 650   fNEventsRead = 0;
 651   if (fEventList) fEventList->Reset();
 652   if (fTrackCounter) fTrackCounter->Reset();
 653 }
 654 /**********************************************************/
 655
 656 TFile* AliReaderESD::OpenFile(Int_t n)
 657 {
 658 //opens fFile with  tree
 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
 675  const TString& dirname = GetDirName(n);
 676  if (dirname == "")
 677   {
 678    Error("OpenFiles","Can not get directory name");
 679    return 0x0;
 680   }
 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());
 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   {
 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
 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();
 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
 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
 749 Int_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 /********************************************************************/
 777 Bool_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
 814 void 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
 822 void 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
 830 void 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
 838 void 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
 847 void 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
 856 void 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
 865 void 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
 874 void 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 }