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

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// AliFemtoEventReaderESDChainKine - the reader class for the Alice ESD and   //
// the model Kinematics information tailored for the Task framework and the   //
// Reads in AliESDfriend to create shared hit/quality information             //
// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                             //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
#include "AliFemtoEventReaderESDChainKine.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"

#include "TParticle.h"
#include "AliFemtoModelHiddenInfo.h"
#include "AliGenHijingEventHeader.h"

ClassImp(AliFemtoEventReaderESDChainKine)

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

using namespace std;
//____________________________
AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
  fFileName(" "),
  fConstrained(true),
  fUseTPCOnly(false),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fEvent(0x0),
  fStack(0x0),
  fGenHeader(0x0)
{
  //constructor with 0 parameters , look at default settings 
}

//__________________
AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
  AliFemtoEventReader(aReader),
  fFileName(" "),
  fConstrained(true),
  fUseTPCOnly(false),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fEvent(0x0),
  fStack(0x0),
  fGenHeader(0x0)
{
  // Copy constructor
  fConstrained = aReader.fConstrained;
  fUseTPCOnly = aReader.fUseTPCOnly;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  fEvent = new AliESDEvent();
  fStack = aReader.fStack;
}
//__________________
AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
{
  //Destructor
  delete fEvent;
}

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

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

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

//__________________
void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
{
  // Select whether to read constrained or not constrained momentum
  fConstrained=constrained;
}
//__________________
bool AliFemtoEventReaderESDChainKine::GetConstrained() const
{
  // Check whether we read constrained or not constrained momentum
  return fConstrained;
}
//__________________
AliFemtoEvent* AliFemtoEventReaderESDChainKine::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];
  double fVCov[6];
  if (fUseTPCOnly) {
    fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
    fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
    if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
      fVCov[4] = -1001.0;
  }
  else {
    fEvent->GetPrimaryVertex()->GetXYZ(fV1);
    fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
    if (!fEvent->GetPrimaryVertex()->GetStatus())
      fVCov[4] = -1001.0;
  }

  AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
  
  hbtEvent->SetPrimVertPos(vertex);
  hbtEvent->SetPrimVertCov(fVCov);

  AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
	
  Double_t tReactionPlane = 0;
  if (hdh)
    {
      tReactionPlane = hdh->ReactionPlaneAngle();
    }
  //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

  Int_t *motherids;
  motherids = new Int_t[fStack->GetNtrack()];
  for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;

  // Read in mother ids
  TParticle *motherpart;
  for (int ip=0; ip<fStack->GetNtrack(); ip++) {
    motherpart = fStack->Particle(ip);
    if (motherpart->GetDaughter(0) > 0)
      motherids[motherpart->GetDaughter(0)] = ip;
    if (motherpart->GetDaughter(1) > 0)
      motherids[motherpart->GetDaughter(1)] = ip;

//     if (motherpart->GetPdgCode() == 211) {
//       cout << "Mother " << ip << " has daughters " 
// 	   << motherpart->GetDaughter(0) << " " 
// 	   << motherpart->GetDaughter(1) << " " 
// 	   << motherpart->Vx() << " " 
// 	   << motherpart->Vy() << " " 
// 	   << motherpart->Vz() << " " 
// 	   << 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];
      double rxyz[3];
      double impact[2];
      double covimpact[3];
      
      if (fUseTPCOnly) {
	if (!esdtrack->GetTPCInnerParam()) {
	  delete trackCopy;
	  continue;
	}

	AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
	param->GetXYZ(rxyz);
	param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
	param->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);

	//some stuff which could be useful 
	trackCopy->SetImpactD(impact[0]);
	trackCopy->SetImpactZ(impact[1]);
	trackCopy->SetCdd(covimpact[0]);
	trackCopy->SetCdz(covimpact[1]);
	trackCopy->SetCzz(covimpact[2]);
	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));	

	delete param;
      }
      else {
	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);
	
	//some stuff which could be useful 
	float imp[2];
	float cim[3];
	esdtrack->GetImpactParameters(imp,cim);

	impact[0] = imp[0];
	impact[1] = imp[1];
	covimpact[0] = cim[0];
	covimpact[1] = cim[1];
	covimpact[2] = cim[2];

	trackCopy->SetImpactD(impact[0]);
	trackCopy->SetImpactZ(impact[1]);
	trackCopy->SetCdd(covimpact[0]);
	trackCopy->SetCdz(covimpact[1]);
	trackCopy->SetCzz(covimpact[2]);
	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
      }
	    	
      trackCopy->SetTrackId(esdtrack->GetID());
      trackCopy->SetFlags(esdtrack->GetStatus());
      trackCopy->SetLabel(esdtrack->GetLabel());
		
      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 xtpc[3];
      esdtrack->GetInnerXYZ(xtpc);
      xtpc[2] -= fV1[2];
      trackCopy->SetNominalTPCEntrancePoint(xtpc);

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

      int indexes[3];
      for (int ik=0; ik<3; ik++) {
	indexes[ik] = esdtrack->GetKinkIndex(ik);
      }
      trackCopy->SetKinkIndexes(indexes);

      // Fill the hidden information with the simulated data
      TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));

      // Check the mother information

      // Using the new way of storing the freeze-out information
      // Final state particle is stored twice on the stack
      // one copy (mother) is stored with original freeze-out information
      //   and is not tracked
      // the other one (daughter) is stored with primary vertex position
      //   and is tracked
	
      // Freeze-out coordinates
      double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
      fpx = tPart->Vx() - fV1[0];
      fpy = tPart->Vy() - fV1[1];
      fpz = tPart->Vz() - fV1[2];
      fpt = tPart->T();

      if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
 	TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
 	// Check if this is the same particle stored twice on the stack
  	if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
	  // It is the same particle
	  // Read in the original freeze-out information
	  // and convert it from to [fm]
 	  fpx = mother->Vx()*1e13;
 	  fpy = mother->Vy()*1e13;
 	  fpz = mother->Vz()*1e13;
 	  fpt = mother->T()*1e13*3e10;
	  
	}
      }

      AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
      tInfo->SetPDGPid(tPart->GetPdgCode());
      tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
      Double_t mass2 = (tPart->Energy() *tPart->Energy() -
			tPart->Px()*tPart->Px() -
			tPart->Py()*tPart->Py() -
			tPart->Pz()*tPart->Pz());
      if (mass2>0.0)
	tInfo->SetMass(TMath::Sqrt(mass2));
      else 
	tInfo->SetMass(0.0);

      tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
      trackCopy->SetHiddenInfo(tInfo);
      
      //decision if we want this track
      //if we using diffrent labels we want that this label was use for first time 
      //if we use hidden info we want to have match between sim data and ESD
      if (tGoodMomentum==true)
	{
	  hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	  realnofTracks++;//real number of tracks
	  //	  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 AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
{
  // The chain loads the ESD for us
  // You must provide the address where it can be found
  fEvent = aESD;
}
//___________________
void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
{
  // The chain loads the stack for us
  // You must provide the address where it can be found
  fStack = aStack;
}
//___________________
void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
{
  // The chain loads the generator event header for us
  // You must provide the address where it can be found
  fGenHeader = aGenHeader;
}

//__________________
void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
{
  fUseTPCOnly=usetpconly;
}

bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
{
  return fUseTPCOnly;
}
Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
{
  // Calculates the number of sigma to the vertex.

  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];

  b[0] = impact[0];
  b[1] = impact[1];
  bCov[0] = covar[0];
  bCov[1] = covar[1];
  bCov[2] = covar[2];

  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
0b3bd1ac	1	////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoEventReaderESDChainKine - the reader class for the Alice ESD and //
	4	// the model Kinematics information tailored for the Task framework and the //
	5	// Reads in AliESDfriend to create shared hit/quality information //
	6	// Authors: Adam Kisiel kisiel@mps.ohio-state.edu //
	7	// //
	8	////////////////////////////////////////////////////////////////////////////////
	9	#include "AliFemtoEventReaderESDChainKine.h"
	10
	11	#include "TFile.h"
	12	#include "TTree.h"
	13	#include "AliESDEvent.h"
	14	#include "AliESDtrack.h"
	15	#include "AliESDVertex.h"
	16
	17	#include "AliFmPhysicalHelixD.h"
	18	#include "AliFmThreeVectorF.h"
	19
	20	#include "SystemOfUnits.h"
	21
	22	#include "AliFemtoEvent.h"
	23
0b3bd1ac	24	#include "TParticle.h"
0b3bd1ac	25	#include "AliFemtoModelHiddenInfo.h"
aad44575	26	#include "AliGenHijingEventHeader.h"
0b3bd1ac	27
	28	ClassImp(AliFemtoEventReaderESDChainKine)
	29
	30	#if !(ST_NO_NAMESPACES)
	31	using namespace units;
	32	#endif
	33
	34	using namespace std;
	35	//____________________________
	36	AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
	37	fFileName(" "),
	38	fConstrained(true),
613e9805	39	fUseTPCOnly(false),
0b3bd1ac	40	fNumberofEvent(0),
	41	fCurEvent(0),
	42	fCurFile(0),
	43	fEvent(0x0),
aad44575	44	fStack(0x0),
aad44575	45	fGenHeader(0x0)
0b3bd1ac	46	{
	47	//constructor with 0 parameters , look at default settings
	48	}
	49
	50	//__________________
	51	AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
fcda1d4e	52	AliFemtoEventReader(aReader),
0b3bd1ac	53	fFileName(" "),
0b3bd1ac	54	fConstrained(true),
613e9805	55	fUseTPCOnly(false),
0b3bd1ac	56	fNumberofEvent(0),
	57	fCurEvent(0),
	58	fCurFile(0),
	59	fEvent(0x0),
aad44575	60	fStack(0x0),
aad44575	61	fGenHeader(0x0)
0b3bd1ac	62	{
	63	// Copy constructor
	64	fConstrained = aReader.fConstrained;
613e9805	65	fUseTPCOnly = aReader.fUseTPCOnly;
0b3bd1ac	66	fNumberofEvent = aReader.fNumberofEvent;
	67	fCurEvent = aReader.fCurEvent;
	68	fCurFile = aReader.fCurFile;
	69	fEvent = new AliESDEvent();
	70	fStack = aReader.fStack;
	71	}
	72	//__________________
	73	AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
	74	{
	75	//Destructor
	76	delete fEvent;
	77	}
	78
	79	//__________________
	80	AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
	81	{
	82	// Assignment operator
	83	if (this == &aReader)
	84	return *this;
	85
	86	fConstrained = aReader.fConstrained;
613e9805	87	fUseTPCOnly = aReader.fUseTPCOnly;
0b3bd1ac	88	fNumberofEvent = aReader.fNumberofEvent;
	89	fCurEvent = aReader.fCurEvent;
	90	fCurFile = aReader.fCurFile;
	91	if (fEvent) delete fEvent;
	92	fEvent = new AliESDEvent();
	93	fStack = aReader.fStack;
aad44575	94	fGenHeader = aReader.fGenHeader;
aad44575	95
0b3bd1ac	96	return *this;
	97	}
	98	//__________________
	99	// Simple report
	100	AliFemtoString AliFemtoEventReaderESDChainKine::Report()
	101	{
	102	AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
	103	return temp;
	104	}
	105
	106	//__________________
	107	void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
	108	{
	109	// Select whether to read constrained or not constrained momentum
	110	fConstrained=constrained;
	111	}
	112	//__________________
	113	bool AliFemtoEventReaderESDChainKine::GetConstrained() const
	114	{
	115	// Check whether we read constrained or not constrained momentum
	116	return fConstrained;
	117	}
	118	//__________________
	119	AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
	120	{
	121	// Get the event, read all the relevant information
	122	// and fill the AliFemtoEvent class
	123	// Returns a valid AliFemtoEvent
	124	AliFemtoEvent *hbtEvent = 0;
	125	string tFriendFileName;
	126
	127	// Get the friend information
	128	cout<<"starting to read event "<<fCurEvent<<endl;
	129	// fEvent->SetESDfriend(fEventFriend);
	130
	131	hbtEvent = new AliFemtoEvent;
	132	//setting basic things
	133	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
	134	hbtEvent->SetRunNumber(fEvent->GetRunNumber());
	135	//hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
	136	hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
	137	hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
	138	hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
	139	hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
	140	hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
	141	hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
	142	hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
	143	hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
	144	hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
	145
	146	//Vertex
	147	double fV1[3];
63a5982a	148	double fVCov[6];
613e9805	149	if (fUseTPCOnly) {
613e9805	150	fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
63a5982a	151	fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
	152	if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
	153	fVCov[4] = -1001.0;
613e9805	154	}
	155	else {
	156	fEvent->GetPrimaryVertex()->GetXYZ(fV1);
63a5982a	157	fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
	158	if (!fEvent->GetPrimaryVertex()->GetStatus())
	159	fVCov[4] = -1001.0;
613e9805	160	}
0b3bd1ac	161
0b3bd1ac	162	AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
63a5982a	163
0b3bd1ac	164	hbtEvent->SetPrimVertPos(vertex);
63a5982a	165	hbtEvent->SetPrimVertCov(fVCov);
aad44575	166
aad44575	167	AliGenHijingEventHeader hdh = dynamic_cast<AliGenHijingEventHeader > (fGenHeader);
0b3bd1ac	168
aad44575	169	Double_t tReactionPlane = 0;
	170	if (hdh)
	171	{
	172	tReactionPlane = hdh->ReactionPlaneAngle();
	173	}
0b3bd1ac	174	//starting to reading tracks
	175	int nofTracks=0; //number of reconstructed tracks in event
	176	nofTracks=fEvent->GetNumberOfTracks();
	177	int realnofTracks=0;//number of track which we use ina analysis
	178
613e9805	179	Int_t *motherids;
	180	motherids = new Int_t[fStack->GetNtrack()];
	181	for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
	182
	183	// Read in mother ids
	184	TParticle *motherpart;
	185	for (int ip=0; ip<fStack->GetNtrack(); ip++) {
	186	motherpart = fStack->Particle(ip);
	187	if (motherpart->GetDaughter(0) > 0)
	188	motherids[motherpart->GetDaughter(0)] = ip;
	189	if (motherpart->GetDaughter(1) > 0)
	190	motherids[motherpart->GetDaughter(1)] = ip;
	191
	192	// if (motherpart->GetPdgCode() == 211) {
	193	// cout << "Mother " << ip << " has daughters "
	194	// << motherpart->GetDaughter(0) << " "
	195	// << motherpart->GetDaughter(1) << " "
	196	// << motherpart->Vx() << " "
	197	// << motherpart->Vy() << " "
	198	// << motherpart->Vz() << " "
	199	// << endl;
	200
	201	// }
	202	}
	203
0b3bd1ac	204	for (int i=0;i<nofTracks;i++)
	205	{
	206	bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
	207
	208	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	209	const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
	210	// const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
	211
	212	trackCopy->SetCharge((short)esdtrack->GetSign());
	213
	214	//in aliroot we have AliPID
	215	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	216	//we use only 5 first
	217	double esdpid[5];
	218	esdtrack->GetESDpid(esdpid);
	219	trackCopy->SetPidProbElectron(esdpid[0]);
	220	trackCopy->SetPidProbMuon(esdpid[1]);
	221	trackCopy->SetPidProbPion(esdpid[2]);
	222	trackCopy->SetPidProbKaon(esdpid[3]);
	223	trackCopy->SetPidProbProton(esdpid[4]);
	224
	225	double pxyz[3];
613e9805	226	double rxyz[3];
	227	double impact[2];
	228	double covimpact[3];
	229
	230	if (fUseTPCOnly) {
	231	if (!esdtrack->GetTPCInnerParam()) {
	232	delete trackCopy;
	233	continue;
	234	}
	235
	236	AliExternalTrackParam param = new AliExternalTrackParam(esdtrack->GetTPCInnerParam());
	237	param->GetXYZ(rxyz);
	238	param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
	239	param->GetPxPyPz(pxyz);//reading noconstarined momentum
0b3bd1ac	240
613e9805	241	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	242	if (v.mag() < 0.0001) {
	243	// cout << "Found 0 momentum ???? " <<endl;
	244	delete trackCopy;
	245	continue;
	246	}
	247	trackCopy->SetP(v);//setting momentum
	248	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
	249
	250	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	251	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	252	//setting helix I do not if it is ok
	253	AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
	254	trackCopy->SetHelix(helix);
	255
	256	//some stuff which could be useful
	257	trackCopy->SetImpactD(impact[0]);
	258	trackCopy->SetImpactZ(impact[1]);
	259	trackCopy->SetCdd(covimpact[0]);
	260	trackCopy->SetCdz(covimpact[1]);
	261	trackCopy->SetCzz(covimpact[2]);
	262	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
	263
	264	delete param;
	265	}
	266	else {
	267	if (fConstrained==true)
	268	tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
	269	else
	270	tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
	271
	272	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	273	if (v.mag() < 0.0001) {
	274	// cout << "Found 0 momentum ???? " <<endl;
	275	delete trackCopy;
	276	continue;
	277	}
	278	trackCopy->SetP(v);//setting momentum
	279	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
	280	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	281	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	282	//setting helix I do not if it is ok
	283	AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
	284	trackCopy->SetHelix(helix);
	285
	286	//some stuff which could be useful
	287	float imp[2];
	288	float cim[3];
	289	esdtrack->GetImpactParameters(imp,cim);
	290
	291	impact[0] = imp[0];
	292	impact[1] = imp[1];
	293	covimpact[0] = cim[0];
	294	covimpact[1] = cim[1];
	295	covimpact[2] = cim[2];
	296
	297	trackCopy->SetImpactD(impact[0]);
	298	trackCopy->SetImpactZ(impact[1]);
	299	trackCopy->SetCdd(covimpact[0]);
	300	trackCopy->SetCdz(covimpact[1]);
	301	trackCopy->SetCzz(covimpact[2]);
	302	trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
0b3bd1ac	303	}
0b3bd1ac	304
	305	trackCopy->SetTrackId(esdtrack->GetID());
	306	trackCopy->SetFlags(esdtrack->GetStatus());
	307	trackCopy->SetLabel(esdtrack->GetLabel());
	308
0b3bd1ac	309	trackCopy->SetITSchi2(esdtrack->GetITSchi2());
	310	trackCopy->SetITSncls(esdtrack->GetNcls(0));
	311	trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
	312	trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
	313	trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
	314	trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
	315	trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
	316
	317
	318	trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
	319	trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
	320
0b3bd1ac	321	double xtpc[3];
0b3bd1ac	322	esdtrack->GetInnerXYZ(xtpc);
613e9805	323	xtpc[2] -= fV1[2];
0b3bd1ac	324	trackCopy->SetNominalTPCEntrancePoint(xtpc);
	325
	326	esdtrack->GetOuterXYZ(xtpc);
613e9805	327	xtpc[2] -= fV1[2];
0b3bd1ac	328	trackCopy->SetNominalTPCExitPoint(xtpc);
	329
	330	int indexes[3];
	331	for (int ik=0; ik<3; ik++) {
	332	indexes[ik] = esdtrack->GetKinkIndex(ik);
	333	}
	334	trackCopy->SetKinkIndexes(indexes);
	335
	336	// Fill the hidden information with the simulated data
	337	TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
613e9805	338
	339	// Check the mother information
	340
	341	// Using the new way of storing the freeze-out information
	342	// Final state particle is stored twice on the stack
	343	// one copy (mother) is stored with original freeze-out information
	344	// and is not tracked
	345	// the other one (daughter) is stored with primary vertex position
	346	// and is tracked
	347
	348	// Freeze-out coordinates
	349	double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
1f445b96	350	fpx = tPart->Vx() - fV1[0];
	351	fpy = tPart->Vy() - fV1[1];
	352	fpz = tPart->Vz() - fV1[2];
613e9805	353	fpt = tPart->T();
	354
	355	if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
	356	TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
	357	// Check if this is the same particle stored twice on the stack
	358	if ((mother->GetPdgCode() == tPart->GetPdgCode() \|\| (mother->Px() == tPart->Px()))) {
	359	// It is the same particle
	360	// Read in the original freeze-out information
	361	// and convert it from to [fm]
	362	fpx = mother->Vx()*1e13;
	363	fpy = mother->Vy()*1e13;
	364	fpz = mother->Vz()*1e13;
	365	fpt = mother->T()1e133e10;
	366
	367	}
	368	}
	369
0b3bd1ac	370	AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
aad44575	371	tInfo->SetPDGPid(tPart->GetPdgCode());
	372	tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
	373	Double_t mass2 = (tPart->Energy() *tPart->Energy() -
	374	tPart->Px()*tPart->Px() -
	375	tPart->Py()*tPart->Py() -
	376	tPart->Pz()*tPart->Pz());
0b3bd1ac	377	if (mass2>0.0)
	378	tInfo->SetMass(TMath::Sqrt(mass2));
	379	else
	380	tInfo->SetMass(0.0);
aad44575	381
613e9805	382	tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
0b3bd1ac	383	trackCopy->SetHiddenInfo(tInfo);
613e9805	384
0b3bd1ac	385	//decision if we want this track
	386	//if we using diffrent labels we want that this label was use for first time
	387	//if we use hidden info we want to have match between sim data and ESD
	388	if (tGoodMomentum==true)
	389	{
	390	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	391	realnofTracks++;//real number of tracks
	392	// delete trackCopy;
	393	}
	394	else
	395	{
	396	delete trackCopy;
	397	}
	398
	399	}
613e9805	400
0b3bd1ac	401	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
	402	fCurEvent++;
	403	cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
	404	return hbtEvent;
	405	}
	406	//___________________
	407	void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
	408	{
	409	// The chain loads the ESD for us
	410	// You must provide the address where it can be found
	411	fEvent = aESD;
	412	}
	413	//___________________
	414	void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
	415	{
	416	// The chain loads the stack for us
	417	// You must provide the address where it can be found
	418	fStack = aStack;
	419	}
aad44575	420	//___________________
	421	void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
	422	{
	423	// The chain loads the generator event header for us
	424	// You must provide the address where it can be found
	425	fGenHeader = aGenHeader;
	426	}
0b3bd1ac	427
613e9805	428	//__________________
	429	void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
	430	{
	431	fUseTPCOnly=usetpconly;
	432	}
0b3bd1ac	433
613e9805	434	bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
	435	{
	436	return fUseTPCOnly;
	437	}
	438	Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double impact, double covar)
	439	{
	440	// Calculates the number of sigma to the vertex.
	441
	442	Float_t b[2];
	443	Float_t bRes[2];
	444	Float_t bCov[3];
	445
	446	b[0] = impact[0];
	447	b[1] = impact[1];
	448	bCov[0] = covar[0];
	449	bCov[1] = covar[1];
	450	bCov[2] = covar[2];
	451
	452	bRes[0] = TMath::Sqrt(bCov[0]);
	453	bRes[1] = TMath::Sqrt(bCov[2]);
	454
	455	// -----------------------------------
	456	// How to get to a n-sigma cut?
	457	//
	458	// The accumulated statistics from 0 to d is
	459	//
	460	// -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
	461	// -> 1 - Exp(-d*2) for a 2-dim gauss (dd = dxdx + dydy != n_sigma)
	462	//
	463	// It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)ErfInv(1 - Exp((-x*2)/2)
	464	// Can this be expressed in a different way?
	465
	466	if (bRes[0] == 0 \|\| bRes[1] ==0)
	467	return -1;
	468
	469	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	470
	471	// stupid rounding problem screws up everything:
	472	// if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
	473	if (TMath::Exp(-d * d / 2) < 1e-10)
	474	return 1000;
	475
	476	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	477	return d;
	478	}