[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"

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->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
      trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());

      // 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(tParticle->GetMostProbable());
      tInfo->SetTrueMomentum(tParticle->Px(), tParticle->Py(), tParticle->Pz());
      Double_t mass2 = (tParticle->E()*tParticle->E() -
			 tParticle->Px()*tParticle->Px() -
			 tParticle->Py()*tParticle->Py() -
			 tParticle->Pz()*tParticle->Pz());
      if (mass2>0.0)
	tInfo->SetMass(TMath::Sqrt(mass2));
      else 
	tInfo->SetMass(0.0);
      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; 
}
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"
	45	#include "AliAODParticle.h"
	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
	57	ClassImp(AliFemtoEventReaderESDKine)
	58
	59	#if !(ST_NO_NAMESPACES)
	60	using namespace units;
	61	#endif
	62
	63	using namespace std;
	64	//____________________________
	65	//constructor with 0 parameters , look at default settings
	66	AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine():
	67	fInputFile(" "),
	68	fFileName(" "),
	69	fConstrained(true),
	70	fNumberofEvent(0),
	71	fCurEvent(0),
cc5faabc	72	fCurRLEvent(0),
f007751a	73	fTree(0x0),
f007751a	74	fEvent(0x0),
cc5faabc	75	fRunLoader(0x0)
f007751a	76	{
f007751a	77	}
	78
	79	AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
	80	fInputFile(" "),
	81	fFileName(" "),
	82	fConstrained(true),
	83	fNumberofEvent(0),
	84	fCurEvent(0),
cc5faabc	85	fCurRLEvent(0),
f007751a	86	fTree(0x0),
f007751a	87	fEvent(0x0),
cc5faabc	88	fRunLoader(0x0)
f007751a	89	{
	90	// copy constructor
	91	fInputFile = aReader.fInputFile;
	92	fFileName = aReader.fFileName;
	93	fConstrained = aReader.fConstrained;
	94	fNumberofEvent = aReader.fNumberofEvent;
	95	fCurEvent = aReader.fCurEvent;
cc5faabc	96	fEvent = new AliESDEvent();
f007751a	97	}
	98	//__________________
	99	//Destructor
	100	AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
	101	{
	102	// destructor
	103	//delete fListOfFiles;
	104	delete fTree;
	105	delete fEvent;
f007751a	106	if (fRunLoader) delete fRunLoader;
f007751a	107	}
	108
	109	//__________________
	110	AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemtoEventReaderESDKine& aReader)
	111	{
	112	// assignment operator
	113	if (this == &aReader)
	114	return *this;
	115
	116	fInputFile = aReader.fInputFile;
	117	fFileName = aReader.fFileName;
	118	fConstrained = aReader.fConstrained;
	119	fNumberofEvent = aReader.fNumberofEvent;
	120	fCurEvent = aReader.fCurEvent;
cc5faabc	121	fCurRLEvent = aReader.fCurRLEvent;
f007751a	122	if (fTree) delete fTree;
cc5faabc	123	// fTree = aReader.fTree->CloneTree();
f007751a	124	if (fEvent) delete fEvent;
cc5faabc	125	fEvent = new AliESDEvent();
f007751a	126	if (fRunLoader) delete fRunLoader;
	127	fRunLoader = new AliRunLoader(*aReader.fRunLoader);
	128
f007751a	129	return *this;
	130	}
	131	//__________________
cc5faabc	132	AliFemtoString AliFemtoEventReaderESDKine::Report() const
f007751a	133	{
	134	// create reader report
	135	AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
	136	return temp;
	137	}
	138
	139	//__________________
	140	void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
	141	{
	142	//setting the name of file where names of ESD file are written
	143	//it takes only this files which have good trees
	144	char buffer[256];
	145	fInputFile=string(inputFile);
	146	cout<<"Input File set on "<<fInputFile<<endl;
	147	ifstream infile(inputFile);
cc5faabc	148
	149	fTree = new TChain("esdTree");
	150
f007751a	151	if(infile.good()==true)
	152	{
	153	//checking if all give files have good tree inside
	154	while (infile.eof()==false)
	155	{
	156	infile.getline(buffer,256);
	157	//ifstream test_file(buffer);
	158	TFile *esdFile=TFile::Open(buffer,"READ");
	159	if (esdFile!=0x0)
	160	{
	161	TTree* tree = (TTree*) esdFile->Get("esdTree");
	162	if (tree!=0x0)
	163	{
	164	cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
cc5faabc	165	fTree->AddFile(buffer);
f007751a	166	delete tree;
	167	}
	168	esdFile->Close();
	169	}
	170	delete esdFile;
	171	}
	172	}
	173	}
	174
f007751a	175	void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
	176	{
	177	fConstrained=constrained;
	178	}
	179
	180	bool AliFemtoEventReaderESDKine::GetConstrained() const
	181	{
	182	return fConstrained;
	183	}
	184
	185	AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
	186	{
	187	// read in a next hbt event from the chain
	188	// convert it to AliFemtoEvent and return
	189	// for further analysis
	190	AliFemtoEvent *hbtEvent = 0;
f007751a	191	TString tGAliceFilename;
	192
	193	if (fCurEvent==fNumberofEvent)//open next file
	194	{
cc5faabc	195	if (fNumberofEvent == 0) {
cc5faabc	196	fEvent=new AliESDEvent();
f007751a	197
f007751a	198	//ESD data
cc5faabc	199	// fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	200	// fTree = (TTree*) fEsdFile->Get("esdTree");
	201
	202	fTree->SetBranchStatus("MuonTracks*",0);
	203	fTree->SetBranchStatus("PmdTracks*",0);
	204	fTree->SetBranchStatus("TrdTracks*",0);
	205	fTree->SetBranchStatus("V0s*",0);
	206	fTree->SetBranchStatus("Cascades*",0);
	207	fTree->SetBranchStatus("Kinks*",0);
	208	fTree->SetBranchStatus("CaloClusters*",0);
	209	fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
	210	fTree->SetBranchStatus("ESDfriend*",0);
	211	fEvent->ReadFromTree(fTree);
f007751a	212
	213	// chain->SetBranchStatus("*",0);
	214	// chain->SetBranchStatus("fUniqueID",1);
	215	// chain->SetBranchStatus("fTracks",1);
	216	// chain->SetBranchStatus("fTracks.*",1);
	217	// chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
cc5faabc	218	// fTree->SetBranchStatus("fTracks.fCalibContainer",0);
	219
	220
	221	fNumberofEvent=fTree->GetEntries();
	222
	223	if (fNumberofEvent == 0) {
	224	cout<<"no event in input "<<endl;
	225	fReaderStatus=1;
	226	return hbtEvent;
f007751a	227	}
cc5faabc	228
	229	cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
	230	fCurEvent=0;
	231	// simulation data reading setup
	232
	233	}
f007751a	234	else //no more data to read
	235	{
	236	cout<<"no more files "<<hbtEvent<<endl;
	237	fReaderStatus=1;
	238	return hbtEvent;
	239	}
	240	}
	241	cout<<"starting to read event "<<fCurEvent<<endl;
	242	fTree->GetEvent(fCurEvent);//getting next event
f007751a	243	// vector<int> tLabelTable;//to check labels
f007751a	244
cc5faabc	245	cout << "fFileName is " << fFileName.Data() << endl;
	246	cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl;
	247	if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) {
	248	fFileName = fTree->GetCurrentFile()->GetName();
	249	tGAliceFilename = fFileName;
	250	tGAliceFilename.ReplaceAll("AliESDs","galice");
	251	cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl;
	252	if (fRunLoader) delete fRunLoader;
	253	fRunLoader = AliRunLoader::Open(tGAliceFilename.Data());
	254	if (fRunLoader==0x0)
	255	{
	256	cout << "No Kine tree in file " << tGAliceFilename.Data() << endl;
	257	exit(0);
	258	}
	259	if(fRunLoader->LoadHeader())
	260	{
	261	cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl;
	262	exit(0);
	263	}
	264	fRunLoader->LoadKinematics();
	265	fCurRLEvent = 0;
	266	}
	267
	268	fRunLoader->GetEvent(fCurRLEvent);
f007751a	269	AliStack* tStack = 0x0;
	270	tStack = fRunLoader->Stack();
	271
	272	hbtEvent = new AliFemtoEvent;
	273	//setting basic things
	274	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
	275	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	276	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	277	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	278	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	279	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	280	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	281	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	282	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	283	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	284	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	285	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	286
	287	//Vertex
	288	double fV1[3];
	289	fEvent->GetVertex()->GetXYZ(fV1);
	290
	291	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	292	hbtEvent->SetPrimVertPos(vertex);
	293
	294	//starting to reading tracks
	295	int nofTracks=0; //number of reconstructed tracks in event
	296	nofTracks=fEvent->GetNumberOfTracks();
	297	int realnofTracks=0;//number of track which we use ina analysis
	298
f007751a	299	for (int i=0;i<nofTracks;i++)
	300	{
	301	bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
	302
	303	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	304	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	305	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
	306
	307	trackCopy->SetCharge((short)esdtrack->GetSign());
	308
	309	//in aliroot we have AliPID
	310	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	311	//we use only 5 first
	312	double esdpid[5];
	313	esdtrack->GetESDpid(esdpid);
	314	trackCopy->SetPidProbElectron(esdpid[0]);
	315	trackCopy->SetPidProbMuon(esdpid[1]);
	316	trackCopy->SetPidProbPion(esdpid[2]);
	317	trackCopy->SetPidProbKaon(esdpid[3]);
	318	trackCopy->SetPidProbProton(esdpid[4]);
	319
	320	double pxyz[3];
	321	if (fConstrained==true)
	322	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
	323	else
	324	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
	325	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	326	trackCopy->SetP(v);//setting momentum
	327	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
	328	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
	329	if (ktP.mag() == 0) {
	330	delete trackCopy;
	331	continue;
	332	}
cc5faabc	333	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
f007751a	334	//setting helix I do not if it is ok
cc5faabc	335	AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
f007751a	336	trackCopy->SetHelix(helix);
	337
	338	trackCopy->SetTrackId(esdtrack->GetID());
	339	trackCopy->SetFlags(esdtrack->GetStatus());
cc5faabc	340	trackCopy->SetLabel(esdtrack->GetLabel());
f007751a	341
	342	//some stuff which could be useful
	343	float impact[2];
	344	float covimpact[3];
	345	esdtrack->GetImpactParameters(impact,covimpact);
	346	trackCopy->SetImpactD(impact[0]);
	347	trackCopy->SetImpactZ(impact[1]);
	348	trackCopy->SetCdd(covimpact[0]);
	349	trackCopy->SetCdz(covimpact[1]);
	350	trackCopy->SetCzz(covimpact[2]);
	351	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
	352	trackCopy->SetITSncls(esdtrack->GetNcls(0));
	353	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
	354	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
	355	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
	356	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
	357	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
	358
cc5faabc	359	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
	360	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
	361
f007751a	362	// Fill the hidden information with the simulated data
	363	TParticle *tPart = tStack->Particle(TMath::Abs(esdtrack->GetLabel()));
	364	AliAODParticle* tParticle= new AliAODParticle(*tPart,i);
	365	AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
	366	tInfo->SetPDGPid(tParticle->GetMostProbable());
	367	tInfo->SetTrueMomentum(tParticle->Px(), tParticle->Py(), tParticle->Pz());
	368	Double_t mass2 = (tParticle->E()*tParticle->E() -
	369	tParticle->Px()*tParticle->Px() -
	370	tParticle->Py()*tParticle->Py() -
	371	tParticle->Pz()*tParticle->Pz());
	372	if (mass2>0.0)
	373	tInfo->SetMass(TMath::Sqrt(mass2));
	374	else
	375	tInfo->SetMass(0.0);
	376	trackCopy->SetHiddenInfo(tInfo);
	377
	378	//decision if we want this track
	379	//if we using diffrent labels we want that this label was use for first time
	380	//if we use hidden info we want to have match between sim data and ESD
	381	if (tGoodMomentum==true)
	382	{
	383	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	384	realnofTracks++;//real number of tracks
	385	}
	386	else
	387	{
	388	delete trackCopy;
	389	}
	390
	391	}
	392
	393	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
	394	fCurEvent++;
cc5faabc	395	fCurRLEvent++;
f007751a	396	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
	397	if (fCurEvent== fNumberofEvent)//if end of current file close all
	398	{
	399	fTree->Reset();
	400	delete fTree;
f007751a	401	}
	402	return hbtEvent;
	403	}