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