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