[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 <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
  //****** 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
  
  if (AliHBTParticle::GetDebug())
    Info("ReadNext","Entered");
    
  TDatabasePDG* pdgdb = TDatabasePDG::Instance();
  if (pdgdb == 0x0)
   {
     Error("Next","Can not get PDG Database Instance.");
     return 1;
   }
  
  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;
      }


     Float_t mf = esd->GetMagneticField(); 
     
     if ( (mf == 0.0) && (fNTrackPoints > 0) )
      {
         Error("ReadNext","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
         fCurrentEvent++;
         continue;
      }
     
     AliStack* stack = 0x0;
     if (fReadParticles && fRunLoader)
      {
        fRunLoader->GetEvent(fCurrentEvent);
        stack = fRunLoader->Stack();
      }
     Info("ReadNext","Reading Event %d",fCurrentEvent);
  
     Int_t ntr = esd->GetNumberOfTracks();
     Info("ReadNext","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;
         }
       
        esdtrack->GetESDpid(pidtable);
        //esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]); 
        esdtrack->GetPxPyPz(mom);
        //esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]);
        esdtrack->GetXYZ(pos);
        
        Double_t extx;
        Double_t extp[5];
        esdtrack->GetExternalParameters(extx,extp);
        
        Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative
        
        //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++)
         {
           Float_t pp = w[s];
           if (pp == 0.0) continue;

           Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
           if(Pass(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
            { 
              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();
            }
         }//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);
      
     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
}
/**********************************************************/
void AliHBTReaderESD::Rewind()
{
  //rewinds reading 
  delete fKeyIterator;
  delete fFile;
  fFile = 0x0;
  fKeyIterator = 0x0;
  delete fRunLoader;
  fRunLoader = 0x0;
  fCurrentDir = 0;
  fNEventsRead = 0;
  fCurrentEvent++;
}
/**********************************************************/

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