[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

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

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

//______________________________________________________________________________
void AliAODDimuon::BookP(){
  //
  // build the TLorentz vector
  //
    fP=new TLorentzVector(Px(),Py(),Pz(),E());
    fP->SetPxPyPzE(Px(),Py(),Pz(),E());
}

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