Bug in the trigger definition corrected. (R. Arnaldi)
[u/mrichter/AliRoot.git] / STEER / 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. An instance of this class has also to be added to the AliAODEvent 
11 // structure to provide additional information that is specific to MUON and 
12 // therefore has not been included into the AOD header.
13 // Two transient data members are not stored on file as they can be recomputed
14 // at runtime.
15 //
16
17 #include "AliAODDimuon.h"
18 #include "TLorentzVector.h"
19 #define AliAODDimuon_CXX
20
21 ClassImp(AliAODDimuon)
22
23 //______________________________________________________________________________
24 AliAODDimuon::AliAODDimuon():AliVParticle(),fP(0),fMProton(0.93827231)
25 {
26   // default constructor
27   fMu[0]=0;
28   fMu[1]=0;
29 }
30
31 //______________________________________________________________________________
32 AliAODDimuon::AliAODDimuon(const AliAODDimuon& dimu):
33   AliVParticle(dimu),
34   fP(0),fMProton(0.93827231)
35 {
36   // copy constructor
37   fMu[0]=dimu.Mu(0);
38   fMu[1]=dimu.Mu(1);
39 }
40
41 //______________________________________________________________________________
42 AliAODDimuon &AliAODDimuon::operator=(const AliAODDimuon& dimu)
43 {
44   // assignment operator
45   if(&dimu != this){
46     fP=0;
47     fMProton=0.93827231;
48     fMu[0]=dimu.Mu(0);
49     fMu[1]=dimu.Mu(1);
50   }
51   return *this;
52 }
53
54 //______________________________________________________________________________
55 AliAODDimuon::AliAODDimuon(TObject *mu0, TObject *mu1):
56   fP(0),fMProton(0.93827231)
57 {
58   // Creates a dimuon pair from two tracks
59
60   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(){
74   //
75   // build the TLorentz vector
76   //
77     fP=new TLorentzVector(Px(),Py(),Pz(),E());
78     fP->SetPxPyPzE(Px(),Py(),Pz(),E());
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(px*px+py*py);
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\n");
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(ebeam*ebeam-mDimu*mDimu);
270   return Pz()/pMax;
271 }
272
273 //______________________________________________________________________________
274 Double_t AliAODDimuon::CostCS(){
275   // Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
276   if(CheckPointers())return -999999999;
277   Double_t ebeam=3500.;  //temporary
278   if(ebeam<=0){
279     printf("Can not compute costCS with EBeam=%f\n",ebeam);
280     return -999999999;
281   }
282   Double_t mp=fMProton;
283   Double_t pbeam=TMath::Sqrt(ebeam*ebeam-mp*mp);
284   Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
285   Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
286   Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
287   Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
288   Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
289   Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
290   Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
291   Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
292   Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
293   Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
294
295   // Fill the Lorentz vector for projectile and target
296   // For the moment we do not consider the crossing angle
297   // Projectile runs towards the MUON arm
298   TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
299   TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
300   //
301   // --- Get the muons parameters in the CM frame
302   //
303   TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
304   TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
305   //
306   // --- Obtain the dimuon parameters in the CM frame
307   //
308   TLorentzVector pDimuCM=pMu1CM+pMu2CM;
309   //
310   // --- Translate the dimuon parameters in the dimuon rest frame
311   //
312   TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
313   TLorentzVector pMu1Dimu=pMu1CM;
314   TLorentzVector pMu2Dimu=pMu2CM;
315   TLorentzVector pProjDimu=pProjCM;
316   TLorentzVector pTargDimu=pTargCM;
317   pMu1Dimu.Boost(beta);
318   pMu2Dimu.Boost(beta);
319   pProjDimu.Boost(beta);
320   pTargDimu.Boost(beta);
321   //
322   // --- Determine the z axis for the CS angle 
323   //
324   TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
325   //
326   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
327   //
328   Double_t cost;
329   if(mu1Charge > 0) {
330     cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
331     // Theta CS is not properly defined for Like-Sign muons
332     if(mu2Charge > 0 && cost<0) cost=-cost;
333   } else { 
334     // Theta CS is not properly defined for Like-Sign muons
335     cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
336     if(mu2Charge < 0 && cost<0) cost=-cost;
337   }
338   return cost;
339 }
340
341 //______________________________________________________________________________
342 Double_t AliAODDimuon::CostHe(){
343   // Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
344   if(CheckPointers())return -999999999;
345   Double_t ebeam=3500; //temporary
346   if(ebeam<=0){
347     printf("Can not compute costCS with EBeam=%f\n",ebeam);
348     return -999999999;
349   }
350   Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
351   Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
352   Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
353   Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
354   Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
355   Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
356   Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
357   Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
358   Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
359   Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
360   //
361   // --- Get the muons parameters in the CM frame
362   //
363   TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
364   TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
365   //
366   // --- Obtain the dimuon parameters in the CM frame
367   //
368   TLorentzVector pDimuCM=pMu1CM+pMu2CM;
369   //
370   // --- Translate the dimuon parameters in the dimuon rest frame
371   //
372   TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
373   TLorentzVector pMu1Dimu=pMu1CM;
374   TLorentzVector pMu2Dimu=pMu2CM;
375   pMu1Dimu.Boost(beta);
376   pMu2Dimu.Boost(beta);
377   //
378   // --- Determine the z axis for the calculation of the polarization angle
379   // (i.e. the direction of the dimuon in the CM system)
380   //
381   TVector3 zaxis;
382   zaxis=(pDimuCM.Vect()).Unit();
383   //
384   // --- Calculation of the polarization angle (Helicity)
385   // (angle between mu+ and the z axis defined above)
386   //
387   Double_t cost;
388   if(mu1Charge > 0) {
389     cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
390     // Theta Helicity is not properly defined for Like-Sign muons
391     if(mu2Charge > 0 && cost<0) cost=-cost;
392   } else { 
393     cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
394     // Theta Helicity is not properly defined for Like-Sign muons
395     if(mu2Charge < 0 && cost<0) cost=-cost;
396   }  
397   return cost;
398 }
399
400 //________________________________________________________________________
401 Double_t AliAODDimuon::PhiCS(){
402   // Cosinus of the Collins-Soper polar decay angle
403   if(CheckPointers())return -999999999;
404   Double_t ebeam=3500.;  //temporary
405   if(ebeam<=0){
406     printf("Can not compute phiCS with EBeam=%f\n",ebeam);
407     return -999999999;
408   }
409   Double_t mp=fMProton;
410   Double_t pbeam=TMath::Sqrt(ebeam*ebeam-mp*mp);
411   Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
412   Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
413   Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
414   Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
415   Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
416   Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
417   Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
418   Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
419   Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
420
421   // Fill the Lorentz vector for projectile and target
422   // For the moment we do not consider the crossing angle
423   // Projectile runs towards the MUON arm
424   TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
425   TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
426   //
427   // --- Get the muons parameters in the CM frame
428   //
429   TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
430   TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
431   //
432   // --- Obtain the dimuon parameters in the CM frame
433   //
434   TLorentzVector pDimuCM=pMu1CM+pMu2CM;
435   //
436   // --- Translate the dimuon parameters in the dimuon rest frame
437   //
438   TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
439   TLorentzVector pMu1Dimu=pMu1CM;
440   TLorentzVector pMu2Dimu=pMu2CM;
441   TLorentzVector pProjDimu=pProjCM;
442   TLorentzVector pTargDimu=pTargCM;
443   pMu1Dimu.Boost(beta);
444   pMu2Dimu.Boost(beta);
445   pProjDimu.Boost(beta);
446   pTargDimu.Boost(beta);
447   //
448   // --- Determine the z axis for the CS angle 
449   //
450   TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
451   //
452   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
453   //
454    TVector3 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
455    TVector3 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
456  
457    Double_t phi;
458    if(mu1Charge>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
459    else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
460      
461    return phi;
462 }
463
464 //______________________________________________________________________________
465 Double_t AliAODDimuon::PhiHe(){
466   // Calculation the Helicity aimuthal angle (adapted from code by R. Arnaldi)
467   if(CheckPointers())return -999999999;
468   Double_t ebeam=3500; //temporary
469   if(ebeam<=0){
470     printf("Can not compute phiHE with EBeam=%f\n",ebeam);
471     return -999999999;
472   }
473   Double_t pbeam=TMath::Sqrt(ebeam*ebeam-fMProton*fMProton);
474   Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
475   Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
476   Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
477   Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
478   Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
479   Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
480   Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
481   Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
482   Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
483
484   // Fill the Lorentz vector for projectile and target
485   // For the moment we consider no crossing angle
486   // Projectile runs towards the MUON arm
487   TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
488   TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
489   //
490   // --- Get the muons parameters in the CM frame
491   //
492   TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
493   TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
494   //
495   // --- Obtain the dimuon parameters in the CM frame
496   //
497   TLorentzVector pDimuCM=pMu1CM+pMu2CM;
498   //
499   // --- Translate the muon parameters in the dimuon rest frame
500   // 
501   TVector3 zaxis=(pDimuCM.Vect()).Unit(); 
502   //
503   // --- Translate the dimuon parameters in the dimuon rest frame
504   //
505   TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
506   TLorentzVector pMu1Dimu=pMu1CM;
507   TLorentzVector pMu2Dimu=pMu2CM;
508   pMu1Dimu.Boost(beta);
509   pMu2Dimu.Boost(beta);
510   
511   TLorentzVector pProjDimu=pProjCM;
512   TLorentzVector pTargDimu=pTargCM; 
513   pProjDimu.Boost(beta);
514   pTargDimu.Boost(beta);
515
516   TVector3 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
517   TVector3 xaxis=(yaxis.Cross(zaxis)).Unit();
518   //
519   // --- Calculation of the azimuthal angle (Helicity)
520   //
521    Double_t phi;
522    if(mu1Charge>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
523    else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
524    
525    return phi;
526 }
527
528 //______________________________________________________________________________
529 Int_t AliAODDimuon::AnyPt(){
530   // Test if the two muons match two trigger tracks
531   if(this->CheckPointers())return 0;
532   return (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger())&&
533          (((AliAODTrack*)fMu[1].GetObject())->MatchTrigger());
534 }
535
536 //______________________________________________________________________________
537 Int_t AliAODDimuon::LowPt(){
538   // Test if the two muons match two trigger tracks with a "Low Pt" cut
539   if(this->CheckPointers())return 0;
540   return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt())&&
541          (((AliAODTrack*)fMu[1].GetObject())->MatchTriggerLowPt());
542 }
543
544 //______________________________________________________________________________
545 Int_t AliAODDimuon::HighPt(){
546   // Test if the two muons match two trigger tracks with a "High Pt" cut
547   if(this->CheckPointers())return 0;
548   return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt())&&
549          (((AliAODTrack*)fMu[1].GetObject())->MatchTriggerHighPt());
550 }
551
552 //______________________________________________________________________________
553 Double_t AliAODDimuon::MaxChi2Match(){
554   // Maximum matching Chi2 between track and trigger track
555   if(this->CheckPointers())return -999999999;
556   return TMath::Max((((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()),
557                     (((AliAODTrack*)fMu[1].GetObject())->GetChi2MatchTrigger()));
558 }