[u/mrichter/AliRoot.git] / PWG3 / 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 "AliDielectronPair.h"
#include "AliVTrack.h"
#include "AliPID.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()
{
  //
  // 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
  //
  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::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]={0,0,0};
  Double_t pxyz2[3]={0,0,0};
  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(d1->Charge()>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 = dynamic_cast<AliVParticle*>(fRefD1.GetObject());
  AliVParticle *d2 = dynamic_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));
  }
}
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
	23	#include "AliDielectronPair.h"
	24	#include "AliVTrack.h"
8df8e382	25	#include "AliPID.h"
b2a297fa	26
	27	ClassImp(AliDielectronPair)
	28
	29	AliDielectronPair::AliDielectronPair() :
b2a297fa	30	fType(-1),
a655b716	31	fLabel(-1),
b2a297fa	32	fPair(),
572b0139	33	fD1(),
572b0139	34	fD2(),
b2a297fa	35	fRefD1(),
	36	fRefD2()
	37	{
	38	//
	39	// Default Constructor
	40	//
	41
	42	}
	43
	44	//______________________________________________
	45	AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
	46	AliVTrack * const particle2, Int_t pid2, Char_t type) :
b2a297fa	47	fType(type),
a655b716	48	fLabel(-1),
b2a297fa	49	fPair(),
572b0139	50	fD1(),
572b0139	51	fD2(),
b2a297fa	52	fRefD1(),
	53	fRefD2()
	54	{
	55	//
	56	// Constructor with tracks
	57	//
	58	SetTracks(particle1, pid1, particle2, pid2);
	59	}
	60
	61	//______________________________________________
	62	AliDielectronPair::~AliDielectronPair()
	63	{
	64	//
	65	// Default Destructor
	66	//
	67
	68	}
	69
	70	//______________________________________________
	71	void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
	72	AliVTrack * const particle2, Int_t pid2)
	73	{
	74	//
572b0139	75	// Sort particles by pt, first particle larget Pt
572b0139	76	// set AliKF daughters and pair
b2a297fa	77	//
b2a297fa	78	fPair.Initialize();
572b0139	79	fD1.Initialize();
572b0139	80	fD2.Initialize();
8df8e382	81
b2a297fa	82	AliKFParticle kf1(*particle1,pid1);
b2a297fa	83	AliKFParticle kf2(*particle2,pid2);
572b0139	84
b2a297fa	85	fPair.AddDaughter(kf1);
b2a297fa	86	fPair.AddDaughter(kf2);
8df8e382	87
a655b716	88	if (particle1->Pt()>particle2->Pt()){
	89	fRefD1 = particle1;
	90	fRefD2 = particle2;
572b0139	91	fD1+=kf1;
572b0139	92	fD2+=kf2;
a655b716	93	} else {
	94	fRefD1 = particle2;
	95	fRefD2 = particle1;
572b0139	96	fD1+=kf2;
572b0139	97	fD2+=kf1;
a655b716	98	}
b2a297fa	99	}
b2a297fa	100
61d106d3	101	//______________________________________________
	102	void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
	103	{
	104	//
	105	// Calculate theta and phi in helicity and Collins-Soper coordinate frame
	106	//
	107	const Double_t kBeamEnergy = 3500.;
	108	Double_t pxyz1[3]={0,0,0};
	109	Double_t pxyz2[3]={0,0,0};
	110	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	111	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	112
	113	AliVParticle d1 = static_cast<AliVParticle>(fRefD1.GetObject());
	114	AliVParticle d2 = static_cast<AliVParticle>(fRefD2.GetObject());
	115
	116	d1->PxPyPz(pxyz1);
	117	d2->PxPyPz(pxyz2);
	118
	119	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	120	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	121
	122	// first & second daughter 4-mom
	123	TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
	124	TMath::Sqrt(pxyz1[0]pxyz1[0]+pxyz1[1]pxyz1[1]+pxyz1[2]pxyz1[2]+eleMasseleMass));
	125	TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
	126	TMath::Sqrt(pxyz2[0]pxyz2[0]+pxyz2[1]pxyz2[1]+pxyz2[2]pxyz2[2]+eleMasseleMass));
	127	// J/Psi 4-momentum vector
	128	TLorentzVector motherMom=p1Mom+p2Mom;
	129
	130	// boost all the 4-mom vectors to the mother rest frame
	131	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	132	p1Mom.Boost(beta);
	133	p2Mom.Boost(beta);
	134	projMom.Boost(beta);
	135	targMom.Boost(beta);
	136
	137	// x,y,z axes
	138	TVector3 zAxisHE = (motherMom.Vect()).Unit();
	139	TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	140	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	141	TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
	142	TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
	143
	144	// fill theta and phi
	145	if(d1->Charge()>0){
	146	thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
	147	thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
	148	phiHE = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
	149	phiCS = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
	150	} else {
	151	thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
	152	thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
	153	phiHE = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
	154	phiCS = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
	155	}
	156	}
	157
8df8e382	158	//______________________________________________
8df8e382	159	Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
61d106d3	160	const Bool_t isHE, const Bool_t isTheta)
	161	{
	162	// The function calculates theta and phi in the mother rest frame with
8df8e382	163	// respect to the helicity coordinate system and Collins-Soper coordinate system
	164	// TO DO: generalize for different decays (only J/Psi->e+e- now)
	165
	166	// Laboratory frame 4-vectors:
	167	// projectile beam & target beam 4-mom
61d106d3	168	// TODO: need to retrieve the beam energy from somewhere
	169	const Double_t kBeamEnergy = 3500.;
	170	Double_t px1=d1->Px();
	171	Double_t py1=d1->Py();
	172	Double_t pz1=d1->Pz();
	173	Double_t px2=d2->Px();
	174	Double_t py2=d2->Py();
	175	Double_t pz2=d2->Pz();
	176	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	177	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	178
	179	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	180	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
8df8e382	181
8df8e382	182	// first & second daughter 4-mom
61d106d3	183	TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1px1+py1py1+pz1pz1+eleMasseleMass));
61d106d3	184	TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2px2+py2py2+pz2pz2+eleMasseleMass));
8df8e382	185	// J/Psi 4-momentum vector
	186	TLorentzVector motherMom=p1Mom+p2Mom;
	187
	188	// boost all the 4-mom vectors to the mother rest frame
	189	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	190	p1Mom.Boost(beta);
	191	p2Mom.Boost(beta);
	192	projMom.Boost(beta);
	193	targMom.Boost(beta);
	194
	195	// x,y,z axes
	196	TVector3 zAxis;
	197	if(isHE) zAxis = (motherMom.Vect()).Unit();
	198	else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	199	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	200	TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
	201
	202	// return either theta or phi
	203	if(isTheta) {
	204	if(d1->Charge()>0)
	205	return zAxis.Dot((p1Mom.Vect()).Unit());
	206	else
	207	return zAxis.Dot((p2Mom.Vect()).Unit());
	208
	209	}
	210	else {
	211	if(d1->Charge()>0)
	212	return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
	213	else
	214	return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
	215	}
	216	}
	217
	218	//______________________________________________
	219	Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
	220	// The function calculates theta and phi in the mother rest frame with
	221	// respect to the helicity coordinate system and Collins-Soper coordinate system
	222	// TO DO: generalize for different decays (only J/Psi->e+e- now)
	223
	224	// Laboratory frame 4-vectors:
	225	// projectile beam & target beam 4-mom
8df8e382	226	AliVParticle d1 = dynamic_cast<AliVParticle>(fRefD1.GetObject());
8df8e382	227	AliVParticle d2 = dynamic_cast<AliVParticle>(fRefD2.GetObject());
61d106d3	228
	229	const Double_t kBeamEnergy = 3500.;
	230	Double_t px1=d1->Px();
	231	Double_t py1=d1->Py();
	232	Double_t pz1=d1->Pz();
	233	Double_t px2=d2->Px();
	234	Double_t py2=d2->Py();
	235	Double_t pz2=d2->Pz();
	236	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	237	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	238
	239	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	240	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	241
	242	// first & second daughter 4-mom
	243	// first & second daughter 4-mom
	244	TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1px1+py1py1+pz1pz1+eleMasseleMass));
	245	TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2px2+py2py2+pz2pz2+eleMasseleMass));
8df8e382	246	// J/Psi 4-momentum vector
	247	TLorentzVector motherMom=p1Mom+p2Mom;
	248
	249	// boost all the 4-mom vectors to the mother rest frame
	250	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	251	p1Mom.Boost(beta);
	252	p2Mom.Boost(beta);
	253	projMom.Boost(beta);
	254	targMom.Boost(beta);
	255
	256	// x,y,z axes
	257	TVector3 zAxis;
	258	if(isHE) zAxis = (motherMom.Vect()).Unit();
	259	else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	260	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	261	TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
	262
	263	// return either theta or phi
	264	if(isTheta) {
	265	if(fD1.GetQ()>0)
	266	return zAxis.Dot((p1Mom.Vect()).Unit());
	267	else
	268	return zAxis.Dot((p2Mom.Vect()).Unit());
	269	}
	270	else {
	271	if(fD1.GetQ()>0)
	272	return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
	273	else
	274	return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
	275	}
	276	}