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.fMu0),
  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 pdg = GetCauseOfCorrelation();
  91   if(pdg < 10) return kFALSE;
  92   Int_t id=pdg%100000;
  93   if(((id-id%10)%110)) return kFALSE;
  94   else return kTRUE;
  95   printf("<AliMUONPairLight::IsAResonance> arriving after this piece of code\n");
  96
  97 }
  98
  99 //====================================
 100
 101 AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index)  {
 102   /// return muon 0 or 1
 103    if (index==0) return &fMu0;
 104    else if (index==1) return &fMu1;
 105    else{ printf ("Index can be either 0 or 1\n"); return 0;}
 106    //   else return &fMu1;
 107 }
 108
 109 //====================================
 110
 111 Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
 112   /// return muon mother pdg code
 113   if (imuon==0) return fMu0.GetParentPDGCode(mother);
 114   else if (imuon==1) return fMu1.GetParentPDGCode(mother);
 115   else { printf ("Index must be only 0 or 1\n"); return -999; }
 116 }
 117
 118 //====================================
 119 void AliMUONPairLight::SetProcess(){
 120   /// finds the process related to the muon pair (open charm/beauty, resonance,
 121   /// uncorrelated...)
 122   AliMUONTrackLight *mu1 = &fMu0;
 123   AliMUONTrackLight *mu2 = &fMu1;
 124
 125   //1.) check if one of the tracks is not a muon
 126   if(IsOneTrackNotAMuon()) {
 127     this->SetCorrelated(kFALSE);
 128     return;
 129   }
 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 = 0; imoth1<npar1; imoth1++) {
 136   for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) {
 137     Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
 138     //    for (Int_t imoth2 = 0; imoth2<npar2; imoth2++) {
 139     for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) {
 140       Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
 141       if(lineMo1 == lineMo2) {
 142         this->SetCorrelated(kTRUE);
 143         this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
 144         if(!IsAResonance()) fCreationProcess = 3;
 145         else fCreationProcess = -1;
 146         return;   // <<<<---------------- RETURN?
 147       }
 148     }
 149   }
 150
 151   // if Open Beauty/Charm we can have 3 creation processes
 152   // (pair creation [0], gluon splitting [1] or flavour excitation [2])
 153   //
 154   // 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same
 155   if (mu1->GetQuarkPythiaLine(2) == mu2->GetQuarkPythiaLine(2) && mu1->GetQuarkPDGCode(2) == 21) {
 156     this->SetCorrelated(kTRUE);
 157     if(GetCauseOfCorrelation() == -1){
 158       this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(2));
 159     }
 160     fCreationProcess = 1;
 161     return;
 162   }
 163
 164   // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
 165   // are filled with a Q and Qbar
 166   Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
 167   Int_t line2 = mu2->GetQuarkPythiaLine(2);
 168   Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(2)); //[2] ... very first quark
 169   Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(2));
 170   if ((line1 == 6 || line1 == 7) && (line2 == 6 || line2 == 7)) {
 171     if((flavour1 == 4 || flavour1 == 5) && (flavour2 == 4 || flavour2 == 5)){
 172       this->SetCorrelated(kTRUE);
 173       fCreationProcess = 0;
 174       return;
 175     }
 176   }
 177
 178   // 3.)flavour excitation:
 179   if((((line1 == 6 || line1 == 7) && (flavour1 == 4 || flavour1 == 5)) && //first track has Q in line 6 or 7
 180       ((line2 == 2 || line2 == 3) && (flavour2 == 21 || flavour2 < 10))) || //second track has a g/q in line 2 or 3
 181       (((line2 == 6 || line2 == 7) && (flavour2 == 4 || flavour2 == 5)) && //or the same,
 182       ((line1 == 2 || line1 == 3) && (flavour1 == 21 || flavour1 < 10)))){ // swapping the track's indices
 183
 184 //     printf("candidate for flavour excitation\n");
 185
 186     //Hermine: is it needed to check the lines 4-7???
 187
 188     //now, we have a candidate for flavour excitation
 189     //we must also verify if in the Pythia listing
 190     //the "incoming" lines 4 and 5 and "outgoing" lines 6 and 7
 191     //are filled with g Q each: e.g. 4(g),5(Q),5(g),6(Q)
 192 //     Int_t pdg4 = TMath::Abs(stack()->Particle(4)->GetPdgCode());
 193 //     Int_t pdg5 = TMath::Abs(stack()->Particle(5)->GetPdgCode());
 194 //     Int_t pdg6 = TMath::Abs(stack()->Particle(6)->GetPdgCode());
 195 //     Int_t pdg7 = TMath::Abs(stack()->Particle(7)->GetPdgCode());
 196 //     if((pdg4 == 21 && pdg5 < 10 || pdg5 == 21 && pdg4 < 10 ) &&
 197 //        (pdg6 == 21 && pdg7 < 10 || pdg7 == 21 && pdg6 < 10 )){
 198 //       this->PrintInfo("H");
 199       this->SetCorrelated(kTRUE);
 200       fCreationProcess = 2;
 201       return;
 202 //     }
 203   }
 204 }
 205
 206 //====================================
 207 void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
 208   /// set the two muons
 209   fMu0 = mu0;
 210   fMu1 = mu1;
 211   this->SetProcess();
 212 }
 213
 214 //====================================
 215 void AliMUONPairLight::PrintInfo(Option_t* opt){
 216   /// print information about muon pairs
 217   /// Options:
 218   /// - "H" single muons' decay histories
 219   /// - "K" dimuon kinematics
 220   /// - "F" dimuon flags
 221   /// - "A" all variables
 222   TString options(opt);
 223   options.ToUpper();
 224
 225   if(options.Contains("H") || options.Contains("A")){//muon decay histories
 226
 227     AliMUONTrackLight *mu1 = &fMu0;
 228     AliMUONTrackLight *mu2 = &fMu1;
 229
 230     printf("========= History =======================\n");
 231     printf("first muon");
 232     mu1->PrintInfo("H");
 233     printf("second muon");
 234     mu2->PrintInfo("H");
 235     printf("=========================================\n");
 236   }
 237   if(options.Contains("F") || options.Contains("A")){//flags
 238     printf("the flags set for this muon pair are:\n");
 239     printf("=====================================\n");
 240     if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
 241     fIsCorrelated ? printf("(*) it is a correlated pair\n") : printf("(*) it is not a correlated pair\n");
 242     if(IsOpenCharm()) printf("(*) correlated open charm: ");
 243     if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
 244     if(IsOpenCharm() || IsOpenBeauty()){
 245       switch(fCreationProcess){
 246       case 0:
 247         printf("pair creation");
 248         break;
 249       case 1:
 250         printf("gluon splitting");
 251         break;
 252       case 2:
 253         printf("flavour excitation");
 254         break;
 255       case 3:
 256         printf("both muons come from same fragmented mother");
 257         break;
 258       }
 259       if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation())
 260         printf("... where oscillation occured\n");
 261       else{
 262         if(IsOpenBeauty())
 263           printf(" (no oscillation)\n");
 264         else
 265           printf("\n");
 266       }
 267     }
 268     IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(fIsFeedDown)) : printf("(*) it is not a resonance\n");
 269     fIsFeedDown ? printf("(*) mother has feed-down: %d --> %d\n", this->GetMuonMotherPDG(1), this->GetMuonMotherPDG(0)) : printf("(*) no feed-down\n");
 270     printf("=====================================\n");
 271   }
 272   if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
 273     Double_t *vtx = this->GetMuon(0)->GetVertex();
 274     TLorentzVector momRec = this->GetPRec();
 275     TLorentzVector momGen = this->GetPGen();
 276     printf("the dimuon charge is %d\n", this->GetCharge());
 277     printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
 278     printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
 279     printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
 280     //rapidity, pT, angles, ...
 281     printf("Rec. variables: 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",
 282            momRec.Pt(), momRec.Eta(),
 283            TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(),
 284            momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
 285            momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
 286   }
 287 }
 288
 289 Double_t AliMUONPairLight::GetOpeningAngle() {
 290   // opening angle between the two muons in the lab frame (in degrees)
 291   TLorentzVector pRecMu0 =  fMu0.GetPRec();
 292   TLorentzVector pRecMu1 =  fMu1.GetPRec();
 293   TVector3 pRecMu03 = pRecMu0.Vect();
 294   TVector3 pRecMu13 = pRecMu1.Vect();
 295   Double_t scalar = pRecMu03.Dot(pRecMu13);
 296   Double_t modMu0 = pRecMu03.Mag();
 297   Double_t modMu1 = pRecMu13.Mag();
 298   Double_t theta = (TMath::ACos(scalar/(modMu0*modMu1)))*(180./TMath::Pi());
 299   return theta;
 300 }