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