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