[u/mrichter/AliRoot.git] / ANALYSIS / AliReaderESD.cxx

#include "AliReaderESD.h"
//____________________________________________________________________
//////////////////////////////////////////////////////////////////////
//                                                                  //
// class AliReaderESD                                               //
//                                                                  //
// reader for ALICE Event Summary Data (ESD).                       //
//                                                                  //
// Piotr.Skowronski@cern.ch                                         //
//                                                                  //
//////////////////////////////////////////////////////////////////////

#include <TPDGCode.h>
#include <TString.h>
#include <TObjString.h>
#include <TTree.h>
#include <TFile.h>
#include <TKey.h>
#include <TParticle.h>
#include <TH1.h>

#include <AliRun.h>
#include <AliRunLoader.h>
#include <AliStack.h>
#include <AliESDtrack.h>
#include <AliESD.h>

#include "AliAnalysis.h"
#include "AliAODRun.h"
#include "AliAOD.h"
#include "AliAODParticle.h"
#include "AliAODParticleCut.h"
#include "AliTrackPoints.h"
#include "AliClusterMap.h"

ClassImp(AliReaderESD)

AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
 fESDFileName(esdfilename),
 fGAlFileName(galfilename),
 fFile(0x0),
 fRunLoader(0x0),
 fKeyIterator(0x0),
 fReadSim(kFALSE),
 fCheckParticlePID(kFALSE),
 fReadMostProbableOnly(kFALSE),
 fNTrackPoints(0),
 fdR(0.0),
 fClusterMap(kFALSE),
 fNTPCClustMin(0),
 fNTPCClustMax(150),
 fTPCChi2PerClustMin(0.0),
 fTPCChi2PerClustMax(10e5),
 fChi2Min(0.0),
 fChi2Max(10e5),
 fC00Min(0.0),
 fC00Max(10e5),
 fC11Min(0.0),
 fC11Max(10e5),
 fC22Min(0.0),
 fC22Max(10e5),
 fC33Min(0.0),
 fC33Max(10e5),
 fC44Min(0.0),
 fC44Max(10e5),
 fTPCC00Min(0.0),
 fTPCC00Max(10e5),
 fTPCC11Min(0.0),
 fTPCC11Max(10e5),
 fTPCC22Min(0.0),
 fTPCC22Max(10e5),
 fTPCC33Min(0.0),
 fTPCC33Max(10e5),
 fTPCC44Min(0.0),
 fTPCC44Max(10e5)

{
  //cosntructor
  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
    Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/
  
AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
 AliReader(dirs), 
 fESDFileName(esdfilename),
 fGAlFileName(galfilename),
 fFile(0x0),
 fRunLoader(0x0),
 fKeyIterator(0x0),
 fReadSim(kFALSE),
 fCheckParticlePID(kFALSE),
 fReadMostProbableOnly(kFALSE),
 fNTrackPoints(0),
 fdR(0.0),
 fClusterMap(kFALSE),
 fNTPCClustMin(0),
 fNTPCClustMax(150),
 fTPCChi2PerClustMin(0.0),
 fTPCChi2PerClustMax(10e5),
 fChi2Min(0.0),
 fChi2Max(10e5),
 fC00Min(0.0),
 fC00Max(10e5),
 fC11Min(0.0),
 fC11Max(10e5),
 fC22Min(0.0),
 fC22Max(10e5),
 fC33Min(0.0),
 fC33Max(10e5),
 fC44Min(0.0),
 fC44Max(10e5),
 fTPCC00Min(0.0),
 fTPCC00Max(10e5),
 fTPCC11Min(0.0),
 fTPCC11Max(10e5),
 fTPCC22Min(0.0),
 fTPCC22Max(10e5),
 fTPCC33Min(0.0),
 fTPCC33Max(10e5),
 fTPCC44Min(0.0),
 fTPCC44Max(10e5)
{
  //cosntructor
  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
    Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/

AliReaderESD::~AliReaderESD()
{
 //desctructor
  delete fRunLoader;
  delete fKeyIterator;
  delete fFile;
}
/**********************************************************/
Int_t AliReaderESD::ReadNext()
{
//reads next event from fFile
  //fRunLoader is for reading Kine
  
  if (AliVAODParticle::GetDebug())
    Info("ReadNext","Entered");
    
  if (fEventSim == 0x0)  fEventSim = new AliAOD();
  if (fEventRec == 0x0)  fEventRec = new AliAOD();
  
  fEventSim->Reset();
  fEventRec->Reset();
        
  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
   {
     if (fFile == 0x0)
      {
       fFile = OpenFile(fCurrentDir);//rl is opened here
       if (fFile == 0x0)
         {
           Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
           fCurrentDir++;
           continue;
         }
       fCurrentEvent = 0;
       fKeyIterator = new TIter(fFile->GetListOfKeys());  
//       fFile->Dump();
//       fFile->GetListOfKeys()->Print();
      } 
     TKey* key = (TKey*)fKeyIterator->Next();
     if (key == 0x0)
      {
        if (AliVAODParticle::GetDebug() > 2 )
          {
            Info("ReadNext","No more keys.");
          }
        fCurrentDir++;
        delete fKeyIterator;
        fKeyIterator = 0x0;
        delete fFile;//we have to assume there is no more ESD objects in the fFile
        fFile = 0x0;
        delete fRunLoader;
        fRunLoader = 0x0;
        continue;
      }
     //try to read
     
     
//     TObject* esdobj = key->ReadObj();
//     if (esdobj == 0x0)
//      {
//        if (AliVAODParticle::GetDebug() > 2 )
//          {
//            Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
//            key->Dump();
//          }
//        continue;
//      }
//     esdobj->Dump();
//     AliESD* esd = dynamic_cast<AliESD*>(esdobj);
     
     TString esdname = "ESD";
     esdname+=fCurrentEvent;
     AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
     if (esd == 0x0)
      {
//        if (AliVAODParticle::GetDebug() > 2 )
//          {
//            Info("ReadNext","This key is not an AliESD object %s",key->GetName());
//          }
        if (AliVAODParticle::GetDebug() > 2 )
          {
            Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
          }
        fCurrentDir++;
        delete fKeyIterator;
        fKeyIterator = 0x0;
        delete fFile;//we have to assume there is no more ESD objects in the fFile
        fFile = 0x0;
        delete fRunLoader;
        fRunLoader = 0x0;
        continue;
      }
     
     ReadESD(esd);
      
     fCurrentEvent++;
     fNEventsRead++;
     delete esd;
     return 0;//success -> read one event
   }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
   
  return 1; //no more directories to read
}
/**********************************************************/

Int_t AliReaderESD::ReadESD(AliESD* esd)
{
  //****** Tentative particle type "concentrations"
  static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
  
  Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
  Double_t w[kNSpecies];
  Double_t mom[3];//momentum
  Double_t pos[3];//position
  Double_t vertexpos[3];//vertex position
  //Reads one ESD
  if (esd == 0x0)
   {
     Error("ReadESD","ESD is NULL");
     return 1;
   }

  TDatabasePDG* pdgdb = TDatabasePDG::Instance();
  if (pdgdb == 0x0)
   {
     Error("ReadESD","Can not get PDG Database Instance.");
     return 1;
   }
   
  Float_t mf = esd->GetMagneticField(); 

  if ( (mf == 0.0) && (fNTrackPoints > 0) )
   {
      Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
      return 1;
   }

  AliStack* stack = 0x0;
  if (fReadSim && fRunLoader)
   {
     fRunLoader->GetEvent(fCurrentEvent);
     stack = fRunLoader->Stack();
   }

  const AliESDVertex* vertex = esd->GetVertex();
  if (vertex == 0x0)
   {
     Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
     vertexpos[0] = 0.0;
     vertexpos[1] = 0.0;
     vertexpos[2] = 0.0;
   }
  else
   {
     vertex->GetXYZ(vertexpos);
   }
   
  if (AliVAODParticle::GetDebug() > 0)
   {
     Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]);
   }
   
  Info("ReadESD","Reading Event %d",fCurrentEvent);

  Int_t ntr = esd->GetNumberOfTracks();
  Info("ReadESD","Found %d tracks.",ntr);
  for (Int_t i = 0;i<ntr; i++)
   {
     AliESDtrack *esdtrack = esd->GetTrack(i);
     if (esdtrack == 0x0)
      {
        Error("Next","Can not get track %d", i);
        continue;
      }

     //if (esdtrack->HasVertexParameters() == kFALSE) 
     if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
      {
        if (AliVAODParticle::GetDebug() > 2) 
          Info("ReadNext","Particle skipped: Data at vertex not available.");
        continue;
      }

     if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE) 
      {
        if (AliVAODParticle::GetDebug() > 2) 
          Info("ReadNext","Particle skipped: PID BIT is not set.");
        continue;
      }


     Double_t extx;
     Double_t extp[5];
     esdtrack->GetConstrainedExternalParameters(extx,extp);
     if (extp[4] == 0.0)
      {
        if (AliVAODParticle::GetDebug() > 2) 
          Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
        continue;
      } 
     esdtrack->GetESDpid(pidtable);
     esdtrack->GetConstrainedPxPyPz(mom);
     esdtrack->GetConstrainedXYZ(pos);
     pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
     pos[1] -= vertexpos[1];
     pos[2] -= vertexpos[2];

     Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive

     //Particle from kinematics
     AliAODParticle* particle = 0;
     Bool_t keeppart = kFALSE;
     if ( fReadSim && stack  )
      {
        if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track 
        TParticle *p = stack->Particle(esdtrack->GetLabel());
        if (p==0x0) 
         {
           Error("ReadNext","Can not find track with such label.");
           continue;
         }
        if(Pass(p->GetPdgCode())) 
         {
           if ( AliVAODParticle::GetDebug() > 5 )
             Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode());
           continue; //check if we are intersted with particles of this type 
         }
//           if(p->GetPdgCode()<0) charge = -1;
        particle = new AliAODParticle(*p,i);

      }
      
     if(CheckTrack(esdtrack)) continue;
      
     //Here we apply Bayes' formula
     Double_t rc=0.;
     for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
     if (rc==0.0) 
      {
        if (AliVAODParticle::GetDebug() > 2) 
          Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
        continue;
      }

     for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;

     if (AliVAODParticle::GetDebug() > 4)
      { 
        Info("ReadNext","###########################################################################");
        Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
        Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
        TString msg("Pid list got from track:");
        for (Int_t s = 0;s<kNSpecies;s++)
         {
           msg+="\n    ";
           msg+=s;
           msg+="(";
           msg+=charge*GetSpeciesPdgCode((ESpecies)s);
           msg+="): ";
           msg+=w[s];
           msg+=" (";
           msg+=pidtable[s];
           msg+=")";
         }
        Info("ReadNext","%s",msg.Data());
      }//if (AliVAODParticle::GetDebug()>4)

      AliTrackPoints* tpts = 0x0;
      if (fNTrackPoints > 0) 
       {
         tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
       }

      AliClusterMap* cmap = 0x0; 
      if (  fClusterMap ) 
       {
         cmap = new AliClusterMap(esdtrack);
       }
     
     //If this flag fReadMostProbableOnly is false the 
     //loop over species (see "LOOPPIDS") is over all possible PIDs
     //in other case the most probablle PID is searched
     //and the loop is limited to that PID only

     Int_t firstspecie = 0;
     Int_t lastspecie = kNSpecies;
     
     if (fReadMostProbableOnly)
      { 
        //find the most probable PID
        Int_t spec = 0;
        Float_t maxprob = w[0];
        for (Int_t s=1; s<AliESDtrack::kSPECIES; s++) 
         {
           if (w[s]>maxprob)
            {
              maxprob = w[s];
              spec = s;
            }
         }
        firstspecie = spec;
        lastspecie = spec + 1;
      }
     
     for (Int_t s = firstspecie; s<lastspecie; s++)//LOOPPIDS
      {
        Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
        Float_t pp = w[s];
        if (pp == 0.0) 
         {
           if ( AliVAODParticle::GetDebug() > 5 )
             Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
           continue;
         }

        if(Pass(pdgcode)) 
         {
           if ( AliVAODParticle::GetDebug() > 5 )
             Info("ReadNext","PID (%d) did not pass the cut.",pdgcode);
           continue; //check if we are intersted with particles of this type 
         }

        Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
        Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track

        AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i, 
                                                   mom[0], mom[1], mom[2], tEtot,
                                                   pos[0], pos[1], pos[2], 0.);
        //copy probabilitis of other species (if not zero)
        for (Int_t k = 0; k<kNSpecies; k++)
         {
           if (k == s) continue;
           if (w[k] == 0.0) continue;
           track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
         }

        if(Pass(track))//check if meets all criteria of any of our cuts
                       //if it does not delete it and take next good track
         { 
           if ( AliVAODParticle::GetDebug() > 4 )
             Info("ReadNext","Track did not pass the cut");
           delete track;
           continue;
         }

        //Single Particle cuts on cluster map and track points rather do not have sense
        if (tpts)
         {
           track->SetTrackPoints(tpts); 
         }

        if (cmap) 
         { 
           track->SetClusterMap(cmap);
         }

        fEventRec->AddParticle(track);
        if (particle) fEventSim->AddParticle(particle);
        keeppart = kTRUE;

        if (AliVAODParticle::GetDebug() > 4 )
         {
           Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
           track->Print();
           if (particle) particle->Print();
           Info("ReadNext","\n----------------------------------------------\n");
         }
      }//for (Int_t s = 0; s<kNSpecies; s++)

     if (keeppart == kFALSE) 
      {
        delete particle;//particle was not stored in event
        delete tpts;
        delete cmap;
      }

   }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks

  Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
         fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(),
         fNEventsRead,fCurrentEvent,fCurrentDir);
  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
  return 0;
}

/**********************************************************/

void AliReaderESD::Rewind()
{
  //rewinds reading 
  delete fKeyIterator;
  delete fFile;
  fFile = 0x0;
  fKeyIterator = 0x0;
  delete fRunLoader;
  fRunLoader = 0x0;
  fCurrentDir = 0;
  fNEventsRead = 0;
  if (fTrackCounter) fTrackCounter->Reset();
}
/**********************************************************/

TFile* AliReaderESD::OpenFile(Int_t n)
{
//opens fFile with kine tree

 const TString& dirname = GetDirName(n);
 if (dirname == "")
  {
   Error("OpenFiles","Can not get directory name");
   return 0x0;
  }
 TString filename = dirname +"/"+ fESDFileName;
 TFile *ret = TFile::Open(filename.Data()); 

 if ( ret == 0x0)
  {
    Error("OpenFiles","Can't open fFile %s",filename.Data());
    return 0x0;
  }
 if (!ret->IsOpen())
  {
    Error("OpenFiles","Can't open fFile  %s",filename.Data());
    return 0x0;
  }
 
 if (fReadSim )
  {
   fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
   if (fRunLoader == 0x0)
    {
      Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
      delete ret;
      return 0x0;
    }
    
   fRunLoader->LoadHeader();
   if (fRunLoader->LoadKinematics())
    {
      Error("Next","Error occured while loading kinematics.");
      delete fRunLoader;
      delete ret;
      return 0x0;
    }
  }
   
 return ret;
}
/**********************************************************/

Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
{
  //returns pdg code from the PID index
  //ask jura about charge
  switch (spec)
   {
     case kESDElectron:
       return kPositron;
       break;
     case kESDMuon:
       return kMuonPlus;
       break;
     case kESDPion:
       return kPiPlus;
       break;
     case kESDKaon:
       return kKPlus;
       break;
     case kESDProton:
       return kProton;
       break;
     default:
       ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
       break;
   }
  return 0;
}
/********************************************************************/
Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const
{
  //Performs check of the track
  
  if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE;
  
  if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE;

  if (t->GetTPCclusters(0x0) > 0)
   {
     Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0));
     if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE;
   }

  Double_t cc[15];
  t->GetConstrainedExternalCovariance(cc);

  if ( (cc[0]  < fC00Min) || (cc[0]  > fC00Max) ) return kTRUE;
  if ( (cc[2]  < fC11Min) || (cc[2]  > fC11Max) ) return kTRUE;
  if ( (cc[5]  < fC22Min) || (cc[5]  > fC22Max) ) return kTRUE;
  if ( (cc[9]  < fC33Min) || (cc[9]  > fC33Max) ) return kTRUE;
  if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;


  t->GetInnerExternalCovariance(cc);

  if ( (cc[0]  < fTPCC00Min) || (cc[0]  > fTPCC00Max) ) return kTRUE;
  if ( (cc[2]  < fTPCC11Min) || (cc[2]  > fTPCC11Max) ) return kTRUE;
  if ( (cc[5]  < fTPCC22Min) || (cc[5]  > fTPCC22Max) ) return kTRUE;
  if ( (cc[9]  < fTPCC33Min) || (cc[9]  > fTPCC33Max) ) return kTRUE;
  if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE;

  return kFALSE;

}
/********************************************************************/

void AliReaderESD::SetChi2Range(Float_t min, Float_t max)
{
  //sets range of Chi2 per Cluster
  fChi2Min = min;
  fChi2Max = max;
}
/********************************************************************/

void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
{
 //sets range of Number Of Clusters that tracks have to have
 fNTPCClustMin = min;
 fNTPCClustMax = max;
}
/********************************************************************/

void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
{
  //sets range of Chi2 per Cluster
  fTPCChi2PerClustMin = min;
  fTPCChi2PerClustMax = max;
}
/********************************************************************/

void AliReaderESD::SetC00Range(Float_t min, Float_t max)
{
 //Sets range of C00 parameter of covariance matrix of the track
 //it defines uncertainty of the momentum
  fC00Min = min;
  fC00Max = max;
}
/********************************************************************/

void AliReaderESD::SetC11Range(Float_t min, Float_t max)
{
 //Sets range of C11 parameter of covariance matrix of the track
 //it defines uncertainty of the momentum
  fC11Min = min;
  fC11Max = max;
}
/********************************************************************/

void AliReaderESD::SetC22Range(Float_t min, Float_t max)
{
 //Sets range of C22 parameter of covariance matrix of the track
 //it defines uncertainty of the momentum
  fC22Min = min;
  fC22Max = max;
}
/********************************************************************/

void AliReaderESD::SetC33Range(Float_t min, Float_t max)
{
 //Sets range of C33 parameter of covariance matrix of the track
 //it defines uncertainty of the momentum
  fC33Min = min;
  fC33Max = max;
}
/********************************************************************/

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