[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronPair.cxx

/*************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                                                                       //
//  Dielectron Pair class. Internally it makes use of AliKFParticle.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////


#include <TDatabasePDG.h>
#include <AliVTrack.h>
#include <AliVVertex.h>
#include <AliPID.h>

#include "AliDielectronPair.h"

ClassImp(AliDielectronPair)

AliDielectronPair::AliDielectronPair() :
  fType(-1),
  fLabel(-1),
  fPair(),
  fD1(),
  fD2(),
  fRefD1(),
  fRefD2()
{
  //
  // Default Constructor
  //
  
}

//______________________________________________
AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
                                     AliVTrack * const particle2, Int_t pid2, Char_t type) :
  fType(type),
  fLabel(-1),
  fPair(),
  fD1(),
  fD2(),
  fRefD1(),
  fRefD2()
{
  //
  // Constructor with tracks
  //
  SetTracks(particle1, pid1, particle2, pid2);
}

//______________________________________________
AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
                                     const AliKFParticle * const particle2,
                                     AliVTrack * const refParticle1,
                                     AliVTrack * const refParticle2, Char_t type) :
  fType(type),
  fLabel(-1),
  fPair(),
  fD1(),
  fD2(),
  fRefD1(),
  fRefD2()
{
  //
  // Constructor with tracks
  //
  SetTracks(particle1, particle2,refParticle1,refParticle2);
}

//______________________________________________
AliDielectronPair::~AliDielectronPair()
{
  //
  // Default Destructor
  //
  
}

//______________________________________________
void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
                                  AliVTrack * const particle2, Int_t pid2)
{
  //
  // Sort particles by pt, first particle larget Pt
  // set AliKF daughters and pair
  // refParticle1 and 2 are the original tracks. In the case of track rotation
  // they are needed in the framework
  //
  fPair.Initialize();
  fD1.Initialize();
  fD2.Initialize();

  AliKFParticle kf1(*particle1,pid1);
  AliKFParticle kf2(*particle2,pid2);

  fPair.AddDaughter(kf1);
  fPair.AddDaughter(kf2);

  if (particle1->Pt()>particle2->Pt()){
    fRefD1 = particle1;
    fRefD2 = particle2;
    fD1+=kf1;
    fD2+=kf2;
  } else {
    fRefD1 = particle2;
    fRefD2 = particle1;
    fD1+=kf2;
    fD2+=kf1;
  }
}

//______________________________________________
void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
                                  const AliKFParticle * const particle2,
                                  AliVTrack * const refParticle1,
                                  AliVTrack * const refParticle2)
{
  //
  // Sort particles by pt, first particle larget Pt
  // set AliKF daughters and pair
  // refParticle1 and 2 are the original tracks. In the case of track rotation
  // they are needed in the framework
  //
  fPair.Initialize();
  fD1.Initialize();
  fD2.Initialize();
  
  AliKFParticle kf1(*particle1);
  AliKFParticle kf2(*particle2);
  
  fPair.AddDaughter(kf1);
  fPair.AddDaughter(kf2);
  
  if (kf1.GetPt()>kf2.GetPt()){
    fRefD1 = refParticle1;
    fRefD2 = refParticle2;
    fD1+=kf1;
    fD2+=kf2;
  } else {
    fRefD1 = refParticle2;
    fRefD2 = refParticle1;
    fD1+=kf2;
    fD2+=kf1;
  }
}

//______________________________________________
void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
{
  //
  // Calculate theta and phi in helicity and Collins-Soper coordinate frame
  //
  const Double_t kBeamEnergy   = 3500.;
  Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
  Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
  Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
  Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
  
//   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
//   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());

//   d1->PxPyPz(pxyz1);
//   d2->PxPyPz(pxyz2);
  
  TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
  TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
  
  // first & second daughter 4-mom
  TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
                       TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
  TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
                       TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
  // J/Psi 4-momentum vector
  TLorentzVector motherMom=p1Mom+p2Mom;
  
  // boost all the 4-mom vectors to the mother rest frame
  TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
  p1Mom.Boost(beta);
  p2Mom.Boost(beta);
  projMom.Boost(beta);
  targMom.Boost(beta);

    // x,y,z axes
  TVector3 zAxisHE = (motherMom.Vect()).Unit();
  TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
  TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
  TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
  TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
  
  // fill theta and phi
  if(fD1.GetQ()>0){
    thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
    thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
    phiHE   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
    phiCS   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
  } else {
    thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
    thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
    phiHE   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
    phiCS   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
  }
}

//______________________________________________
Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2, 
                                       const Bool_t isHE, const Bool_t isTheta)
{
  // The function calculates theta and phi in the mother rest frame with
  // respect to the helicity coordinate system and Collins-Soper coordinate system
  // TO DO: generalize for different decays (only J/Psi->e+e- now)

  // Laboratory frame 4-vectors:
  // projectile beam & target beam 4-mom
  // TODO: need to retrieve the beam energy from somewhere
  const Double_t kBeamEnergy   = 3500.;
  Double_t px1=d1->Px();
  Double_t py1=d1->Py();
  Double_t pz1=d1->Pz();
  Double_t px2=d2->Px();
  Double_t py2=d2->Py();
  Double_t pz2=d2->Pz();
  Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
  Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
  
  TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
  TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
  
  // first & second daughter 4-mom
  TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
  TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
  // J/Psi 4-momentum vector
  TLorentzVector motherMom=p1Mom+p2Mom;
  
  // boost all the 4-mom vectors to the mother rest frame
  TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
  p1Mom.Boost(beta);
  p2Mom.Boost(beta);
  projMom.Boost(beta);
  targMom.Boost(beta);
  
  // x,y,z axes 
  TVector3 zAxis;
  if(isHE) zAxis = (motherMom.Vect()).Unit();
  else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
  TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
  TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
  
  // return either theta or phi
  if(isTheta) {
    if(d1->Charge()>0)
      return zAxis.Dot((p1Mom.Vect()).Unit());
    else 
      return zAxis.Dot((p2Mom.Vect()).Unit());

  }
  else {
    if(d1->Charge()>0)
      return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
    else
      return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
  }
}

//______________________________________________
Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
  // The function calculates theta and phi in the mother rest frame with 
  // respect to the helicity coordinate system and Collins-Soper coordinate system
  // TO DO: generalize for different decays (only J/Psi->e+e- now)

  // Laboratory frame 4-vectors:
  // projectile beam & target beam 4-mom
  AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
  AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
  
  const Double_t kBeamEnergy   = 3500.;
  Double_t px1=d1->Px();
  Double_t py1=d1->Py();
  Double_t pz1=d1->Pz();
  Double_t px2=d2->Px();
  Double_t py2=d2->Py();
  Double_t pz2=d2->Pz();
  Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
  Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
  
  TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
  TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
  
  // first & second daughter 4-mom
  // first & second daughter 4-mom
  TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
  TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
  // J/Psi 4-momentum vector
  TLorentzVector motherMom=p1Mom+p2Mom;

  // boost all the 4-mom vectors to the mother rest frame
  TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
  p1Mom.Boost(beta);
  p2Mom.Boost(beta);
  projMom.Boost(beta);
  targMom.Boost(beta);

  // x,y,z axes 
  TVector3 zAxis;
  if(isHE) zAxis = (motherMom.Vect()).Unit();
  else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
  TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
  TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();

  // return either theta or phi
  if(isTheta) {
    if(fD1.GetQ()>0) 
      return zAxis.Dot((p1Mom.Vect()).Unit());
    else
      return zAxis.Dot((p2Mom.Vect()).Unit());
  }
  else {
    if(fD1.GetQ()>0)
      return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
    else
      return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
  }
}

//______________________________________________
Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
{
  //
  // Calculate the decay length in XY taking into account the primary vertex position
  //
  if(!vtx) return 0;
  return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt()  ;
}

//______________________________________________
Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
{
  //
  // Calculate the pseudo proper time
  //
  Double_t lxy=GetLXY(vtx);
  Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
  return psProperDecayLength;
}
Commit	Line	Data
b2a297fa	1	/*************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	///////////////////////////////////////////////////////////////////////////
	17	// //
	18	// Dielectron Pair class. Internally it makes use of AliKFParticle. //
	19	// //
	20	///////////////////////////////////////////////////////////////////////////
	21
	22
ba15fdfb	23	#include <TDatabasePDG.h>
	24	#include <AliVTrack.h>
	25	#include <AliVVertex.h>
	26	#include <AliPID.h>
	27
b2a297fa	28	#include "AliDielectronPair.h"
b2a297fa	29
	30	ClassImp(AliDielectronPair)
	31
	32	AliDielectronPair::AliDielectronPair() :
b2a297fa	33	fType(-1),
a655b716	34	fLabel(-1),
b2a297fa	35	fPair(),
572b0139	36	fD1(),
572b0139	37	fD2(),
b2a297fa	38	fRefD1(),
	39	fRefD2()
	40	{
	41	//
	42	// Default Constructor
	43	//
	44
	45	}
	46
	47	//______________________________________________
	48	AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
	49	AliVTrack * const particle2, Int_t pid2, Char_t type) :
b2a297fa	50	fType(type),
a655b716	51	fLabel(-1),
b2a297fa	52	fPair(),
572b0139	53	fD1(),
572b0139	54	fD2(),
b2a297fa	55	fRefD1(),
	56	fRefD2()
	57	{
	58	//
	59	// Constructor with tracks
	60	//
	61	SetTracks(particle1, pid1, particle2, pid2);
	62	}
	63
1201a1a9	64	//______________________________________________
	65	AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
	66	const AliKFParticle * const particle2,
	67	AliVTrack * const refParticle1,
	68	AliVTrack * const refParticle2, Char_t type) :
	69	fType(type),
	70	fLabel(-1),
	71	fPair(),
	72	fD1(),
	73	fD2(),
	74	fRefD1(),
	75	fRefD2()
	76	{
	77	//
	78	// Constructor with tracks
	79	//
	80	SetTracks(particle1, particle2,refParticle1,refParticle2);
	81	}
	82
b2a297fa	83	//______________________________________________
	84	AliDielectronPair::~AliDielectronPair()
	85	{
	86	//
	87	// Default Destructor
	88	//
	89
	90	}
	91
	92	//______________________________________________
	93	void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
	94	AliVTrack * const particle2, Int_t pid2)
	95	{
	96	//
572b0139	97	// Sort particles by pt, first particle larget Pt
572b0139	98	// set AliKF daughters and pair
1201a1a9	99	// refParticle1 and 2 are the original tracks. In the case of track rotation
1201a1a9	100	// they are needed in the framework
b2a297fa	101	//
b2a297fa	102	fPair.Initialize();
572b0139	103	fD1.Initialize();
572b0139	104	fD2.Initialize();
8df8e382	105
b2a297fa	106	AliKFParticle kf1(*particle1,pid1);
b2a297fa	107	AliKFParticle kf2(*particle2,pid2);
572b0139	108
b2a297fa	109	fPair.AddDaughter(kf1);
b2a297fa	110	fPair.AddDaughter(kf2);
8df8e382	111
a655b716	112	if (particle1->Pt()>particle2->Pt()){
	113	fRefD1 = particle1;
	114	fRefD2 = particle2;
572b0139	115	fD1+=kf1;
572b0139	116	fD2+=kf2;
a655b716	117	} else {
	118	fRefD1 = particle2;
	119	fRefD2 = particle1;
572b0139	120	fD1+=kf2;
572b0139	121	fD2+=kf1;
a655b716	122	}
b2a297fa	123	}
b2a297fa	124
1201a1a9	125	//______________________________________________
	126	void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
	127	const AliKFParticle * const particle2,
	128	AliVTrack * const refParticle1,
	129	AliVTrack * const refParticle2)
	130	{
	131	//
	132	// Sort particles by pt, first particle larget Pt
	133	// set AliKF daughters and pair
	134	// refParticle1 and 2 are the original tracks. In the case of track rotation
	135	// they are needed in the framework
	136	//
	137	fPair.Initialize();
	138	fD1.Initialize();
	139	fD2.Initialize();
	140
	141	AliKFParticle kf1(*particle1);
	142	AliKFParticle kf2(*particle2);
	143
	144	fPair.AddDaughter(kf1);
	145	fPair.AddDaughter(kf2);
	146
	147	if (kf1.GetPt()>kf2.GetPt()){
	148	fRefD1 = refParticle1;
	149	fRefD2 = refParticle2;
	150	fD1+=kf1;
	151	fD2+=kf2;
	152	} else {
	153	fRefD1 = refParticle2;
	154	fRefD2 = refParticle1;
	155	fD1+=kf2;
	156	fD2+=kf1;
	157	}
	158	}
	159
61d106d3	160	//______________________________________________
	161	void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
	162	{
	163	//
	164	// Calculate theta and phi in helicity and Collins-Soper coordinate frame
	165	//
	166	const Double_t kBeamEnergy = 3500.;
1201a1a9	167	Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
1201a1a9	168	Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
61d106d3	169	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	170	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	171
1201a1a9	172	// AliVParticle d1 = static_cast<AliVParticle>(fRefD1.GetObject());
1201a1a9	173	// AliVParticle d2 = static_cast<AliVParticle>(fRefD2.GetObject());
61d106d3	174
1201a1a9	175	// d1->PxPyPz(pxyz1);
1201a1a9	176	// d2->PxPyPz(pxyz2);
61d106d3	177
	178	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	179	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	180
	181	// first & second daughter 4-mom
	182	TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
	183	TMath::Sqrt(pxyz1[0]pxyz1[0]+pxyz1[1]pxyz1[1]+pxyz1[2]pxyz1[2]+eleMasseleMass));
	184	TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
	185	TMath::Sqrt(pxyz2[0]pxyz2[0]+pxyz2[1]pxyz2[1]+pxyz2[2]pxyz2[2]+eleMasseleMass));
	186	// J/Psi 4-momentum vector
	187	TLorentzVector motherMom=p1Mom+p2Mom;
	188
	189	// boost all the 4-mom vectors to the mother rest frame
	190	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	191	p1Mom.Boost(beta);
	192	p2Mom.Boost(beta);
	193	projMom.Boost(beta);
	194	targMom.Boost(beta);
	195
	196	// x,y,z axes
	197	TVector3 zAxisHE = (motherMom.Vect()).Unit();
	198	TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	199	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	200	TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
	201	TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
	202
	203	// fill theta and phi
1201a1a9	204	if(fD1.GetQ()>0){
61d106d3	205	thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
	206	thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
	207	phiHE = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
	208	phiCS = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
	209	} else {
	210	thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
	211	thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
	212	phiHE = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
	213	phiCS = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
	214	}
	215	}
	216
8df8e382	217	//______________________________________________
8df8e382	218	Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
61d106d3	219	const Bool_t isHE, const Bool_t isTheta)
	220	{
	221	// The function calculates theta and phi in the mother rest frame with
8df8e382	222	// respect to the helicity coordinate system and Collins-Soper coordinate system
	223	// TO DO: generalize for different decays (only J/Psi->e+e- now)
	224
	225	// Laboratory frame 4-vectors:
	226	// projectile beam & target beam 4-mom
61d106d3	227	// TODO: need to retrieve the beam energy from somewhere
	228	const Double_t kBeamEnergy = 3500.;
	229	Double_t px1=d1->Px();
	230	Double_t py1=d1->Py();
	231	Double_t pz1=d1->Pz();
	232	Double_t px2=d2->Px();
	233	Double_t py2=d2->Py();
	234	Double_t pz2=d2->Pz();
	235	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	236	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	237
	238	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	239	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
8df8e382	240
8df8e382	241	// first & second daughter 4-mom
61d106d3	242	TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1px1+py1py1+pz1pz1+eleMasseleMass));
61d106d3	243	TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2px2+py2py2+pz2pz2+eleMasseleMass));
8df8e382	244	// J/Psi 4-momentum vector
	245	TLorentzVector motherMom=p1Mom+p2Mom;
	246
	247	// boost all the 4-mom vectors to the mother rest frame
	248	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	249	p1Mom.Boost(beta);
	250	p2Mom.Boost(beta);
	251	projMom.Boost(beta);
	252	targMom.Boost(beta);
	253
	254	// x,y,z axes
	255	TVector3 zAxis;
	256	if(isHE) zAxis = (motherMom.Vect()).Unit();
	257	else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	258	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	259	TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
	260
	261	// return either theta or phi
	262	if(isTheta) {
	263	if(d1->Charge()>0)
	264	return zAxis.Dot((p1Mom.Vect()).Unit());
	265	else
	266	return zAxis.Dot((p2Mom.Vect()).Unit());
	267
	268	}
	269	else {
	270	if(d1->Charge()>0)
	271	return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
	272	else
	273	return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
	274	}
	275	}
	276
	277	//______________________________________________
	278	Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
	279	// The function calculates theta and phi in the mother rest frame with
	280	// respect to the helicity coordinate system and Collins-Soper coordinate system
	281	// TO DO: generalize for different decays (only J/Psi->e+e- now)
	282
	283	// Laboratory frame 4-vectors:
	284	// projectile beam & target beam 4-mom
45b2b1b8	285	AliVParticle d1 = static_cast<AliVParticle>(fRefD1.GetObject());
45b2b1b8	286	AliVParticle d2 = static_cast<AliVParticle>(fRefD2.GetObject());
61d106d3	287
	288	const Double_t kBeamEnergy = 3500.;
	289	Double_t px1=d1->Px();
	290	Double_t py1=d1->Py();
	291	Double_t pz1=d1->Pz();
	292	Double_t px2=d2->Px();
	293	Double_t py2=d2->Py();
	294	Double_t pz2=d2->Pz();
	295	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	296	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	297
	298	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	299	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	300
	301	// first & second daughter 4-mom
	302	// first & second daughter 4-mom
	303	TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1px1+py1py1+pz1pz1+eleMasseleMass));
	304	TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2px2+py2py2+pz2pz2+eleMasseleMass));
8df8e382	305	// J/Psi 4-momentum vector
	306	TLorentzVector motherMom=p1Mom+p2Mom;
	307
	308	// boost all the 4-mom vectors to the mother rest frame
	309	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	310	p1Mom.Boost(beta);
	311	p2Mom.Boost(beta);
	312	projMom.Boost(beta);
	313	targMom.Boost(beta);
	314
	315	// x,y,z axes
	316	TVector3 zAxis;
	317	if(isHE) zAxis = (motherMom.Vect()).Unit();
	318	else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	319	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	320	TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
	321
	322	// return either theta or phi
	323	if(isTheta) {
	324	if(fD1.GetQ()>0)
	325	return zAxis.Dot((p1Mom.Vect()).Unit());
	326	else
	327	return zAxis.Dot((p2Mom.Vect()).Unit());
	328	}
	329	else {
	330	if(fD1.GetQ()>0)
	331	return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
	332	else
	333	return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
	334	}
	335	}
ba15fdfb	336
	337	//______________________________________________
	338	Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
	339	{
	340	//
	341	// Calculate the decay length in XY taking into account the primary vertex position
	342	//
	343	if(!vtx) return 0;
	344	return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ;
	345	}
	346
	347	//______________________________________________
	348	Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
	349	{
	350	//
	351	// Calculate the pseudo proper time
	352	//
	353	Double_t lxy=GetLXY(vtx);
	354	Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
	355	return psProperDecayLength;
	356	}