Bug in the trigger definition corrected. (R. Arnaldi)
[u/mrichter/AliRoot.git] / STEER / AliAODDimuon.cxx
CommitLineData
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
21ClassImp(AliAODDimuon)
22
23//______________________________________________________________________________
24AliAODDimuon::AliAODDimuon():AliVParticle(),fP(0),fMProton(0.93827231)
25{
26 // default constructor
27 fMu[0]=0;
28 fMu[1]=0;
29}
30
31//______________________________________________________________________________
32AliAODDimuon::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//______________________________________________________________________________
42AliAODDimuon &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//______________________________________________________________________________
55AliAODDimuon::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//______________________________________________________________________________
65AliAODDimuon::~AliAODDimuon()
66{
67 // destructor
68 if(fP)delete fP;
69 fP=0;
70}
71
72//______________________________________________________________________________
73void 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//______________________________________________________________________________
82Double_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//______________________________________________________________________________
90Double_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//______________________________________________________________________________
98Double_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//______________________________________________________________________________
106Double_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//______________________________________________________________________________
116Double_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//______________________________________________________________________________
124Double_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//______________________________________________________________________________
131Double_t AliAODDimuon::P() {
132 // Dimuon momentum
133 if(this->CheckPointers())return -999999999;
134 BookP();
135 return fP->P();
136}
137
138//______________________________________________________________________________
139Double_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//______________________________________________________________________________
146Double_t AliAODDimuon::M() {
147 // Dimuon invariant mass
148 if(this->CheckPointers())return -999999999;
149 BookP();
150 return fP->M();
151}
152
153//______________________________________________________________________________
154Double_t AliAODDimuon::Mass() {
155 // Dimuon invariant mass
156 if(this->CheckPointers())return -999999999;
157 BookP();
158 return fP->M();
159}
160
161//______________________________________________________________________________
162Double_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//______________________________________________________________________________
169Double_t AliAODDimuon::Eta() {
170 // Dimuon pseudorapidity
171 if(this->CheckPointers())return -999999999;
172 BookP();
173 return fP->Eta();
174}
175
176//______________________________________________________________________________
177Double_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//______________________________________________________________________________
184Double_t AliAODDimuon::Phi() {
185 // Dimuon asimuthal angle
186 if(this->CheckPointers())return -999999999;
187 BookP();
188 return fP->Phi();
189}
190//______________________________________________________________________________
191Double_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//______________________________________________________________________________
198Double_t AliAODDimuon::Theta() {
199 // Dimuon polar angle
200 if(this->CheckPointers())return -999999999;
201 BookP();
202 return fP->Theta();
203}
204
205//______________________________________________________________________________
206Double_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//______________________________________________________________________________
213Double_t AliAODDimuon::Y() {
214 // Dimuon rapidity
215 if(this->CheckPointers())return -999999999;
216 BookP();
217 return fP->Rapidity();
218}
219
220//______________________________________________________________________________
221Short_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//______________________________________________________________________________
229Int_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){
48453d08 236 printf("Can not get objects\n");
866d8d78 237 return -999;
238 }
239 return 0;
240}
241
242//______________________________________________________________________________
243void 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//______________________________________________________________________________
251void AliAODDimuon::SetMuons(AliAODTrack *mu0, AliAODTrack *mu1){
252 // Assign the track pointers
253 fMu[0]=mu0;
254 fMu[1]=mu1;
255}
256
257//______________________________________________________________________________
258Double_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//______________________________________________________________________________
866d8d78 274Double_t AliAODDimuon::CostCS(){
48453d08 275 // Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
866d8d78 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
48453d08 298 TLorentzVector pProjCM(0.,0.,-pbeam,ebeam); // projectile
299 TLorentzVector pTargCM(0.,0., pbeam,ebeam); // target
866d8d78 300 //
48453d08 301 // --- Get the muons parameters in the CM frame
866d8d78 302 //
48453d08 303 TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
304 TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
866d8d78 305 //
48453d08 306 // --- Obtain the dimuon parameters in the CM frame
866d8d78 307 //
48453d08 308 TLorentzVector pDimuCM=pMu1CM+pMu2CM;
866d8d78 309 //
310 // --- Translate the dimuon parameters in the dimuon rest frame
311 //
48453d08 312 TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
313 TLorentzVector pMu1Dimu=pMu1CM;
314 TLorentzVector pMu2Dimu=pMu2CM;
315 TLorentzVector pProjDimu=pProjCM;
316 TLorentzVector pTargDimu=pTargCM;
866d8d78 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//______________________________________________________________________________
866d8d78 342Double_t AliAODDimuon::CostHe(){
48453d08 343 // Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
866d8d78 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 }
866d8d78 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();
866d8d78 360 //
48453d08 361 // --- Get the muons parameters in the CM frame
866d8d78 362 //
48453d08 363 TLorentzVector pMu1CM(pla10,pla11,pla12,e1);
364 TLorentzVector pMu2CM(pla20,pla21,pla22,e2);
866d8d78 365 //
48453d08 366 // --- Obtain the dimuon parameters in the CM frame
866d8d78 367 //
48453d08 368 TLorentzVector pDimuCM=pMu1CM+pMu2CM;
866d8d78 369 //
370 // --- Translate the dimuon parameters in the dimuon rest frame
371 //
48453d08 372 TVector3 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
373 TLorentzVector pMu1Dimu=pMu1CM;
374 TLorentzVector pMu2Dimu=pMu2CM;
866d8d78 375 pMu1Dimu.Boost(beta);
376 pMu2Dimu.Boost(beta);
377 //
866d8d78 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
48453d08 400//________________________________________________________________________
401Double_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//______________________________________________________________________________
465Double_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
866d8d78 528//______________________________________________________________________________
529Int_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())&&
8f6c57a4 533 (((AliAODTrack*)fMu[1].GetObject())->MatchTrigger());
866d8d78 534}
535
536//______________________________________________________________________________
537Int_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())&&
8f6c57a4 541 (((AliAODTrack*)fMu[1].GetObject())->MatchTriggerLowPt());
866d8d78 542}
543
544//______________________________________________________________________________
545Int_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())&&
8f6c57a4 549 (((AliAODTrack*)fMu[1].GetObject())->MatchTriggerHighPt());
866d8d78 550}
551
552//______________________________________________________________________________
553Double_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()),
8f6c57a4 557 (((AliAODTrack*)fMu[1].GetObject())->GetChi2MatchTrigger()));
866d8d78 558}