[u/mrichter/AliRoot.git] / STEER / 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){
    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
  
  //printf("Creating dimuon from %p %p\n",mu0,mu1);
  fMu[0]=mu0;
  fMu[1]=mu1;
}

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

//______________________________________________________________________________
void AliAODDimuon::BookP(){
  // Fills the dimuon momentum if not filled yet
  static UInt_t unID[2]={0,0};
  if(!fP){
    fP=new TLorentzVector(Px(),Py(),Pz(),E());
    unID[0]=fMu[0].GetUniqueID();
    unID[1]=fMu[1].GetUniqueID();
  }
  // For efficiency reasons
  if((unID[0]!=fMu[0].GetUniqueID())||(unID[1]!=fMu[1].GetUniqueID())){
    fP->SetPxPyPzE(Px(),Py(),Pz(),E());
    unID[0]=fMu[0].GetUniqueID();
    unID[1]=fMu[1].GetUniqueID();
  }
}

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

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

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

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

//______________________________________________________________________________
Double_t AliAODDimuon::E() const {
  // Dimuon energy
  if(this->CheckPointers())return -999999999;
  return ((AliAODTrack*)fMu[0].GetObject())->E()+
         ((AliAODTrack*)fMu[1].GetObject())->E();
}

//______________________________________________________________________________
Double_t AliAODDimuon::P() const {
  // This is just to override the virtual function
  printf("You should never call: Double_t AliAODDimuon::P() const\n");
  return -999999999;
}

//______________________________________________________________________________
Double_t AliAODDimuon::P() {
  // Dimuon momentum
  if(this->CheckPointers())return -999999999;
  BookP();
  return fP->P();
}

//______________________________________________________________________________
Double_t AliAODDimuon::M() const {
  // This is just to override the virtual function
  printf("You should never call: Double_t AliAODDimuon::M() const\n");
  return -999999999;
}

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

//______________________________________________________________________________
Double_t AliAODDimuon::Mass() {
  // Dimuon invariant mass
  if(this->CheckPointers())return -999999999;
  BookP();
  return fP->M();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Eta() const {
  // This is just to override the virtual function
  printf("You should never call: Double_t AliAODDimuon::Eta() const\n");
  return -999999999;
}

//______________________________________________________________________________
Double_t AliAODDimuon::Eta() {
  // Dimuon pseudorapidity
  if(this->CheckPointers())return -999999999;
  BookP();
  return fP->Eta();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Phi() const {
  // This is just to override the virtual function
  printf("You should never call: Double_t AliAODDimuon::Phi() const\n");
  return -999999999;
}

//______________________________________________________________________________
Double_t AliAODDimuon::Phi() {
  // Dimuon asimuthal angle
  if(this->CheckPointers())return -999999999;
  BookP();
  return fP->Phi();
}
//______________________________________________________________________________
Double_t AliAODDimuon::Theta() const {
  // This is just to override the virtual function
  printf("You should never call: Double_t AliAODDimuon::Theta() const\n");
  return -999999999;
}

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

//______________________________________________________________________________
Double_t AliAODDimuon::Y() const {
  // This is just to override the virtual function
  printf("You should never call: Double_t AliAODDimuon::Y() const\n");
  return -999999999;
}

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

//______________________________________________________________________________
Short_t AliAODDimuon::Charge() const {
  // Dimuon charge
  if(this->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. Got: %p %p\n",fMu[0].GetObject(),fMu[1].GetObject());
    return -999;
  }
  return 0;
}

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

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

//______________________________________________________________________________
Double_t AliAODDimuon::XF() {
  // Dimuon Feynman x
  //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;
  }
  if(this->CheckPointers())return -999999999;
  BookP();
  Double_t mDimu=M();
  Double_t pMax=TMath::Sqrt(ebeam*ebeam-mDimu*mDimu);
  return Pz()/pMax;
}

//______________________________________________________________________________
// Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
Double_t AliAODDimuon::CostCS(){
  // Cosinus of the Collins-Soper polar decay angle
  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 pProjLab(0.,0.,-pbeam,ebeam); // projectile
  TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
  //
  // --- Get the muons parameters in the LAB frame
  //
  TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
  TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
  //
  // --- Obtain the dimuon parameters in the LAB frame
  //
  TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
  //
  // --- Translate the dimuon parameters in the dimuon rest frame
  //
  TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
  TLorentzVector pMu1Dimu=pMu1Lab;
  TLorentzVector pMu2Dimu=pMu2Lab;
  TLorentzVector pProjDimu=pProjLab;
  TLorentzVector pTargDimu=pTargLab;
  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;
}

//______________________________________________________________________________
// Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
Double_t AliAODDimuon::CostHe(){
  // Cosinus of the polar decay angle in the Helicity reference frame
  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 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();
  Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();

  // Fill the Lorentz vector for projectile and target
  // For the moment we consider no crossing angle
  // Projectile runs towards the MUON arm
  TLorentzVector pProjLab(0.,0.,-pbeam,ebeam); // projectile
  TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
  //
  // --- Get the muons parameters in the LAB frame
  //
  TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
  TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
  //
  // --- Obtain the dimuon parameters in the LAB frame
  //
  TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
  //
  // --- Translate the dimuon parameters in the dimuon rest frame
  //
  TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
  TLorentzVector pMu1Dimu=pMu1Lab;
  TLorentzVector pMu2Dimu=pMu2Lab;
  pMu1Dimu.Boost(beta);
  pMu2Dimu.Boost(beta);
  //
  // --- Translate the dimuon parameters in the CM frame
  //
  TLorentzVector pDimuCM; //CM frame
  TVector3 beta2;
  beta2=(-1./(fMProton+pProjLab.E()))*pProjLab.Vect();
  pDimuCM=pDimuLab;
  pDimuCM.Boost(beta2);
  //
  // --- 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;
}

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

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

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

//______________________________________________________________________________
Double_t AliAODDimuon::MaxChi2Match(){
  // Maximum matching Chi2 between track and trigger track
  if(this->CheckPointers())return -999999999;
  return TMath::Max((((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()),
                    (((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()));
}
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){
	46	fP=0;
	47	fMProton=0.93827231;
	48	fMu[0]=dimu.Mu(0);
	49	fMu[1]=dimu.Mu(1);
	50	}
	51	return *this;
	52	}
	53
	54	//______________________________________________________________________________
	55	AliAODDimuon::AliAODDimuon(TObject mu0, TObject mu1):
	56	fP(0),fMProton(0.93827231)
	57	{
	58	// Creates a dimuon pair from two tracks
	59
	60	//printf("Creating dimuon from %p %p\n",mu0,mu1);
	61	fMu[0]=mu0;
	62	fMu[1]=mu1;
	63	}
	64
65	//______________________________________________________________________________
66	AliAODDimuon::~AliAODDimuon()
67	{
68	// destructor
69	if(fP)delete fP;
70	fP=0;
71	}
72
73	//______________________________________________________________________________
74	void AliAODDimuon::BookP(){
75	// Fills the dimuon momentum if not filled yet
76	static UInt_t unID[2]={0,0};
77	if(!fP){
78	fP=new TLorentzVector(Px(),Py(),Pz(),E());
79	unID[0]=fMu[0].GetUniqueID();
80	unID[1]=fMu[1].GetUniqueID();
81	}
82	// For efficiency reasons
83	if((unID[0]!=fMu[0].GetUniqueID())\|\|(unID[1]!=fMu[1].GetUniqueID())){
84	fP->SetPxPyPzE(Px(),Py(),Pz(),E());
85	unID[0]=fMu[0].GetUniqueID();
86	unID[1]=fMu[1].GetUniqueID();
87	}
88	}
89
90	//______________________________________________________________________________
91	Double_t AliAODDimuon::Px() const {
92	// Px of the dimuon
93	if(this->CheckPointers())return -999999999;
94	return ((AliAODTrack*)fMu[0].GetObject())->Px()+
95	((AliAODTrack*)fMu[1].GetObject())->Px();
96	}
97
98	//______________________________________________________________________________
99	Double_t AliAODDimuon::Py() const {
100	// Py of the dimuon
101	if(this->CheckPointers())return -999999999;
102	return ((AliAODTrack*)fMu[0].GetObject())->Py()+
103	((AliAODTrack*)fMu[1].GetObject())->Py();
104	}
105
106	//______________________________________________________________________________
107	Double_t AliAODDimuon::Pz() const {
108	// Pz of the dimuon
109	if(this->CheckPointers())return -999999999;
110	return ((AliAODTrack*)fMu[0].GetObject())->Pz()+
111	((AliAODTrack*)fMu[1].GetObject())->Pz();
112	}
113
114	//______________________________________________________________________________
115	Double_t AliAODDimuon::Pt() const {
116	// Pt of the dimuon
117	if(this->CheckPointers())return -999999999;
118	Double_t px=Px();
119	Double_t py=Py();
120	return TMath::Sqrt(pxpx+pypy);
121	return -999999999;
122	}
123
124	//______________________________________________________________________________
125	Double_t AliAODDimuon::E() const {
126	// Dimuon energy
127	if(this->CheckPointers())return -999999999;
128	return ((AliAODTrack*)fMu[0].GetObject())->E()+
129	((AliAODTrack*)fMu[1].GetObject())->E();
130	}
131
132	//______________________________________________________________________________
133	Double_t AliAODDimuon::P() const {
134	// This is just to override the virtual function
135	printf("You should never call: Double_t AliAODDimuon::P() const\n");
136	return -999999999;
137	}
138
139	//______________________________________________________________________________
140	Double_t AliAODDimuon::P() {
141	// Dimuon momentum
142	if(this->CheckPointers())return -999999999;
143	BookP();
144	return fP->P();
145	}
146
147	//______________________________________________________________________________
148	Double_t AliAODDimuon::M() const {
149	// This is just to override the virtual function
150	printf("You should never call: Double_t AliAODDimuon::M() const\n");
151	return -999999999;
152	}
153
154	//______________________________________________________________________________
155	Double_t AliAODDimuon::M() {
156	// Dimuon invariant mass
157	if(this->CheckPointers())return -999999999;
158	BookP();
159	return fP->M();
160	}
161
162	//______________________________________________________________________________
163	Double_t AliAODDimuon::Mass() {
164	// Dimuon invariant mass
165	if(this->CheckPointers())return -999999999;
166	BookP();
167	return fP->M();
168	}
169
170	//______________________________________________________________________________
171	Double_t AliAODDimuon::Eta() const {
172	// This is just to override the virtual function
173	printf("You should never call: Double_t AliAODDimuon::Eta() const\n");
174	return -999999999;
175	}
176
177	//______________________________________________________________________________
178	Double_t AliAODDimuon::Eta() {
179	// Dimuon pseudorapidity
180	if(this->CheckPointers())return -999999999;
181	BookP();
182	return fP->Eta();
183	}
184
185	//______________________________________________________________________________
186	Double_t AliAODDimuon::Phi() const {
187	// This is just to override the virtual function
188	printf("You should never call: Double_t AliAODDimuon::Phi() const\n");
189	return -999999999;
190	}
191
192	//______________________________________________________________________________
193	Double_t AliAODDimuon::Phi() {
194	// Dimuon asimuthal angle
195	if(this->CheckPointers())return -999999999;
196	BookP();
197	return fP->Phi();
198	}
199	//______________________________________________________________________________
200	Double_t AliAODDimuon::Theta() const {
201	// This is just to override the virtual function
202	printf("You should never call: Double_t AliAODDimuon::Theta() const\n");
203	return -999999999;
204	}
205
206	//______________________________________________________________________________
207	Double_t AliAODDimuon::Theta() {
208	// Dimuon polar angle
209	if(this->CheckPointers())return -999999999;
210	BookP();
211	return fP->Theta();
212	}
213
214	//______________________________________________________________________________
215	Double_t AliAODDimuon::Y() const {
216	// This is just to override the virtual function
217	printf("You should never call: Double_t AliAODDimuon::Y() const\n");
218	return -999999999;
219	}
220
221	//______________________________________________________________________________
222	Double_t AliAODDimuon::Y() {
223	// Dimuon rapidity
224	if(this->CheckPointers())return -999999999;
225	BookP();
226	return fP->Rapidity();
227	}
228
229	//______________________________________________________________________________
230	Short_t AliAODDimuon::Charge() const {
231	// Dimuon charge
232	if(this->CheckPointers())return -999;
233	return ((AliAODTrack*)fMu[0].GetObject())->Charge()+
234	((AliAODTrack*)fMu[1].GetObject())->Charge();
235	}
236
237	//______________________________________________________________________________
238	Int_t AliAODDimuon::CheckPointers() const{
239	// Checks if the track pointers have been initialized
240	if(fMu[0]==0\|\|fMu[1]==0){
241	printf("Dimuon not initialized\n");
242	return -999;
243	}
244	if((fMu[0].GetObject())==0\|\|(fMu[1].GetObject())==0){
245	printf("Can not get objects. Got: %p %p\n",fMu[0].GetObject(),fMu[1].GetObject());
246	return -999;
247	}
248	return 0;
249	}
250
251	//______________________________________________________________________________
252	void AliAODDimuon::SetMu(Int_t imu, AliAODTrack *mu){
253	// Assign a track pointer
254	if (imu==0\|\|imu==1){
255	fMu[imu]=mu;
256	}
257	}
258
259	//______________________________________________________________________________
260	void AliAODDimuon::SetMuons(AliAODTrack mu0, AliAODTrack mu1){
261	// Assign the track pointers
262	fMu[0]=mu0;
263	fMu[1]=mu1;
264	}
265
266	//______________________________________________________________________________
267	Double_t AliAODDimuon::XF() {
268	// Dimuon Feynman x
269	//Double_t ebeam=((AliAODEventInfo*)fEi.GetObject())->EBeam();
270	Double_t ebeam = 3500.; // temporary
271	if(ebeam<=0){
272	printf("AliAODDimuon::xf: can not compute xf with EBeam=%f\n",ebeam);
273	return -999999999;
274	}
275	if(this->CheckPointers())return -999999999;
276	BookP();
277	Double_t mDimu=M();
278	Double_t pMax=TMath::Sqrt(ebeamebeam-mDimumDimu);
279	return Pz()/pMax;
280	}
281
282	//______________________________________________________________________________
283	// Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
284	Double_t AliAODDimuon::CostCS(){
285	// Cosinus of the Collins-Soper polar decay angle
286	if(CheckPointers())return -999999999;
287	Double_t ebeam=3500.; //temporary
288	if(ebeam<=0){
289	printf("Can not compute costCS with EBeam=%f\n",ebeam);
290	return -999999999;
291	}
292	Double_t mp=fMProton;
293	Double_t pbeam=TMath::Sqrt(ebeamebeam-mpmp);
294	Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
295	Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
296	Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
297	Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
298	Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
299	Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
300	Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
301	Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
302	Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
303	Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
304
305	// Fill the Lorentz vector for projectile and target
306	// For the moment we do not consider the crossing angle
307	// Projectile runs towards the MUON arm
308	TLorentzVector pProjLab(0.,0.,-pbeam,ebeam); // projectile
309	TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
310	//
311	// --- Get the muons parameters in the LAB frame
312	//
313	TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
314	TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
315	//
316	// --- Obtain the dimuon parameters in the LAB frame
317	//
318	TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
319	//
320	// --- Translate the dimuon parameters in the dimuon rest frame
321	//
322	TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
323	TLorentzVector pMu1Dimu=pMu1Lab;
324	TLorentzVector pMu2Dimu=pMu2Lab;
325	TLorentzVector pProjDimu=pProjLab;
326	TLorentzVector pTargDimu=pTargLab;
327	pMu1Dimu.Boost(beta);
328	pMu2Dimu.Boost(beta);
329	pProjDimu.Boost(beta);
330	pTargDimu.Boost(beta);
331	//
332	// --- Determine the z axis for the CS angle
333	//
334	TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
335	//
336	// --- Determine the CS angle (angle between mu+ and the z axis defined above)
337	//
338	Double_t cost;
339	if(mu1Charge > 0) {
340	cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
341	// Theta CS is not properly defined for Like-Sign muons
342	if(mu2Charge > 0 && cost<0) cost=-cost;
343	} else {
344	// Theta CS is not properly defined for Like-Sign muons
345	cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
346	if(mu2Charge < 0 && cost<0) cost=-cost;
347	}
348	return cost;
349	}
350
351	//______________________________________________________________________________
352	// Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
353	Double_t AliAODDimuon::CostHe(){
354	// Cosinus of the polar decay angle in the Helicity reference frame
355	if(CheckPointers())return -999999999;
356	Double_t ebeam=3500; //temporary
357	if(ebeam<=0){
358	printf("Can not compute costCS with EBeam=%f\n",ebeam);
359	return -999999999;
360	}
361	Double_t pbeam=TMath::Sqrt(ebeamebeam-fMProtonfMProton);
362	Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
363	Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
364	Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
365	Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
366	Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
367	Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
368	Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
369	Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
370	Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
371	Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
372
373	// Fill the Lorentz vector for projectile and target
374	// For the moment we consider no crossing angle
375	// Projectile runs towards the MUON arm
376	TLorentzVector pProjLab(0.,0.,-pbeam,ebeam); // projectile
377	TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
378	//
379	// --- Get the muons parameters in the LAB frame
380	//
381	TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
382	TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
383	//
384	// --- Obtain the dimuon parameters in the LAB frame
385	//
386	TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
387	//
388	// --- Translate the dimuon parameters in the dimuon rest frame
389	//
390	TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
391	TLorentzVector pMu1Dimu=pMu1Lab;
392	TLorentzVector pMu2Dimu=pMu2Lab;
393	pMu1Dimu.Boost(beta);
394	pMu2Dimu.Boost(beta);
395	//
396	// --- Translate the dimuon parameters in the CM frame
397	//
398	TLorentzVector pDimuCM; //CM frame
399	TVector3 beta2;
400	beta2=(-1./(fMProton+pProjLab.E()))*pProjLab.Vect();
401	pDimuCM=pDimuLab;
402	pDimuCM.Boost(beta2);
403	//
404	// --- Determine the z axis for the calculation of the polarization angle
405	// (i.e. the direction of the dimuon in the CM system)
406	//
407	TVector3 zaxis;
408	zaxis=(pDimuCM.Vect()).Unit();
409	//
410	// --- Calculation of the polarization angle (Helicity)
411	// (angle between mu+ and the z axis defined above)
412	//
413	Double_t cost;
414	if(mu1Charge > 0) {
415	cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
416	// Theta Helicity is not properly defined for Like-Sign muons
417	if(mu2Charge > 0 && cost<0) cost=-cost;
418	} else {
419	cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
420	// Theta Helicity is not properly defined for Like-Sign muons
421	if(mu2Charge < 0 && cost<0) cost=-cost;
422	}
423	return cost;
424	}
425
426	//______________________________________________________________________________
427	Int_t AliAODDimuon::AnyPt(){
428	// Test if the two muons match two trigger tracks
429	if(this->CheckPointers())return 0;
430	return (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger())&&
431	(((AliAODTrack*)fMu[0].GetObject())->MatchTrigger());
432	}
433
434	//______________________________________________________________________________
435	Int_t AliAODDimuon::LowPt(){
436	// Test if the two muons match two trigger tracks with a "Low Pt" cut
437	if(this->CheckPointers())return 0;
438	return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt())&&
439	(((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt());
440	}
441
442	//______________________________________________________________________________
443	Int_t AliAODDimuon::HighPt(){
444	// Test if the two muons match two trigger tracks with a "High Pt" cut
445	if(this->CheckPointers())return 0;
446	return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt())&&
447	(((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt());
448	}
449
450	//______________________________________________________________________________
451	Double_t AliAODDimuon::MaxChi2Match(){
452	// Maximum matching Chi2 between track and trigger track
453	if(this->CheckPointers())return -999999999;
454	return TMath::Max((((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()),
455	(((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()));
456	}