[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 <AliExternalTrackParam.h>

#include "AliDielectronPair.h"

ClassImp(AliDielectronPair)

AliDielectronPair::AliDielectronPair() :
  fType(-1),
  fLabel(-1),
  fPdgCode(0),
  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),
  fPdgCode(0),
  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),
  fPdgCode(0),
  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::PsiPair(Double_t MagField) const
{
  //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
  //to ID conversions. Adapted from AliTRDv0Info class
  Double_t x, y, z;
  x = fPair.GetX();
  y = fPair.GetY();
  z = fPair.GetZ();

  Double_t m1[3] = {0,0,0};
  Double_t m2[3] = {0,0,0};

  m1[0] = fD1.GetPx();
  m1[1] = fD1.GetPy();
  m1[2] = fD1.GetPz();  

  m2[0] = fD2.GetPx();
  m2[1] = fD2.GetPy();
  m2[2] = fD2.GetPz();

  Double_t deltat = 1.;
  deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
	TMath::ATan(m1[2]/(TMath::Sqrt(m1[0]*m1[0] + m1[1]*m1[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis

  Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated

  Double_t mom1Prop[3];
  Double_t mom2Prop[3];

  AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject());
  AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject());

  AliExternalTrackParam nt(*d1), pt(*d2);

  Double_t fPsiPair = 4.;
  if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
	fPsiPair =  -5.;
  if(pt.PropagateTo(radiussum,MagField) == 0)
	fPsiPair = -5.;
  pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
  nt.GetPxPyPz(mom2Prop);


  Double_t pEle =
	TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
  Double_t pPos =
	TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val

  Double_t scalarproduct =
	mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit

  Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks

  fPsiPair =  TMath::Abs(TMath::ASin(deltat/chipair));

  return fPsiPair;

}

//______________________________________________
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>
236e1bda	27	#include <AliExternalTrackParam.h>
ba15fdfb	28
b2a297fa	29	#include "AliDielectronPair.h"
b2a297fa	30
	31	ClassImp(AliDielectronPair)
	32
	33	AliDielectronPair::AliDielectronPair() :
b2a297fa	34	fType(-1),
a655b716	35	fLabel(-1),
e4339752	36	fPdgCode(0),
b2a297fa	37	fPair(),
572b0139	38	fD1(),
572b0139	39	fD2(),
b2a297fa	40	fRefD1(),
	41	fRefD2()
	42	{
	43	//
	44	// Default Constructor
	45	//
	46
	47	}
	48
	49	//______________________________________________
	50	AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
	51	AliVTrack * const particle2, Int_t pid2, Char_t type) :
b2a297fa	52	fType(type),
a655b716	53	fLabel(-1),
e4339752	54	fPdgCode(0),
b2a297fa	55	fPair(),
572b0139	56	fD1(),
572b0139	57	fD2(),
b2a297fa	58	fRefD1(),
	59	fRefD2()
	60	{
	61	//
	62	// Constructor with tracks
	63	//
	64	SetTracks(particle1, pid1, particle2, pid2);
	65	}
	66
1201a1a9	67	//______________________________________________
	68	AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
	69	const AliKFParticle * const particle2,
	70	AliVTrack * const refParticle1,
	71	AliVTrack * const refParticle2, Char_t type) :
	72	fType(type),
	73	fLabel(-1),
e4339752	74	fPdgCode(0),
1201a1a9	75	fPair(),
	76	fD1(),
	77	fD2(),
	78	fRefD1(),
	79	fRefD2()
	80	{
	81	//
	82	// Constructor with tracks
	83	//
	84	SetTracks(particle1, particle2,refParticle1,refParticle2);
	85	}
	86
b2a297fa	87	//______________________________________________
	88	AliDielectronPair::~AliDielectronPair()
	89	{
	90	//
	91	// Default Destructor
	92	//
	93
	94	}
	95
	96	//______________________________________________
	97	void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
	98	AliVTrack * const particle2, Int_t pid2)
	99	{
	100	//
572b0139	101	// Sort particles by pt, first particle larget Pt
572b0139	102	// set AliKF daughters and pair
1201a1a9	103	// refParticle1 and 2 are the original tracks. In the case of track rotation
1201a1a9	104	// they are needed in the framework
b2a297fa	105	//
b2a297fa	106	fPair.Initialize();
572b0139	107	fD1.Initialize();
572b0139	108	fD2.Initialize();
8df8e382	109
b2a297fa	110	AliKFParticle kf1(*particle1,pid1);
b2a297fa	111	AliKFParticle kf2(*particle2,pid2);
572b0139	112
b2a297fa	113	fPair.AddDaughter(kf1);
b2a297fa	114	fPair.AddDaughter(kf2);
8df8e382	115
a655b716	116	if (particle1->Pt()>particle2->Pt()){
	117	fRefD1 = particle1;
	118	fRefD2 = particle2;
572b0139	119	fD1+=kf1;
572b0139	120	fD2+=kf2;
a655b716	121	} else {
	122	fRefD1 = particle2;
	123	fRefD2 = particle1;
572b0139	124	fD1+=kf2;
572b0139	125	fD2+=kf1;
a655b716	126	}
b2a297fa	127	}
b2a297fa	128
1201a1a9	129	//______________________________________________
	130	void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
	131	const AliKFParticle * const particle2,
	132	AliVTrack * const refParticle1,
	133	AliVTrack * const refParticle2)
	134	{
	135	//
	136	// Sort particles by pt, first particle larget Pt
	137	// set AliKF daughters and pair
	138	// refParticle1 and 2 are the original tracks. In the case of track rotation
	139	// they are needed in the framework
	140	//
	141	fPair.Initialize();
	142	fD1.Initialize();
	143	fD2.Initialize();
	144
	145	AliKFParticle kf1(*particle1);
	146	AliKFParticle kf2(*particle2);
	147
	148	fPair.AddDaughter(kf1);
	149	fPair.AddDaughter(kf2);
	150
	151	if (kf1.GetPt()>kf2.GetPt()){
	152	fRefD1 = refParticle1;
	153	fRefD2 = refParticle2;
	154	fD1+=kf1;
	155	fD2+=kf2;
	156	} else {
	157	fRefD1 = refParticle2;
	158	fRefD2 = refParticle1;
	159	fD1+=kf2;
	160	fD2+=kf1;
	161	}
	162	}
	163
61d106d3	164	//______________________________________________
	165	void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
	166	{
	167	//
	168	// Calculate theta and phi in helicity and Collins-Soper coordinate frame
	169	//
	170	const Double_t kBeamEnergy = 3500.;
1201a1a9	171	Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
1201a1a9	172	Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
61d106d3	173	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	174	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	175
1201a1a9	176	// AliVParticle d1 = static_cast<AliVParticle>(fRefD1.GetObject());
1201a1a9	177	// AliVParticle d2 = static_cast<AliVParticle>(fRefD2.GetObject());
61d106d3	178
1201a1a9	179	// d1->PxPyPz(pxyz1);
1201a1a9	180	// d2->PxPyPz(pxyz2);
61d106d3	181
	182	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	183	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	184
	185	// first & second daughter 4-mom
	186	TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
	187	TMath::Sqrt(pxyz1[0]pxyz1[0]+pxyz1[1]pxyz1[1]+pxyz1[2]pxyz1[2]+eleMasseleMass));
	188	TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
	189	TMath::Sqrt(pxyz2[0]pxyz2[0]+pxyz2[1]pxyz2[1]+pxyz2[2]pxyz2[2]+eleMasseleMass));
	190	// J/Psi 4-momentum vector
	191	TLorentzVector motherMom=p1Mom+p2Mom;
	192
	193	// boost all the 4-mom vectors to the mother rest frame
	194	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	195	p1Mom.Boost(beta);
	196	p2Mom.Boost(beta);
	197	projMom.Boost(beta);
	198	targMom.Boost(beta);
	199
	200	// x,y,z axes
	201	TVector3 zAxisHE = (motherMom.Vect()).Unit();
	202	TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	203	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	204	TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
	205	TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
	206
	207	// fill theta and phi
1201a1a9	208	if(fD1.GetQ()>0){
61d106d3	209	thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
	210	thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
	211	phiHE = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
	212	phiCS = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
	213	} else {
	214	thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
	215	thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
	216	phiHE = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
	217	phiCS = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
	218	}
	219	}
	220
236e1bda	221	//______________________________________________
	222	Double_t AliDielectronPair::PsiPair(Double_t MagField) const
	223	{
	224	//Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
	225	//to ID conversions. Adapted from AliTRDv0Info class
	226	Double_t x, y, z;
	227	x = fPair.GetX();
	228	y = fPair.GetY();
	229	z = fPair.GetZ();
	230
	231	Double_t m1[3] = {0,0,0};
	232	Double_t m2[3] = {0,0,0};
	233
	234	m1[0] = fD1.GetPx();
	235	m1[1] = fD1.GetPy();
	236	m1[2] = fD1.GetPz();
	237
	238	m2[0] = fD2.GetPx();
	239	m2[1] = fD2.GetPy();
	240	m2[2] = fD2.GetPz();
	241
	242	Double_t deltat = 1.;
	243	deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]m2[0] + m2[1]m2[1])+1.e-13))-
	244	TMath::ATan(m1[2]/(TMath::Sqrt(m1[0]m1[0] + m1[1]m1[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
	245
	246	Double_t radiussum = TMath::Sqrt(xx + yy) + 50;//radius to which tracks shall be propagated
	247
	248	Double_t mom1Prop[3];
	249	Double_t mom2Prop[3];
	250
	251	AliExternalTrackParam d1 = static_cast<AliExternalTrackParam>(fRefD1.GetObject());
	252	AliExternalTrackParam d2 = static_cast<AliExternalTrackParam>(fRefD2.GetObject());
	253
	254	AliExternalTrackParam nt(d1), pt(d2);
236e1bda	255
	256	Double_t fPsiPair = 4.;
	257	if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
	258	fPsiPair = -5.;
	259	if(pt.PropagateTo(radiussum,MagField) == 0)
	260	fPsiPair = -5.;
	261	pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
	262	nt.GetPxPyPz(mom2Prop);
	263
	264
	265
	266	Double_t pEle =
	267	TMath::Sqrt(mom2Prop[0]mom2Prop[0]+mom2Prop[1]mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
	268	Double_t pPos =
	269	TMath::Sqrt(mom1Prop[0]mom1Prop[0]+mom1Prop[1]mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
	270
	271	Double_t scalarproduct =
	272	mom1Prop[0]mom2Prop[0]+mom1Prop[1]mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
	273
	274	Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
	275
	276	fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
	277
	278	return fPsiPair;
	279
	280	}
	281
8df8e382	282	//______________________________________________
8df8e382	283	Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
61d106d3	284	const Bool_t isHE, const Bool_t isTheta)
	285	{
	286	// The function calculates theta and phi in the mother rest frame with
8df8e382	287	// respect to the helicity coordinate system and Collins-Soper coordinate system
	288	// TO DO: generalize for different decays (only J/Psi->e+e- now)
	289
	290	// Laboratory frame 4-vectors:
	291	// projectile beam & target beam 4-mom
61d106d3	292	// TODO: need to retrieve the beam energy from somewhere
	293	const Double_t kBeamEnergy = 3500.;
	294	Double_t px1=d1->Px();
	295	Double_t py1=d1->Py();
	296	Double_t pz1=d1->Pz();
	297	Double_t px2=d2->Px();
	298	Double_t py2=d2->Py();
	299	Double_t pz2=d2->Pz();
	300	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	301	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	302
	303	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	304	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
8df8e382	305
8df8e382	306	// first & second daughter 4-mom
61d106d3	307	TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1px1+py1py1+pz1pz1+eleMasseleMass));
61d106d3	308	TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2px2+py2py2+pz2pz2+eleMasseleMass));
8df8e382	309	// J/Psi 4-momentum vector
	310	TLorentzVector motherMom=p1Mom+p2Mom;
	311
	312	// boost all the 4-mom vectors to the mother rest frame
	313	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	314	p1Mom.Boost(beta);
	315	p2Mom.Boost(beta);
	316	projMom.Boost(beta);
	317	targMom.Boost(beta);
	318
	319	// x,y,z axes
	320	TVector3 zAxis;
	321	if(isHE) zAxis = (motherMom.Vect()).Unit();
	322	else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	323	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	324	TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
	325
	326	// return either theta or phi
	327	if(isTheta) {
	328	if(d1->Charge()>0)
	329	return zAxis.Dot((p1Mom.Vect()).Unit());
	330	else
	331	return zAxis.Dot((p2Mom.Vect()).Unit());
	332
	333	}
	334	else {
	335	if(d1->Charge()>0)
	336	return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
	337	else
	338	return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
	339	}
	340	}
	341
	342	//______________________________________________
	343	Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
	344	// The function calculates theta and phi in the mother rest frame with
	345	// respect to the helicity coordinate system and Collins-Soper coordinate system
	346	// TO DO: generalize for different decays (only J/Psi->e+e- now)
	347
	348	// Laboratory frame 4-vectors:
	349	// projectile beam & target beam 4-mom
45b2b1b8	350	AliVParticle d1 = static_cast<AliVParticle>(fRefD1.GetObject());
45b2b1b8	351	AliVParticle d2 = static_cast<AliVParticle>(fRefD2.GetObject());
61d106d3	352
	353	const Double_t kBeamEnergy = 3500.;
	354	Double_t px1=d1->Px();
	355	Double_t py1=d1->Py();
	356	Double_t pz1=d1->Pz();
	357	Double_t px2=d2->Px();
	358	Double_t py2=d2->Py();
	359	Double_t pz2=d2->Pz();
	360	Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
	361	Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
	362
	363	TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	364	TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergykBeamEnergy+proMassproMass));
	365
	366	// first & second daughter 4-mom
	367	// first & second daughter 4-mom
	368	TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1px1+py1py1+pz1pz1+eleMasseleMass));
	369	TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2px2+py2py2+pz2pz2+eleMasseleMass));
8df8e382	370	// J/Psi 4-momentum vector
	371	TLorentzVector motherMom=p1Mom+p2Mom;
	372
	373	// boost all the 4-mom vectors to the mother rest frame
	374	TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
	375	p1Mom.Boost(beta);
	376	p2Mom.Boost(beta);
	377	projMom.Boost(beta);
	378	targMom.Boost(beta);
	379
	380	// x,y,z axes
	381	TVector3 zAxis;
	382	if(isHE) zAxis = (motherMom.Vect()).Unit();
	383	else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
	384	TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
	385	TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
	386
	387	// return either theta or phi
	388	if(isTheta) {
	389	if(fD1.GetQ()>0)
	390	return zAxis.Dot((p1Mom.Vect()).Unit());
	391	else
	392	return zAxis.Dot((p2Mom.Vect()).Unit());
	393	}
	394	else {
	395	if(fD1.GetQ()>0)
	396	return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
	397	else
	398	return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
	399	}
	400	}
ba15fdfb	401
5720c765	402	// //______________________________________________
	403	// Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
	404	// {
	405	// //
	406	// // Calculate the decay length in XY taking into account the primary vertex position
	407	// //
	408	// if(!vtx) return 0;
	409	// return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt() ;
	410	// }
ba15fdfb	411
5720c765	412	// //______________________________________________
	413	// Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
	414	// {
	415	// //
	416	// // Calculate the pseudo proper time
	417	// //
	418	// Double_t lxy=GetLXY(vtx);
	419	// Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
	420	// return psProperDecayLength;
	421	// }