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():AliVParticle(),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):
  34   AliVParticle(dimu),
  35   fEi(dimu.fEi),fP(0),fMProton(0.93827231)
  36 {
  37   // copy constructor
  38   fMu[0]=dimu.Mu(0);
  39   fMu[1]=dimu.Mu(1);
  40 }
  41
  42 //______________________________________________________________________________
  43 AliAODDimuon &AliAODDimuon::operator=(const AliAODDimuon& dimu)
  44 {
  45   // assignment operator
  46   if(&dimu != this){
  47     fP=0;
  48     fMProton=0.93827231;
  49     fMu[0]=dimu.Mu(0);
  50     fMu[1]=dimu.Mu(1);
  51     fEi=dimu.Ei();
  52   }
  53   return *this;
  54 }
  55
  56 //______________________________________________________________________________
  57 AliAODDimuon::AliAODDimuon(TObject *mu0, TObject *mu1, TObject *ei):
  58   fEi(ei),fP(0),fMProton(0.93827231)
  59 {
  60   // Creates a dimuon pair from two tracks and the EventInfo
  61
  62   //printf("Creating dimuon from %p %p\n",mu0,mu1);
  63   fMu[0]=mu0;
  64   fMu[1]=mu1;
  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 }