[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoEventReaderESDKine.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoEventReaderESDKine - the reader class for the Alice ESD              ///
/// Reads in ESD information and converts it into internal AliFemtoEvent     ///
/// Reads in AliESDfriend to create shared hit/quality information           ///
/// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl                        ///
///          Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////

/*
 *$Id$
 *$Log$
 *Revision 1.1  2007/05/25 12:42:54  akisiel
 *Adding a reader for the Kine information
 *
 *Revision 1.2  2007/05/22 09:01:42  akisiel
 *Add the possibiloity to save cut settings in the ROOT file
 *
 *Revision 1.1  2007/05/16 10:22:11  akisiel
 *Making the directory structure of AliFemto flat. All files go into one common directory
 *
 *Revision 1.5  2007/05/03 09:45:20  akisiel
 *Fixing Effective C++ warnings
 *
 *Revision 1.4  2007/04/27 07:28:34  akisiel
 *Remove event number reading due to interface changes
 *
 *Revision 1.3  2007/04/27 07:25:16  akisiel
 *Make revisions needed for compilation from the main AliRoot tree
 *
 *Revision 1.1.1.1  2007/04/25 15:38:41  panos
 *Importing the HBT code dir
 *
 */

#include "AliFemtoEventReaderESDKine.h"

#include "TFile.h"
#include "TChain.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliStack.h"
//#include "AliAODParticle.h"
#include "TParticle.h"

//#include "TSystem.h"

#include "AliFmPhysicalHelixD.h"
#include "AliFmThreeVectorF.h"

#include "SystemOfUnits.h"

#include "AliFemtoEvent.h"

#include "TMath.h"

ClassImp(AliFemtoEventReaderESDKine)

#if !(ST_NO_NAMESPACES)
  using namespace units;
#endif

using namespace std;
//____________________________
//constructor with 0 parameters , look at default settings 
AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine():
  fInputFile(" "),
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurRLEvent(0),
  fTree(0x0),
  fEvent(0x0),
  fRunLoader(0x0)
{
}

AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
  fInputFile(" "),
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurRLEvent(0),
  fTree(0x0),
  fEvent(0x0),
  fRunLoader(0x0)
{
  // copy constructor
  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fEvent = new AliESDEvent();
}
//__________________
//Destructor
AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
{
  // destructor
  //delete fListOfFiles;
  delete fTree;
  delete fEvent;
  if (fRunLoader) delete fRunLoader;
}

//__________________
AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemtoEventReaderESDKine& aReader)
{
  // assignment operator
  if (this == &aReader)
    return *this;

  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurRLEvent = aReader.fCurRLEvent;
  if (fTree) delete fTree;
  //  fTree = aReader.fTree->CloneTree();
  if (fEvent) delete fEvent;
  fEvent = new AliESDEvent();
  if (fRunLoader) delete fRunLoader;
  fRunLoader = new AliRunLoader(*aReader.fRunLoader);

  return *this;
}
//__________________
AliFemtoString AliFemtoEventReaderESDKine::Report() const
{
  // create reader report
  AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
  return temp;
}

//__________________
void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
{
  //setting the name of file where names of ESD file are written 
  //it takes only this files which have good trees
  char buffer[256];
  fInputFile=string(inputFile);
  cout<<"Input File set on "<<fInputFile<<endl;
  ifstream infile(inputFile);

  fTree = new TChain("esdTree");

  if(infile.good()==true)
    { 
      //checking if all give files have good tree inside
      while (infile.eof()==false)
	{
	  infile.getline(buffer,256);
	  //ifstream test_file(buffer);
	  TFile *esdFile=TFile::Open(buffer,"READ");
	  if (esdFile!=0x0)
	    {	
	      TTree* tree = (TTree*) esdFile->Get("esdTree");
	      if (tree!=0x0)
		{
		  cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
		  fTree->AddFile(buffer);
		  delete tree;
		}
	      esdFile->Close();	
	    }
	  delete esdFile;
	}
    }
}

void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
{
  fConstrained=constrained;
}

bool AliFemtoEventReaderESDKine::GetConstrained() const
{
  return fConstrained;
}

AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
{
  // read in a next hbt event from the chain
  // convert it to AliFemtoEvent and return
  // for further analysis
  AliFemtoEvent *hbtEvent = 0;
  TString tGAliceFilename;

  if (fCurEvent==fNumberofEvent)//open next file  
    {
      if (fNumberofEvent == 0) {
	fEvent=new AliESDEvent();
		
	  //ESD data
// 	  fEsdFile=TFile::Open(fFileName.c_str(),"READ");
// 	  fTree = (TTree*) fEsdFile->Get("esdTree");			

	  fTree->SetBranchStatus("MuonTracks*",0);
	  fTree->SetBranchStatus("PmdTracks*",0);
	  fTree->SetBranchStatus("TrdTracks*",0);
	  fTree->SetBranchStatus("V0s*",0);
	  fTree->SetBranchStatus("Cascades*",0);
	  fTree->SetBranchStatus("Kinks*",0);
	  fTree->SetBranchStatus("CaloClusters*",0);
	  fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
	  fTree->SetBranchStatus("ESDfriend*",0);
	  fEvent->ReadFromTree(fTree);

// 	  chain->SetBranchStatus("*",0);
// 	  chain->SetBranchStatus("fUniqueID",1);
// 	  chain->SetBranchStatus("fTracks",1);
// 	  chain->SetBranchStatus("fTracks.*",1);
// 	  chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
//	  fTree->SetBranchStatus("fTracks.fCalibContainer",0);


	fNumberofEvent=fTree->GetEntries();

	if (fNumberofEvent == 0) {
	  cout<<"no event in input "<<endl;
	  fReaderStatus=1;
	  return hbtEvent; 
	}

	cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
	fCurEvent=0;
	// simulation data reading setup
	
      }
      else //no more data to read
	{
	  cout<<"no more files "<<hbtEvent<<endl;
	  fReaderStatus=1;
	  return hbtEvent; 
	}
    }		
  cout<<"starting to read event "<<fCurEvent<<endl;
  fTree->GetEvent(fCurEvent);//getting next event
  //  vector<int> tLabelTable;//to check labels
  
  cout << "fFileName is " << fFileName.Data() << endl;
  cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
  if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
    fFileName = fTree->GetCurrentFile()->GetName();
    tGAliceFilename = fFileName;
    tGAliceFilename.ReplaceAll("AliESDs","galice");
    cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
    if (fRunLoader) delete fRunLoader;
    fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
    if (fRunLoader==0x0)
      {
	cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
	exit(0);
      }
    if(fRunLoader->LoadHeader())
      {
	cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
	exit(0);
      }
    fRunLoader->LoadKinematics();
    fCurRLEvent = 0;
  }

  fRunLoader->GetEvent(fCurRLEvent);
  AliStack* tStack = 0x0;
  tStack = fRunLoader->Stack();
	
  hbtEvent = new AliFemtoEvent;
  //setting basic things
  //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
  hbtEvent->SetRunNumber(fEvent->GetRunNumber());
  //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
  hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
  hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
  hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
  hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
  hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
  hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
  hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
  hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
  hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	
  //Vertex
  double fV1[3];
  fEvent->GetVertex()->GetXYZ(fV1);

  AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
  hbtEvent->SetPrimVertPos(vertex);
	
  //starting to reading tracks
  int nofTracks=0;  //number of reconstructed tracks in event
  nofTracks=fEvent->GetNumberOfTracks();
  int realnofTracks=0;//number of track which we use ina analysis

  for (int i=0;i<nofTracks;i++)
    {
      bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
		
      AliFemtoTrack* trackCopy = new AliFemtoTrack();	
      const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
      //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();

      trackCopy->SetCharge((short)esdtrack->GetSign());

      //in aliroot we have AliPID 
      //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
      //we use only 5 first
      double esdpid[5];
      esdtrack->GetESDpid(esdpid);
      trackCopy->SetPidProbElectron(esdpid[0]);
      trackCopy->SetPidProbMuon(esdpid[1]);
      trackCopy->SetPidProbPion(esdpid[2]);
      trackCopy->SetPidProbKaon(esdpid[3]);
      trackCopy->SetPidProbProton(esdpid[4]);
						
      double pxyz[3];
      if (fConstrained==true)		    
	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
      else
	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
      AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
      trackCopy->SetP(v);//setting momentum
      trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
      const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
      if (ktP.mag() == 0) {
	delete trackCopy;
	continue;
      }
      const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
      //setting helix I do not if it is ok
      AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
      trackCopy->SetHelix(helix);
	    	
      trackCopy->SetTrackId(esdtrack->GetID());
      trackCopy->SetFlags(esdtrack->GetStatus());
      trackCopy->SetLabel(esdtrack->GetLabel());
		
      //some stuff which could be useful 
      float impact[2];
      float covimpact[3];
      esdtrack->GetImpactParameters(impact,covimpact);
      trackCopy->SetImpactD(impact[0]);
      trackCopy->SetImpactZ(impact[1]);
      trackCopy->SetCdd(covimpact[0]);
      trackCopy->SetCdz(covimpact[1]);
      trackCopy->SetCzz(covimpact[2]);
      trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
      trackCopy->SetITSncls(esdtrack->GetNcls(0));     
      trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
      trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
      trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
      trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
      trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 
      trackCopy->SetSigmaToVertex(GetSigmaToVertex(esdtrack));

      trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
      trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());

      double pvrt[3];
      fEvent->GetPrimaryVertex()->GetXYZ(pvrt);

      double xtpc[3];
      esdtrack->GetInnerXYZ(xtpc);
      xtpc[2] -= pvrt[2];
      trackCopy->SetNominalTPCEntrancePoint(xtpc);

      esdtrack->GetOuterXYZ(xtpc);
      xtpc[2] -= pvrt[2];
      trackCopy->SetNominalTPCExitPoint(xtpc);

      int indexes[3];
      for (int ik=0; ik<3; ik++) {
	indexes[ik] = esdtrack->GetKinkIndex(ik);
      }
      trackCopy->SetKinkIndexes(indexes);

      // Fill the hidden information with the simulated data
      TParticle *tPart = tStack->Particle(TMath::Abs(esdtrack->GetLabel()));
      //      AliAODParticle* tParticle= new AliAODParticle(*tPart,i);
      AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
      tInfo->SetPDGPid(tPart->GetPdgCode());
      tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
      Double_t mass2 = (tPart->Energy() *tPart->Energy() -
			tPart->Px()*tPart->Px() -
			tPart->Py()*tPart->Py() -
			tPart->Pz()*tPart->Pz());
      if (mass2>0.0)
	tInfo->SetMass(TMath::Sqrt(mass2));
      else 
	tInfo->SetMass(0.0);
      tInfo->SetEmissionPoint(tPart->Vx()*1e13, tPart->Vy()*1e13, tPart->Vz()*1e13, tPart->T()*1e13*300000);
      trackCopy->SetHiddenInfo(tInfo);

      //decision if we want this track
      //if we using diffrent labels we want that this label was use for first time 
      //if we use hidden info we want to have match between sim data and ESD
      if (tGoodMomentum==true)
	{
	  hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	  realnofTracks++;//real number of tracks
	}
      else
	{
	  delete  trackCopy;
	}
		
    }

  hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event	
  fCurEvent++;	
  fCurRLEvent++;
  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
  if (fCurEvent== fNumberofEvent)//if end of current file close all
    {   
      fTree->Reset(); 
      delete fTree;
    }
  return hbtEvent; 
}
//____________________________________________________________________
Float_t AliFemtoEventReaderESDKine::GetSigmaToVertex(const AliESDtrack* esdTrack)
{
  // Calculates the number of sigma to the vertex.

  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];
  esdTrack->GetImpactParameters(b,bCov);
//   if (bCov[0]<=0 || bCov[2]<=0) {
//     AliDebug(1, "Estimated b resolution lower or equal zero!");
//     bCov[0]=0; bCov[2]=0;
//   }
  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);

  // -----------------------------------
  // How to get to a n-sigma cut?
  //
  // The accumulated statistics from 0 to d is
  //
  // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
  // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
  //
  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
  // Can this be expressed in a different way?

  if (bRes[0] == 0 || bRes[1] ==0)
    return -1;

  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));

  // stupid rounding problem screws up everything:
  // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
  if (TMath::Exp(-d * d / 2) < 1e-10)
    return 1000;

  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  return d;
}
Commit	Line	Data
f007751a	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoEventReaderESDKine - the reader class for the Alice ESD ///
	4	/// Reads in ESD information and converts it into internal AliFemtoEvent ///
	5	/// Reads in AliESDfriend to create shared hit/quality information ///
	6	/// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl ///
	7	/// Adam Kisiel kisiel@mps.ohio-state.edu ///
	8	/// ///
	9	////////////////////////////////////////////////////////////////////////////////
	10
	11	/*
	12	*$Id$
	13	*$Log$
cc5faabc	14	*Revision 1.1 2007/05/25 12:42:54 akisiel
	15	*Adding a reader for the Kine information
	16	*
f007751a	17	*Revision 1.2 2007/05/22 09:01:42 akisiel
	18	*Add the possibiloity to save cut settings in the ROOT file
	19	*
	20	*Revision 1.1 2007/05/16 10:22:11 akisiel
	21	*Making the directory structure of AliFemto flat. All files go into one common directory
	22	*
	23	*Revision 1.5 2007/05/03 09:45:20 akisiel
	24	*Fixing Effective C++ warnings
	25	*
	26	*Revision 1.4 2007/04/27 07:28:34 akisiel
	27	*Remove event number reading due to interface changes
	28	*
	29	*Revision 1.3 2007/04/27 07:25:16 akisiel
	30	*Make revisions needed for compilation from the main AliRoot tree
	31	*
	32	*Revision 1.1.1.1 2007/04/25 15:38:41 panos
	33	*Importing the HBT code dir
	34	*
	35	*/
	36
	37	#include "AliFemtoEventReaderESDKine.h"
	38
	39	#include "TFile.h"
cc5faabc	40	#include "TChain.h"
cc5faabc	41	#include "AliESDEvent.h"
f007751a	42	#include "AliESDtrack.h"
cc5faabc	43	#include "AliESDVertex.h"
f007751a	44	#include "AliStack.h"
0b3bd1ac	45	//#include "AliAODParticle.h"
f007751a	46	#include "TParticle.h"
	47
	48	//#include "TSystem.h"
	49
	50	#include "AliFmPhysicalHelixD.h"
	51	#include "AliFmThreeVectorF.h"
	52
	53	#include "SystemOfUnits.h"
	54
	55	#include "AliFemtoEvent.h"
	56
0b3bd1ac	57	#include "TMath.h"
0b3bd1ac	58
f007751a	59	ClassImp(AliFemtoEventReaderESDKine)
	60
	61	#if !(ST_NO_NAMESPACES)
	62	using namespace units;
	63	#endif
	64
	65	using namespace std;
	66	//____________________________
	67	//constructor with 0 parameters , look at default settings
	68	AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine():
	69	fInputFile(" "),
	70	fFileName(" "),
	71	fConstrained(true),
	72	fNumberofEvent(0),
	73	fCurEvent(0),
cc5faabc	74	fCurRLEvent(0),
f007751a	75	fTree(0x0),
f007751a	76	fEvent(0x0),
cc5faabc	77	fRunLoader(0x0)
f007751a	78	{
f007751a	79	}
	80
	81	AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
	82	fInputFile(" "),
	83	fFileName(" "),
	84	fConstrained(true),
	85	fNumberofEvent(0),
	86	fCurEvent(0),
cc5faabc	87	fCurRLEvent(0),
f007751a	88	fTree(0x0),
f007751a	89	fEvent(0x0),
cc5faabc	90	fRunLoader(0x0)
f007751a	91	{
	92	// copy constructor
	93	fInputFile = aReader.fInputFile;
	94	fFileName = aReader.fFileName;
	95	fConstrained = aReader.fConstrained;
	96	fNumberofEvent = aReader.fNumberofEvent;
	97	fCurEvent = aReader.fCurEvent;
cc5faabc	98	fEvent = new AliESDEvent();
f007751a	99	}
	100	//__________________
	101	//Destructor
	102	AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
	103	{
	104	// destructor
	105	//delete fListOfFiles;
	106	delete fTree;
	107	delete fEvent;
f007751a	108	if (fRunLoader) delete fRunLoader;
f007751a	109	}
	110
	111	//__________________
	112	AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemtoEventReaderESDKine& aReader)
	113	{
	114	// assignment operator
	115	if (this == &aReader)
	116	return *this;
	117
	118	fInputFile = aReader.fInputFile;
	119	fFileName = aReader.fFileName;
	120	fConstrained = aReader.fConstrained;
	121	fNumberofEvent = aReader.fNumberofEvent;
	122	fCurEvent = aReader.fCurEvent;
cc5faabc	123	fCurRLEvent = aReader.fCurRLEvent;
f007751a	124	if (fTree) delete fTree;
cc5faabc	125	// fTree = aReader.fTree->CloneTree();
f007751a	126	if (fEvent) delete fEvent;
cc5faabc	127	fEvent = new AliESDEvent();
f007751a	128	if (fRunLoader) delete fRunLoader;
	129	fRunLoader = new AliRunLoader(*aReader.fRunLoader);
	130
f007751a	131	return *this;
	132	}
	133	//__________________
cc5faabc	134	AliFemtoString AliFemtoEventReaderESDKine::Report() const
f007751a	135	{
	136	// create reader report
	137	AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
	138	return temp;
	139	}
	140
	141	//__________________
	142	void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
	143	{
	144	//setting the name of file where names of ESD file are written
	145	//it takes only this files which have good trees
	146	char buffer[256];
	147	fInputFile=string(inputFile);
	148	cout<<"Input File set on "<<fInputFile<<endl;
	149	ifstream infile(inputFile);
cc5faabc	150
	151	fTree = new TChain("esdTree");
	152
f007751a	153	if(infile.good()==true)
	154	{
	155	//checking if all give files have good tree inside
	156	while (infile.eof()==false)
	157	{
	158	infile.getline(buffer,256);
	159	//ifstream test_file(buffer);
	160	TFile *esdFile=TFile::Open(buffer,"READ");
	161	if (esdFile!=0x0)
	162	{
	163	TTree* tree = (TTree*) esdFile->Get("esdTree");
	164	if (tree!=0x0)
	165	{
	166	cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
cc5faabc	167	fTree->AddFile(buffer);
f007751a	168	delete tree;
	169	}
	170	esdFile->Close();
	171	}
	172	delete esdFile;
	173	}
	174	}
	175	}
	176
f007751a	177	void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
	178	{
	179	fConstrained=constrained;
	180	}
	181
	182	bool AliFemtoEventReaderESDKine::GetConstrained() const
	183	{
	184	return fConstrained;
	185	}
	186
	187	AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
	188	{
	189	// read in a next hbt event from the chain
	190	// convert it to AliFemtoEvent and return
	191	// for further analysis
	192	AliFemtoEvent *hbtEvent = 0;
f007751a	193	TString tGAliceFilename;
	194
	195	if (fCurEvent==fNumberofEvent)//open next file
	196	{
cc5faabc	197	if (fNumberofEvent == 0) {
cc5faabc	198	fEvent=new AliESDEvent();
f007751a	199
f007751a	200	//ESD data
cc5faabc	201	// fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	202	// fTree = (TTree*) fEsdFile->Get("esdTree");
	203
	204	fTree->SetBranchStatus("MuonTracks*",0);
	205	fTree->SetBranchStatus("PmdTracks*",0);
	206	fTree->SetBranchStatus("TrdTracks*",0);
	207	fTree->SetBranchStatus("V0s*",0);
	208	fTree->SetBranchStatus("Cascades*",0);
	209	fTree->SetBranchStatus("Kinks*",0);
	210	fTree->SetBranchStatus("CaloClusters*",0);
	211	fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
	212	fTree->SetBranchStatus("ESDfriend*",0);
	213	fEvent->ReadFromTree(fTree);
f007751a	214
	215	// chain->SetBranchStatus("*",0);
	216	// chain->SetBranchStatus("fUniqueID",1);
	217	// chain->SetBranchStatus("fTracks",1);
	218	// chain->SetBranchStatus("fTracks.*",1);
	219	// chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
cc5faabc	220	// fTree->SetBranchStatus("fTracks.fCalibContainer",0);
	221
	222
	223	fNumberofEvent=fTree->GetEntries();
	224
	225	if (fNumberofEvent == 0) {
	226	cout<<"no event in input "<<endl;
	227	fReaderStatus=1;
	228	return hbtEvent;
f007751a	229	}
cc5faabc	230
	231	cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
	232	fCurEvent=0;
	233	// simulation data reading setup
	234
	235	}
f007751a	236	else //no more data to read
	237	{
	238	cout<<"no more files "<<hbtEvent<<endl;
	239	fReaderStatus=1;
	240	return hbtEvent;
	241	}
	242	}
	243	cout<<"starting to read event "<<fCurEvent<<endl;
	244	fTree->GetEvent(fCurEvent);//getting next event
f007751a	245	// vector<int> tLabelTable;//to check labels
f007751a	246
cc5faabc	247	cout << "fFileName is " << fFileName.Data() << endl;
	248	cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
	249	if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
	250	fFileName = fTree->GetCurrentFile()->GetName();
	251	tGAliceFilename = fFileName;
	252	tGAliceFilename.ReplaceAll("AliESDs","galice");
	253	cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
	254	if (fRunLoader) delete fRunLoader;
	255	fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
	256	if (fRunLoader==0x0)
	257	{
	258	cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
	259	exit(0);
	260	}
	261	if(fRunLoader->LoadHeader())
	262	{
	263	cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
	264	exit(0);
	265	}
	266	fRunLoader->LoadKinematics();
	267	fCurRLEvent = 0;
	268	}
	269
	270	fRunLoader->GetEvent(fCurRLEvent);
f007751a	271	AliStack* tStack = 0x0;
	272	tStack = fRunLoader->Stack();
	273
	274	hbtEvent = new AliFemtoEvent;
	275	//setting basic things
	276	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
	277	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	278	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	279	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	280	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	281	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	282	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	283	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	284	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	285	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	286	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	287	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	288
	289	//Vertex
	290	double fV1[3];
	291	fEvent->GetVertex()->GetXYZ(fV1);
	292
	293	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	294	hbtEvent->SetPrimVertPos(vertex);
	295
	296	//starting to reading tracks
	297	int nofTracks=0; //number of reconstructed tracks in event
	298	nofTracks=fEvent->GetNumberOfTracks();
	299	int realnofTracks=0;//number of track which we use ina analysis
	300
f007751a	301	for (int i=0;i<nofTracks;i++)
	302	{
	303	bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
	304
	305	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	306	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	307	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
	308
	309	trackCopy->SetCharge((short)esdtrack->GetSign());
	310
	311	//in aliroot we have AliPID
	312	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	313	//we use only 5 first
	314	double esdpid[5];
	315	esdtrack->GetESDpid(esdpid);
	316	trackCopy->SetPidProbElectron(esdpid[0]);
	317	trackCopy->SetPidProbMuon(esdpid[1]);
	318	trackCopy->SetPidProbPion(esdpid[2]);
	319	trackCopy->SetPidProbKaon(esdpid[3]);
	320	trackCopy->SetPidProbProton(esdpid[4]);
	321
	322	double pxyz[3];
	323	if (fConstrained==true)
	324	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
	325	else
	326	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
	327	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	328	trackCopy->SetP(v);//setting momentum
	329	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
	330	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
	331	if (ktP.mag() == 0) {
	332	delete trackCopy;
	333	continue;
	334	}
cc5faabc	335	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
f007751a	336	//setting helix I do not if it is ok
cc5faabc	337	AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
f007751a	338	trackCopy->SetHelix(helix);
	339
	340	trackCopy->SetTrackId(esdtrack->GetID());
	341	trackCopy->SetFlags(esdtrack->GetStatus());
cc5faabc	342	trackCopy->SetLabel(esdtrack->GetLabel());
f007751a	343
	344	//some stuff which could be useful
	345	float impact[2];
	346	float covimpact[3];
	347	esdtrack->GetImpactParameters(impact,covimpact);
	348	trackCopy->SetImpactD(impact[0]);
	349	trackCopy->SetImpactZ(impact[1]);
	350	trackCopy->SetCdd(covimpact[0]);
	351	trackCopy->SetCdz(covimpact[1]);
	352	trackCopy->SetCzz(covimpact[2]);
	353	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
	354	trackCopy->SetITSncls(esdtrack->GetNcls(0));
	355	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
	356	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
	357	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
	358	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
	359	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
0b3bd1ac	360	trackCopy->SetSigmaToVertex(GetSigmaToVertex(esdtrack));
f007751a	361
cc5faabc	362	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
	363	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
	364
0b3bd1ac	365	double pvrt[3];
	366	fEvent->GetPrimaryVertex()->GetXYZ(pvrt);
	367
	368	double xtpc[3];
	369	esdtrack->GetInnerXYZ(xtpc);
	370	xtpc[2] -= pvrt[2];
	371	trackCopy->SetNominalTPCEntrancePoint(xtpc);
	372
	373	esdtrack->GetOuterXYZ(xtpc);
	374	xtpc[2] -= pvrt[2];
	375	trackCopy->SetNominalTPCExitPoint(xtpc);
	376
	377	int indexes[3];
	378	for (int ik=0; ik<3; ik++) {
	379	indexes[ik] = esdtrack->GetKinkIndex(ik);
	380	}
	381	trackCopy->SetKinkIndexes(indexes);
	382
f007751a	383	// Fill the hidden information with the simulated data
f007751a	384	TParticle *tPart = tStack->Particle(TMath::Abs(esdtrack->GetLabel()));
0b3bd1ac	385	// AliAODParticle* tParticle= new AliAODParticle(*tPart,i);
f007751a	386	AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
0b3bd1ac	387	tInfo->SetPDGPid(tPart->GetPdgCode());
	388	tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
	389	Double_t mass2 = (tPart->Energy() *tPart->Energy() -
	390	tPart->Px()*tPart->Px() -
	391	tPart->Py()*tPart->Py() -
	392	tPart->Pz()*tPart->Pz());
f007751a	393	if (mass2>0.0)
	394	tInfo->SetMass(TMath::Sqrt(mass2));
	395	else
	396	tInfo->SetMass(0.0);
0b3bd1ac	397	tInfo->SetEmissionPoint(tPart->Vx()1e13, tPart->Vy()1e13, tPart->Vz()1e13, tPart->T()1e13*300000);
f007751a	398	trackCopy->SetHiddenInfo(tInfo);
	399
	400	//decision if we want this track
	401	//if we using diffrent labels we want that this label was use for first time
	402	//if we use hidden info we want to have match between sim data and ESD
	403	if (tGoodMomentum==true)
	404	{
	405	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	406	realnofTracks++;//real number of tracks
	407	}
	408	else
	409	{
	410	delete trackCopy;
	411	}
	412
	413	}
	414
	415	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
	416	fCurEvent++;
cc5faabc	417	fCurRLEvent++;
f007751a	418	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
	419	if (fCurEvent== fNumberofEvent)//if end of current file close all
	420	{
	421	fTree->Reset();
	422	delete fTree;
f007751a	423	}
	424	return hbtEvent;
	425	}
0b3bd1ac	426	//____________________________________________________________________
	427	Float_t AliFemtoEventReaderESDKine::GetSigmaToVertex(const AliESDtrack* esdTrack)
	428	{
	429	// Calculates the number of sigma to the vertex.
	430
	431	Float_t b[2];
	432	Float_t bRes[2];
	433	Float_t bCov[3];
	434	esdTrack->GetImpactParameters(b,bCov);
	435	// if (bCov[0]<=0 \|\| bCov[2]<=0) {
	436	// AliDebug(1, "Estimated b resolution lower or equal zero!");
	437	// bCov[0]=0; bCov[2]=0;
	438	// }
	439	bRes[0] = TMath::Sqrt(bCov[0]);
	440	bRes[1] = TMath::Sqrt(bCov[2]);
	441
	442	// -----------------------------------
	443	// How to get to a n-sigma cut?
	444	//
	445	// The accumulated statistics from 0 to d is
	446	//
	447	// -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
	448	// -> 1 - Exp(-d*2) for a 2-dim gauss (dd = dxdx + dydy != n_sigma)
	449	//
	450	// It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)ErfInv(1 - Exp((-x*2)/2)
	451	// Can this be expressed in a different way?
	452
	453	if (bRes[0] == 0 \|\| bRes[1] ==0)
	454	return -1;
	455
	456	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	457
	458	// stupid rounding problem screws up everything:
	459	// if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
	460	if (TMath::Exp(-d * d / 2) < 1e-10)
	461	return 1000;
	462
	463	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	464	return d;
	465	}