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

#include "AliFemtoEventReaderESDChain.h"

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

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

#include "Base/SystemOfUnits.h"

#include "Infrastructure/AliFemtoEvent.h"

ClassImp(AliFemtoEventReaderESDChain)

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

using namespace std;
//____________________________
//constructor with 0 parameters , look at default settings 
AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fEvent(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
AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
{
  delete fEvent;

  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;
}

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

//__________________
// Select whether to read constrained or not constrained momentum
void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
{
  fConstrained=constrained;
}
//__________________
// Check whether we read constrained or not constrained momentum
bool AliFemtoEventReaderESDChain::GetConstrained() const
{
  return fConstrained;
}
//__________________
AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
{
  // Get the event, read all the relevant information
  // and fill the AliFemtoEvent class
  // Returns a valid AliFemtoEvent
  AliFemtoEvent *hbtEvent = 0;
  string tFriendFileName;

  // Get the friend information
  cout<<"starting to read event "<<fCurEvent<<endl;
  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]);
      if (v.mag() == 0) {
	//	cout << "Found 0 momentum ???? " <<endl;
	delete trackCopy;
	continue;
      }
      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]);
      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);
	  }
	}
      }

      if (good_momentum==true)
	{
	  hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	  realnofTracks++;//real number of tracks
	  //	  delete trackCopy;
	}
      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;
  return hbtEvent; 
}
//___________________
// The chain loads the ESD for us
// You must provide the address where it can be found
void AliFemtoEventReaderESDChain::SetESDSource(AliESD *aESD)
{
  fEvent = aESD;
}
//___________________
// We need the ESD tree to obtain
// information about the friend file location
void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
{
  fEventFriend = aFriend;
}
Commit	Line	Data
67427ff7	1	#include "AliFemtoEventReaderESDChain.h"
	2
	3	#include "TFile.h"
	4	#include "TTree.h"
	5	#include "AliESD.h"
	6	#include "AliESDtrack.h"
	7
b2f60a91	8	#include "Infrastructure/AliFmPhysicalHelixD.h"
b2f60a91	9	#include "Infrastructure/AliFmThreeVectorF.h"
67427ff7	10
b2f60a91	11	#include "Base/SystemOfUnits.h"
67427ff7	12
b2f60a91	13	#include "Infrastructure/AliFemtoEvent.h"
67427ff7	14
	15	ClassImp(AliFemtoEventReaderESDChain)
	16
	17	#if !(ST_NO_NAMESPACES)
	18	using namespace units;
	19	#endif
	20
	21	using namespace std;
	22	//____________________________
	23	//constructor with 0 parameters , look at default settings
	24	AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
	25	fFileName(" "),
	26	fConstrained(true),
	27	fNumberofEvent(0),
	28	fCurEvent(0),
	29	fCurFile(0),
	30	fEvent(0x0),
	31	fEventFriend(0)
	32	{
	33	fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	34	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	35	fClusterPerPadrow[tPad] = new list<Int_t>();
	36	}
	37	fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	38	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	39	fSharedList[tPad] = new list<Int_t>();
	40	}
	41	}
	42
	43	//__________________
	44	//Destructor
	45	AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
	46	{
	47	delete fEvent;
	48
	49	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	50	fClusterPerPadrow[tPad]->clear();
	51	delete fClusterPerPadrow[tPad];
	52	}
	53	delete [] fClusterPerPadrow;
	54	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	55	fSharedList[tPad]->clear();
	56	delete fSharedList[tPad];
	57	}
	58	delete [] fSharedList;
	59	}
	60
	61	//__________________
	62	// Simple report
	63	AliFemtoString AliFemtoEventReaderESDChain::Report()
	64	{
	65	AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
	66	return temp;
	67	}
	68
	69	//__________________
	70	// Select whether to read constrained or not constrained momentum
	71	void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
	72	{
	73	fConstrained=constrained;
	74	}
	75	//__________________
	76	// Check whether we read constrained or not constrained momentum
	77	bool AliFemtoEventReaderESDChain::GetConstrained() const
78	{
79	return fConstrained;
80	}
81	//__________________
82	AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
83	{
84	// Get the event, read all the relevant information
85	// and fill the AliFemtoEvent class
86	// Returns a valid AliFemtoEvent
87	AliFemtoEvent *hbtEvent = 0;
88	string tFriendFileName;
89
90	// Get the friend information
91	cout<<"starting to read event "<<fCurEvent<<endl;
92	fEvent->SetESDfriend(fEventFriend);
93	vector<int> label_table;//to check labels
94
95	hbtEvent = new AliFemtoEvent;
96	//setting basic things
a531c6fc	97	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7	98	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	99	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	100	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	101	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	102	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	103	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	104	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	105	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	106	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	107	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	108	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	109
	110	//Vertex
	111	double fV1[3];
	112	fEvent->GetVertex()->GetXYZ(fV1);
	113
	114	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	115	hbtEvent->SetPrimVertPos(vertex);
	116
	117	//starting to reading tracks
	118	int nofTracks=0; //number of reconstructed tracks in event
	119	nofTracks=fEvent->GetNumberOfTracks();
	120	int realnofTracks=0;//number of track which we use ina analysis
	121
	122	// Clear the shared cluster list
	123	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	124	fClusterPerPadrow[tPad]->clear();
	125	}
	126	for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	127	fSharedList[tPad]->clear();
	128	}
	129
	130
	131	for (int i=0;i<nofTracks;i++) {
	132	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	133
	134	list<Int_t>::iterator tClustIter;
	135
	136	Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
	137	Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
	138	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	139	if (tTrackIndices[tNcl] >= 0) {
	140	tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
	141	if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
	142	fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
	143	}
	144	else {
	145	fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
	146	}
	147	}
	148	}
	149
	150	}
	151
	152	for (int i=0;i<nofTracks;i++)
	153	{
	154	bool good_momentum=true; //flaga to chcek if we can read momentum of this track
	155
	156	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	157	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	158	const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
	159
	160	trackCopy->SetCharge((short)esdtrack->GetSign());
	161
162	//in aliroot we have AliPID
163	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
164	//we use only 5 first
165	double esdpid[5];
166	esdtrack->GetESDpid(esdpid);
167	trackCopy->SetPidProbElectron(esdpid[0]);
168	trackCopy->SetPidProbMuon(esdpid[1]);
169	trackCopy->SetPidProbPion(esdpid[2]);
170	trackCopy->SetPidProbKaon(esdpid[3]);
171	trackCopy->SetPidProbProton(esdpid[4]);
172
173	double pxyz[3];
174	if (fConstrained==true)
175	good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
176	else
177	good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
178
179	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
180	if (v.mag() == 0) {
181	// cout << "Found 0 momentum ???? " <<endl;
182	delete trackCopy;
183	continue;
184	}
185	trackCopy->SetP(v);//setting momentum
186	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
187	const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
188	const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
189	//setting helix I do not if it is ok
190	AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
191	trackCopy->SetHelix(helix);
192
193	trackCopy->SetTrackId(esdtrack->GetID());
194	trackCopy->SetFlags(esdtrack->GetStatus());
195	//trackCopy->SetLabel(esdtrack->GetLabel());
196
197	//some stuff which could be useful
198	float impact[2];
199	float covimpact[3];
200	esdtrack->GetImpactParameters(impact,covimpact);
201	trackCopy->SetImpactD(impact[0]);
202	trackCopy->SetImpactZ(impact[1]);
203	trackCopy->SetCdd(covimpact[0]);
204	trackCopy->SetCdz(covimpact[1]);
205	trackCopy->SetCzz(covimpact[2]);
206	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
207	trackCopy->SetITSncls(esdtrack->GetNcls(0));
208	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
209	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
210	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
211	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
212	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
213
214	// Fill cluster per padrow information
215	Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
216	Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
217	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
218	if (tTrackIndices[tNcl] > 0)
219	trackCopy->SetTPCcluster(tNcl, 1);
220	else
221	trackCopy->SetTPCcluster(tNcl, 0);
222	}
223
224	// Fill shared cluster information
225	list<Int_t>::iterator tClustIter;
226
227	for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
228	if (tTrackIndices[tNcl] > 0) {
229	tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
230	if (tClustIter != fSharedList[tNcl]->end()) {
231	trackCopy->SetTPCshared(tNcl, 1);
232	cout << "Event next" << endl;
233	cout << "Track: " << i << endl;
234	cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
235	}
236	else {
237	trackCopy->SetTPCshared(tNcl, 0);
238	}
239	}
240	}
241
242	if (good_momentum==true)
243	{
244	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
245	realnofTracks++;//real number of tracks
246	// delete trackCopy;
247	}
248	else
249	{
250	delete trackCopy;
251	}
252
253	}
254
255	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
256	fCurEvent++;
257	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
258	return hbtEvent;
259	}
260	//___________________
261	// The chain loads the ESD for us
262	// You must provide the address where it can be found
263	void AliFemtoEventReaderESDChain::SetESDSource(AliESD *aESD)
264	{
265	fEvent = aESD;
266	}
267	//___________________
268	// We need the ESD tree to obtain
269	// information about the friend file location
270	void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
271	{
272	fEventFriend = aFriend;
273	}
274
275
276
277
278
279
280
281
282