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