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

////////////////////////////////////////////////////////////////////////////////
///                                                                          ///
/// AliFemtoEventReaderESDChain - the reader class for the Alice ESD         ///
/// tailored for the Task framework                                 ///
/// Reads in AliESDfriend to create shared hit/quality information           ///
/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
///                                                                          ///
////////////////////////////////////////////////////////////////////////////////
#include "AliFemtoEventReaderESDChain.h"

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

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

#include "SystemOfUnits.h"

#include "AliFemtoEvent.h"

ClassImp(AliFemtoEventReaderESDChain)

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

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

//__________________
AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fEvent(0x0)
{
  // Copy constructor
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  //  fEvent = new AliESD(*aReader.fEvent);
  fEvent = new AliESDEvent();
//   fEventFriend = aReader.fEventFriend;
//   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>();
//     list<Int_t>::iterator iter;
//     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
//       fClusterPerPadrow[tPad]->push_back(*iter);
//     }
//   }
//   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>();
//     list<Int_t>::iterator iter;
//     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
//       fSharedList[tPad]->push_back(*iter);
//     }
//   }
}
//__________________
AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
{
  //Destructor
  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;
}

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

  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  if (fEvent) delete fEvent;
  fEvent = new AliESDEvent();

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

//   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>();
//     list<Int_t>::iterator iter;
//     for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
//       fClusterPerPadrow[tPad]->push_back(*iter);
//     }
//   }
//   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>();
//     list<Int_t>::iterator iter;
//     for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
//       fSharedList[tPad]->push_back(*iter);
//     }
//   }
  
  return *this;
}
//__________________
// Simple report
AliFemtoString AliFemtoEventReaderESDChain::Report()
{
  AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
  return temp;
}

//__________________
void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
{
  // Select whether to read constrained or not constrained momentum
  fConstrained=constrained;
}
//__________________
bool AliFemtoEventReaderESDChain::GetConstrained() const
{
  // Check whether we read constrained or not constrained momentum
  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);
	
  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  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]);
      if (v.mag() < 0.0001) {
	//	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 kP(pxyz[0],pxyz[1],pxyz[2]);
      const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
      //setting helix I do not if it is ok
      AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
      trackCopy->SetHelix(helix);
	    	
      trackCopy->SetTrackId(esdtrack->GetID());
      trackCopy->SetFlags(esdtrack->GetStatus());
      trackCopy->SetLabel(esdtrack->GetLabel());
		
      //some stuff which could be useful 
      float impact[2];
      float covimpact[3];
      esdtrack->GetImpactParameters(impact,covimpact);
      trackCopy->SetImpactD(impact[0]);
      trackCopy->SetImpactZ(impact[1]);
      trackCopy->SetCdd(covimpact[0]);
      trackCopy->SetCdz(covimpact[1]);
      trackCopy->SetCzz(covimpact[2]);
      trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
      trackCopy->SetITSncls(esdtrack->GetNcls(0));     
      trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
      trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
      trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
      trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
      trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 
      trackCopy->SetSigmaToVertex(GetSigmaToVertex(esdtrack));

      trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
      trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());

      double pvrt[3];
      fEvent->GetPrimaryVertex()->GetXYZ(pvrt);

      double xtpc[3];
      esdtrack->GetInnerXYZ(xtpc);
      xtpc[2] -= pvrt[2];
      trackCopy->SetNominalTPCEntrancePoint(xtpc);

      esdtrack->GetOuterXYZ(xtpc);
      xtpc[2] -= pvrt[2];
      trackCopy->SetNominalTPCExitPoint(xtpc);

      int indexes[3];
      for (int ik=0; ik<3; ik++) {
	indexes[ik] = esdtrack->GetKinkIndex(ik);
      }
      trackCopy->SetKinkIndexes(indexes);
      //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
	  //	  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; 
}
//___________________
void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
{
  // The chain loads the ESD for us
  // You must provide the address where it can be found
  fEvent = aESD;
}
//___________________
// void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
// {
//   // We need the ESD tree to obtain
//   // information about the friend file location
//   fEventFriend = aFriend;
// }

//____________________________________________________________________
Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(const AliESDtrack* esdTrack)
{
  // Calculates the number of sigma to the vertex.

  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];
  esdTrack->GetImpactParameters(b,bCov);
//   if (bCov[0]<=0 || bCov[2]<=0) {
//     AliDebug(1, "Estimated b resolution lower or equal zero!");
//     bCov[0]=0; bCov[2]=0;
//   }
  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);

  // -----------------------------------
  // How to get to a n-sigma cut?
  //
  // The accumulated statistics from 0 to d is
  //
  // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
  // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
  //
  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
  // Can this be expressed in a different way?

  if (bRes[0] == 0 || bRes[1] ==0)
    return -1;

  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));

  // stupid rounding problem screws up everything:
  // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
  if (TMath::Exp(-d * d / 2) < 1e-10)
    return 1000;

  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  return d;
}
Commit	Line	Data
d0e92d9a	1	////////////////////////////////////////////////////////////////////////////////
	2	/// ///
	3	/// AliFemtoEventReaderESDChain - the reader class for the Alice ESD ///
	4	/// tailored for the Task framework ///
	5	/// Reads in AliESDfriend to create shared hit/quality information ///
	6	/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
	7	/// ///
	8	////////////////////////////////////////////////////////////////////////////////
67427ff7	9	#include "AliFemtoEventReaderESDChain.h"
	10
	11	#include "TFile.h"
	12	#include "TTree.h"
ea77036b	13	#include "AliESDEvent.h"
67427ff7	14	#include "AliESDtrack.h"
ea77036b	15	#include "AliESDVertex.h"
67427ff7	16
d0e92d9a	17	#include "AliFmPhysicalHelixD.h"
d0e92d9a	18	#include "AliFmThreeVectorF.h"
67427ff7	19
d0e92d9a	20	#include "SystemOfUnits.h"
67427ff7	21
d0e92d9a	22	#include "AliFemtoEvent.h"
67427ff7	23
	24	ClassImp(AliFemtoEventReaderESDChain)
	25
	26	#if !(ST_NO_NAMESPACES)
	27	using namespace units;
	28	#endif
	29
	30	using namespace std;
	31	//____________________________
67427ff7	32	AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
	33	fFileName(" "),
	34	fConstrained(true),
	35	fNumberofEvent(0),
	36	fCurEvent(0),
	37	fCurFile(0),
ea77036b	38	fEvent(0x0)
67427ff7	39	{
d0e92d9a	40	//constructor with 0 parameters , look at default settings
ea77036b	41	// fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	42	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	43	// fClusterPerPadrow[tPad] = new list<Int_t>();
	44	// }
	45	// fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	46	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	47	// fSharedList[tPad] = new list<Int_t>();
	48	// }
67427ff7	49	}
67427ff7	50
0215f606	51	//__________________
0215f606	52	AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
	53	fFileName(" "),
	54	fConstrained(true),
	55	fNumberofEvent(0),
	56	fCurEvent(0),
	57	fCurFile(0),
ea77036b	58	fEvent(0x0)
0215f606	59	{
d0e92d9a	60	// Copy constructor
0215f606	61	fConstrained = aReader.fConstrained;
	62	fNumberofEvent = aReader.fNumberofEvent;
	63	fCurEvent = aReader.fCurEvent;
	64	fCurFile = aReader.fCurFile;
d0e92d9a	65	// fEvent = new AliESD(*aReader.fEvent);
ea77036b	66	fEvent = new AliESDEvent();
	67	// fEventFriend = aReader.fEventFriend;
	68	// fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	69	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	70	// fClusterPerPadrow[tPad] = new list<Int_t>();
	71	// list<Int_t>::iterator iter;
	72	// for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
	73	// fClusterPerPadrow[tPad]->push_back(*iter);
	74	// }
	75	// }
	76	// fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	77	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	78	// fSharedList[tPad] = new list<Int_t>();
	79	// list<Int_t>::iterator iter;
	80	// for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
	81	// fSharedList[tPad]->push_back(*iter);
	82	// }
	83	// }
0215f606	84	}
67427ff7	85	//__________________
67427ff7	86	AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
67427ff7	87	{
d0e92d9a	88	//Destructor
67427ff7	89	delete fEvent;
67427ff7	90
ea77036b	91	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	92	// fClusterPerPadrow[tPad]->clear();
	93	// delete fClusterPerPadrow[tPad];
	94	// }
	95	// delete [] fClusterPerPadrow;
	96	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	97	// fSharedList[tPad]->clear();
	98	// delete fSharedList[tPad];
	99	// }
	100	// delete [] fSharedList;
67427ff7	101	}
67427ff7	102
0215f606	103	//__________________
0215f606	104	AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
0215f606	105	{
d0e92d9a	106	// Assignment operator
0215f606	107	if (this == &aReader)
	108	return *this;
	109
	110	fConstrained = aReader.fConstrained;
	111	fNumberofEvent = aReader.fNumberofEvent;
	112	fCurEvent = aReader.fCurEvent;
	113	fCurFile = aReader.fCurFile;
	114	if (fEvent) delete fEvent;
ea77036b	115	fEvent = new AliESDEvent();
0215f606	116
ea77036b	117	// fEventFriend = aReader.fEventFriend;
0215f606	118
ea77036b	119	// if (fClusterPerPadrow) {
	120	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	121	// fClusterPerPadrow[tPad]->clear();
	122	// delete fClusterPerPadrow[tPad];
	123	// }
	124	// delete [] fClusterPerPadrow;
	125	// }
0215f606	126
ea77036b	127	// if (fSharedList) {
	128	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	129	// fSharedList[tPad]->clear();
	130	// delete fSharedList[tPad];
	131	// }
	132	// delete [] fSharedList;
	133	// }
	134
	135	// fClusterPerPadrow = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	136	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	137	// fClusterPerPadrow[tPad] = new list<Int_t>();
	138	// list<Int_t>::iterator iter;
	139	// for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
	140	// fClusterPerPadrow[tPad]->push_back(*iter);
	141	// }
	142	// }
	143	// fSharedList = (list<Int_t> *) malloc(sizeof(list<Int_t> ) * AliESDfriendTrack::kMaxTPCcluster);
	144	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	145	// fSharedList[tPad] = new list<Int_t>();
	146	// list<Int_t>::iterator iter;
	147	// for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
	148	// fSharedList[tPad]->push_back(*iter);
	149	// }
	150	// }
0215f606	151
	152	return *this;
	153	}
67427ff7	154	//__________________
	155	// Simple report
	156	AliFemtoString AliFemtoEventReaderESDChain::Report()
	157	{
	158	AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
	159	return temp;
	160	}
	161
	162	//__________________
67427ff7	163	void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
67427ff7	164	{
d0e92d9a	165	// Select whether to read constrained or not constrained momentum
67427ff7	166	fConstrained=constrained;
	167	}
	168	//__________________
67427ff7	169	bool AliFemtoEventReaderESDChain::GetConstrained() const
67427ff7	170	{
d0e92d9a	171	// Check whether we read constrained or not constrained momentum
67427ff7	172	return fConstrained;
	173	}
	174	//__________________
	175	AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
	176	{
	177	// Get the event, read all the relevant information
	178	// and fill the AliFemtoEvent class
	179	// Returns a valid AliFemtoEvent
	180	AliFemtoEvent *hbtEvent = 0;
	181	string tFriendFileName;
	182
	183	// Get the friend information
	184	cout<<"starting to read event "<<fCurEvent<<endl;
ea77036b	185	// fEvent->SetESDfriend(fEventFriend);
67427ff7	186
	187	hbtEvent = new AliFemtoEvent;
	188	//setting basic things
a531c6fc	189	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7	190	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	191	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	192	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	193	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	194	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	195	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	196	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	197	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	198	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	199	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	200	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	201
	202	//Vertex
	203	double fV1[3];
	204	fEvent->GetVertex()->GetXYZ(fV1);
	205
	206	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
	207	hbtEvent->SetPrimVertPos(vertex);
	208
	209	//starting to reading tracks
	210	int nofTracks=0; //number of reconstructed tracks in event
	211	nofTracks=fEvent->GetNumberOfTracks();
	212	int realnofTracks=0;//number of track which we use ina analysis
	213
ea77036b	214	// // Clear the shared cluster list
	215	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	216	// fClusterPerPadrow[tPad]->clear();
	217	// }
	218	// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
	219	// fSharedList[tPad]->clear();
	220	// }
	221
	222
	223	// for (int i=0;i<nofTracks;i++) {
	224	// const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	225
	226	// list<Int_t>::iterator tClustIter;
	227
	228	// Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
	229	// Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
	230	// for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
	231	// if (tTrackIndices[tNcl] >= 0) {
	232	// tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
	233	// if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
	234	// fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
	235	// }
	236	// else {
	237	// fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
	238	// }
	239	// }
	240	// }
67427ff7	241
ea77036b	242	// }
67427ff7	243
	244	for (int i=0;i<nofTracks;i++)
	245	{
d0e92d9a	246	bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
67427ff7	247
	248	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	249	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
0215f606	250	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
67427ff7	251
	252	trackCopy->SetCharge((short)esdtrack->GetSign());
	253
	254	//in aliroot we have AliPID
	255	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	256	//we use only 5 first
	257	double esdpid[5];
	258	esdtrack->GetESDpid(esdpid);
	259	trackCopy->SetPidProbElectron(esdpid[0]);
	260	trackCopy->SetPidProbMuon(esdpid[1]);
	261	trackCopy->SetPidProbPion(esdpid[2]);
	262	trackCopy->SetPidProbKaon(esdpid[3]);
	263	trackCopy->SetPidProbProton(esdpid[4]);
	264
	265	double pxyz[3];
	266	if (fConstrained==true)
d0e92d9a	267	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
67427ff7	268	else
d0e92d9a	269	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
67427ff7	270
67427ff7	271	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
0b3bd1ac	272	if (v.mag() < 0.0001) {
67427ff7	273	// cout << "Found 0 momentum ???? " <<endl;
	274	delete trackCopy;
	275	continue;
	276	}
	277	trackCopy->SetP(v);//setting momentum
	278	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
d0e92d9a	279	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
d0e92d9a	280	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
67427ff7	281	//setting helix I do not if it is ok
d0e92d9a	282	AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
67427ff7	283	trackCopy->SetHelix(helix);
	284
	285	trackCopy->SetTrackId(esdtrack->GetID());
	286	trackCopy->SetFlags(esdtrack->GetStatus());
d9241b1b	287	trackCopy->SetLabel(esdtrack->GetLabel());
67427ff7	288
	289	//some stuff which could be useful
	290	float impact[2];
	291	float covimpact[3];
	292	esdtrack->GetImpactParameters(impact,covimpact);
	293	trackCopy->SetImpactD(impact[0]);
	294	trackCopy->SetImpactZ(impact[1]);
	295	trackCopy->SetCdd(covimpact[0]);
	296	trackCopy->SetCdz(covimpact[1]);
	297	trackCopy->SetCzz(covimpact[2]);
	298	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
	299	trackCopy->SetITSncls(esdtrack->GetNcls(0));
	300	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
	301	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
	302	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
	303	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
	304	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
0b3bd1ac	305	trackCopy->SetSigmaToVertex(GetSigmaToVertex(esdtrack));
ea77036b	306
	307	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
	308	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
67427ff7	309
0b3bd1ac	310	double pvrt[3];
	311	fEvent->GetPrimaryVertex()->GetXYZ(pvrt);
	312
	313	double xtpc[3];
	314	esdtrack->GetInnerXYZ(xtpc);
	315	xtpc[2] -= pvrt[2];
	316	trackCopy->SetNominalTPCEntrancePoint(xtpc);
	317
	318	esdtrack->GetOuterXYZ(xtpc);
	319	xtpc[2] -= pvrt[2];
	320	trackCopy->SetNominalTPCExitPoint(xtpc);
	321
	322	int indexes[3];
	323	for (int ik=0; ik<3; ik++) {
	324	indexes[ik] = esdtrack->GetKinkIndex(ik);
	325	}
	326	trackCopy->SetKinkIndexes(indexes);
	327	//decision if we want this track
	328	//if we using diffrent labels we want that this label was use for first time
	329	//if we use hidden info we want to have match between sim data and ESD
d0e92d9a	330	if (tGoodMomentum==true)
67427ff7	331	{
	332	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	333	realnofTracks++;//real number of tracks
	334	// delete trackCopy;
	335	}
	336	else
	337	{
	338	delete trackCopy;
	339	}
	340
	341	}
	342
	343	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
	344	fCurEvent++;
	345	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
	346	return hbtEvent;
	347	}
	348	//___________________
ea77036b	349	void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
67427ff7	350	{
d0e92d9a	351	// The chain loads the ESD for us
d0e92d9a	352	// You must provide the address where it can be found
67427ff7	353	fEvent = aESD;
	354	}
	355	//___________________
ea77036b	356	// void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
	357	// {
	358	// // We need the ESD tree to obtain
	359	// // information about the friend file location
	360	// fEventFriend = aFriend;
	361	// }
67427ff7	362
0b3bd1ac	363	//____________________________________________________________________
	364	Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(const AliESDtrack* esdTrack)
	365	{
	366	// Calculates the number of sigma to the vertex.
	367
	368	Float_t b[2];
	369	Float_t bRes[2];
	370	Float_t bCov[3];
	371	esdTrack->GetImpactParameters(b,bCov);
	372	// if (bCov[0]<=0 \|\| bCov[2]<=0) {
	373	// AliDebug(1, "Estimated b resolution lower or equal zero!");
	374	// bCov[0]=0; bCov[2]=0;
	375	// }
	376	bRes[0] = TMath::Sqrt(bCov[0]);
	377	bRes[1] = TMath::Sqrt(bCov[2]);
	378
	379	// -----------------------------------
	380	// How to get to a n-sigma cut?
	381	//
	382	// The accumulated statistics from 0 to d is
	383	//
	384	// -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
	385	// -> 1 - Exp(-d*2) for a 2-dim gauss (dd = dxdx + dydy != n_sigma)
	386	//
	387	// It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)ErfInv(1 - Exp((-x*2)/2)
	388	// Can this be expressed in a different way?
	389
	390	if (bRes[0] == 0 \|\| bRes[1] ==0)
	391	return -1;
	392
	393	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	394
	395	// stupid rounding problem screws up everything:
	396	// if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
	397	if (TMath::Exp(-d * d / 2) < 1e-10)
	398	return 1000;
	399
	400	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	401	return d;
	402	}
67427ff7	403
	404
	405
	406
	407
	408
	409
	410