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