MUON/AliMUONPairLight.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *      SigmaEffect_thetadegrees                                                                  *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpeateose. It is      *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------------
  19 //This class was prepared by INFN Cagliari, July 2006
  20 //(authors: H.Woehri, A.de Falco)
  21 //
  22 // Compact information for the generated muon pairs in the MUON arm
  23 // useful at the last stage of the analysis chain
  24 // Pairs are built with two AliMUONTrackLight objects
  25 // Using the class AliMUONTrackLight this class combines the decay
  26 // information ("history") of the reconstructed tracks and fills
  27 // a series of flags for the formed reconstructed dimuon:
  28 // fIsCorrelated, fCreationProcess, fIsFeedDown, ...
  29 // for information about the dimuon, use PrintInfo with the appropriate
  30 // printflag
  31 // To be used together with AliMUONTrackLight
  32 //-----------------------------------------------------------------------------
  33
  34
  35 //MUON classes
  36 #include "AliMUONPairLight.h"
  37 //Root classes
  38 #include "TString.h"
  39
  40 ClassImp(AliMUONPairLight)
  41
  42 //====================================
  43 AliMUONPairLight::AliMUONPairLight() :
  44   TObject(),
  45   fMu0(),
  46   fMu1(),
  47   fCreationProcess(-999),
  48   fIsCorrelated(kFALSE),
  49   fCauseOfCorrelation (-1),
  50   fIsFeedDown(kFALSE)
  51 {
  52   /// default constructor
  53   ;
  54 }
  55
  56 //====================================
  57
  58 AliMUONPairLight::AliMUONPairLight(AliMUONPairLight &dimuCopy)
  59   : TObject(dimuCopy),
  60     fMu0(dimuCopy.fMu0),
  61     fMu1(dimuCopy.fMu1),
  62     fCreationProcess(dimuCopy.fCreationProcess),
  63     fIsCorrelated(dimuCopy.fIsCorrelated),
  64     fCauseOfCorrelation (dimuCopy.fCauseOfCorrelation),
  65     fIsFeedDown(dimuCopy.fIsFeedDown)
  66 {
  67 /// copy constructor
  68 ///   fMu0 = AliMUONTrackLight(dimuCopy.fMu0);
  69 ///   fMu1 = AliMUONTrackLight(dimuCopy.fMu1);
  70 ///   fIsCorrelated = dimuCopy.fIsCorrelated;
  71 ///   fCauseOfCorrelation = dimuCopy.fCauseOfCorrelation;
  72 ///   fCreationProcess = dimuCopy.fCreationProcess;
  73 ///   fIsFeedDown = dimuCopy.fIsFeedDown;
  74   ;
  75 }
  76
  77 //====================================
  78
  79 AliMUONPairLight::~AliMUONPairLight(){
  80   /// destructor
  81 }
  82
  83 //====================================
  84
  85 Bool_t AliMUONPairLight::IsAResonance(){
  86   /// checks if muon pair comes from a resonance decay
  87   if (!fIsCorrelated) return kFALSE;   //if muons not correlated, cannot be a resonance
  88   //if muons are correlated, check if the PDG of the
  89   //common mother is a resonance
  90   Int_t nparents0 = fMu0.GetNParents();
  91   Int_t nparents1 = fMu1.GetNParents();
  92
  93   Int_t minP = TMath::Min(nparents0, nparents1);
  94   for (Int_t i = 0 ; i < minP; i++) {
  95     if (fMu0.IsMotherAResonance(nparents0-1-i) && fMu1.IsMotherAResonance(nparents1-1-i) &&
  96         fMu0.GetParentPythiaLine(nparents0-1-i)==fMu1.GetParentPythiaLine(nparents1-1-i)) {
  97       if (nparents0-1-i) SetFeedDown(nparents0-1-i);
  98       return kTRUE;
  99     }
 100   }
 101   return kFALSE;
 102 }
 103
 104 //====================================
 105
 106 AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index)  {
 107   /// return muon 0 or 1
 108    if (index==0) return &fMu0;
 109    else if (index==1) return &fMu1;
 110    else{ printf ("Index can be either 0 or 1\n"); return 0;}
 111    //   else return &fMu1;
 112 }
 113
 114 //====================================
 115
 116 Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
 117   /// return muon mother pdg code
 118   if (imuon==0) return fMu0.GetParentPDGCode(mother);
 119   else if (imuon==1) return fMu1.GetParentPDGCode(mother);
 120   else { printf ("Index must be only 0 or 1\n"); return -999; }
 121 }
 122
 123 //====================================
 124 void AliMUONPairLight::SetProcess(){
 125   /// finds the process related to the muon pair (open charm/beauty, resonance,
 126   /// uncorrelated...)
 127
 128   AliMUONTrackLight *mu1 = &fMu0;
 129   AliMUONTrackLight *mu2 = &fMu1;
 130
 131   // check if the two muons are correlated
 132   // first check if they come from the same hadron (resonance or beauty/charm meson)
 133   Int_t npar1 = mu1->GetNParents();
 134   Int_t npar2 = mu2->GetNParents();
 135   for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) {
 136     Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
 137     for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) {
 138       Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
 139       if(lineMo1 == lineMo2) {
 140         //reject "diquark" mothers
 141         if(mu1->IsDiquark(mu1->GetParentPDGCode(imoth1)))return;
 142 //      if(IsDiquark(mu1->GetParentPDGCode(imoth1))) return;
 143         this->SetCorrelated(kTRUE);
 144         this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
 145         if(!IsAResonance()) fCreationProcess = 3;
 146         else fCreationProcess = -1;
 147         return;
 148       }
 149     }
 150   }
 151
 152   //now, check if we have a correlated pi/K:
 153   if(this->IsDimuonFromCorrPiK()){
 154     this->SetCorrelated(kTRUE);
 155     this->SetCauseOfCorrelation(mu1->GetParentPDGCode(0));
 156     fCreationProcess = -1;
 157   }
 158
 159   // if Open Beauty/Charm we can have 3 creation processes
 160   // (pair creation [0], gluon splitting [1] or flavour excitation [2])
 161   // 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same
 162   Int_t flavPar1 = mu1->GetParentFlavour(0);
 163   Int_t flavPar2 = mu2->GetParentFlavour(0);
 164   for (Int_t imoth1 = 0; imoth1 < 4; imoth1++) {
 165     Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
 166     for (Int_t imoth2 = 0; imoth2 < 4; imoth2++) {
 167       Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
 168       if(lineMo1 == lineMo2 && mu1->GetQuarkPDGCode(imoth1) == 21) {
 169         //now, check also that the string fragmented into two hadrons
 170         //of the same flavour (string usually splits into many hadrons
 171         //among which there are mostly soft particles)
 172         if(flavPar1 == flavPar2){
 173           this->SetCorrelated(kTRUE);
 174           if(GetCauseOfCorrelation() == -1)
 175             this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(imoth1));
 176
 177           fCreationProcess = 1;
 178           return;
 179         }
 180       }
 181     }
 182   }
 183
 184   Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
 185   Int_t line2 = mu2->GetQuarkPythiaLine(2);
 186
 187   Int_t line6or7[2] = {-1, -1}; //holds the index of quark in line 6 or 7
 188   Int_t flavourLine6or7[2] = {-1, -1};
 189   // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
 190   // are filled with a Q and Qbar
 191   for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
 192     Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
 193     Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
 194     if(lineMo1 == 6 || lineMo1 == 7){ //track 0 has a mother in line 6 or 7
 195       line6or7[0] = imoth1;
 196       flavourLine6or7[0] = flavour1;
 197     }
 198     for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
 199       Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
 200       Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
 201       if(lineMo2 == 6 || lineMo2 == 7){ //track 1 has a mother in line 6 or 7
 202         line6or7[1] = imoth2;
 203         flavourLine6or7[1] = flavour2;
 204       }
 205       if((line6or7[0] > 0 && line6or7[1] > 0) && //both tracks must have an entry in line 6 or 7
 206          (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5) && //this entry must be a c or b quark
 207          (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5) && // == " ==
 208          (flavPar1 == flavPar2)){ //make sure that the first hadronised parents of the 2 tracks are of the same flavour
 209         this->SetCorrelated(kTRUE);
 210         fCreationProcess = 0;
 211         return;
 212       }
 213     }
 214   }
 215
 216   // 3.)flavour excitation: if pythia line 6 of one track *and* pythia line 7 of second track
 217   // are filled with a Q and Qbar and if in addition there is another heavy quark in line(s) 4 and/or 5
 218   Int_t line2or3[2] = {-1, -1}; //holds the index of g/q in line 2 or 3
 219   Int_t flavourLine2or3[2] = {-1, -1};
 220   for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
 221     Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
 222     Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
 223     if(lineMo1 == 2 || lineMo1 == 3){ //track 0 has a mother in line 2 or 3
 224       line2or3[0] = imoth1;
 225       flavourLine2or3[0] = flavour1;
 226     }
 227     for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
 228       Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
 229       Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
 230       if(lineMo2 == 2 || lineMo2 == 3){ //track 1 has a mother in line 2 or 3
 231         line2or3[1] = imoth2;
 232         flavourLine2or3[1] = flavour2;
 233       }
 234       if(((line6or7[0] > 0 && (flavourLine6or7[0] == 4  || flavourLine6or7[0] == 5)) && //first track has Q in line 6 or 7
 235           (line2or3[1] > 0 && (flavourLine2or3[1] == 21 || flavourLine2or3[1] < 10))) || //second track has a g/q in line 2 or 3
 236          ((line6or7[1] > 0 && (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5)) &&  //or the same,
 237           (line2or3[0] > 0 && (flavourLine2or3[0] == 21 || flavourLine2or3[0] < 10)))){ // swapping the track's indices
 238         //now, check also that the string fragmented into two hadrons
 239         //of the same flavour (string usually splits into many hadrons
 240         //among which there are mostly soft particles)
 241         if(flavPar1 == flavPar2){
 242           this->SetCorrelated(kTRUE);
 243           fCreationProcess = 2;
 244           return;
 245         }
 246       }
 247     }
 248   }
 249
 250   //now flag (rare) processes in which only the incoming parton in line 2 or 3
 251   //radiates a gluon which produces a QQbar pair:
 252   //exclude the light quarks
 253   if(line1 == line2 && (line1 == 2 || line1 == 3)){
 254     if((TMath::Abs(mu1->GetQuarkPDGCode(1)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 4) ||
 255        (TMath::Abs(mu1->GetQuarkPDGCode(1)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 5)){
 256
 257       //now, check also that the string fragmented into two hadrons
 258       //of the same flavour (string usually splits into many hadrons
 259       //among which there are mostly soft particles)
 260       if(flavPar1 == flavPar2){
 261
 262         this->SetCorrelated(kTRUE);
 263         fCreationProcess = 1;
 264         if(GetCauseOfCorrelation() == -1){
 265           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
 266         }
 267         return;
 268       }
 269     }
 270   }
 271
 272   //in initial-state-radiation produced QQbar events the "mother quark"
 273   //is acknowledged as the second quark [1] and sits in line 2 or 3
 274   //is part of gluon splitting
 275   line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark of outgoing quark in [0]
 276   line2 = mu2->GetQuarkPythiaLine(1);
 277   if(line1 == line2 && (line1 == 2 || line1 == 3)){
 278     if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
 279        (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
 280
 281       //now, check also that the string fragmented into two hadrons
 282       //of the same flavour (string usually splits into many hadrons
 283       //among which there are mostly soft particles)
 284       if(flavPar1 == flavPar2){
 285
 286         this->SetCorrelated(kTRUE);
 287         fCreationProcess = 1;
 288         if(GetCauseOfCorrelation() == -1){
 289           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1)); //should be flagged as initial state radiation?
 290         }
 291         return;
 292       }
 293     }
 294   }
 295
 296   //in final-state-radiation produced QQbar events the "mother quark"
 297   //is acknowledged as the first quark [1] and sits in line 6 or 7
 298   //is part of gluon splitting
 299   line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark
 300   line2 = mu2->GetQuarkPythiaLine(1);
 301   if(line1 == line2 && (line1 == 6 || line1 == 7)){
 302     if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
 303        (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
 304
 305       //now, check also that the string fragmented into two hadrons
 306       //of the same flavour (string usually splits into many hadrons
 307       //among which there are mostly soft particles)
 308       if(flavPar1 == flavPar2){
 309
 310         this->SetCorrelated(kTRUE);
 311         fCreationProcess = 1;
 312         if(GetCauseOfCorrelation() == -1){
 313           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
 314         }
 315         return;
 316       }
 317     }
 318   }
 319 }
 320
 321 //====================================
 322 void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
 323   /// set the two muons
 324   fMu0 = mu0;
 325   fMu1 = mu1;
 326   this->SetProcess();
 327 }
 328
 329 //====================================
 330 void AliMUONPairLight::PrintInfo(Option_t* opt){
 331   /// print information about muon pairs
 332   /// Options:
 333   /// - "H" single muons' decay histories
 334   /// - "K" dimuon kinematics
 335   /// - "F" dimuon flags
 336   /// - "A" all variables
 337   TString options(opt);
 338   options.ToUpper();
 339
 340   if(options.Contains("H") || options.Contains("A")){//muon decay histories
 341
 342     AliMUONTrackLight *mu1 = &fMu0;
 343     AliMUONTrackLight *mu2 = &fMu1;
 344
 345     printf("========= History =======================\n");
 346     printf("first muon");
 347     mu1->PrintInfo("H");
 348     printf("second muon");
 349     mu2->PrintInfo("H");
 350     printf("=========================================\n");
 351   }
 352   if(options.Contains("F") || options.Contains("A")){//flags
 353     printf("the flags set for this muon pair are:\n");
 354     printf("=====================================\n");
 355     if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
 356     fIsCorrelated ? printf("(*) it is a correlated pair\n") : printf("(*) it is not a correlated pair\n");
 357     if(IsOpenCharm()) printf("(*) correlated open charm: ");
 358     if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
 359     if(IsOpenCharm() || IsOpenBeauty()){
 360       switch(fCreationProcess){
 361       case 0:
 362         printf("pair creation");
 363         break;
 364       case 1:
 365         printf("gluon splitting");
 366         break;
 367       case 2:
 368         printf("flavour excitation");
 369         break;
 370       case 3:
 371         printf("both muons come from same fragmented mother");
 372         break;
 373       }
 374       if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation())
 375         printf("... where oscillation occured\n");
 376       else{
 377         if(IsOpenBeauty())
 378           printf(" (no oscillation)\n");
 379         else
 380           printf("\n");
 381       }
 382     }
 383     IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(0, fIsFeedDown)) : printf("(*) it is not a resonance\n");
 384     fIsFeedDown ? printf("(*) mother has feed-down: %d --> %d\n", this->GetMuonMotherPDG(0,fMu0.GetNParents()-2), this->GetMuonMotherPDG(0,fMu0.GetNParents()-1)) : printf("(*) no feed-down\n");
 385     printf("=====================================\n");
 386   }
 387   if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
 388     Double_t *vtx = this->GetMuon(0)->GetVertex();
 389     TLorentzVector momRec = this->GetPRec();
 390     TLorentzVector momGen = this->GetPGen();
 391     printf("the dimuon charge is %d\n", this->GetCharge());
 392     printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
 393     printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
 394     printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
 395     //rapidity, pT, angles, ...
 396     printf("Rec. variables: mass %1.3f, pT %1.3f, pseudo-rapidity %1.3f, openingAngle %1.3f (%1.3f degree), theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
 397            momRec.M(), momRec.Pt(), momRec.Eta(),
 398            TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(),
 399            momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
 400            momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
 401   }
 402 }
 403
 404 //====================================
 405 Double_t AliMUONPairLight::GetOpeningAngle() {
 406   /// opening angle between the two muons in the lab frame (in degrees)
 407   TLorentzVector pRecMu0 =  fMu0.GetPRec();
 408   TLorentzVector pRecMu1 =  fMu1.GetPRec();
 409   TVector3 pRecMu03 = pRecMu0.Vect();
 410   TVector3 pRecMu13 = pRecMu1.Vect();
 411   Double_t scalar = pRecMu03.Dot(pRecMu13);
 412   Double_t modMu0 = pRecMu03.Mag();
 413   Double_t modMu1 = pRecMu13.Mag();
 414   Double_t theta = (TMath::ACos(scalar/(modMu0*modMu1)))*(180./TMath::Pi());
 415   return theta;
 416 }
 417 //================================================
 418 Bool_t AliMUONPairLight::IsDimuonFromCorrPiK(){
 419   ///check if we have a correlated pi/K
 420
 421   AliMUONTrackLight *mu0 = this->GetMuon(0), *mu1 = this->GetMuon(1);
 422   Bool_t fromSameLine = kFALSE;
 423   if (mu0->IsParentPionOrKaon() &&
 424       mu1->IsParentPionOrKaon() &&
 425       mu1->GetQuarkPythiaLine() == mu0->GetQuarkPythiaLine()
 426       ) fromSameLine = kTRUE;
 427
 428   return fromSameLine;
 429 }