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. Got: %p %p\n",fMu[0].GetObject(),fMu[1].GetObject());
 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 // Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
 275 Double_t AliAODDimuon::CostCS(){
 276   // Cosinus of the Collins-Soper polar decay angle
 277   if(CheckPointers())return -999999999;
 278   Double_t ebeam=3500.;  //temporary
 279   if(ebeam<=0){
 280     printf("Can not compute costCS with EBeam=%f\n",ebeam);
 281     return -999999999;
 282   }
 283   Double_t mp=fMProton;
 284   Double_t pbeam=TMath::Sqrt(ebeam*ebeam-mp*mp);
 285   Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
 286   Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
 287   Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
 288   Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
 289   Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
 290   Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
 291   Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
 292   Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
 293   Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
 294   Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
 295
 296   // Fill the Lorentz vector for projectile and target
 297   // For the moment we do not consider the crossing angle
 298   // Projectile runs towards the MUON arm
 299   TLorentzVector pProjLab(0.,0.,-pbeam,ebeam); // projectile
 300   TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
 301   //
 302   // --- Get the muons parameters in the LAB frame
 303   //
 304   TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
 305   TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
 306   //
 307   // --- Obtain the dimuon parameters in the LAB frame
 308   //
 309   TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
 310   //
 311   // --- Translate the dimuon parameters in the dimuon rest frame
 312   //
 313   TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
 314   TLorentzVector pMu1Dimu=pMu1Lab;
 315   TLorentzVector pMu2Dimu=pMu2Lab;
 316   TLorentzVector pProjDimu=pProjLab;
 317   TLorentzVector pTargDimu=pTargLab;
 318   pMu1Dimu.Boost(beta);
 319   pMu2Dimu.Boost(beta);
 320   pProjDimu.Boost(beta);
 321   pTargDimu.Boost(beta);
 322   //
 323   // --- Determine the z axis for the CS angle
 324   //
 325   TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 326   //
 327   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 328   //
 329   Double_t cost;
 330   if(mu1Charge > 0) {
 331     cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
 332     // Theta CS is not properly defined for Like-Sign muons
 333     if(mu2Charge > 0 && cost<0) cost=-cost;
 334   } else {
 335     // Theta CS is not properly defined for Like-Sign muons
 336     cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
 337     if(mu2Charge < 0 && cost<0) cost=-cost;
 338   }
 339   return cost;
 340 }
 341
 342 //______________________________________________________________________________
 343 // Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
 344 Double_t AliAODDimuon::CostHe(){
 345   // Cosinus of the polar decay angle in the Helicity reference frame
 346   if(CheckPointers())return -999999999;
 347   Double_t ebeam=3500; //temporary
 348   if(ebeam<=0){
 349     printf("Can not compute costCS with EBeam=%f\n",ebeam);
 350     return -999999999;
 351   }
 352   Double_t pbeam=TMath::Sqrt(ebeam*ebeam-fMProton*fMProton);
 353   Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
 354   Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
 355   Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
 356   Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
 357   Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
 358   Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
 359   Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
 360   Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
 361   Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
 362   Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
 363
 364   // Fill the Lorentz vector for projectile and target
 365   // For the moment we consider no crossing angle
 366   // Projectile runs towards the MUON arm
 367   TLorentzVector pProjLab(0.,0.,-pbeam,ebeam); // projectile
 368   TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
 369   //
 370   // --- Get the muons parameters in the LAB frame
 371   //
 372   TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
 373   TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
 374   //
 375   // --- Obtain the dimuon parameters in the LAB frame
 376   //
 377   TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
 378   //
 379   // --- Translate the dimuon parameters in the dimuon rest frame
 380   //
 381   TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
 382   TLorentzVector pMu1Dimu=pMu1Lab;
 383   TLorentzVector pMu2Dimu=pMu2Lab;
 384   pMu1Dimu.Boost(beta);
 385   pMu2Dimu.Boost(beta);
 386   //
 387   // --- Translate the dimuon parameters in the CM frame
 388   //
 389   TLorentzVector pDimuCM; //CM frame
 390   TVector3 beta2;
 391   beta2=(-1./(fMProton+pProjLab.E()))*pProjLab.Vect();
 392   pDimuCM=pDimuLab;
 393   pDimuCM.Boost(beta2);
 394   //
 395   // --- Determine the z axis for the calculation of the polarization angle
 396   // (i.e. the direction of the dimuon in the CM system)
 397   //
 398   TVector3 zaxis;
 399   zaxis=(pDimuCM.Vect()).Unit();
 400   //
 401   // --- Calculation of the polarization angle (Helicity)
 402   // (angle between mu+ and the z axis defined above)
 403   //
 404   Double_t cost;
 405   if(mu1Charge > 0) {
 406     cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
 407     // Theta Helicity is not properly defined for Like-Sign muons
 408     if(mu2Charge > 0 && cost<0) cost=-cost;
 409   } else {
 410     cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
 411     // Theta Helicity is not properly defined for Like-Sign muons
 412     if(mu2Charge < 0 && cost<0) cost=-cost;
 413   }
 414   return cost;
 415 }
 416
 417 //______________________________________________________________________________
 418 Int_t AliAODDimuon::AnyPt(){
 419   // Test if the two muons match two trigger tracks
 420   if(this->CheckPointers())return 0;
 421   return (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger())&&
 422          (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger());
 423 }
 424
 425 //______________________________________________________________________________
 426 Int_t AliAODDimuon::LowPt(){
 427   // Test if the two muons match two trigger tracks with a "Low Pt" cut
 428   if(this->CheckPointers())return 0;
 429   return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt())&&
 430          (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt());
 431 }
 432
 433 //______________________________________________________________________________
 434 Int_t AliAODDimuon::HighPt(){
 435   // Test if the two muons match two trigger tracks with a "High Pt" cut
 436   if(this->CheckPointers())return 0;
 437   return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt())&&
 438          (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt());
 439 }
 440
 441 //______________________________________________________________________________
 442 Double_t AliAODDimuon::MaxChi2Match(){
 443   // Maximum matching Chi2 between track and trigger track
 444   if(this->CheckPointers())return -999999999;
 445   return TMath::Max((((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()),
 446                     (((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()));
 447 }