[u/mrichter/AliRoot.git] / STEER / AOD / AliAODDimuon.cxx

// AliAODDimuon: a class for AODs for the MUON Arm of the ALICE Experiment
// Author: P. Cortese, Universita' del Piemonte Orientale in Alessandria and
// INFN of Torino - Italy
//
// The class defines a dimuon pair object from two AliAODTrack objects.
// AliAODDimuon objects are supposed to be added to the AliAODEvent structure
// during analysis. They would then allow to calculate the dimuon-related
// kinematic variables with a minimal disk occupancy.
// The payload of the class has been reduced to two pointers to the two
// tracks. An instance of this class has also to be added to the AliAODEvent 
// structure to provide additional information that is specific to MUON and 
// therefore has not been included into the AOD header.
// Two transient data members are not stored on file as they can be recomputed
// at runtime.
//

#include "AliAODDimuon.h"
#include "TLorentzVector.h"
#define AliAODDimuon_CXX

ClassImp(AliAODDimuon)

//______________________________________________________________________________
AliAODDimuon::AliAODDimuon():AliVParticle(),fP(0),fMProton(0.93827231)
{
  // default constructor
  fMu[0]=0;
  fMu[1]=0;
}

//______________________________________________________________________________
AliAODDimuon::AliAODDimuon(const AliAODDimuon& dimu):
  AliVParticle(dimu),
  fP(0),fMProton(0.93827231)
{
  // copy constructor
  fMu[0]=dimu.Mu(0);
  fMu[1]=dimu.Mu(1);
}

//______________________________________________________________________________
AliAODDimuon &AliAODDimuon::operator=(const AliAODDimuon& dimu)
{
  // assignment operator
  if(&dimu != this){
    delete fP;
    fP=0;
    fMProton=0.93827231;
    fMu[0]=dimu.Mu(0);
    fMu[1]=dimu.Mu(1);
  }
  return *this;
}

//______________________________________________________________________________
AliAODDimuon::AliAODDimuon(TObject *mu0, TObject *mu1):
  fP(0),fMProton(0.93827231)
{
  // Creates a dimuon pair from two tracks

  fMu[0]=mu0;
  fMu[1]=mu1;
}

//______________________________________________________________________________
AliAODDimuon::~AliAODDimuon()
{
  // destructor
  delete fP;
}

//______________________________________________________________________________
Double_t AliAODDimuon::Px() const {
  // Px of the dimuon
  if(CheckPointers())return -999999999;
  return ((AliAODTrack*)fMu[0].GetObject())->Px()+
         ((AliAODTrack*)fMu[1].GetObject())->Px();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Py() const {
  // Py of the dimuon
  if(CheckPointers())return -999999999;
  return ((AliAODTrack*)fMu[0].GetObject())->Py()+
         ((AliAODTrack*)fMu[1].GetObject())->Py();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Pz() const {
  // Pz of the dimuon
  if(CheckPointers())return -999999999;
  return ((AliAODTrack*)fMu[0].GetObject())->Pz()+
         ((AliAODTrack*)fMu[1].GetObject())->Pz();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Pt() const {
  // Pt of the dimuon
  if(CheckPointers())return -999999999;
  Double_t px=Px();
  Double_t py=Py();
  return TMath::Sqrt(px*px+py*py);
}

//______________________________________________________________________________
Double_t AliAODDimuon::E() const 
{
    // Dimuon energy
  
  if(CheckPointers())return -999999999;

  return ((AliAODTrack*)fMu[0].GetObject())->E()+ ((AliAODTrack*)fMu[1].GetObject())->E();
}

//______________________________________________________________________________
Double_t AliAODDimuon::P() const {
  // Dimuon momentum
  if(CheckPointers())return -999999999;
  return TLV()->P();
}

//______________________________________________________________________________
Double_t AliAODDimuon::M() const {
  // Dimuon invariant mass
  if(CheckPointers())return -999999999;
  return TLV()->M();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Eta() const {
  // Dimuon pseudorapidity
  if(CheckPointers())return -999999999;
  return TLV()->Eta();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Phi() const {
  // Dimuon asimuthal angle
  if(CheckPointers())return -999999999;
  return TLV()->Phi();
}


//______________________________________________________________________________
Double_t AliAODDimuon::Theta() const {
  // Dimuon polar angle
  if(CheckPointers())return -999999999;
  return TLV()->Theta();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Y() const {
  // Dimuon rapidity
  if(CheckPointers())return -999999999;
  return TLV()->Rapidity();
}

//______________________________________________________________________________
Short_t AliAODDimuon::Charge() const {
  // Dimuon charge
  if(CheckPointers())return -999;
  return ((AliAODTrack*)fMu[0].GetObject())->Charge()+
         ((AliAODTrack*)fMu[1].GetObject())->Charge();
}

//______________________________________________________________________________
Int_t AliAODDimuon::CheckPointers() const{
  // Checks if the track pointers have been initialized
  if(fMu[0]==0||fMu[1]==0){
    printf("Dimuon not initialized\n");
    return -999;
  }
  if((fMu[0].GetObject())==0||(fMu[1].GetObject())==0){
    printf("Can not get objects\n");
    return -999;
  }
  return 0;
}

//______________________________________________________________________________
void AliAODDimuon::SetMu(Int_t imu, AliAODTrack *mu){
  // Assign a track pointer
  if (imu==0||imu==1){
    delete fP;
    fP=0;
    fMu[imu]=mu;
  }
}

//______________________________________________________________________________
void AliAODDimuon::SetMuons(AliAODTrack *mu0, AliAODTrack *mu1){
  // Assign the track pointers
  fMu[0]=mu0;
  fMu[1]=mu1;
  delete fP;
  fP=0;
}

//______________________________________________________________________________
Double_t AliAODDimuon::XF() {
  // Dimuon Feynman x

  if(CheckPointers())return -999999999;

  //Double_t ebeam=((AliAODEventInfo*)fEi.GetObject())->EBeam();
  Double_t ebeam = 3500.; // temporary
  if(ebeam<=0){
    printf("AliAODDimuon::xf: can not compute xf with EBeam=%f\n",ebeam);
    return -999999999;
  }
  Double_t mDimu=M();
  Double_t pMax=TMath::Sqrt(ebeam*ebeam-mDimu*mDimu);
  return Pz()/pMax;
}

//______________________________________________________________________________
Double_t AliAODDimuon::CostCS(){
  // Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
  if(CheckPointers())return -999999999;
  Double_t ebeam=3500.;  //temporary
  if(ebeam<=0){
    printf("Can not compute costCS with EBeam=%f\n",ebeam);
    return -999999999;
  }
  Double_t mp=fMProton;
  Double_t pbeam=TMath::Sqrt(ebeam*ebeam-mp*mp);
  Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
  Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
  Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
  Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
  Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
  Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
  Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
  Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
  Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
  Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();

  // Fill the Lorentz vector for projectile and target
  // For the moment we do not consider the crossing angle
  // Projectile runs towards the MUON arm
  TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
  TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
  //
  // --- Get the muons parameters in the CM frame
  //
  TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
  TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
  //
  // --- Obtain the dimuon parameters in the CM frame
  //
  TLorentzVector pDimuCM=pMu1CM+pMu2CM;
  //
  // --- Translate the dimuon parameters in the dimuon rest frame
  //
  TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
  TLorentzVector pMu1Dimu=pMu1CM;
  TLorentzVector pMu2Dimu=pMu2CM;
  TLorentzVector pProjDimu=pProjCM;
  TLorentzVector pTargDimu=pTargCM;
  pMu1Dimu.Boost(beta);
  pMu2Dimu.Boost(beta);
  pProjDimu.Boost(beta);
  pTargDimu.Boost(beta);
  //
  // --- Determine the z axis for the CS angle 
  //
  TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
  //
  // --- Determine the CS angle (angle between mu+ and the z axis defined above)
  //
  Double_t cost;
  if(mu1Charge > 0) {
    cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
    // Theta CS is not properly defined for Like-Sign muons
    if(mu2Charge > 0 && cost<0) cost=-cost;
  } else { 
    // Theta CS is not properly defined for Like-Sign muons
    cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
    if(mu2Charge < 0 && cost<0) cost=-cost;
  }
  return cost;
}

//______________________________________________________________________________
Double_t AliAODDimuon::CostHe(){
  // Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
  if(CheckPointers())return -999999999;
  Double_t ebeam=3500; //temporary
  if(ebeam<=0){
    printf("Can not compute costCS with EBeam=%f\n",ebeam);
    return -999999999;
  }
  Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
  Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
  Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
  Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
  Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
  Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
  Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
  Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
  Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
  Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
  //
  // --- Get the muons parameters in the CM frame
  //
  TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
  TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
  //
  // --- Obtain the dimuon parameters in the CM frame
  //
  TLorentzVector pDimuCM=pMu1CM+pMu2CM;
  //
  // --- Translate the dimuon parameters in the dimuon rest frame
  //
  TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
  TLorentzVector pMu1Dimu=pMu1CM;
  TLorentzVector pMu2Dimu=pMu2CM;
  pMu1Dimu.Boost(beta);
  pMu2Dimu.Boost(beta);
  //
  // --- Determine the z axis for the calculation of the polarization angle
  // (i.e. the direction of the dimuon in the CM system)
  //
  TVector3 zaxis;
  zaxis=(pDimuCM.Vect()).Unit();
  //
  // --- Calculation of the polarization angle (Helicity)
  // (angle between mu+ and the z axis defined above)
  //
  Double_t cost;
  if(mu1Charge > 0) {
    cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
    // Theta Helicity is not properly defined for Like-Sign muons
    if(mu2Charge > 0 && cost<0) cost=-cost;
  } else { 
    cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
    // Theta Helicity is not properly defined for Like-Sign muons
    if(mu2Charge < 0 && cost<0) cost=-cost;
  }  
  return cost;
}

//________________________________________________________________________
Double_t AliAODDimuon::PhiCS(){
  // Cosinus of the Collins-Soper polar decay angle
  if(CheckPointers())return -999999999;
  Double_t ebeam=3500.;  //temporary
  if(ebeam<=0){
    printf("Can not compute phiCS with EBeam=%f\n",ebeam);
    return -999999999;
  }
  Double_t mp=fMProton;
  Double_t pbeam=TMath::Sqrt(ebeam*ebeam-mp*mp);
  Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
  Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
  Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
  Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
  Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
  Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
  Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
  Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
  Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();

  // Fill the Lorentz vector for projectile and target
  // For the moment we do not consider the crossing angle
  // Projectile runs towards the MUON arm
  TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
  TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
  //
  // --- Get the muons parameters in the CM frame
  //
  TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
  TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
  //
  // --- Obtain the dimuon parameters in the CM frame
  //
  TLorentzVector pDimuCM=pMu1CM+pMu2CM;
  //
  // --- Translate the dimuon parameters in the dimuon rest frame
  //
  TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
  TLorentzVector pMu1Dimu=pMu1CM;
  TLorentzVector pMu2Dimu=pMu2CM;
  TLorentzVector pProjDimu=pProjCM;
  TLorentzVector pTargDimu=pTargCM;
  pMu1Dimu.Boost(beta);
  pMu2Dimu.Boost(beta);
  pProjDimu.Boost(beta);
  pTargDimu.Boost(beta);
  //
  // --- Determine the z axis for the CS angle 
  //
  TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
  //
  // --- Determine the CS angle (angle between mu+ and the z axis defined above)
  //
   TVector3 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
   TVector3 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 
   Double_t phi;
   if(mu1Charge>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
   else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
     
   return phi;
}

//______________________________________________________________________________
Double_t AliAODDimuon::PhiHe(){
  // Calculation the Helicity aimuthal angle (adapted from code by R. Arnaldi)
  if(CheckPointers())return -999999999;
  Double_t ebeam=3500; //temporary
  if(ebeam<=0){
    printf("Can not compute phiHE with EBeam=%f\n",ebeam);
    return -999999999;
  }
  Double_t pbeam=TMath::Sqrt(ebeam*ebeam-fMProton*fMProton);
  Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
  Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
  Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
  Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
  Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
  Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
  Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
  Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
  Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();

  // Fill the Lorentz vector for projectile and target
  // For the moment we consider no crossing angle
  // Projectile runs towards the MUON arm
  TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
  TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
  //
  // --- Get the muons parameters in the CM frame
  //
  TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
  TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
  //
  // --- Obtain the dimuon parameters in the CM frame
  //
  TLorentzVector pDimuCM=pMu1CM+pMu2CM;
  //
  // --- Translate the muon parameters in the dimuon rest frame
  // 
  TVector3 zaxis=(pDimuCM.Vect()).Unit(); 
  //
  // --- Translate the dimuon parameters in the dimuon rest frame
  //
  TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
  TLorentzVector pMu1Dimu=pMu1CM;
  TLorentzVector pMu2Dimu=pMu2CM;
  pMu1Dimu.Boost(beta);
  pMu2Dimu.Boost(beta);
  
  TLorentzVector pProjDimu=pProjCM;
  TLorentzVector pTargDimu=pTargCM; 
  pProjDimu.Boost(beta);
  pTargDimu.Boost(beta);

  TVector3 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
  TVector3 xaxis=(yaxis.Cross(zaxis)).Unit();
  //
  // --- Calculation of the azimuthal angle (Helicity)
  //
   Double_t phi;
   if(mu1Charge>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
   else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
   
   return phi;
}

//______________________________________________________________________________
Int_t AliAODDimuon::AnyPt(){
  // Test if the two muons match two trigger tracks
  if(CheckPointers())return 0;
  return (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger())&&
         (((AliAODTrack*)fMu[1].GetObject())->MatchTrigger());
}

//______________________________________________________________________________
Int_t AliAODDimuon::LowPt(){
  // Test if the two muons match two trigger tracks with a "Low Pt" cut
  if(CheckPointers())return 0;
  return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt())&&
         (((AliAODTrack*)fMu[1].GetObject())->MatchTriggerLowPt());
}

//______________________________________________________________________________
Int_t AliAODDimuon::HighPt(){
  // Test if the two muons match two trigger tracks with a "High Pt" cut
  if(CheckPointers())return 0;
  return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt())&&
         (((AliAODTrack*)fMu[1].GetObject())->MatchTriggerHighPt());
}

//______________________________________________________________________________
Double_t AliAODDimuon::MaxChi2Match(){
  // Maximum matching Chi2 between track and trigger track
  if(CheckPointers())return -999999999;
  return TMath::Max((((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()),
                    (((AliAODTrack*)fMu[1].GetObject())->GetChi2MatchTrigger()));
}

//______________________________________________________________________________
TLorentzVector* AliAODDimuon::TLV() const
{
  delete fP;
  fP = new TLorentzVector(Px(),Py(),Pz(),E());
  return fP;  
}
Commit	Line	Data
866d8d78	1	// AliAODDimuon: a class for AODs for the MUON Arm of the ALICE Experiment
	2	// Author: P. Cortese, Universita' del Piemonte Orientale in Alessandria and
	3	// INFN of Torino - Italy
	4	//
	5	// The class defines a dimuon pair object from two AliAODTrack objects.
	6	// AliAODDimuon objects are supposed to be added to the AliAODEvent structure
	7	// during analysis. They would then allow to calculate the dimuon-related
	8	// kinematic variables with a minimal disk occupancy.
	9	// The payload of the class has been reduced to two pointers to the two
	10	// tracks. An instance of this class has also to be added to the AliAODEvent
	11	// structure to provide additional information that is specific to MUON and
	12	// therefore has not been included into the AOD header.
	13	// Two transient data members are not stored on file as they can be recomputed
	14	// at runtime.
	15	//
	16
	17	#include "AliAODDimuon.h"
	18	#include "TLorentzVector.h"
	19	#define AliAODDimuon_CXX
	20
	21	ClassImp(AliAODDimuon)
	22
	23	//______________________________________________________________________________
	24	AliAODDimuon::AliAODDimuon():AliVParticle(),fP(0),fMProton(0.93827231)
	25	{
	26	// default constructor
	27	fMu[0]=0;
	28	fMu[1]=0;
	29	}
	30
	31	//______________________________________________________________________________
	32	AliAODDimuon::AliAODDimuon(const AliAODDimuon& dimu):
	33	AliVParticle(dimu),
	34	fP(0),fMProton(0.93827231)
	35	{
	36	// copy constructor
	37	fMu[0]=dimu.Mu(0);
	38	fMu[1]=dimu.Mu(1);
	39	}
	40
	41	//______________________________________________________________________________
	42	AliAODDimuon &AliAODDimuon::operator=(const AliAODDimuon& dimu)
	43	{
	44	// assignment operator
	45	if(&dimu != this){
e7329cf1	46	delete fP;
866d8d78	47	fP=0;
	48	fMProton=0.93827231;
	49	fMu[0]=dimu.Mu(0);
	50	fMu[1]=dimu.Mu(1);
	51	}
	52	return *this;
	53	}
	54
	55	//______________________________________________________________________________
	56	AliAODDimuon::AliAODDimuon(TObject mu0, TObject mu1):
	57	fP(0),fMProton(0.93827231)
	58	{
	59	// Creates a dimuon pair from two tracks
8f6e7d10	60
866d8d78	61	fMu[0]=mu0;
	62	fMu[1]=mu1;
	63	}
	64
	65	//______________________________________________________________________________
	66	AliAODDimuon::~AliAODDimuon()
	67	{
	68	// destructor
e7329cf1	69	delete fP;
866d8d78	70	}
	71
	72	//______________________________________________________________________________
	73	Double_t AliAODDimuon::Px() const {
	74	// Px of the dimuon
e7329cf1	75	if(CheckPointers())return -999999999;
866d8d78	76	return ((AliAODTrack*)fMu[0].GetObject())->Px()+
	77	((AliAODTrack*)fMu[1].GetObject())->Px();
	78	}
	79
	80	//______________________________________________________________________________
	81	Double_t AliAODDimuon::Py() const {
	82	// Py of the dimuon
e7329cf1	83	if(CheckPointers())return -999999999;
866d8d78	84	return ((AliAODTrack*)fMu[0].GetObject())->Py()+
	85	((AliAODTrack*)fMu[1].GetObject())->Py();
	86	}
	87
	88	//______________________________________________________________________________
	89	Double_t AliAODDimuon::Pz() const {
	90	// Pz of the dimuon
e7329cf1	91	if(CheckPointers())return -999999999;
866d8d78	92	return ((AliAODTrack*)fMu[0].GetObject())->Pz()+
	93	((AliAODTrack*)fMu[1].GetObject())->Pz();
	94	}
	95
	96	//______________________________________________________________________________
	97	Double_t AliAODDimuon::Pt() const {
	98	// Pt of the dimuon
e7329cf1	99	if(CheckPointers())return -999999999;
866d8d78	100	Double_t px=Px();
	101	Double_t py=Py();
	102	return TMath::Sqrt(pxpx+pypy);
866d8d78	103	}
	104
	105	//______________________________________________________________________________
e7329cf1	106	Double_t AliAODDimuon::E() const
	107	{
	108	// Dimuon energy
	109
	110	if(CheckPointers())return -999999999;
866d8d78	111
e7329cf1	112	return ((AliAODTrack)fMu[0].GetObject())->E()+ ((AliAODTrack)fMu[1].GetObject())->E();
866d8d78	113	}
	114
	115	//______________________________________________________________________________
e7329cf1	116	Double_t AliAODDimuon::P() const {
866d8d78	117	// Dimuon momentum
e7329cf1	118	if(CheckPointers())return -999999999;
e7329cf1	119	return TLV()->P();
866d8d78	120	}
	121
	122	//______________________________________________________________________________
	123	Double_t AliAODDimuon::M() const {
866d8d78	124	// Dimuon invariant mass
e7329cf1	125	if(CheckPointers())return -999999999;
e7329cf1	126	return TLV()->M();
866d8d78	127	}
	128
	129	//______________________________________________________________________________
	130	Double_t AliAODDimuon::Eta() const {
866d8d78	131	// Dimuon pseudorapidity
e7329cf1	132	if(CheckPointers())return -999999999;
e7329cf1	133	return TLV()->Eta();
866d8d78	134	}
	135
	136	//______________________________________________________________________________
	137	Double_t AliAODDimuon::Phi() const {
866d8d78	138	// Dimuon asimuthal angle
e7329cf1	139	if(CheckPointers())return -999999999;
e7329cf1	140	return TLV()->Phi();
866d8d78	141	}
866d8d78	142
e7329cf1	143
866d8d78	144	//______________________________________________________________________________
e7329cf1	145	Double_t AliAODDimuon::Theta() const {
866d8d78	146	// Dimuon polar angle
e7329cf1	147	if(CheckPointers())return -999999999;
e7329cf1	148	return TLV()->Theta();
866d8d78	149	}
	150
	151	//______________________________________________________________________________
	152	Double_t AliAODDimuon::Y() const {
866d8d78	153	// Dimuon rapidity
e7329cf1	154	if(CheckPointers())return -999999999;
e7329cf1	155	return TLV()->Rapidity();
866d8d78	156	}
	157
	158	//______________________________________________________________________________
	159	Short_t AliAODDimuon::Charge() const {
	160	// Dimuon charge
e7329cf1	161	if(CheckPointers())return -999;
866d8d78	162	return ((AliAODTrack*)fMu[0].GetObject())->Charge()+
	163	((AliAODTrack*)fMu[1].GetObject())->Charge();
	164	}
	165
	166	//______________________________________________________________________________
	167	Int_t AliAODDimuon::CheckPointers() const{
	168	// Checks if the track pointers have been initialized
	169	if(fMu[0]==0\|\|fMu[1]==0){
	170	printf("Dimuon not initialized\n");
	171	return -999;
	172	}
	173	if((fMu[0].GetObject())==0\|\|(fMu[1].GetObject())==0){
48453d08	174	printf("Can not get objects\n");
866d8d78	175	return -999;
	176	}
	177	return 0;
	178	}
	179
	180	//______________________________________________________________________________
	181	void AliAODDimuon::SetMu(Int_t imu, AliAODTrack *mu){
	182	// Assign a track pointer
	183	if (imu==0\|\|imu==1){
e7329cf1	184	delete fP;
e7329cf1	185	fP=0;
866d8d78	186	fMu[imu]=mu;
	187	}
	188	}
	189
	190	//______________________________________________________________________________
	191	void AliAODDimuon::SetMuons(AliAODTrack mu0, AliAODTrack mu1){
	192	// Assign the track pointers
	193	fMu[0]=mu0;
	194	fMu[1]=mu1;
e7329cf1	195	delete fP;
e7329cf1	196	fP=0;
866d8d78	197	}
	198
	199	//______________________________________________________________________________
	200	Double_t AliAODDimuon::XF() {
	201	// Dimuon Feynman x
e7329cf1	202
	203	if(CheckPointers())return -999999999;
	204
866d8d78	205	//Double_t ebeam=((AliAODEventInfo*)fEi.GetObject())->EBeam();
	206	Double_t ebeam = 3500.; // temporary
	207	if(ebeam<=0){
	208	printf("AliAODDimuon::xf: can not compute xf with EBeam=%f\n",ebeam);
	209	return -999999999;
	210	}
866d8d78	211	Double_t mDimu=M();
	212	Double_t pMax=TMath::Sqrt(ebeamebeam-mDimumDimu);
	213	return Pz()/pMax;
	214	}
	215
	216	//______________________________________________________________________________
866d8d78	217	Double_t AliAODDimuon::CostCS(){
48453d08	218	// Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
866d8d78	219	if(CheckPointers())return -999999999;
	220	Double_t ebeam=3500.; //temporary
	221	if(ebeam<=0){
	222	printf("Can not compute costCS with EBeam=%f\n",ebeam);
	223	return -999999999;
	224	}
	225	Double_t mp=fMProton;
	226	Double_t pbeam=TMath::Sqrt(ebeamebeam-mpmp);
	227	Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
	228	Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
	229	Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
	230	Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
	231	Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
	232	Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
	233	Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
	234	Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
	235	Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
	236	Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
	237
	238	// Fill the Lorentz vector for projectile and target
	239	// For the moment we do not consider the crossing angle
	240	// Projectile runs towards the MUON arm
48453d08	241	TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
48453d08	242	TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
866d8d78	243	//
48453d08	244	// --- Get the muons parameters in the CM frame
866d8d78	245	//
48453d08	246	TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
48453d08	247	TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
866d8d78	248	//
48453d08	249	// --- Obtain the dimuon parameters in the CM frame
866d8d78	250	//
48453d08	251	TLorentzVector pDimuCM=pMu1CM+pMu2CM;
866d8d78	252	//
	253	// --- Translate the dimuon parameters in the dimuon rest frame
	254	//
48453d08	255	TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
	256	TLorentzVector pMu1Dimu=pMu1CM;
	257	TLorentzVector pMu2Dimu=pMu2CM;
	258	TLorentzVector pProjDimu=pProjCM;
	259	TLorentzVector pTargDimu=pTargCM;
866d8d78	260	pMu1Dimu.Boost(beta);
	261	pMu2Dimu.Boost(beta);
	262	pProjDimu.Boost(beta);
	263	pTargDimu.Boost(beta);
	264	//
	265	// --- Determine the z axis for the CS angle
	266	//
	267	TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
	268	//
	269	// --- Determine the CS angle (angle between mu+ and the z axis defined above)
	270	//
	271	Double_t cost;
	272	if(mu1Charge > 0) {
	273	cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
	274	// Theta CS is not properly defined for Like-Sign muons
	275	if(mu2Charge > 0 && cost<0) cost=-cost;
	276	} else {
	277	// Theta CS is not properly defined for Like-Sign muons
	278	cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
	279	if(mu2Charge < 0 && cost<0) cost=-cost;
	280	}
	281	return cost;
	282	}
	283
	284	//______________________________________________________________________________
866d8d78	285	Double_t AliAODDimuon::CostHe(){
48453d08	286	// Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
866d8d78	287	if(CheckPointers())return -999999999;
	288	Double_t ebeam=3500; //temporary
	289	if(ebeam<=0){
	290	printf("Can not compute costCS with EBeam=%f\n",ebeam);
	291	return -999999999;
	292	}
866d8d78	293	Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
	294	Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
	295	Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
	296	Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
	297	Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
	298	Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
	299	Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
	300	Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
	301	Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
	302	Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
866d8d78	303	//
48453d08	304	// --- Get the muons parameters in the CM frame
866d8d78	305	//
48453d08	306	TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
48453d08	307	TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
866d8d78	308	//
48453d08	309	// --- Obtain the dimuon parameters in the CM frame
866d8d78	310	//
48453d08	311	TLorentzVector pDimuCM=pMu1CM+pMu2CM;
866d8d78	312	//
	313	// --- Translate the dimuon parameters in the dimuon rest frame
	314	//
48453d08	315	TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
	316	TLorentzVector pMu1Dimu=pMu1CM;
	317	TLorentzVector pMu2Dimu=pMu2CM;
866d8d78	318	pMu1Dimu.Boost(beta);
	319	pMu2Dimu.Boost(beta);
	320	//
866d8d78	321	// --- Determine the z axis for the calculation of the polarization angle
	322	// (i.e. the direction of the dimuon in the CM system)
	323	//
	324	TVector3 zaxis;
	325	zaxis=(pDimuCM.Vect()).Unit();
	326	//
	327	// --- Calculation of the polarization angle (Helicity)
	328	// (angle between mu+ and the z axis defined above)
	329	//
	330	Double_t cost;
	331	if(mu1Charge > 0) {
	332	cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
	333	// Theta Helicity is not properly defined for Like-Sign muons
	334	if(mu2Charge > 0 && cost<0) cost=-cost;
	335	} else {
	336	cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
	337	// Theta Helicity is not properly defined for Like-Sign muons
	338	if(mu2Charge < 0 && cost<0) cost=-cost;
	339	}
	340	return cost;
	341	}
	342
48453d08	343	//________________________________________________________________________
	344	Double_t AliAODDimuon::PhiCS(){
	345	// Cosinus of the Collins-Soper polar decay angle
	346	if(CheckPointers())return -999999999;
	347	Double_t ebeam=3500.; //temporary
	348	if(ebeam<=0){
	349	printf("Can not compute phiCS with EBeam=%f\n",ebeam);
	350	return -999999999;
	351	}
	352	Double_t mp=fMProton;
	353	Double_t pbeam=TMath::Sqrt(ebeamebeam-mpmp);
	354	Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
	355	Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
	356	Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
	357	Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
	358	Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
	359	Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
	360	Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
	361	Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
	362	Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
	363
	364	// Fill the Lorentz vector for projectile and target
	365	// For the moment we do not consider the crossing angle
	366	// Projectile runs towards the MUON arm
	367	TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
	368	TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
	369	//
	370	// --- Get the muons parameters in the CM frame
	371	//
	372	TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
	373	TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
	374	//
	375	// --- Obtain the dimuon parameters in the CM frame
	376	//
	377	TLorentzVector pDimuCM=pMu1CM+pMu2CM;
	378	//
	379	// --- Translate the dimuon parameters in the dimuon rest frame
	380	//
	381	TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
	382	TLorentzVector pMu1Dimu=pMu1CM;
	383	TLorentzVector pMu2Dimu=pMu2CM;
	384	TLorentzVector pProjDimu=pProjCM;
	385	TLorentzVector pTargDimu=pTargCM;
	386	pMu1Dimu.Boost(beta);
	387	pMu2Dimu.Boost(beta);
	388	pProjDimu.Boost(beta);
	389	pTargDimu.Boost(beta);
	390	//
	391	// --- Determine the z axis for the CS angle
	392	//
	393	TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
	394	//
	395	// --- Determine the CS angle (angle between mu+ and the z axis defined above)
	396	//
	397	TVector3 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
	398	TVector3 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
	399
	400	Double_t phi;
	401	if(mu1Charge>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
	402	else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
	403
	404	return phi;
	405	}
	406
407	//______________________________________________________________________________
408	Double_t AliAODDimuon::PhiHe(){
409	// Calculation the Helicity aimuthal angle (adapted from code by R. Arnaldi)
410	if(CheckPointers())return -999999999;
411	Double_t ebeam=3500; //temporary
412	if(ebeam<=0){
413	printf("Can not compute phiHE with EBeam=%f\n",ebeam);
414	return -999999999;
415	}
416	Double_t pbeam=TMath::Sqrt(ebeamebeam-fMProtonfMProton);
417	Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
418	Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
419	Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
420	Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
421	Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
422	Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
423	Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
424	Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
425	Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
426
427	// Fill the Lorentz vector for projectile and target
428	// For the moment we consider no crossing angle
429	// Projectile runs towards the MUON arm
430	TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
431	TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
432	//
433	// --- Get the muons parameters in the CM frame
434	//
435	TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
436	TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
437	//
438	// --- Obtain the dimuon parameters in the CM frame
439	//
440	TLorentzVector pDimuCM=pMu1CM+pMu2CM;
441	//
442	// --- Translate the muon parameters in the dimuon rest frame
443	//
444	TVector3 zaxis=(pDimuCM.Vect()).Unit();
445	//
446	// --- Translate the dimuon parameters in the dimuon rest frame
447	//
448	TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
449	TLorentzVector pMu1Dimu=pMu1CM;
450	TLorentzVector pMu2Dimu=pMu2CM;
451	pMu1Dimu.Boost(beta);
452	pMu2Dimu.Boost(beta);
453
454	TLorentzVector pProjDimu=pProjCM;
455	TLorentzVector pTargDimu=pTargCM;
456	pProjDimu.Boost(beta);
457	pTargDimu.Boost(beta);
458
459	TVector3 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
460	TVector3 xaxis=(yaxis.Cross(zaxis)).Unit();
461	//
462	// --- Calculation of the azimuthal angle (Helicity)
463	//
464	Double_t phi;
465	if(mu1Charge>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
466	else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
467
468	return phi;
469	}
470
866d8d78	471	//______________________________________________________________________________
	472	Int_t AliAODDimuon::AnyPt(){
	473	// Test if the two muons match two trigger tracks
e7329cf1	474	if(CheckPointers())return 0;
866d8d78	475	return (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger())&&
8f6c57a4	476	(((AliAODTrack*)fMu[1].GetObject())->MatchTrigger());
866d8d78	477	}
	478
	479	//______________________________________________________________________________
	480	Int_t AliAODDimuon::LowPt(){
	481	// Test if the two muons match two trigger tracks with a "Low Pt" cut
e7329cf1	482	if(CheckPointers())return 0;
866d8d78	483	return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt())&&
8f6c57a4	484	(((AliAODTrack*)fMu[1].GetObject())->MatchTriggerLowPt());
866d8d78	485	}
	486
	487	//______________________________________________________________________________
	488	Int_t AliAODDimuon::HighPt(){
	489	// Test if the two muons match two trigger tracks with a "High Pt" cut
e7329cf1	490	if(CheckPointers())return 0;
866d8d78	491	return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt())&&
8f6c57a4	492	(((AliAODTrack*)fMu[1].GetObject())->MatchTriggerHighPt());
866d8d78	493	}
	494
	495	//______________________________________________________________________________
	496	Double_t AliAODDimuon::MaxChi2Match(){
	497	// Maximum matching Chi2 between track and trigger track
e7329cf1	498	if(CheckPointers())return -999999999;
866d8d78	499	return TMath::Max((((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()),
8f6c57a4	500	(((AliAODTrack*)fMu[1].GetObject())->GetChi2MatchTrigger()));
866d8d78	501	}
e7329cf1	502
	503	//______________________________________________________________________________
	504	TLorentzVector* AliAODDimuon::TLV() const
	505	{
	506	delete fP;
	507	fP = new TLorentzVector(Px(),Py(),Pz(),E());
	508	return fP;
	509	}