[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderESD.cxx

#include "AliHBTReaderESD.h"
//____________________________________________________________________
//////////////////////////////////////////////////////////////////////
//                                                                  //
// class AliHBTReaderESD                                            //
//                                                                  //
// 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 "AliHBTRun.h"
#include "AliHBTEvent.h"
#include "AliHBTParticle.h"
#include "AliHBTParticleCut.h"
#include "AliHBTTrackPoints.h"
#include "AliHBTClusterMap.h"

ClassImp(AliHBTReaderESD)

AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
 fESDFileName(esdfilename),
 fGAlFileName(galfilename),
 fFile(0x0),
 fRunLoader(0x0),
 fKeyIterator(0x0),
 fReadParticles(kFALSE),
 fNTrackPoints(0),
 fdR(0.0),
 fClusterMap(kFALSE),
 fNTPCClustMin(0),
 fNTPCClustMax(150),
 fTPCChi2PerClustMin(0.0),
 fTPCChi2PerClustMax(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("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/
  
AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
 AliHBTReader(dirs), 
 fESDFileName(esdfilename),
 fGAlFileName(galfilename),
 fFile(0x0),
 fRunLoader(0x0),
 fKeyIterator(0x0),
 fReadParticles(kFALSE),
 fNTrackPoints(0),
 fdR(0.0),
 fClusterMap(kFALSE),
 fNTPCClustMin(0),
 fNTPCClustMax(150),
 fTPCChi2PerClustMin(0.0),
 fTPCChi2PerClustMax(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("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/

AliHBTReaderESD::~AliHBTReaderESD()
{
 //desctructor
  delete fRunLoader;
  delete fKeyIterator;
  delete fFile;
}
/**********************************************************/
Int_t AliHBTReaderESD::ReadNext()
{
//reads next event from fFile
  //fRunLoader is for reading Kine
  
  if (AliHBTParticle::GetDebug())
    Info("ReadNext","Entered");
    
  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
  
  fParticlesEvent->Reset();
  fTracksEvent->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 (AliHBTParticle::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 (AliHBTParticle::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 (AliHBTParticle::GetDebug() > 2 )
//          {
//            Info("ReadNext","This key is not an AliESD object %s",key->GetName());
//          }
        if (AliHBTParticle::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 AliHBTReaderESD::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 (fReadParticles && 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 (AliHBTParticle::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 (AliHBTParticle::GetDebug() > 2) 
          Info("ReadNext","Particle skipped: Data at vertex not available.");
        continue;
      }

     if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE) 
      {
        if (AliHBTParticle::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 (AliHBTParticle::GetDebug() > 2) 
          Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
        continue;
      } 
     esdtrack->GetESDpid(pidtable);
     //esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]); 
     esdtrack->GetConstrainedPxPyPz(mom);
     //esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]);
     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
     AliHBTParticle* particle = 0;
     Bool_t keeppart = kFALSE;
     if ( fReadParticles && 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(p->GetPdgCode()<0) charge = -1;
        particle = new AliHBTParticle(*p,i);

      }

     //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 (AliHBTParticle::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 (AliHBTParticle::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 (AliHBTParticle::GetDebug()>4)

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

      AliHBTClusterMap* cmap = 0x0; 
      if (  fClusterMap ) 
       {
         cmap = new AliHBTClusterMap(esdtrack);
       }

     for (Int_t s = 0; s<kNSpecies; s++)
      {
        Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
        Float_t pp = w[s];
        if (pp == 0.0) 
         {
           if ( AliHBTParticle::GetDebug() > 5 )
             Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
           continue;
         }

        if(Pass(pdgcode)) 
         {
           if ( AliHBTParticle::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

        AliHBTParticle* track = new AliHBTParticle(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 ( AliHBTParticle::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);
         }

        fTracksEvent->AddParticle(track);
        if (particle) fParticlesEvent->AddParticle(particle);
        keeppart = kTRUE;

        if (AliHBTParticle::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).",
         fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
         fNEventsRead,fCurrentEvent,fCurrentDir);
  fTrackCounter->Fill(fTracksEvent->GetNumberOfParticles());
  return 0;
}

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

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

TFile* AliHBTReaderESD::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 (fReadParticles )
  {
   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 AliHBTReaderESD::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;
}
/********************************************************************/

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

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

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