[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 "TList.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 "AliFemtoModelGlobalHiddenInfo.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenCocktailEventHeader.h"

#include "AliVertexerTracks.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),
  fRotateToEventPlane(0)
{
  //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),
  fRotateToEventPlane(0)
{
  // Copy constructor
  fConstrained = aReader.fConstrained;
  fUseTPCOnly = aReader.fUseTPCOnly;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  fEvent = new AliESDEvent();
  fStack = aReader.fStack;
  fRotateToEventPlane = aReader.fRotateToEventPlane;
}
//__________________
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;
  fRotateToEventPlane = aReader.fRotateToEventPlane;
  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 << "AliFemtoEventReaderESDChainKine::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 {
    if (fEvent->GetPrimaryVertex()) {
      fEvent->GetPrimaryVertex()->GetXYZ(fV1);
      fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);

      if (!fEvent->GetPrimaryVertex()->GetStatus()) {
	// Get the vertex from SPD
	fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
	fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
	
	
	if (!fEvent->GetPrimaryVertexSPD()->GetStatus())
	  fVCov[4] = -1001.0;
	else {
	  fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
	  fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
	}
      }
    }
    else {
      if (fEvent->GetPrimaryVertexSPD()) {
	fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
	fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
      }
    }
    if ((!fEvent->GetPrimaryVertex()) && (!fEvent->GetPrimaryVertexSPD()))
      {
	cout << "No vertex found !!!" << endl;
	fV1[0] = 10000.0;
	fV1[1] = 10000.0;
	fV1[2] = 10000.0;
	fVCov[4] = -1001.0;
      }
  }

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

  Double_t tReactionPlane = 0;

  AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);	
  if (!hdh) {
    AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
    if (cdh) {
      TList *tGenHeaders = cdh->GetHeaders();
      for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
	hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);	
	if (hdh) break;
      }
    }
  }

  if (hdh)
    {
      tReactionPlane = hdh->ReactionPlaneAngle();
      cout << "Got reaction plane " << tReactionPlane << endl;
    }

  hbtEvent->SetReactionPlaneAngle(tReactionPlane);

  //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++)
    {
      //      cout << "Reading track " << i << endl;
      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()) {
	  cout << "No TPC inner param !" << endl;
	  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

	if (fRotateToEventPlane) {
	  double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
	  double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
	  
	  pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
	  pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
	}

	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	if (v.mag() < 0.0001) {
	  //	  cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << 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
	

	if (fRotateToEventPlane) {
	  double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
	  double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
	  
	  pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
	  pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
	}

	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	if (v.mag() < 0.0001) {
	  //	  cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << 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
      if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
	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();
	
	AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
	tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
	
	fpx *= 1e13;
	fpy *= 1e13;
	fpz *= 1e13;
	fpt *= 1e13;
	
	//      cout << "Looking for mother ids " << endl;
	if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
	  //	cout << "Got mother id" << endl;
	  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]
	    
	    // EPOS style 
	    fpx = mother->Vx()*1e13*0.197327;
	    fpy = mother->Vy()*1e13*0.197327;
	    fpz = mother->Vz()*1e13*0.197327;
	    fpt = mother->T() *1e13*0.197327;
	    
	    
	    // Therminator style 
// 	    fpx = mother->Vx()*1e13;
// 	    fpy = mother->Vy()*1e13;
// 	    fpz = mother->Vz()*1e13;
// 	    fpt = mother->T() *1e13*3e10;
	    
	  }
	}
	
	if (fRotateToEventPlane) {
	  double tPhi = TMath::ATan2(fpy, fpx);
	  double tRad = TMath::Hypot(fpx, fpy);
	  
	  fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
	  fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
	}
	
	tInfo->SetPDGPid(tPart->GetPdgCode());
	
	if (fRotateToEventPlane) {
	  double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
	  double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
	  
	  tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
				 tRad*TMath::Sin(tPhi - tReactionPlane),
				 tPart->Pz());
	}
	else
	  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);
      }
      else {
	AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
	tInfo->SetMass(0.0);
	double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
	fpx = fV1[0]*1e13;
	fpy = fV1[1]*1e13;
	fpz = fV1[2]*1e13;
	fpt = 0.0;
	tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);

	tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
	
	trackCopy->SetHiddenInfo(tInfo);
      }
      //      cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " <<  mass2 << " " << tPart->GetPdgCode() << endl; 
      
      //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;
	}
		
    }

  delete motherids;
  
  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;
}
void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
{
  fRotateToEventPlane=dorotate;
}

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