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

/*
 *$Id$
 *$Log$
 *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
 *
 *Revision 1.5  2007-04-03 16:00:08  mchojnacki
 *Changes to iprove memory managing
 *
 *Revision 1.4  2007/03/13 15:30:03  mchojnacki
 *adding reader for simulated data
 *
 *Revision 1.3  2007/03/08 14:58:03  mchojnacki
 *adding some alice stuff
 *
 *Revision 1.2  2007/03/07 13:36:16  mchojnacki
 *Add some comments
 *
 *Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
 *First version on CVS
 *
 */

#include "AliFemtoEventReaderESD.h"

#include "TFile.h"
#include "TTree.h"
#include "AliESD.h"
#include "AliESDtrack.h"

//#include "TSystem.h"

#include "Infrastructure/AliFmPhysicalHelixD.h"
#include "Infrastructure/AliFmThreeVectorF.h"

#include "Base/SystemOfUnits.h"

#include "Infrastructure/AliFemtoEvent.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),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fTree(0x0),
  fEvent(0x0),
  fEsdFile(0x0),
  fEventFriend(0)
{
  fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad] = new list<Int_t>();
  }
  fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad] = new list<Int_t>();
  }
}

//__________________
//Destructor
AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
{
  //delete fListOfFiles;
  delete fTree;
  delete fEvent;
  delete fEsdFile;

  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad]->clear();
    delete fClusterPerPadrow[tPad];
  }
  delete [] fClusterPerPadrow;
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad]->clear();
    delete fSharedList[tPad];
  }
  delete [] fSharedList;
}

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

//__________________
//setting the name of file where names of ESD file are written 
//it takes only this files which have good trees
void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
{
  char buffer[256];
  fInputFile=string(inputFile);
  cout<<"Input File set on "<<fInputFile<<endl;
  ifstream infile(inputFile);
  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;
		  fListOfFiles.push_back(string(buffer));
		  delete tree;
		}
	      esdFile->Close();	
	    }
	  delete esdFile;
	}
    }
}

//setting the next file to read	
bool AliFemtoEventReaderESD::GetNextFile()
{ 	
  if (fCurFile>=fListOfFiles.size())
    return false;
  fFileName=fListOfFiles.at(fCurFile);	
  cout<<"FileName set on "<<fFileName<<" "<<fCurFile<<endl;

  fCurFile++;
  return true;
}
void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
{
  fConstrained=constrained;
}

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

AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
{
  AliFemtoEvent *hbtEvent = 0;
  string tFriendFileName;

  if (fCurEvent==fNumberofEvent)//open next file  
    {
      cout<<"next file"<<endl;
      if(GetNextFile())	
	{
	  delete fEventFriend;
	  fEventFriend = 0;
	  delete fEvent;//added 1.04.2007
	  fEvent=new AliESD();
	  //	  delete fTree;
	  //fTree=0;
	  delete fEsdFile;
		
	  //ESD data
	  fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	  fTree = (TTree*) fEsdFile->Get("esdTree");			
	  fTree->SetBranchAddress("ESD", &fEvent);			

	  // Attach the friend tree with additional information
 	  tFriendFileName = fFileName;
 	  tFriendFileName.insert(tFriendFileName.find("s.root"),"friend");
 	  cout << "Reading friend " << tFriendFileName.c_str() << endl;;
  	  fTree->AddFriend("esdFriendTree",tFriendFileName.c_str());
  	  fTree->SetBranchAddress("ESDfriend",&fEventFriend);

// 	  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();
	  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
  fEvent->SetESDfriend(fEventFriend);
  vector<int> label_table;//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

  // Clear the shared cluster list
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fClusterPerPadrow[tPad]->clear();
  }
  for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
    fSharedList[tPad]->clear();
  }


  for (int i=0;i<nofTracks;i++) {
    const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track

    list<Int_t>::iterator tClustIter;

    Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
    Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
    for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
      if (tTrackIndices[tNcl] >= 0) {
	tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
	if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
	  fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
	}
	else {
	  fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
	}
      }
    }
      
  }

  for (int i=0;i<nofTracks;i++)
    {
      bool  good_momentum=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)		    
	good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
      else
	good_momentum=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 p(pxyz[0],pxyz[1],pxyz[2]);
      if (p.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(p,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()); 

      // Fill cluster per padrow information
      Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
      Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
      for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	if (tTrackIndices[tNcl] > 0)
	  trackCopy->SetTPCcluster(tNcl, 1);
	else
	  trackCopy->SetTPCcluster(tNcl, 0);
      }
      
      // Fill shared cluster information
      list<Int_t>::iterator tClustIter;

      for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	if (tTrackIndices[tNcl] > 0) {
	  tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
	  if (tClustIter != fSharedList[tNcl]->end()) {
	    trackCopy->SetTPCshared(tNcl, 1);
	    cout << "Event next" <<  endl;
	    cout << "Track: " << i << endl;
	    cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
	  }
	  else {
	    trackCopy->SetTPCshared(tNcl, 0);
	  }
	}
      }

      //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 (good_momentum==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
67427ff7	1	/*
	2	*$Id$
	3	*$Log$
a531c6fc	4	*Revision 1.3 2007/04/27 07:25:16 akisiel
	5	*Make revisions needed for compilation from the main AliRoot tree
	6	*
b2f60a91	7	*Revision 1.1.1.1 2007/04/25 15:38:41 panos
	8	*Importing the HBT code dir
	9	*
67427ff7	10	*Revision 1.5 2007-04-03 16:00:08 mchojnacki
	11	*Changes to iprove memory managing
	12	*
	13	*Revision 1.4 2007/03/13 15:30:03 mchojnacki
	14	*adding reader for simulated data
	15	*
	16	*Revision 1.3 2007/03/08 14:58:03 mchojnacki
	17	*adding some alice stuff
	18	*
	19	*Revision 1.2 2007/03/07 13:36:16 mchojnacki
	20	*Add some comments
	21	*
	22	*Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
	23	*First version on CVS
	24	*
	25	*/
	26
	27	#include "AliFemtoEventReaderESD.h"
	28
	29	#include "TFile.h"
	30	#include "TTree.h"
	31	#include "AliESD.h"
	32	#include "AliESDtrack.h"
	33
	34	//#include "TSystem.h"
	35
b2f60a91	36	#include "Infrastructure/AliFmPhysicalHelixD.h"
b2f60a91	37	#include "Infrastructure/AliFmThreeVectorF.h"
67427ff7	38
b2f60a91	39	#include "Base/SystemOfUnits.h"
67427ff7	40
b2f60a91	41	#include "Infrastructure/AliFemtoEvent.h"
67427ff7	42
	43	ClassImp(AliFemtoEventReaderESD)
	44
	45	#if !(ST_NO_NAMESPACES)
	46	using namespace units;
	47	#endif
	48
	49	using namespace std;
	50	//____________________________
	51	//constructor with 0 parameters , look at default settings
	52	AliFemtoEventReaderESD::AliFemtoEventReaderESD():
	53	fInputFile(" "),
	54	fFileName(" "),
	55	fConstrained(true),
	56	fNumberofEvent(0),
	57	fCurEvent(0),
	58	fCurFile(0),
	59	fTree(0x0),
	60	fEvent(0x0),
	61	fEsdFile(0x0),
	62	fEventFriend(0)
	63	{
	64	fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	65	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	66	fClusterPerPadrow[tPad] = new list<Int_t>();
	67	}
	68	fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	69	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	70	fSharedList[tPad] = new list<Int_t>();
	71	}
	72	}
	73
	74	//__________________
	75	//Destructor
	76	AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
	77	{
	78	//delete fListOfFiles;
	79	delete fTree;
	80	delete fEvent;
	81	delete fEsdFile;
	82
	83	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	84	fClusterPerPadrow[tPad]->clear();
	85	delete fClusterPerPadrow[tPad];
	86	}
	87	delete [] fClusterPerPadrow;
	88	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	89	fSharedList[tPad]->clear();
	90	delete fSharedList[tPad];
	91	}
	92	delete [] fSharedList;
	93	}
	94
	95	//__________________
	96	AliFemtoString AliFemtoEventReaderESD::Report()
	97	{
	98	AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
	99	return temp;
	100	}
	101
	102	//__________________
	103	//setting the name of file where names of ESD file are written
	104	//it takes only this files which have good trees
	105	void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
106	{
107	char buffer[256];
108	fInputFile=string(inputFile);
109	cout<<"Input File set on "<<fInputFile<<endl;
110	ifstream infile(inputFile);
111	if(infile.good()==true)
112	{
113	//checking if all give files have good tree inside
114	while (infile.eof()==false)
115	{
116	infile.getline(buffer,256);
117	//ifstream test_file(buffer);
118	TFile *esdFile=TFile::Open(buffer,"READ");
119	if (esdFile!=0x0)
120	{
121	TTree* tree = (TTree*) esdFile->Get("esdTree");
122	if (tree!=0x0)
123	{
124	cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
125	fListOfFiles.push_back(string(buffer));
126	delete tree;
127	}
128	esdFile->Close();
129	}
130	delete esdFile;
131	}
132	}
133	}
134
135	//setting the next file to read
136	bool AliFemtoEventReaderESD::GetNextFile()
137	{
138	if (fCurFile>=fListOfFiles.size())
139	return false;
140	fFileName=fListOfFiles.at(fCurFile);
141	cout<<"FileName set on "<<fFileName<<" "<<fCurFile<<endl;
142
143	fCurFile++;
144	return true;
145	}
146	void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
147	{
148	fConstrained=constrained;
149	}
150
151	bool AliFemtoEventReaderESD::GetConstrained() const
152	{
153	return fConstrained;
154	}
155
156	AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
157	{
158	AliFemtoEvent *hbtEvent = 0;
159	string tFriendFileName;
160
161	if (fCurEvent==fNumberofEvent)//open next file
162	{
163	cout<<"next file"<<endl;
164	if(GetNextFile())
165	{
166	delete fEventFriend;
167	fEventFriend = 0;
168	delete fEvent;//added 1.04.2007
169	fEvent=new AliESD();
1e1d2996	170	// delete fTree;
67427ff7	171	//fTree=0;
	172	delete fEsdFile;
	173
	174	//ESD data
	175	fEsdFile=TFile::Open(fFileName.c_str(),"READ");
	176	fTree = (TTree*) fEsdFile->Get("esdTree");
	177	fTree->SetBranchAddress("ESD", &fEvent);
	178
	179	// Attach the friend tree with additional information
	180	tFriendFileName = fFileName;
	181	tFriendFileName.insert(tFriendFileName.find("s.root"),"friend");
	182	cout << "Reading friend " << tFriendFileName.c_str() << endl;;
	183	fTree->AddFriend("esdFriendTree",tFriendFileName.c_str());
	184	fTree->SetBranchAddress("ESDfriend",&fEventFriend);
	185
1e1d2996	186	// chain->SetBranchStatus("*",0);
	187	// chain->SetBranchStatus("fUniqueID",1);
	188	// chain->SetBranchStatus("fTracks",1);
	189	// chain->SetBranchStatus("fTracks.*",1);
	190	// chain->SetBranchStatus("fTracks.fTPCindex[160]",1);
	191	fTree->SetBranchStatus("fTracks.fCalibContainer",0);
	192
	193
67427ff7	194	fNumberofEvent=fTree->GetEntries();
	195	cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
	196	fCurEvent=0;
	197	//sim data
	198	}
	199	else //no more data to read
	200	{
	201	cout<<"no more files "<<hbtEvent<<endl;
	202	fReaderStatus=1;
	203	return hbtEvent;
	204	}
	205	}
	206	cout<<"starting to read event "<<fCurEvent<<endl;
	207	fTree->GetEvent(fCurEvent);//getting next event
	208	fEvent->SetESDfriend(fEventFriend);
	209	vector<int> label_table;//to check labels
	210
	211	hbtEvent = new AliFemtoEvent;
	212	//setting basic things
a531c6fc	213	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7	214	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	215	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	216	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	217	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	218	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	219	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	220	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	221	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	222	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	223	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	224	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	225
	226	//Vertex
	227	double fV1[3];
	228	fEvent->GetVertex()->GetXYZ(fV1);
	229
	230	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	231	hbtEvent->SetPrimVertPos(vertex);
	232
	233	//starting to reading tracks
	234	int nofTracks=0; //number of reconstructed tracks in event
	235	nofTracks=fEvent->GetNumberOfTracks();
	236	int realnofTracks=0;//number of track which we use ina analysis
	237
	238	// Clear the shared cluster list
	239	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	240	fClusterPerPadrow[tPad]->clear();
	241	}
	242	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	243	fSharedList[tPad]->clear();
	244	}
	245
	246
	247	for (int i=0;i<nofTracks;i++) {
	248	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	249
	250	list<Int_t>::iterator tClustIter;
	251
	252	Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
	253	Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
	254	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	255	if (tTrackIndices[tNcl] >= 0) {
	256	tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
	257	if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
	258	fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
	259	}
	260	else {
	261	fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
	262	}
	263	}
	264	}
	265
	266	}
	267
	268	for (int i=0;i<nofTracks;i++)
	269	{
	270	bool good_momentum=true; //flaga to chcek if we can read momentum of this track
	271
	272	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	273	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	274	const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
	275
	276	trackCopy->SetCharge((short)esdtrack->GetSign());
	277
278	//in aliroot we have AliPID
279	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
280	//we use only 5 first
281	double esdpid[5];
282	esdtrack->GetESDpid(esdpid);
283	trackCopy->SetPidProbElectron(esdpid[0]);
284	trackCopy->SetPidProbMuon(esdpid[1]);
285	trackCopy->SetPidProbPion(esdpid[2]);
286	trackCopy->SetPidProbKaon(esdpid[3]);
287	trackCopy->SetPidProbProton(esdpid[4]);
288
289	double pxyz[3];
290	if (fConstrained==true)
291	good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
292	else
293	good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
294	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
295	trackCopy->SetP(v);//setting momentum
296	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
297	const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
298	if (p.mag() == 0) {
299	delete trackCopy;
300	continue;
301	}
302	const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
303	//setting helix I do not if it is ok
304	AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
305	trackCopy->SetHelix(helix);
306
307	trackCopy->SetTrackId(esdtrack->GetID());
308	trackCopy->SetFlags(esdtrack->GetStatus());
309	//trackCopy->SetLabel(esdtrack->GetLabel());
310
311	//some stuff which could be useful
312	float impact[2];
313	float covimpact[3];
314	esdtrack->GetImpactParameters(impact,covimpact);
315	trackCopy->SetImpactD(impact[0]);
316	trackCopy->SetImpactZ(impact[1]);
317	trackCopy->SetCdd(covimpact[0]);
318	trackCopy->SetCdz(covimpact[1]);
319	trackCopy->SetCzz(covimpact[2]);
320	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
321	trackCopy->SetITSncls(esdtrack->GetNcls(0));
322	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
323	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
324	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
325	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
326	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
327
328	// Fill cluster per padrow information
329	Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
330	Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
331	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
332	if (tTrackIndices[tNcl] > 0)
333	trackCopy->SetTPCcluster(tNcl, 1);
334	else
335	trackCopy->SetTPCcluster(tNcl, 0);
336	}
337
338	// Fill shared cluster information
339	list<Int_t>::iterator tClustIter;
340
341	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
342	if (tTrackIndices[tNcl] > 0) {
343	tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
344	if (tClustIter != fSharedList[tNcl]->end()) {
345	trackCopy->SetTPCshared(tNcl, 1);
346	cout << "Event next" << endl;
347	cout << "Track: " << i << endl;
348	cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
349	}
350	else {
351	trackCopy->SetTPCshared(tNcl, 0);
352	}
353	}
354	}
355
356	//decision if we want this track
357	//if we using diffrent labels we want that this label was use for first time
358	//if we use hidden info we want to have match between sim data and ESD
359	if (good_momentum==true)
360	{
361	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
362	realnofTracks++;//real number of tracks
363	}
364	else
365	{
366	delete trackCopy;
367	}
368
369	}
370
371	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
372	fCurEvent++;
373	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
374	if (fCurEvent== fNumberofEvent)//if end of current file close all
375	{
376	fTree->Reset();
377	delete fTree;
378	fEsdFile->Close();
379	}
380	return hbtEvent;
381	}
382
383
384
385
386
387
388
389
390