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