[u/mrichter/AliRoot.git] / PWG3 / 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
 */

#include "AliAODDimuon.h"
#define AliAODDimuon_CXX

ClassImp(AliAODDimuon)

AliAODDimuon::AliAODDimuon():ei(0),p(0),MProton(0.93827231)
{
  // default constructor
  //mu[0]=0;
  //mu[1]=0;
}

AliAODDimuon::AliAODDimuon(const AliAODDimuon& dimu):p(0),MProton(0.93827231)
{
  // default constructor
  mu[0]=dimu.mu[0];
  mu[1]=dimu.mu[1];
  ei=dimu.ei;
}

AliAODDimuon::AliAODDimuon(TObject *mu0, TObject *mu1, TObject *eipoint):p(0),MProton(0.93827231)
{
  ///printf("Creating dimuon from %p %p\n",mu0,mu1);
  mu[0]=mu0;
  mu[1]=mu1;
  ei=eipoint;
}

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

void AliAODDimuon::BookP(){
  static UInt_t UnID[2]={0,0};
  if(!p){
    p=new TLorentzVector(Px(),Py(),Pz(),E());
    UnID[0]=mu[0].GetUniqueID();
    UnID[1]=mu[1].GetUniqueID();
  }
  // For efficiency reasons
  if((UnID[0]!=mu[0].GetUniqueID())||(UnID[1]!=mu[1].GetUniqueID())){
    p->SetPxPyPzE(Px(),Py(),Pz(),E());
    UnID[0]=mu[0].GetUniqueID();
    UnID[1]=mu[1].GetUniqueID();
  }
}

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

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

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

//______________________________________________________________________________
Double_t AliAODDimuon::Pt() const {
  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 {
  if(this->CheckPointers())return -999999999;
  return ((AliAODTrack*)mu[0].GetObject())->E()+((AliAODTrack*)mu[1].GetObject())->E();
}

//______________________________________________________________________________
Double_t AliAODDimuon::P() const {
  printf("You should never call: Double_t AliAODDimuon::P() const\n");
  return -999999999;
}

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

//______________________________________________________________________________
Double_t AliAODDimuon::M() const {
  printf("You should never call: Double_t AliAODDimuon::M() const\n");
  return -999999999;
}

//______________________________________________________________________________
Double_t AliAODDimuon::M() {
  if(this->CheckPointers())return -999999999;
  BookP();
  return p->M();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Eta() const {
  printf("You should never call: Double_t AliAODDimuon::Eta() const\n");
  return -999999999;
}

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

//______________________________________________________________________________
Double_t AliAODDimuon::Phi() const {
  printf("You should never call: Double_t AliAODDimuon::Phi() const\n");
  return -999999999;
}

//______________________________________________________________________________
Double_t AliAODDimuon::Phi() {
  if(this->CheckPointers())return -999999999;
  BookP();
  return p->Phi();
}
//______________________________________________________________________________
Double_t AliAODDimuon::Theta() const {
  printf("You should never call: Double_t AliAODDimuon::Theta() const\n");
  return -999999999;
}

//______________________________________________________________________________
Double_t AliAODDimuon::Theta() {
  if(this->CheckPointers())return -999999999;
  BookP();
  return p->Theta();
}

//______________________________________________________________________________
Double_t AliAODDimuon::Y() const {
  printf("You should never call: Double_t AliAODDimuon::Y() const\n");
  return -999999999;
}

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

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

//______________________________________________________________________________
Int_t AliAODDimuon::CheckPointers() const{
  if(mu[0]==0||mu[1]==0){
    printf("Dimuon not initialized\n");
    return -999;
  }
  if((mu[0].GetObject())==0||(mu[1].GetObject())==0){
    printf("Can not get objects. Got: %p %p\n",mu[0].GetObject(),mu[1].GetObject());
    return -999;
  }
  return 0;
}

//______________________________________________________________________________
Double_t AliAODDimuon::xf() {
  Double_t EBeam=((AliAODEventInfo*)ei.GetObject())->EBeam();
  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(){
  if(CheckPointers())return -999999999;
  if(ei==0){
    printf("Pointer to MuonHeader not initialized\n");
    return -999999999;
  }
  if(ei.GetObject()==0){
    printf("Can not get MuonHeader object\n");
    return -999999999;
  }
  Double_t EBeam=((AliAODEventInfo*)ei.GetObject())->EBeam();
  if(EBeam<=0){
    printf("Can not compute costCS with EBeam=%f\n",EBeam);
    return -999999999;
  }
  Double_t mp=MProton;
  Double_t PBeam=TMath::Sqrt(EBeam*EBeam-mp*mp);
  Double_t pla10=((AliAODTrack*)mu[0].GetObject())->Px();
  Double_t pla11=((AliAODTrack*)mu[0].GetObject())->Py();
  Double_t pla12=((AliAODTrack*)mu[0].GetObject())->Pz();
  Double_t e1=((AliAODTrack*)mu[0].GetObject())->E();
  Double_t Mu1Charge=((AliAODTrack*)mu[0].GetObject())->Charge();
  Double_t pla20=((AliAODTrack*)mu[1].GetObject())->Px();
  Double_t pla21=((AliAODTrack*)mu[1].GetObject())->Py();
  Double_t pla22=((AliAODTrack*)mu[1].GetObject())->Pz();
  Double_t e2=((AliAODTrack*)mu[1].GetObject())->E();
  Double_t Mu2Charge=((AliAODTrack*)mu[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;
  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(){
  if(CheckPointers())return -999999999;
  if(ei==0){
    printf("Pointer to MuonHeader not initialized\n");
    return -999999999;
  }
  if(ei.GetObject()==0){
    printf("Can not get MuonHeader object\n");
    return -999999999;
  }
  Double_t EBeam=((AliAODEventInfo*)ei.GetObject())->EBeam();
  if(EBeam<=0){
    printf("Can not compute costCS with EBeam=%f\n",EBeam);
    return -999999999;
  }
  Double_t PBeam=TMath::Sqrt(EBeam*EBeam-MProton*MProton);
  Double_t pla10=((AliAODTrack*)mu[0].GetObject())->Px();
  Double_t pla11=((AliAODTrack*)mu[0].GetObject())->Py();
  Double_t pla12=((AliAODTrack*)mu[0].GetObject())->Pz();
  Double_t e1=((AliAODTrack*)mu[0].GetObject())->E();
  Double_t Mu1Charge=((AliAODTrack*)mu[0].GetObject())->Charge();
  Double_t pla20=((AliAODTrack*)mu[1].GetObject())->Px();
  Double_t pla21=((AliAODTrack*)mu[1].GetObject())->Py();
  Double_t pla22=((AliAODTrack*)mu[1].GetObject())->Pz();
  Double_t e2=((AliAODTrack*)mu[1].GetObject())->E();
  Double_t Mu2Charge=((AliAODTrack*)mu[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./(MProton+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 (Kharzeev) (angle between mu+ and the z axis defined above)
  //
  Double_t cost;
  if(Mu1Charge > 0) {
    cost = zaxis.Dot((PMu1Dimu.Vect()).Unit());
    // Theta Kharzeev is not properly defined for Like-Sign muons
    if(Mu2Charge > 0 && cost<0) cost=-cost;
  } else { 
    cost = zaxis.Dot((PMu2Dimu.Vect()).Unit());
    // Theta Kharzeev is not properly defined for Like-Sign muons
    if(Mu2Charge < 0 && cost<0) cost=-cost;
  }  
  return cost;
}

Int_t AliAODDimuon::AnyPt(){
  if(this->CheckPointers())return 0;
  return (((AliAODTrack*)mu[0].GetObject())->MatchTriggerAnyPt())&&
         (((AliAODTrack*)mu[0].GetObject())->MatchTriggerAnyPt());
}

Int_t AliAODDimuon::LowPt(){
  if(this->CheckPointers())return 0;
  return (((AliAODTrack*)mu[0].GetObject())->MatchTriggerLowPt())&&
         (((AliAODTrack*)mu[0].GetObject())->MatchTriggerLowPt());
}

Int_t AliAODDimuon::HighPt(){
  if(this->CheckPointers())return 0;
  return (((AliAODTrack*)mu[0].GetObject())->MatchTriggerHighPt())&&
         (((AliAODTrack*)mu[0].GetObject())->MatchTriggerHighPt());
}

Double_t AliAODDimuon::MaxChi2Match(){
  if(this->CheckPointers())return -999999999;
  return TMath::Max((((AliAODTrack*)mu[0].GetObject())->GetChi2MatchTrigger()),
                    (((AliAODTrack*)mu[0].GetObject())->GetChi2MatchTrigger()));
}
Commit	Line	Data
71b7d225	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
	6	#include "AliAODDimuon.h"
	7	#define AliAODDimuon_CXX
	8
	9	ClassImp(AliAODDimuon)
	10
	11	AliAODDimuon::AliAODDimuon():ei(0),p(0),MProton(0.93827231)
	12	{
	13	// default constructor
	14	//mu[0]=0;
	15	//mu[1]=0;
	16	}
	17
	18	AliAODDimuon::AliAODDimuon(const AliAODDimuon& dimu):p(0),MProton(0.93827231)
	19	{
	20	// default constructor
	21	mu[0]=dimu.mu[0];
	22	mu[1]=dimu.mu[1];
	23	ei=dimu.ei;
	24	}
	25
	26	AliAODDimuon::AliAODDimuon(TObject mu0, TObject mu1, TObject *eipoint):p(0),MProton(0.93827231)
	27	{
	28	///printf("Creating dimuon from %p %p\n",mu0,mu1);
	29	mu[0]=mu0;
	30	mu[1]=mu1;
	31	ei=eipoint;
	32	}
	33
	34	//______________________________________________________________________________
	35	AliAODDimuon::~AliAODDimuon()
	36	{
	37	// destructor
	38	if(p)delete p;
	39	p=0;
	40	}
	41
	42	void AliAODDimuon::BookP(){
	43	static UInt_t UnID[2]={0,0};
	44	if(!p){
	45	p=new TLorentzVector(Px(),Py(),Pz(),E());
	46	UnID[0]=mu[0].GetUniqueID();
	47	UnID[1]=mu[1].GetUniqueID();
	48	}
	49	// For efficiency reasons
	50	if((UnID[0]!=mu[0].GetUniqueID())\|\|(UnID[1]!=mu[1].GetUniqueID())){
	51	p->SetPxPyPzE(Px(),Py(),Pz(),E());
	52	UnID[0]=mu[0].GetUniqueID();
	53	UnID[1]=mu[1].GetUniqueID();
	54	}
	55	}
	56
	57	//______________________________________________________________________________
	58	Double_t AliAODDimuon::Px() const {
	59	if(this->CheckPointers())return -999999999;
	60	return ((AliAODTrack)mu[0].GetObject())->Px()+((AliAODTrack)mu[1].GetObject())->Px();
	61	}
	62
	63	//______________________________________________________________________________
	64	Double_t AliAODDimuon::Py() const {
65	if(this->CheckPointers())return -999999999;
66	return ((AliAODTrack)mu[0].GetObject())->Py()+((AliAODTrack)mu[1].GetObject())->Py();
67	}
68
69	//______________________________________________________________________________
70	Double_t AliAODDimuon::Pz() const {
71	if(this->CheckPointers())return -999999999;
72	return ((AliAODTrack)mu[0].GetObject())->Pz()+((AliAODTrack)mu[1].GetObject())->Pz();
73	}
74
75	//______________________________________________________________________________
76	Double_t AliAODDimuon::Pt() const {
77	if(this->CheckPointers())return -999999999;
78	Double_t px=Px();
79	Double_t py=Py();
80	return TMath::Sqrt(pxpx+pypy);
81	return -999999999;
82	}
83
84	//______________________________________________________________________________
85	Double_t AliAODDimuon::E() const {
86	if(this->CheckPointers())return -999999999;
87	return ((AliAODTrack)mu[0].GetObject())->E()+((AliAODTrack)mu[1].GetObject())->E();
88	}
89
90	//______________________________________________________________________________
91	Double_t AliAODDimuon::P() const {
92	printf("You should never call: Double_t AliAODDimuon::P() const\n");
93	return -999999999;
94	}
95
96	//______________________________________________________________________________
97	Double_t AliAODDimuon::P() {
98	if(this->CheckPointers())return -999999999;
99	BookP();
100	return p->P();
101	}
102
103	//______________________________________________________________________________
104	Double_t AliAODDimuon::M() const {
105	printf("You should never call: Double_t AliAODDimuon::M() const\n");
106	return -999999999;
107	}
108
109	//______________________________________________________________________________
110	Double_t AliAODDimuon::M() {
111	if(this->CheckPointers())return -999999999;
112	BookP();
113	return p->M();
114	}
115
116	//______________________________________________________________________________
117	Double_t AliAODDimuon::Eta() const {
118	printf("You should never call: Double_t AliAODDimuon::Eta() const\n");
119	return -999999999;
120	}
121
122	//______________________________________________________________________________
123	Double_t AliAODDimuon::Eta() {
124	if(this->CheckPointers())return -999999999;
125	BookP();
126	return p->Eta();
127	}
128
129	//______________________________________________________________________________
130	Double_t AliAODDimuon::Phi() const {
131	printf("You should never call: Double_t AliAODDimuon::Phi() const\n");
132	return -999999999;
133	}
134
135	//______________________________________________________________________________
136	Double_t AliAODDimuon::Phi() {
137	if(this->CheckPointers())return -999999999;
138	BookP();
139	return p->Phi();
140	}
141	//______________________________________________________________________________
142	Double_t AliAODDimuon::Theta() const {
143	printf("You should never call: Double_t AliAODDimuon::Theta() const\n");
144	return -999999999;
145	}
146
147	//______________________________________________________________________________
148	Double_t AliAODDimuon::Theta() {
149	if(this->CheckPointers())return -999999999;
150	BookP();
151	return p->Theta();
152	}
153
154	//______________________________________________________________________________
155	Double_t AliAODDimuon::Y() const {
156	printf("You should never call: Double_t AliAODDimuon::Y() const\n");
157	return -999999999;
158	}
159
160	//______________________________________________________________________________
161	Double_t AliAODDimuon::Y() {
162	if(this->CheckPointers())return -999999999;
163	BookP();
164	return p->Rapidity();
165	}
166
167	//______________________________________________________________________________
168	Short_t AliAODDimuon::Charge() const
169	{
170	if(this->CheckPointers())return -999;
171	return ((AliAODTrack)mu[0].GetObject())->Charge()+((AliAODTrack)mu[1].GetObject())->Charge();
172	}
173
174	//______________________________________________________________________________
175	Int_t AliAODDimuon::CheckPointers() const{
176	if(mu[0]==0\|\|mu[1]==0){
177	printf("Dimuon not initialized\n");
178	return -999;
179	}
180	if((mu[0].GetObject())==0\|\|(mu[1].GetObject())==0){
181	printf("Can not get objects. Got: %p %p\n",mu[0].GetObject(),mu[1].GetObject());
182	return -999;
183	}
184	return 0;
185	}
186
187	//______________________________________________________________________________
188	Double_t AliAODDimuon::xf() {
189	Double_t EBeam=((AliAODEventInfo*)ei.GetObject())->EBeam();
190	if(EBeam<=0){
191	printf("AliAODDimuon::xf: can not compute xf with EBeam=%f\n",EBeam);
192	return -999999999;
193	}
194	if(this->CheckPointers())return -999999999;
195	BookP();
196	Double_t MDimu=M();
197	Double_t PMax=TMath::Sqrt(EBeamEBeam-MDimuMDimu);
198	return Pz()/PMax;
199	}
200
201	//______________________________________________________________________________
202	// Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
203	Double_t AliAODDimuon::CostCS(){
204	if(CheckPointers())return -999999999;
205	if(ei==0){
206	printf("Pointer to MuonHeader not initialized\n");
207	return -999999999;
208	}
209	if(ei.GetObject()==0){
210	printf("Can not get MuonHeader object\n");
211	return -999999999;
212	}
213	Double_t EBeam=((AliAODEventInfo*)ei.GetObject())->EBeam();
214	if(EBeam<=0){
215	printf("Can not compute costCS with EBeam=%f\n",EBeam);
216	return -999999999;
217	}
218	Double_t mp=MProton;
219	Double_t PBeam=TMath::Sqrt(EBeamEBeam-mpmp);
220	Double_t pla10=((AliAODTrack*)mu[0].GetObject())->Px();
221	Double_t pla11=((AliAODTrack*)mu[0].GetObject())->Py();
222	Double_t pla12=((AliAODTrack*)mu[0].GetObject())->Pz();
223	Double_t e1=((AliAODTrack*)mu[0].GetObject())->E();
224	Double_t Mu1Charge=((AliAODTrack*)mu[0].GetObject())->Charge();
225	Double_t pla20=((AliAODTrack*)mu[1].GetObject())->Px();
226	Double_t pla21=((AliAODTrack*)mu[1].GetObject())->Py();
227	Double_t pla22=((AliAODTrack*)mu[1].GetObject())->Pz();
228	Double_t e2=((AliAODTrack*)mu[1].GetObject())->E();
229	Double_t Mu2Charge=((AliAODTrack*)mu[1].GetObject())->Charge();
230
231	// Fill the Lorentz vector for projectile and target
232	// For the moment we consider no crossing angle
233	// Projectile runs towards the MUON arm
234	TLorentzVector PProjLab(0.,0.,-PBeam,EBeam); // projectile
235	TLorentzVector PTargLab(0.,0., PBeam,EBeam); // target
236	//
237	// --- Get the muons parameters in the LAB frame
238	//
239	TLorentzVector PMu1Lab(pla10,pla11,pla12,e1);
240	TLorentzVector PMu2Lab(pla20,pla21,pla22,e2);
241	//
242	// --- Obtain the dimuon parameters in the LAB frame
243	//
244	TLorentzVector PDimuLab=PMu1Lab+PMu2Lab;
245	//
246	// --- Translate the dimuon parameters in the dimuon rest frame
247	//
248	TVector3 beta=(-1./PDimuLab.E())*PDimuLab.Vect();
249	TLorentzVector PMu1Dimu=PMu1Lab;
250	TLorentzVector PMu2Dimu=PMu2Lab;
251	TLorentzVector PProjDimu=PProjLab;
252	TLorentzVector PTargDimu=PTargLab;
253	PMu1Dimu.Boost(beta);
254	PMu2Dimu.Boost(beta);
255	PProjDimu.Boost(beta);
256	PTargDimu.Boost(beta);
257	//
258	// --- Determine the z axis for the CS angle
259	//
260	TVector3 zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
261	//
262	// --- Determine the CS angle (angle between mu+ and the z axis defined above)
263	//
264	Double_t cost;
265	if(Mu1Charge > 0) {
266	cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());
267	// Theta CS is not properly defined for Like-Sign muons
268	if(Mu2Charge > 0 && cost<0) cost=-cost;
269	} else {
270	// Theta CS is not properly defined for Like-Sign muons
271	cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());
272	if(Mu2Charge < 0 && cost<0) cost=-cost;
273	}
274	return cost;
275	}
276
277	//______________________________________________________________________________
278	// Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
279	Double_t AliAODDimuon::CostHe(){
280	if(CheckPointers())return -999999999;
281	if(ei==0){
282	printf("Pointer to MuonHeader not initialized\n");
283	return -999999999;
284	}
285	if(ei.GetObject()==0){
286	printf("Can not get MuonHeader object\n");
287	return -999999999;
288	}
289	Double_t EBeam=((AliAODEventInfo*)ei.GetObject())->EBeam();
290	if(EBeam<=0){
291	printf("Can not compute costCS with EBeam=%f\n",EBeam);
292	return -999999999;
293	}
294	Double_t PBeam=TMath::Sqrt(EBeamEBeam-MProtonMProton);
295	Double_t pla10=((AliAODTrack*)mu[0].GetObject())->Px();
296	Double_t pla11=((AliAODTrack*)mu[0].GetObject())->Py();
297	Double_t pla12=((AliAODTrack*)mu[0].GetObject())->Pz();
298	Double_t e1=((AliAODTrack*)mu[0].GetObject())->E();
299	Double_t Mu1Charge=((AliAODTrack*)mu[0].GetObject())->Charge();
300	Double_t pla20=((AliAODTrack*)mu[1].GetObject())->Px();
301	Double_t pla21=((AliAODTrack*)mu[1].GetObject())->Py();
302	Double_t pla22=((AliAODTrack*)mu[1].GetObject())->Pz();
303	Double_t e2=((AliAODTrack*)mu[1].GetObject())->E();
304	Double_t Mu2Charge=((AliAODTrack*)mu[1].GetObject())->Charge();
305
306	// Fill the Lorentz vector for projectile and target
307	// For the moment we consider no crossing angle
308	// Projectile runs towards the MUON arm
309	TLorentzVector PProjLab(0.,0.,-PBeam,EBeam); // projectile
310	TLorentzVector PTargLab(0.,0., PBeam,EBeam); // target
311	//
312	// --- Get the muons parameters in the LAB frame
313	//
314	TLorentzVector PMu1Lab(pla10,pla11,pla12,e1);
315	TLorentzVector PMu2Lab(pla20,pla21,pla22,e2);
316	//
317	// --- Obtain the dimuon parameters in the LAB frame
318	//
319	TLorentzVector PDimuLab=PMu1Lab+PMu2Lab;
320	//
321	// --- Translate the dimuon parameters in the dimuon rest frame
322	//
323	TVector3 beta=(-1./PDimuLab.E())*PDimuLab.Vect();
324	TLorentzVector PMu1Dimu=PMu1Lab;
325	TLorentzVector PMu2Dimu=PMu2Lab;
326	PMu1Dimu.Boost(beta);
327	PMu2Dimu.Boost(beta);
328	//
329	// --- Translate the dimuon parameters in the CM frame
330	//
331	TLorentzVector PDimuCM; //CM frame
332	TVector3 beta2;
333	beta2=(-1./(MProton+PProjLab.E()))*PProjLab.Vect();
334	PDimuCM=PDimuLab;
335	PDimuCM.Boost(beta2);
336	//
337	// --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
338	//
339	TVector3 zaxis;
340	zaxis=(PDimuCM.Vect()).Unit();
341	//
342	// --- Calculation of the polarization angle (Kharzeev) (angle between mu+ and the z axis defined above)
343	//
344	Double_t cost;
345	if(Mu1Charge > 0) {
346	cost = zaxis.Dot((PMu1Dimu.Vect()).Unit());
347	// Theta Kharzeev is not properly defined for Like-Sign muons
348	if(Mu2Charge > 0 && cost<0) cost=-cost;
349	} else {
350	cost = zaxis.Dot((PMu2Dimu.Vect()).Unit());
351	// Theta Kharzeev is not properly defined for Like-Sign muons
352	if(Mu2Charge < 0 && cost<0) cost=-cost;
353	}
354	return cost;
355	}
356
357	Int_t AliAODDimuon::AnyPt(){
358	if(this->CheckPointers())return 0;
359	return (((AliAODTrack*)mu[0].GetObject())->MatchTriggerAnyPt())&&
360	(((AliAODTrack*)mu[0].GetObject())->MatchTriggerAnyPt());
361	}
362
363	Int_t AliAODDimuon::LowPt(){
364	if(this->CheckPointers())return 0;
365	return (((AliAODTrack*)mu[0].GetObject())->MatchTriggerLowPt())&&
366	(((AliAODTrack*)mu[0].GetObject())->MatchTriggerLowPt());
367	}
368
369	Int_t AliAODDimuon::HighPt(){
370	if(this->CheckPointers())return 0;
371	return (((AliAODTrack*)mu[0].GetObject())->MatchTriggerHighPt())&&
372	(((AliAODTrack*)mu[0].GetObject())->MatchTriggerHighPt());
373	}
374
375	Double_t AliAODDimuon::MaxChi2Match(){
376	if(this->CheckPointers())return -999999999;
377	return TMath::Max((((AliAODTrack*)mu[0].GetObject())->GetChi2MatchTrigger()),
378	(((AliAODTrack*)mu[0].GetObject())->GetChi2MatchTrigger()));
379	}