[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoEventReaderKinematicsChain.cxx

/////////////////////////////////////////////////////////////////////////////////////
//                                                                                 //
// AliFemtoEventReaderKinematicsChain - 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: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch     //
//          Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
//  	       				                                           //
/////////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoEventReaderKinematicsChain.h"

#include "TFile.h"
#include "TTree.h"
#include "TList.h"

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

#include "SystemOfUnits.h"

#include "AliFemtoEvent.h"

#include "TParticle.h"
#include "AliStack.h"
#include "TParticlePDG.h"
#include "AliFemtoModelHiddenInfo.h"
#include "AliFemtoModelGlobalHiddenInfo.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenCocktailEventHeader.h"

#include "AliVertexerTracks.h"

ClassImp(AliFemtoEventReaderKinematicsChain)

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

using namespace std;
//____________________________
AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain():
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fStack(0x0),
  fGenHeader(0x0),
  fEstEventMult(kGlobalCount),
  fRotateToEventPlane(0)
{
  //constructor with 0 parameters , look at default settings 
}

//__________________
AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader):
  AliFemtoEventReader(aReader),
  fFileName(" "),
  fConstrained(true),
  fNumberofEvent(0),
  fCurEvent(0),
  fCurFile(0),
  fStack(0x0),
  fGenHeader(0x0),
  fEstEventMult(kGlobalCount),
  fRotateToEventPlane(0)
{
  // Copy constructor
  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  fStack = aReader.fStack;
  fEstEventMult = aReader.fEstEventMult;
  fRotateToEventPlane = aReader.fRotateToEventPlane;
}
//__________________
AliFemtoEventReaderKinematicsChain::~AliFemtoEventReaderKinematicsChain()
{
  //Destructor
  //delete fEvent;
}

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

  fConstrained = aReader.fConstrained;
  fNumberofEvent = aReader.fNumberofEvent;
  fCurEvent = aReader.fCurEvent;
  fCurFile = aReader.fCurFile;
  fStack = aReader.fStack;
  fGenHeader = aReader.fGenHeader;
  fEstEventMult = aReader.fEstEventMult;
  fRotateToEventPlane = aReader.fRotateToEventPlane;
  return *this;
}
//__________________
// Simple report
AliFemtoString AliFemtoEventReaderKinematicsChain::Report()
{
  AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChain\n";
  return temp;
}

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

  cout << "AliFemtoEventReaderKinematlaicsChain::Starting to read event: "<<fCurEvent<<endl;
	
  hbtEvent = new AliFemtoEvent;
  //setting basic things
  //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
  hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
  hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
  hbtEvent->SetZDCN1Energy(0);
  hbtEvent->SetZDCP1Energy(0);
  hbtEvent->SetZDCN2Energy(0);
  hbtEvent->SetZDCP2Energy(0);
  hbtEvent->SetZDCEMEnergy(0);
  hbtEvent->SetZDCParticipants(0);
  hbtEvent->SetTriggerMask(0);
  hbtEvent->SetTriggerCluster(0);
	
  //Vertex
  double fV1[3] = {0.0,0.0,0.0};
  double fVCov[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};


  AliFmThreeVectorF vertex(0,0,0);
    
  
  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 all tracks in MC event
  nofTracks=fStack->GetNtrack(); 
  int realnofTracks=0;//number of track which we use in analysis


  int tNormMult = 0;
  int tV0direction = 0;
  for (int i=0;i<nofTracks;i++)
    {
      //take only primaries
      if(!fStack->IsPhysicalPrimary(i)) {continue;}
	  	  
      AliFemtoTrack* trackCopy = new AliFemtoTrack();	
	
      	  //getting next track
      TParticle *kinetrack= fStack->Particle(i);

      //setting multiplicity
        realnofTracks++;//real number of tracks (only primary particles)

      //setting normalized multiplicity
      
	if(kinetrack->GetPDG()->Charge()/3!=0)
	  if (kinetrack->Pt() > 0.15 && kinetrack->Pt() < 20) 
	    if (kinetrack->Eta() < 0.8)
	      tNormMult++;
	  
	//counting particles that go into direction of VZERO detector
	if(kinetrack->Eta() > 2.8 && kinetrack->Eta() < 5.1) //VZERO-A
	  tV0direction++;
	if(kinetrack->Eta() > -3.7 && kinetrack->Eta() < -1.7)//VZERO-C
	  tV0direction++;	
 
	  //charge
      trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));


      //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 kinepid[5];
      for(int pid_iter=0;pid_iter<5;pid_iter++)
	  kinepid[pid_iter]=0;

      int pdgcode = kinetrack->GetPdgCode();
      //proton
      if(pdgcode==2212 || pdgcode==-2212)
        kinepid[4]=1000;
      //kaon
      else if(pdgcode==321 || pdgcode==-321 )
        kinepid[3]=1000;
      //pion
      else if( pdgcode==211 || pdgcode==-211)
        kinepid[2]=1000;
      //electron
      else if(pdgcode==11 || pdgcode==-11)
        kinepid[0]=1000;
      //muon
      else if(pdgcode==13 || pdgcode==-13)
        kinepid[1]=1000;
      else
	continue;
      trackCopy->SetPidProbElectron(kinepid[0]);
      trackCopy->SetPidProbMuon(kinepid[1]);
      trackCopy->SetPidProbPion(kinepid[2]);
      trackCopy->SetPidProbKaon(kinepid[3]);
      trackCopy->SetPidProbProton(kinepid[4]);
					
					
	//Momentum
      double pxyz[3];
      // double rxyz[3];
     
	pxyz[0]=kinetrack->Px();
	pxyz[1]=kinetrack->Py();
	pxyz[2]=kinetrack->Pz();

	// rxyz[0]=kinetrack->Vx();
	// rxyz[1]=kinetrack->Vy();
	// rxyz[2]=kinetrack->Vz();

	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]);

	//label
	trackCopy->SetLabel(i);


	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
	//cout<<"Track added: "<<i<<endl;
		
    }
  
  hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
  if (fEstEventMult == kGlobalCount) 
    hbtEvent->SetNormalizedMult(tNormMult);
  else if(fEstEventMult == kVZERO)
    hbtEvent->SetNormalizedMult(tV0direction);
  fCurEvent++;	


  //V0 analysis code - no V0 finder for Kinematics, we can only check if it is primary and if it has at least 2 daughters.

  for (int i=0;i<nofTracks;i++)
    {
      //do not take primaries
      if(!fStack->IsPhysicalPrimary(i)) {continue;}
      //getting next track
      TParticle *kinetrack= fStack->Particle(i);
      if (!kinetrack) continue;
      
      if(kinetrack->GetPDG()->Charge()!=0) continue; //charge - neutral
      //if(kinetrack->GetPDG()->Stable()==1) continue; //particle is not stable
      if(kinetrack->GetDaughter(0)<1) continue; //has 1'st daughter
      if(kinetrack->GetDaughter(1)<1) continue;  //has 2'nd daughter
 
      //we want one positive, one negative particle. Or two neutral.
      // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()>=0)
      // 	if((fStack->Particle(kinetrack->GetDaughter(1)))->GetPDG()->Charge()>0)
      // 	  continue;
      // if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<=0)
      // 	if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<0)
      // 	  continue;
      
      if(kinetrack->Pt()<0.00001)
	continue;
      
      AliFemtoV0* trackCopyV0 = new AliFemtoV0();
      CopyAODtoFemtoV0(kinetrack, trackCopyV0);
      hbtEvent->V0Collection()->push_back(trackCopyV0);
    //cout<<"Pushback v0 to v0collection"<<endl;
    }


  cout<<"Number of tracks: "<<realnofTracks<<endl;

  return hbtEvent; 
}

//___________________
void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
{
  // The chain loads the stack for us
  // You must provide the address where it can be found
  fStack = aStack;
}
//___________________
void AliFemtoEventReaderKinematicsChain::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 AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
{
  fRotateToEventPlane=dorotate;
}

void AliFemtoEventReaderKinematicsChain::SetUseMultiplicity(EstEventMult aType)
{
  fEstEventMult = aType;
}

Float_t AliFemtoEventReaderKinematicsChain::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;
}


 void AliFemtoEventReaderKinematicsChain::CopyAODtoFemtoV0(TParticle *tv0, AliFemtoV0 *tFemtoV0 )
{
  tFemtoV0->SetEtaV0(tv0->Eta());
  tFemtoV0->SetEtaV0(tv0->Phi());
  tFemtoV0->SetptV0(tv0->Pt()); 
  tFemtoV0->SetptotV0(tv0->P());

  tFemtoV0->SetmomV0X(tv0->Px());
  tFemtoV0->SetmomV0Y(tv0->Py());
  tFemtoV0->SetmomV0Z(tv0->Pz());
  AliFemtoThreeVector momv0(tv0->Px(),tv0->Py(),tv0->Pz());
  tFemtoV0->SetmomV0(momv0);


  TParticle *trackpos;
  TParticle *trackneg;

  //daughters
  if(fStack->Particle(tv0->GetDaughter(0))->GetPDG()->Charge()>=0) //first positive, second negative
    {
      trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
      trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
      tFemtoV0->SetidPos(tv0->GetDaughter(0));
      tFemtoV0->SetidNeg(tv0->GetDaughter(1));
    }
  else //first negative, second positive
    {
      trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
      trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
      tFemtoV0->SetidPos(tv0->GetDaughter(1));
      tFemtoV0->SetidNeg(tv0->GetDaughter(0));
    }

  tFemtoV0->SetEtaPos(trackpos->Eta());
  tFemtoV0->SetEtaNeg(trackneg->Eta());
  
  tFemtoV0->SetptPos(trackpos->Pt());
  tFemtoV0->SetptNeg(trackneg->Pt());
   
  tFemtoV0->SetptotPos(trackpos->P());
  tFemtoV0->SetptotNeg(trackneg->P());

  tFemtoV0->SetmomPosX(trackpos->Px());
  tFemtoV0->SetmomPosY(trackpos->Py());
  tFemtoV0->SetmomPosZ(trackpos->Pz());
  AliFemtoThreeVector mompos(trackpos->Px(),trackpos->Py(),trackpos->Pz());
  tFemtoV0->SetmomPos(mompos);

  tFemtoV0->SetmomNegX(trackneg->Px());
  tFemtoV0->SetmomNegY(trackneg->Py());
  tFemtoV0->SetmomNegZ(trackneg->Pz());
  AliFemtoThreeVector momneg(trackneg->Px(),trackneg->Py(),trackneg->Pz());
  tFemtoV0->SetmomNeg(momneg);

    
  tFemtoV0->SetmassLambda(tv0->GetMass());
  tFemtoV0->SetmassAntiLambda(tv0->GetMass());
  tFemtoV0->SetmassK0Short(tv0->GetMass());

  tFemtoV0->SetYV0(tv0->Y());

  tFemtoV0->SetdecayVertexV0X(trackpos->Vx()); //vertex of the decay is set as the vertex of creation of daughters
  tFemtoV0->SetdecayVertexV0Y(trackpos->Vy());
  tFemtoV0->SetdecayVertexV0Z(trackpos->Vz());
  AliFemtoThreeVector decayvertex(trackpos->Vx(),trackpos->Vy(),trackpos->Vz());
  tFemtoV0->SetdecayVertexV0(decayvertex);

  tFemtoV0->SetdcaV0Daughters(0);
  tFemtoV0->SetCosPointingAngle(1);


  tFemtoV0->SetStatusPos(1);
  tFemtoV0->SetStatusNeg(1);

  
  if(trackpos->GetPdgCode()==2212) //proton
    {
      tFemtoV0->SetPosNSigmaTPCK(1000);
      tFemtoV0->SetPosNSigmaTPCPi(1000);
      tFemtoV0->SetPosNSigmaTPCP(0);
    }
  if(trackneg->GetPdgCode()==-2212) //antiproton
    {
      tFemtoV0->SetNegNSigmaTPCK(1000);
      tFemtoV0->SetNegNSigmaTPCPi(1000);
      tFemtoV0->SetNegNSigmaTPCP(0);
    }
  if(trackpos->GetPdgCode()==211) //pion plus
    {
      tFemtoV0->SetPosNSigmaTPCK(1000);
      tFemtoV0->SetPosNSigmaTPCPi(0);
      tFemtoV0->SetPosNSigmaTPCP(1000);
    }
  if(trackneg->GetPdgCode()==-211) //pion minus
    {
      tFemtoV0->SetNegNSigmaTPCK(1000);
      tFemtoV0->SetNegNSigmaTPCPi(0);
      tFemtoV0->SetNegNSigmaTPCP(1000);
    }
  if(trackpos->GetPdgCode()==321) //K+
    {
      tFemtoV0->SetPosNSigmaTPCK(0);
      tFemtoV0->SetPosNSigmaTPCPi(1000);
      tFemtoV0->SetPosNSigmaTPCP(1000);
    }
  if(trackneg->GetPdgCode()==-321) //K-
    {
      tFemtoV0->SetNegNSigmaTPCK(0);
      tFemtoV0->SetNegNSigmaTPCPi(1000);
      tFemtoV0->SetNegNSigmaTPCP(1000);
    }

  
}
Commit	Line	Data
76ce4b5b	1	/////////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoEventReaderKinematicsChain - 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: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch //
	7	// Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
	8	// //
	9	/////////////////////////////////////////////////////////////////////////////////////
	10
	11	#include "AliFemtoEventReaderKinematicsChain.h"
	12
	13	#include "TFile.h"
	14	#include "TTree.h"
	15	#include "TList.h"
	16
	17	#include "AliFmPhysicalHelixD.h"
	18	#include "AliFmThreeVectorF.h"
	19
	20	#include "SystemOfUnits.h"
	21
	22	#include "AliFemtoEvent.h"
	23
	24	#include "TParticle.h"
	25	#include "AliStack.h"
	26	#include "TParticlePDG.h"
	27	#include "AliFemtoModelHiddenInfo.h"
	28	#include "AliFemtoModelGlobalHiddenInfo.h"
	29	#include "AliGenHijingEventHeader.h"
	30	#include "AliGenCocktailEventHeader.h"
	31
	32	#include "AliVertexerTracks.h"
	33
	34	ClassImp(AliFemtoEventReaderKinematicsChain)
	35
	36	#if !(ST_NO_NAMESPACES)
	37	using namespace units;
	38	#endif
	39
	40	using namespace std;
	41	//____________________________
	42	AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain():
	43	fFileName(" "),
	44	fConstrained(true),
	45	fNumberofEvent(0),
	46	fCurEvent(0),
	47	fCurFile(0),
	48	fStack(0x0),
	49	fGenHeader(0x0),
4eac0b05	50	fEstEventMult(kGlobalCount),
76ce4b5b	51	fRotateToEventPlane(0)
	52	{
	53	//constructor with 0 parameters , look at default settings
	54	}
	55
	56	//__________________
	57	AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader):
	58	AliFemtoEventReader(aReader),
	59	fFileName(" "),
	60	fConstrained(true),
	61	fNumberofEvent(0),
	62	fCurEvent(0),
	63	fCurFile(0),
	64	fStack(0x0),
	65	fGenHeader(0x0),
4eac0b05	66	fEstEventMult(kGlobalCount),
76ce4b5b	67	fRotateToEventPlane(0)
	68	{
	69	// Copy constructor
	70	fConstrained = aReader.fConstrained;
	71	fNumberofEvent = aReader.fNumberofEvent;
	72	fCurEvent = aReader.fCurEvent;
	73	fCurFile = aReader.fCurFile;
	74	fStack = aReader.fStack;
4eac0b05	75	fEstEventMult = aReader.fEstEventMult;
76ce4b5b	76	fRotateToEventPlane = aReader.fRotateToEventPlane;
	77	}
	78	//__________________
	79	AliFemtoEventReaderKinematicsChain::~AliFemtoEventReaderKinematicsChain()
	80	{
	81	//Destructor
	82	//delete fEvent;
	83	}
	84
	85	//__________________
	86	AliFemtoEventReaderKinematicsChain& AliFemtoEventReaderKinematicsChain::operator=(const AliFemtoEventReaderKinematicsChain& aReader)
	87	{
	88	// Assignment operator
	89	if (this == &aReader)
	90	return *this;
	91
	92	fConstrained = aReader.fConstrained;
	93	fNumberofEvent = aReader.fNumberofEvent;
	94	fCurEvent = aReader.fCurEvent;
	95	fCurFile = aReader.fCurFile;
	96	fStack = aReader.fStack;
	97	fGenHeader = aReader.fGenHeader;
4eac0b05	98	fEstEventMult = aReader.fEstEventMult;
76ce4b5b	99	fRotateToEventPlane = aReader.fRotateToEventPlane;
	100	return *this;
	101	}
	102	//__________________
	103	// Simple report
	104	AliFemtoString AliFemtoEventReaderKinematicsChain::Report()
	105	{
	106	AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChain\n";
	107	return temp;
	108	}
	109
	110	//__________________
	111	void AliFemtoEventReaderKinematicsChain::SetConstrained(const bool constrained)
	112	{
	113	// Select whether to read constrained or not constrained momentum
	114	fConstrained=constrained;
	115	}
	116	//__________________
	117	bool AliFemtoEventReaderKinematicsChain::GetConstrained() const
	118	{
	119	// Check whether we read constrained or not constrained momentum
	120	return fConstrained;
	121	}
	122	//__________________
	123	AliFemtoEvent* AliFemtoEventReaderKinematicsChain::ReturnHbtEvent()
	124	{
	125	// Get the event, read all the relevant information from the stack
	126	// and fill the AliFemtoEvent class
	127	// Returns a valid AliFemtoEvent
	128	AliFemtoEvent *hbtEvent = 0;
	129	string tFriendFileName;
	130
	131	cout << "AliFemtoEventReaderKinematlaicsChain::Starting to read event: "<<fCurEvent<<endl;
	132
	133	hbtEvent = new AliFemtoEvent;
	134	//setting basic things
	135	// hbtEvent->SetEventNumber(fEvent->GetEventNumber());
	136	hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
	137	hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
	138	hbtEvent->SetZDCN1Energy(0);
	139	hbtEvent->SetZDCP1Energy(0);
	140	hbtEvent->SetZDCN2Energy(0);
	141	hbtEvent->SetZDCP2Energy(0);
	142	hbtEvent->SetZDCEMEnergy(0);
	143	hbtEvent->SetZDCParticipants(0);
	144	hbtEvent->SetTriggerMask(0);
	145	hbtEvent->SetTriggerCluster(0);
	146
	147	//Vertex
	148	double fV1[3] = {0.0,0.0,0.0};
	149	double fVCov[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
	150
	151
	152	AliFmThreeVectorF vertex(0,0,0);
	153
	154
	155	hbtEvent->SetPrimVertPos(vertex);
	156	hbtEvent->SetPrimVertCov(fVCov);
	157
	158	Double_t tReactionPlane = 0;
	159
	160	AliGenHijingEventHeader hdh = dynamic_cast<AliGenHijingEventHeader > (fGenHeader);
	161	if (!hdh) {
	162	AliGenCocktailEventHeader cdh = dynamic_cast<AliGenCocktailEventHeader > (fGenHeader);
163	if (cdh) {
164	TList *tGenHeaders = cdh->GetHeaders();
165	for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
166	hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
167	if (hdh) break;
168	}
169	}
170	}
171
172	if (hdh)
173	{
174	tReactionPlane = hdh->ReactionPlaneAngle();
175	cout << "Got reaction plane " << tReactionPlane << endl;
176	}
177
178	hbtEvent->SetReactionPlaneAngle(tReactionPlane);
179
180	//starting to reading tracks
181	int nofTracks=0; //number of all tracks in MC event
182	nofTracks=fStack->GetNtrack();
183	int realnofTracks=0;//number of track which we use in analysis
184
185
186	int tNormMult = 0;
4eac0b05	187	int tV0direction = 0;
76ce4b5b	188	for (int i=0;i<nofTracks;i++)
76ce4b5b	189	{
76ce4b5b	190	//take only primaries
	191	if(!fStack->IsPhysicalPrimary(i)) {continue;}
	192
	193	AliFemtoTrack* trackCopy = new AliFemtoTrack();
	194
	195	//getting next track
	196	TParticle *kinetrack= fStack->Particle(i);
	197
	198	//setting multiplicity
	199	realnofTracks++;//real number of tracks (only primary particles)
	200
	201	//setting normalized multiplicity
1777446b	202
76ce4b5b	203	if(kinetrack->GetPDG()->Charge()/3!=0)
1777446b	204	if (kinetrack->Pt() > 0.15 && kinetrack->Pt() < 20)
	205	if (kinetrack->Eta() < 0.8)
	206	tNormMult++;
76ce4b5b	207
4eac0b05	208	//counting particles that go into direction of VZERO detector
	209	if(kinetrack->Eta() > 2.8 && kinetrack->Eta() < 5.1) //VZERO-A
	210	tV0direction++;
	211	if(kinetrack->Eta() > -3.7 && kinetrack->Eta() < -1.7)//VZERO-C
	212	tV0direction++;
	213
76ce4b5b	214	//charge
	215	trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));
	216
	217
	218	//in aliroot we have AliPID
	219	//0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
	220	//we use only 5 first
	221	double kinepid[5];
	222	for(int pid_iter=0;pid_iter<5;pid_iter++)
	223	kinepid[pid_iter]=0;
	224
	225	int pdgcode = kinetrack->GetPdgCode();
	226	//proton
	227	if(pdgcode==2212 \|\| pdgcode==-2212)
	228	kinepid[4]=1000;
	229	//kaon
ce7b3d98	230	else if(pdgcode==321 \|\| pdgcode==-321 )
76ce4b5b	231	kinepid[3]=1000;
76ce4b5b	232	//pion
ce7b3d98	233	else if( pdgcode==211 \|\| pdgcode==-211)
76ce4b5b	234	kinepid[2]=1000;
76ce4b5b	235	//electron
ce7b3d98	236	else if(pdgcode==11 \|\| pdgcode==-11)
76ce4b5b	237	kinepid[0]=1000;
76ce4b5b	238	//muon
ce7b3d98	239	else if(pdgcode==13 \|\| pdgcode==-13)
76ce4b5b	240	kinepid[1]=1000;
ce7b3d98	241	else
ce7b3d98	242	continue;
76ce4b5b	243	trackCopy->SetPidProbElectron(kinepid[0]);
	244	trackCopy->SetPidProbMuon(kinepid[1]);
	245	trackCopy->SetPidProbPion(kinepid[2]);
	246	trackCopy->SetPidProbKaon(kinepid[3]);
	247	trackCopy->SetPidProbProton(kinepid[4]);
	248
	249
	250	//Momentum
	251	double pxyz[3];
1777446b	252	// double rxyz[3];
76ce4b5b	253
	254	pxyz[0]=kinetrack->Px();
	255	pxyz[1]=kinetrack->Py();
	256	pxyz[2]=kinetrack->Pz();
	257
1777446b	258	// rxyz[0]=kinetrack->Vx();
	259	// rxyz[1]=kinetrack->Vy();
	260	// rxyz[2]=kinetrack->Vz();
76ce4b5b	261
	262	if (fRotateToEventPlane) {
	263	double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
	264	double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
	265
	266	pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
	267	pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
	268	}
	269
	270	AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
	271	if (v.Mag() < 0.0001) {
973a91f8	272	//cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
76ce4b5b	273	delete trackCopy;
	274	continue;
	275	}
	276
	277	trackCopy->SetP(v);//setting momentum
	278	trackCopy->SetPt(sqrt(pxyz[0]pxyz[0]+pxyz[1]pxyz[1]));
	279	const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
	280	const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
	281
	282	//label
973a91f8	283	trackCopy->SetLabel(i);
76ce4b5b	284
	285
	286	hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
973a91f8	287	//cout<<"Track added: "<<i<<endl;
76ce4b5b	288
	289	}
	290
	291	hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
4eac0b05	292	if (fEstEventMult == kGlobalCount)
	293	hbtEvent->SetNormalizedMult(tNormMult);
	294	else if(fEstEventMult == kVZERO)
	295	hbtEvent->SetNormalizedMult(tV0direction);
76ce4b5b	296	fCurEvent++;
76ce4b5b	297
973a91f8	298
	299	//V0 analysis code - no V0 finder for Kinematics, we can only check if it is primary and if it has at least 2 daughters.
	300
	301	for (int i=0;i<nofTracks;i++)
	302	{
	303	//do not take primaries
	304	if(!fStack->IsPhysicalPrimary(i)) {continue;}
	305	//getting next track
	306	TParticle *kinetrack= fStack->Particle(i);
	307	if (!kinetrack) continue;
	308
	309	if(kinetrack->GetPDG()->Charge()!=0) continue; //charge - neutral
	310	//if(kinetrack->GetPDG()->Stable()==1) continue; //particle is not stable
	311	if(kinetrack->GetDaughter(0)<1) continue; //has 1'st daughter
	312	if(kinetrack->GetDaughter(1)<1) continue; //has 2'nd daughter
	313
	314	//we want one positive, one negative particle. Or two neutral.
	315	// if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()>=0)
	316	// if((fStack->Particle(kinetrack->GetDaughter(1)))->GetPDG()->Charge()>0)
	317	// continue;
	318	// if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<=0)
	319	// if((fStack->Particle(kinetrack->GetDaughter(0)))->GetPDG()->Charge()<0)
	320	// continue;
	321
	322	if(kinetrack->Pt()<0.00001)
	323	continue;
	324
	325	AliFemtoV0* trackCopyV0 = new AliFemtoV0();
	326	CopyAODtoFemtoV0(kinetrack, trackCopyV0);
	327	hbtEvent->V0Collection()->push_back(trackCopyV0);
	328	//cout<<"Pushback v0 to v0collection"<<endl;
	329	}
	330
	331
	332	cout<<"Number of tracks: "<<realnofTracks<<endl;
	333
76ce4b5b	334	return hbtEvent;
	335	}
	336
	337	//___________________
	338	void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
	339	{
	340	// The chain loads the stack for us
	341	// You must provide the address where it can be found
	342	fStack = aStack;
	343	}
	344	//___________________
	345	void AliFemtoEventReaderKinematicsChain::SetGenEventHeader(AliGenEventHeader *aGenHeader)
	346	{
	347	// The chain loads the generator event header for us
	348	// You must provide the address where it can be found
	349	fGenHeader = aGenHeader;
	350	}
	351
	352	//__________________
	353	void AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
	354	{
	355	fRotateToEventPlane=dorotate;
	356	}
	357
4eac0b05	358	void AliFemtoEventReaderKinematicsChain::SetUseMultiplicity(EstEventMult aType)
	359	{
	360	fEstEventMult = aType;
	361	}
	362
76ce4b5b	363	Float_t AliFemtoEventReaderKinematicsChain::GetSigmaToVertex(double impact, double covar)
	364	{
	365	// Calculates the number of sigma to the vertex.
	366
	367	Float_t b[2];
	368	Float_t bRes[2];
	369	Float_t bCov[3];
	370
	371	b[0] = impact[0];
	372	b[1] = impact[1];
	373	bCov[0] = covar[0];
	374	bCov[1] = covar[1];
	375	bCov[2] = covar[2];
	376
	377	bRes[0] = TMath::Sqrt(bCov[0]);
	378	bRes[1] = TMath::Sqrt(bCov[2]);
	379
	380	// -----------------------------------
	381	// How to get to a n-sigma cut?
	382	//
	383	// The accumulated statistics from 0 to d is
	384	//
	385	// -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
	386	// -> 1 - Exp(-d*2) for a 2-dim gauss (dd = dxdx + dydy != n_sigma)
	387	//
	388	// It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)ErfInv(1 - Exp((-x*2)/2)
	389	// Can this be expressed in a different way?
	390
	391	if (bRes[0] == 0 \|\| bRes[1] ==0)
	392	return -1;
	393
	394	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	395
	396	// stupid rounding problem screws up everything:
	397	// if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
	398	if (TMath::Exp(-d * d / 2) < 1e-10)
	399	return 1000;
	400
	401	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	402	return d;
	403	}
973a91f8	404
	405
	406
	407	void AliFemtoEventReaderKinematicsChain::CopyAODtoFemtoV0(TParticle tv0, AliFemtoV0 tFemtoV0 )
	408	{
	409	tFemtoV0->SetEtaV0(tv0->Eta());
	410	tFemtoV0->SetEtaV0(tv0->Phi());
	411	tFemtoV0->SetptV0(tv0->Pt());
	412	tFemtoV0->SetptotV0(tv0->P());
	413
	414	tFemtoV0->SetmomV0X(tv0->Px());
	415	tFemtoV0->SetmomV0Y(tv0->Py());
	416	tFemtoV0->SetmomV0Z(tv0->Pz());
	417	AliFemtoThreeVector momv0(tv0->Px(),tv0->Py(),tv0->Pz());
	418	tFemtoV0->SetmomV0(momv0);
	419
	420
	421	TParticle *trackpos;
	422	TParticle *trackneg;
	423
	424	//daughters
	425	if(fStack->Particle(tv0->GetDaughter(0))->GetPDG()->Charge()>=0) //first positive, second negative
	426	{
	427	trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
	428	trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
	429	tFemtoV0->SetidPos(tv0->GetDaughter(0));
	430	tFemtoV0->SetidNeg(tv0->GetDaughter(1));
	431	}
	432	else //first negative, second positive
	433	{
	434	trackpos = (TParticle*)(fStack->Particle(tv0->GetDaughter(1)));
	435	trackneg = (TParticle*)(fStack->Particle(tv0->GetDaughter(0)));
	436	tFemtoV0->SetidPos(tv0->GetDaughter(1));
	437	tFemtoV0->SetidNeg(tv0->GetDaughter(0));
	438	}
	439
	440	tFemtoV0->SetEtaPos(trackpos->Eta());
	441	tFemtoV0->SetEtaNeg(trackneg->Eta());
	442
	443	tFemtoV0->SetptPos(trackpos->Pt());
	444	tFemtoV0->SetptNeg(trackneg->Pt());
	445
	446	tFemtoV0->SetptotPos(trackpos->P());
	447	tFemtoV0->SetptotNeg(trackneg->P());
	448
	449	tFemtoV0->SetmomPosX(trackpos->Px());
	450	tFemtoV0->SetmomPosY(trackpos->Py());
	451	tFemtoV0->SetmomPosZ(trackpos->Pz());
	452	AliFemtoThreeVector mompos(trackpos->Px(),trackpos->Py(),trackpos->Pz());
	453	tFemtoV0->SetmomPos(mompos);
	454
	455	tFemtoV0->SetmomNegX(trackneg->Px());
	456	tFemtoV0->SetmomNegY(trackneg->Py());
	457	tFemtoV0->SetmomNegZ(trackneg->Pz());
	458	AliFemtoThreeVector momneg(trackneg->Px(),trackneg->Py(),trackneg->Pz());
	459	tFemtoV0->SetmomNeg(momneg);
	460
	461
	462	tFemtoV0->SetmassLambda(tv0->GetMass());
	463	tFemtoV0->SetmassAntiLambda(tv0->GetMass());
	464	tFemtoV0->SetmassK0Short(tv0->GetMass());
	465
	466	tFemtoV0->SetYV0(tv0->Y());
	467
468	tFemtoV0->SetdecayVertexV0X(trackpos->Vx()); //vertex of the decay is set as the vertex of creation of daughters
469	tFemtoV0->SetdecayVertexV0Y(trackpos->Vy());
470	tFemtoV0->SetdecayVertexV0Z(trackpos->Vz());
471	AliFemtoThreeVector decayvertex(trackpos->Vx(),trackpos->Vy(),trackpos->Vz());
472	tFemtoV0->SetdecayVertexV0(decayvertex);
473
474	tFemtoV0->SetdcaV0Daughters(0);
475	tFemtoV0->SetCosPointingAngle(1);
476
477
478	tFemtoV0->SetStatusPos(1);
479	tFemtoV0->SetStatusNeg(1);
480
481
482	if(trackpos->GetPdgCode()==2212) //proton
483	{
484	tFemtoV0->SetPosNSigmaTPCK(1000);
485	tFemtoV0->SetPosNSigmaTPCPi(1000);
486	tFemtoV0->SetPosNSigmaTPCP(0);
487	}
488	if(trackneg->GetPdgCode()==-2212) //antiproton
489	{
490	tFemtoV0->SetNegNSigmaTPCK(1000);
491	tFemtoV0->SetNegNSigmaTPCPi(1000);
492	tFemtoV0->SetNegNSigmaTPCP(0);
493	}
494	if(trackpos->GetPdgCode()==211) //pion plus
495	{
496	tFemtoV0->SetPosNSigmaTPCK(1000);
497	tFemtoV0->SetPosNSigmaTPCPi(0);
498	tFemtoV0->SetPosNSigmaTPCP(1000);
499	}
500	if(trackneg->GetPdgCode()==-211) //pion minus
501	{
502	tFemtoV0->SetNegNSigmaTPCK(1000);
503	tFemtoV0->SetNegNSigmaTPCPi(0);
504	tFemtoV0->SetNegNSigmaTPCP(1000);
505	}
506	if(trackpos->GetPdgCode()==321) //K+
507	{
508	tFemtoV0->SetPosNSigmaTPCK(0);
509	tFemtoV0->SetPosNSigmaTPCPi(1000);
510	tFemtoV0->SetPosNSigmaTPCP(1000);
511	}
512	if(trackneg->GetPdgCode()==-321) //K-
513	{
514	tFemtoV0->SetNegNSigmaTPCK(0);
515	tFemtoV0->SetNegNSigmaTPCPi(1000);
516	tFemtoV0->SetNegNSigmaTPCP(1000);
517	}
518
519
520	}