[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoEventReaderESD.cxx

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoEventReaderESD - 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.2.2.2  2007/10/04 13:10:52  akisiel
 *Add Kink index storageAliFemtoEventReaderESD.cxx AliFemtoTrack.cxx AliFemtoTrack.h
 *
 *Revision 1.2.2.1  2007/09/30 11:38:59  akisiel
 *Adapt the readers to the new AliESDEvent structure
 *
 *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 "AliFemtoEventReaderESD.h"

#include "TFile.h"
#include "TTree.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"

//#include "TSystem.h"

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

#include "SystemOfUnits.h"

#include "AliFemtoEvent.h"
#include "AliFemtoModelHiddenInfo.h"

ClassImp(AliFemtoEventReaderESD)

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

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

AliFemtoEventReaderESD::AliFemtoEventReaderESD(const AliFemtoEventReaderESD &aReader) :
  AliFemtoEventReader(aReader),
  fInputFile(" "),
  fFileName(" "),
  fConstrained(true),
  fReadInner(false),
  fNumberofEvent(0),
  fCurEvent(0),
  fTree(0x0),
  fEsdFile(0x0),
  fEvent(0x0)
{
  // copy constructor
  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fReadInner = aReader.fReadInner;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  //  fTree = aReader.fTree->CloneTree();
  //  fEvent = new AliESD(*aReader.fEvent);
  fEvent = new AliESDEvent();
  fEsdFile = new TFile(aReader.fEsdFile->GetName());
}
//__________________
//Destructor
AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
{
  // destructor
  //delete fListOfFiles;
  delete fTree;
  delete fEvent;
  delete fEsdFile;
}

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

  fInputFile = aReader.fInputFile;
  fFileName  = aReader.fFileName;
  fConstrained = aReader.fConstrained;
  fReadInner = aReader.fReadInner;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  if (fTree) delete fTree;
  //  fTree = aReader.fTree->CloneTree();
  if (fEvent) delete fEvent;
  fEvent = new AliESDEvent();
  if (fEsdFile) delete fEsdFile;
  fEsdFile = new TFile(aReader.fEsdFile->GetName());

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

//__________________
void AliFemtoEventReaderESD::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 AliFemtoEventReaderESD::SetConstrained(const bool constrained)
{
  fConstrained=constrained;
}

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

void AliFemtoEventReaderESD::SetReadTPCInner(const bool readinner)
{
  fReadInner=readinner;
}

bool AliFemtoEventReaderESD::GetReadTPCInner() const
{
  return fReadInner;
}

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

  if (fCurEvent==fNumberofEvent)//open next file  
    {
      if(fNumberofEvent==0)	
	{
	  //	  delete fEvent;//added 1.04.2007
	  fEvent=new AliESDEvent();
	  //	  delete fTree;
	  //fTree=0;
	  //	  delete fEsdFile;
		
	  //ESD data
	  //	  fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	  //	  fTree = (TTree*) fEsdFile->Get("esdTree");			
	  //	  fTree->SetBranchAddress("ESD", &fEvent);			
	  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);

	  fNumberofEvent=fTree->GetEntries();
	  cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
	  fCurEvent=0;
	  //sim data
	}
      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
  cout << "Read event " << fEvent << " from file " << fTree << endl;
  //  vector<int> tLabelTable;//to check labels
	
  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
  cout << "Event has " << nofTracks << " tracks " << endl;

  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 (fReadInner == true) {
	
	if (esdtrack->GetTPCInnerParam()) {
	  AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
	  param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
	  param->GetPxPyPz(pxyz);//reading noconstarined momentum
	  delete param;

	  AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
	  tInfo->SetPDGPid(211);
	  tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
	  tInfo->SetMass(0.13957);
	  trackCopy->SetHiddenInfo(tInfo);
	}
      }
      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 origin(fV1[0],fV1[1],fV1[2]);
      //setting helix I do not if it is ok
      AliFmPhysicalHelixD helix(ktP,origin,(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());

      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);
      //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++;	
  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
//   if (fCurEvent== fNumberofEvent)//if end of current file close all
//     {   
//       fTree->Reset(); 
//       delete fTree;
//       fEsdFile->Close();
//     }
  return hbtEvent; 
}
Commit	Line	Data
76ce4b5b	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoEventReaderESD - 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$
	14	*Revision 1.2.2.2 2007/10/04 13:10:52 akisiel
	15	*Add Kink index storageAliFemtoEventReaderESD.cxx AliFemtoTrack.cxx AliFemtoTrack.h
	16	*
	17	*Revision 1.2.2.1 2007/09/30 11:38:59 akisiel
	18	*Adapt the readers to the new AliESDEvent structure
	19	*
	20	*Revision 1.2 2007/05/22 09:01:42 akisiel
	21	*Add the possibiloity to save cut settings in the ROOT file
	22	*
	23	*Revision 1.1 2007/05/16 10:22:11 akisiel
	24	*Making the directory structure of AliFemto flat. All files go into one common directory
	25	*
	26	*Revision 1.5 2007/05/03 09:45:20 akisiel
	27	*Fixing Effective C++ warnings
	28	*
	29	*Revision 1.4 2007/04/27 07:28:34 akisiel
	30	*Remove event number reading due to interface changes
	31	*
	32	*Revision 1.3 2007/04/27 07:25:16 akisiel
	33	*Make revisions needed for compilation from the main AliRoot tree
	34	*
	35	*Revision 1.1.1.1 2007/04/25 15:38:41 panos
	36	*Importing the HBT code dir
	37	*
	38	*/
	39
	40	#include "AliFemtoEventReaderESD.h"
	41
	42	#include "TFile.h"
	43	#include "TTree.h"
	44	#include "AliESDEvent.h"
	45	#include "AliESDtrack.h"
	46	#include "AliESDVertex.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	#include "AliFemtoModelHiddenInfo.h"
	57
	58	ClassImp(AliFemtoEventReaderESD)
	59
	60	#if !(ST_NO_NAMESPACES)
	61	using namespace units;
	62	#endif
	63
	64	using namespace std;
65	//____________________________
66	//constructor with 0 parameters , look at default settings
67	AliFemtoEventReaderESD::AliFemtoEventReaderESD():
68	fInputFile(" "),
69	fFileName(" "),
70	fConstrained(true),
71	fReadInner(false),
72	fNumberofEvent(0),
73	fCurEvent(0),
74	fTree(0x0),
75	fEsdFile(0x0),
76	fEvent(0x0)
77	{
78	// default constructor
79	}
80
81	AliFemtoEventReaderESD::AliFemtoEventReaderESD(const AliFemtoEventReaderESD &aReader) :
82	AliFemtoEventReader(aReader),
83	fInputFile(" "),
84	fFileName(" "),
85	fConstrained(true),
86	fReadInner(false),
87	fNumberofEvent(0),
88	fCurEvent(0),
89	fTree(0x0),
90	fEsdFile(0x0),
91	fEvent(0x0)
92	{
93	// copy constructor
94	fInputFile = aReader.fInputFile;
95	fFileName = aReader.fFileName;
96	fConstrained = aReader.fConstrained;
97	fReadInner = aReader.fReadInner;
98	fNumberofEvent = aReader.fNumberofEvent;
99	fCurEvent = aReader.fCurEvent;
100	// fTree = aReader.fTree->CloneTree();
101	// fEvent = new AliESD(*aReader.fEvent);
102	fEvent = new AliESDEvent();
103	fEsdFile = new TFile(aReader.fEsdFile->GetName());
104	}
105	//__________________
106	//Destructor
107	AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
108	{
109	// destructor
110	//delete fListOfFiles;
111	delete fTree;
112	delete fEvent;
113	delete fEsdFile;
114	}
115
116	//__________________
117	AliFemtoEventReaderESD& AliFemtoEventReaderESD::operator=(const AliFemtoEventReaderESD& aReader)
118	{
119	// assignment operator
120	if (this == &aReader)
121	return *this;
122
123	fInputFile = aReader.fInputFile;
124	fFileName = aReader.fFileName;
125	fConstrained = aReader.fConstrained;
126	fReadInner = aReader.fReadInner;
127	fNumberofEvent = aReader.fNumberofEvent;
128	fCurEvent = aReader.fCurEvent;
129	if (fTree) delete fTree;
130	// fTree = aReader.fTree->CloneTree();
131	if (fEvent) delete fEvent;
132	fEvent = new AliESDEvent();
133	if (fEsdFile) delete fEsdFile;
134	fEsdFile = new TFile(aReader.fEsdFile->GetName());
135
136	return *this;
137	}
138	//__________________
139	AliFemtoString AliFemtoEventReaderESD::Report()
140	{
141	// create reader report
142	AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
143	return temp;
144	}
145
146	//__________________
147	void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
148	{
149	//setting the name of file where names of ESD file are written
150	//it takes only this files which have good trees
151	char buffer[256];
152	fInputFile=string(inputFile);
153	cout<<"Input File set on "<<fInputFile<<endl;
154	ifstream infile(inputFile);
155
156	fTree = new TChain("esdTree");
157
158	if(infile.good()==true)
159	{
160	//checking if all give files have good tree inside
161	while (infile.eof()==false)
162	{
163	infile.getline(buffer,256);
164	//ifstream test_file(buffer);
165	TFile *esdFile=TFile::Open(buffer,"READ");
166	if (esdFile!=0x0)
167	{
168	TTree* tree = (TTree*) esdFile->Get("esdTree");
169	if (tree!=0x0)
170	{
171	cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
172	fTree->AddFile(buffer);
173	delete tree;
174	}
175	esdFile->Close();
176	}
177	delete esdFile;
178	}
179	}
180	}
181
182	void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
183	{
184	fConstrained=constrained;
185	}
186
187	bool AliFemtoEventReaderESD::GetConstrained() const
188	{
189	return fConstrained;
190	}
191
192	void AliFemtoEventReaderESD::SetReadTPCInner(const bool readinner)
193	{
194	fReadInner=readinner;
195	}
196
197	bool AliFemtoEventReaderESD::GetReadTPCInner() const
198	{
199	return fReadInner;
200	}
201
202	AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
203	{
204	// read in a next hbt event from the chain
205	// convert it to AliFemtoEvent and return
206	// for further analysis
207	AliFemtoEvent *hbtEvent = 0;
208
209	if (fCurEvent==fNumberofEvent)//open next file
210	{
211	if(fNumberofEvent==0)
212	{
213	// delete fEvent;//added 1.04.2007
214	fEvent=new AliESDEvent();
215	// delete fTree;
216	//fTree=0;
217	// delete fEsdFile;
218
219	//ESD data
220	// fEsdFile=TFile::Open(fFileName.c_str(),"READ");
221	// fTree = (TTree*) fEsdFile->Get("esdTree");
222	// fTree->SetBranchAddress("ESD", &fEvent);
223	fTree->SetBranchStatus("MuonTracks*",0);
224	fTree->SetBranchStatus("PmdTracks*",0);
225	fTree->SetBranchStatus("TrdTracks*",0);
226	fTree->SetBranchStatus("V0s*",0);
227	fTree->SetBranchStatus("Cascades*",0);
228	fTree->SetBranchStatus("Kinks*",0);
229	fTree->SetBranchStatus("CaloClusters*",0);
230	fTree->SetBranchStatus("AliRawDataErrorLogs*",0);
231	fTree->SetBranchStatus("ESDfriend*",0);
232	fEvent->ReadFromTree(fTree);
233
234	fNumberofEvent=fTree->GetEntries();
235	cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
236	fCurEvent=0;
237	//sim data
238	}
239	else //no more data to read
240	{
241	cout<<"no more files "<<hbtEvent<<endl;
242	fReaderStatus=1;
243	return hbtEvent;
244	}
245	}
246	cout<<"starting to read event "<<fCurEvent<<endl;
247	fTree->GetEvent(fCurEvent);//getting next event
248	cout << "Read event " << fEvent << " from file " << fTree << endl;
249	// vector<int> tLabelTable;//to check labels
250
251	hbtEvent = new AliFemtoEvent;
252	//setting basic things
253	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
254	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
255	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
256	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
257	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
258	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
259	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
260	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
261	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
262	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
263	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
264	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
265
266	//Vertex
267	double fV1[3];
268	fEvent->GetVertex()->GetXYZ(fV1);
269
270	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
271	hbtEvent->SetPrimVertPos(vertex);
272
273	//starting to reading tracks
274	int nofTracks=0; //number of reconstructed tracks in event
275	nofTracks=fEvent->GetNumberOfTracks();
276	int realnofTracks=0;//number of track which we use ina analysis
277	cout << "Event has " << nofTracks << " tracks " << endl;
278
279	for (int i=0;i<nofTracks;i++)
280	{
281	bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
282
283	AliFemtoTrack* trackCopy = new AliFemtoTrack();
284	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
285	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
286
287	trackCopy->SetCharge((short)esdtrack->GetSign());
288
289	//in aliroot we have AliPID
290	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
291	//we use only 5 first
292	double esdpid[5];
293	esdtrack->GetESDpid(esdpid);
294	trackCopy->SetPidProbElectron(esdpid[0]);
295	trackCopy->SetPidProbMuon(esdpid[1]);
296	trackCopy->SetPidProbPion(esdpid[2]);
297	trackCopy->SetPidProbKaon(esdpid[3]);
298	trackCopy->SetPidProbProton(esdpid[4]);
299
300	double pxyz[3];
301	if (fReadInner == true) {
302
303	if (esdtrack->GetTPCInnerParam()) {
304	AliExternalTrackParam param = new AliExternalTrackParam(esdtrack->GetTPCInnerParam());
305	param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
306	param->GetPxPyPz(pxyz);//reading noconstarined momentum
307	delete param;
308
309	AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
310	tInfo->SetPDGPid(211);
311	tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
312	tInfo->SetMass(0.13957);
313	trackCopy->SetHiddenInfo(tInfo);
314	}
315	}
316	if (fConstrained==true)
317	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
318	else
319	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
320	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
321	trackCopy->SetP(v);//setting momentum
322	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
323	const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
324	if (ktP.Mag() == 0) {
325	delete trackCopy;
326	continue;
327	}
328	const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
329	//setting helix I do not if it is ok
330	AliFmPhysicalHelixD helix(ktP,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
331	trackCopy->SetHelix(helix);
332
333	trackCopy->SetTrackId(esdtrack->GetID());
334	trackCopy->SetFlags(esdtrack->GetStatus());
335	trackCopy->SetLabel(esdtrack->GetLabel());
336
337	//some stuff which could be useful
338	float impact[2];
339	float covimpact[3];
340	esdtrack->GetImpactParameters(impact,covimpact);
341	trackCopy->SetImpactD(impact[0]);
342	trackCopy->SetImpactZ(impact[1]);
343	trackCopy->SetCdd(covimpact[0]);
344	trackCopy->SetCdz(covimpact[1]);
345	trackCopy->SetCzz(covimpact[2]);
346	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
347	trackCopy->SetITSncls(esdtrack->GetNcls(0));
348	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
349	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
350	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
351	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
352	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
353
354	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
355	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
356
357	double pvrt[3];
358	fEvent->GetPrimaryVertex()->GetXYZ(pvrt);
359
360	double xtpc[3];
361	esdtrack->GetInnerXYZ(xtpc);
362	xtpc[2] -= pvrt[2];
363	trackCopy->SetNominalTPCEntrancePoint(xtpc);
364
365	esdtrack->GetOuterXYZ(xtpc);
366	xtpc[2] -= pvrt[2];
367	trackCopy->SetNominalTPCExitPoint(xtpc);
368
369	int indexes[3];
370	for (int ik=0; ik<3; ik++) {
371	indexes[ik] = esdtrack->GetKinkIndex(ik);
372	}
373	trackCopy->SetKinkIndexes(indexes);
374	//decision if we want this track
375	//if we using diffrent labels we want that this label was use for first time
376	//if we use hidden info we want to have match between sim data and ESD
377	if (tGoodMomentum==true)
378	{
379	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
380	realnofTracks++;//real number of tracks
381	}
382	else
383	{
384	delete trackCopy;
385	}
386
387	}
388
389	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
390	fCurEvent++;
391	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
392	// if (fCurEvent== fNumberofEvent)//if end of current file close all
393	// {
394	// fTree->Reset();
395	// delete fTree;
396	// fEsdFile->Close();
397	// }
398	return hbtEvent;
399	}
400
401
402
403
404
405
406
407
408