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 AliMUONPairLight& AliMUONPairLight::operator=(const AliMUONPairLight& dimuCopy)
  86 {
  87   // check assignment to self
  88   if (this == &dimuCopy) return *this;
  89
  90   // base class assignment
  91   TObject::operator=(dimuCopy);
  92
  93   // assignment operator
  94   fMu0 = dimuCopy.fMu0;
  95   fMu1 = dimuCopy.fMu1;
  96   fCreationProcess = dimuCopy.fCreationProcess;
  97   fIsCorrelated = dimuCopy.fIsCorrelated;
  98   fCauseOfCorrelation  = dimuCopy.fCauseOfCorrelation;
  99   fIsFeedDown = dimuCopy.fIsFeedDown;
 100
 101   return *this;
 102 }
 103
 104 //====================================
 105
 106 Bool_t AliMUONPairLight::IsAResonance(){
 107   /// checks if muon pair comes from a resonance decay
 108   if (!fIsCorrelated) return kFALSE;   //if muons not correlated, cannot be a resonance
 109   //if muons are correlated, check if the PDG of the
 110   //common mother is a resonance
 111   Int_t nparents0 = fMu0.GetNParents();
 112   Int_t nparents1 = fMu1.GetNParents();
 113
 114   Int_t minP = TMath::Min(nparents0, nparents1);
 115   for (Int_t i = 0 ; i < minP; i++) {
 116     if (fMu0.IsMotherAResonance(nparents0-1-i) && fMu1.IsMotherAResonance(nparents1-1-i) &&
 117         fMu0.GetParentPythiaLine(nparents0-1-i)==fMu1.GetParentPythiaLine(nparents1-1-i)) {
 118       if (nparents0-1-i) SetFeedDown(nparents0-1-i);
 119       return kTRUE;
 120     }
 121   }
 122   return kFALSE;
 123 }
 124
 125 //====================================
 126
 127 AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index)  {
 128   /// return muon 0 or 1
 129    if (index==0) return &fMu0;
 130    else if (index==1) return &fMu1;
 131    else{ printf ("Index can be either 0 or 1\n"); return 0;}
 132    //   else return &fMu1;
 133 }
 134
 135 //====================================
 136
 137 Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
 138   /// return muon mother pdg code
 139   if (imuon==0) return fMu0.GetParentPDGCode(mother);
 140   else if (imuon==1) return fMu1.GetParentPDGCode(mother);
 141   else { printf ("Index must be only 0 or 1\n"); return -999; }
 142 }
 143
 144 //====================================
 145 void AliMUONPairLight::SetProcess(){
 146   /// finds the process related to the muon pair (open charm/beauty, resonance,
 147   /// uncorrelated...)
 148
 149   AliMUONTrackLight *mu1 = &fMu0;
 150   AliMUONTrackLight *mu2 = &fMu1;
 151
 152   // check if the two muons are correlated
 153   // first check if they come from the same hadron (resonance or beauty/charm meson)
 154   Int_t npar1 = mu1->GetNParents();
 155   Int_t npar2 = mu2->GetNParents();
 156   for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) {
 157     Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
 158     for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) {
 159       Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
 160       if(lineMo1 == lineMo2) {
 161         //reject "diquark" mothers
 162         if(mu1->IsDiquark(mu1->GetParentPDGCode(imoth1)))return;
 163 //      if(IsDiquark(mu1->GetParentPDGCode(imoth1))) return;
 164         this->SetCorrelated(kTRUE);
 165         this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
 166         if(!IsAResonance()) fCreationProcess = 3;
 167         else fCreationProcess = -1;
 168         return;
 169       }
 170     }
 171   }
 172
 173   //now, check if we have a correlated pi/K:
 174   if(this->IsDimuonFromCorrPiK()){
 175     this->SetCorrelated(kTRUE);
 176     this->SetCauseOfCorrelation(mu1->GetParentPDGCode(0));
 177     fCreationProcess = -1;
 178   }
 179
 180   // if Open Beauty/Charm we can have 3 creation processes
 181   // (pair creation [0], gluon splitting [1] or flavour excitation [2])
 182   // 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same
 183   Int_t flavPar1 = mu1->GetParentFlavour(0);
 184   Int_t flavPar2 = mu2->GetParentFlavour(0);
 185   for (Int_t imoth1 = 0; imoth1 < 4; imoth1++) {
 186     Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
 187     for (Int_t imoth2 = 0; imoth2 < 4; imoth2++) {
 188       Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
 189       if(lineMo1 == lineMo2 && mu1->GetQuarkPDGCode(imoth1) == 21) {
 190         //now, check also that the string fragmented into two hadrons
 191         //of the same flavour (string usually splits into many hadrons
 192         //among which there are mostly soft particles)
 193         if(flavPar1 == flavPar2){
 194           this->SetCorrelated(kTRUE);
 195           if(GetCauseOfCorrelation() == -1)
 196             this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(imoth1));
 197
 198           fCreationProcess = 1;
 199           return;
 200         }
 201       }
 202     }
 203   }
 204
 205   Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
 206   Int_t line2 = mu2->GetQuarkPythiaLine(2);
 207
 208   Int_t line6or7[2] = {-1, -1}; //holds the index of quark in line 6 or 7
 209   Int_t flavourLine6or7[2] = {-1, -1};
 210   // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
 211   // are filled with a Q and Qbar
 212   for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
 213     Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
 214     Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
 215     if(lineMo1 == 6 || lineMo1 == 7){ //track 0 has a mother in line 6 or 7
 216       line6or7[0] = imoth1;
 217       flavourLine6or7[0] = flavour1;
 218     }
 219     for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
 220       Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
 221       Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
 222       if(lineMo2 == 6 || lineMo2 == 7){ //track 1 has a mother in line 6 or 7
 223         line6or7[1] = imoth2;
 224         flavourLine6or7[1] = flavour2;
 225       }
 226       if((line6or7[0] > 0 && line6or7[1] > 0) && //both tracks must have an entry in line 6 or 7
 227          (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5) && //this entry must be a c or b quark
 228          (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5) && // == " ==
 229          (flavPar1 == flavPar2)){ //make sure that the first hadronised parents of the 2 tracks are of the same flavour
 230         this->SetCorrelated(kTRUE);
 231         fCreationProcess = 0;
 232         return;
 233       }
 234     }
 235   }
 236
 237   // 3.)flavour excitation: if pythia line 6 of one track *and* pythia line 7 of second track
 238   // are filled with a Q and Qbar and if in addition there is another heavy quark in line(s) 4 and/or 5
 239   Int_t line2or3[2] = {-1, -1}; //holds the index of g/q in line 2 or 3
 240   Int_t flavourLine2or3[2] = {-1, -1};
 241   for (Int_t imoth1 = 3; imoth1>=0; imoth1--) {
 242     Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
 243     Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
 244     if(lineMo1 == 2 || lineMo1 == 3){ //track 0 has a mother in line 2 or 3
 245       line2or3[0] = imoth1;
 246       flavourLine2or3[0] = flavour1;
 247     }
 248     for (Int_t imoth2 = 3; imoth2>=0; imoth2--) {
 249       Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
 250       Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
 251       if(lineMo2 == 2 || lineMo2 == 3){ //track 1 has a mother in line 2 or 3
 252         line2or3[1] = imoth2;
 253         flavourLine2or3[1] = flavour2;
 254       }
 255       if(((line6or7[0] > 0 && (flavourLine6or7[0] == 4  || flavourLine6or7[0] == 5)) && //first track has Q in line 6 or 7
 256           (line2or3[1] > 0 && (flavourLine2or3[1] == 21 || flavourLine2or3[1] < 10))) || //second track has a g/q in line 2 or 3
 257          ((line6or7[1] > 0 && (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5)) &&  //or the same,
 258           (line2or3[0] > 0 && (flavourLine2or3[0] == 21 || flavourLine2or3[0] < 10)))){ // swapping the track's indices
 259         //now, check also that the string fragmented into two hadrons
 260         //of the same flavour (string usually splits into many hadrons
 261         //among which there are mostly soft particles)
 262         if(flavPar1 == flavPar2){
 263           this->SetCorrelated(kTRUE);
 264           fCreationProcess = 2;
 265           return;
 266         }
 267       }
 268     }
 269   }
 270
 271   //now flag (rare) processes in which only the incoming parton in line 2 or 3
 272   //radiates a gluon which produces a QQbar pair:
 273   //exclude the light quarks
 274   if(line1 == line2 && (line1 == 2 || line1 == 3)){
 275     if((TMath::Abs(mu1->GetQuarkPDGCode(1)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 4) ||
 276        (TMath::Abs(mu1->GetQuarkPDGCode(1)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 5)){
 277
 278       //now, check also that the string fragmented into two hadrons
 279       //of the same flavour (string usually splits into many hadrons
 280       //among which there are mostly soft particles)
 281       if(flavPar1 == flavPar2){
 282
 283         this->SetCorrelated(kTRUE);
 284         fCreationProcess = 1;
 285         if(GetCauseOfCorrelation() == -1){
 286           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
 287         }
 288         return;
 289       }
 290     }
 291   }
 292
 293   //in initial-state-radiation produced QQbar events the "mother quark"
 294   //is acknowledged as the second quark [1] and sits in line 2 or 3
 295   //is part of gluon splitting
 296   line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark of outgoing quark in [0]
 297   line2 = mu2->GetQuarkPythiaLine(1);
 298   if(line1 == line2 && (line1 == 2 || line1 == 3)){
 299     if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
 300        (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
 301
 302       //now, check also that the string fragmented into two hadrons
 303       //of the same flavour (string usually splits into many hadrons
 304       //among which there are mostly soft particles)
 305       if(flavPar1 == flavPar2){
 306
 307         this->SetCorrelated(kTRUE);
 308         fCreationProcess = 1;
 309         if(GetCauseOfCorrelation() == -1){
 310           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1)); //should be flagged as initial state radiation?
 311         }
 312         return;
 313       }
 314     }
 315   }
 316
 317   //in final-state-radiation produced QQbar events the "mother quark"
 318   //is acknowledged as the first quark [1] and sits in line 6 or 7
 319   //is part of gluon splitting
 320   line1 = mu1->GetQuarkPythiaLine(1); //[1] ... direct mother quark
 321   line2 = mu2->GetQuarkPythiaLine(1);
 322   if(line1 == line2 && (line1 == 6 || line1 == 7)){
 323     if((TMath::Abs(mu1->GetQuarkPDGCode(0)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 4) ||
 324        (TMath::Abs(mu1->GetQuarkPDGCode(0)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(0)) == 5)){
 325
 326       //now, check also that the string fragmented into two hadrons
 327       //of the same flavour (string usually splits into many hadrons
 328       //among which there are mostly soft particles)
 329       if(flavPar1 == flavPar2){
 330
 331         this->SetCorrelated(kTRUE);
 332         fCreationProcess = 1;
 333         if(GetCauseOfCorrelation() == -1){
 334           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
 335         }
 336         return;
 337       }
 338     }
 339   }
 340 }
 341
 342 //====================================
 343 void AliMUONPairLight::SetMuons(const AliMUONTrackLight& mu0, const AliMUONTrackLight& mu1){
 344   /// set the two muons
 345   fMu0 = mu0;
 346   fMu1 = mu1;
 347   this->SetProcess();
 348 }
 349
 350 //====================================
 351 void AliMUONPairLight::PrintInfo(const Option_t* opt){
 352   /// print information about muon pairs
 353   /// Options:
 354   /// - "H" single muons' decay histories
 355   /// - "K" dimuon kinematics
 356   /// - "F" dimuon flags
 357   /// - "A" all variables
 358   TString options(opt);
 359   options.ToUpper();
 360
 361   if(options.Contains("H") || options.Contains("A")){//muon decay histories
 362
 363     AliMUONTrackLight *mu1 = &fMu0;
 364     AliMUONTrackLight *mu2 = &fMu1;
 365
 366     printf("========= History =======================\n");
 367     printf("first muon");
 368     mu1->PrintInfo("H");
 369     printf("second muon");
 370     mu2->PrintInfo("H");
 371     printf("=========================================\n");
 372   }
 373   if(options.Contains("F") || options.Contains("A")){//flags
 374     printf("the flags set for this muon pair are:\n");
 375     printf("=====================================\n");
 376     if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
 377     fIsCorrelated ? printf("(*) it is a correlated pair\n") : printf("(*) it is not a correlated pair\n");
 378     if(IsOpenCharm()) printf("(*) correlated open charm: ");
 379     if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
 380     if(IsOpenCharm() || IsOpenBeauty()){
 381       switch(fCreationProcess){
 382       case 0:
 383         printf("pair creation");
 384         break;
 385       case 1:
 386         printf("gluon splitting");
 387         break;
 388       case 2:
 389         printf("flavour excitation");
 390         break;
 391       case 3:
 392         printf("both muons come from same fragmented mother");
 393         break;
 394       }
 395       if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation())
 396         printf("... where oscillation occured\n");
 397       else{
 398         if(IsOpenBeauty())
 399           printf(" (no oscillation)\n");
 400         else
 401           printf("\n");
 402       }
 403     }
 404     IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(0, fIsFeedDown)) : printf("(*) it is not a resonance\n");
 405     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");
 406     printf("=====================================\n");
 407   }
 408   if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
 409     Double_t *vtx = this->GetMuon(0)->GetVertex();
 410     TLorentzVector momRec = this->GetPRec();
 411     TLorentzVector momGen = this->GetPGen();
 412     printf("the dimuon charge is %d\n", this->GetCharge());
 413     printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
 414     printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
 415     printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
 416     //rapidity, pT, angles, ...
 417     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",
 418            momRec.M(), momRec.Pt(), momRec.Eta(),
 419            TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(),
 420            momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
 421            momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
 422   }
 423 }
 424
 425 //====================================
 426 Double_t AliMUONPairLight::GetOpeningAngle() {
 427   /// opening angle between the two muons in the lab frame (in degrees)
 428   TLorentzVector pRecMu0 =  fMu0.GetPRec();
 429   TLorentzVector pRecMu1 =  fMu1.GetPRec();
 430   TVector3 pRecMu03 = pRecMu0.Vect();
 431   TVector3 pRecMu13 = pRecMu1.Vect();
 432   Double_t scalar = pRecMu03.Dot(pRecMu13);
 433   Double_t modMu0 = pRecMu03.Mag();
 434   Double_t modMu1 = pRecMu13.Mag();
 435   Double_t theta = (TMath::ACos(scalar/(modMu0*modMu1)))*(180./TMath::Pi());
 436   return theta;
 437 }
 438 //================================================
 439 Bool_t AliMUONPairLight::IsDimuonFromCorrPiK(){
 440   ///check if we have a correlated pi/K
 441
 442   AliMUONTrackLight *mu0 = this->GetMuon(0), *mu1 = this->GetMuon(1);
 443   Bool_t fromSameLine = kFALSE;
 444   if (mu0->IsParentPionOrKaon() &&
 445       mu1->IsParentPionOrKaon() &&
 446       mu1->GetQuarkPythiaLine() == mu0->GetQuarkPythiaLine()
 447       ) fromSameLine = kTRUE;
 448
 449   return fromSameLine;
 450 }