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