- ESD can be used instead of AliMUONTrack objects to access the reconstructed variables.
authorivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 19 Mar 2007 15:26:47 +0000 (15:26 +0000)
committerivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 19 Mar 2007 15:26:47 +0000 (15:26 +0000)
- the zvertex from ITS can be used.
- Improvement in the Pythia decoding
- Cleaning of the code and small bug corrections.
(Alessandro, Hermine)
- Comments for Doxygen (mostly added comments for inline functions)
(Ivana)

MUON/AliMUONPairLight.cxx
MUON/AliMUONPairLight.h
MUON/AliMUONTrackLight.cxx
MUON/AliMUONTrackLight.h

index 1d68580..61ee91e 100644 (file)
@@ -87,13 +87,18 @@ Bool_t AliMUONPairLight::IsAResonance(){
   if (!fIsCorrelated) return kFALSE;   //if muons not correlated, cannot be a resonance
   //if muons are correlated, check if the PDG of the
   //common mother is a resonance
-  Int_t pdg = GetCauseOfCorrelation();
-  if(pdg < 10) return kFALSE;
-  Int_t id=pdg%100000;
-  if(((id-id%10)%110)) return kFALSE;
-  else return kTRUE;
-  printf("<AliMUONPairLight::IsAResonance> arriving after this piece of code\n");
+  Int_t nparents0 = fMu0.GetNParents(); 
+  Int_t nparents1 = fMu1.GetNParents(); 
 
+  Int_t minP = TMath::Min(nparents0, nparents1);
+  for (Int_t i = 0 ; i < minP; i++) { 
+    if (fMu0.IsMotherAResonance(nparents0-1-i) && fMu1.IsMotherAResonance(nparents1-1-i) && 
+       fMu0.GetParentPythiaLine(nparents0-1-i)==fMu1.GetParentPythiaLine(nparents1-1-i)) {
+      if (nparents0-1-i) SetFeedDown(nparents0-1-i);
+      return kTRUE;
+    }
+  }
+  return kFALSE; 
 }
 
 //====================================
@@ -119,87 +124,140 @@ Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) {
 void AliMUONPairLight::SetProcess(){
   /// finds the process related to the muon pair (open charm/beauty, resonance, 
   /// uncorrelated...) 
+
   AliMUONTrackLight *mu1 = &fMu0;
   AliMUONTrackLight *mu2 = &fMu1;
 
-  //1.) check if one of the tracks is not a muon
-  if(IsOneTrackNotAMuon()) { 
-    this->SetCorrelated(kFALSE);
-    return;
-  }
-
   // check if the two muons are correlated
   // first check if they come from the same hadron (resonance or beauty/charm meson)
   Int_t npar1 = mu1->GetNParents(); 
   Int_t npar2 = mu2->GetNParents(); 
-  //  for (Int_t imoth1 = 0; imoth1<npar1; imoth1++) { 
   for (Int_t imoth1 = npar1-1; imoth1>=0; imoth1--) { 
     Int_t lineMo1 = mu1->GetParentPythiaLine(imoth1);
-    //    for (Int_t imoth2 = 0; imoth2<npar2; imoth2++) { 
     for (Int_t imoth2 = npar2-1; imoth2>=0; imoth2--) { 
       Int_t lineMo2 = mu2->GetParentPythiaLine(imoth2);
       if(lineMo1 == lineMo2) { 
+       //reject "diquark" mothers
+       if(mu1->IsDiquark(mu1->GetParentPDGCode(imoth1)))return;
+//     if(IsDiquark(mu1->GetParentPDGCode(imoth1))) return;
        this->SetCorrelated(kTRUE); 
        this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
        if(!IsAResonance()) fCreationProcess = 3; 
        else fCreationProcess = -1;
-       return;   // <<<<---------------- RETURN? 
+       return;
       }
     }
   }  
-  
   // if Open Beauty/Charm we can have 3 creation processes 
   // (pair creation [0], gluon splitting [1] or flavour excitation [2])
-  //
   // 1.) gluon splitting: gluon (stored with index 2, id=21) must be the same 
-  if (mu1->GetQuarkPythiaLine(2) == mu2->GetQuarkPythiaLine(2) && mu1->GetQuarkPDGCode(2) == 21) { 
-    this->SetCorrelated(kTRUE); 
-    if(GetCauseOfCorrelation() == -1){
-      this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(2));
+  Int_t flavPar1 = mu1->GetParentFlavour(0);
+  Int_t flavPar2 = mu2->GetParentFlavour(0);
+  for (Int_t imoth1 = 0; imoth1 < 4; imoth1++) { 
+    Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
+    for (Int_t imoth2 = 0; imoth2 < 4; imoth2++) { 
+      Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
+      if(lineMo1 == lineMo2 && mu1->GetQuarkPDGCode(imoth1) == 21) {
+       //now, check also that the string fragmented into two hadrons
+       //of the same flavour (string usually splits into many hadrons
+       //among which there are mostly soft particles)
+       if(flavPar1 == flavPar2){
+         this->SetCorrelated(kTRUE); 
+         if(GetCauseOfCorrelation() == -1)
+           this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(imoth1));
+
+         fCreationProcess = 1; 
+         return;
+       }
+      }
     }
-    fCreationProcess = 1; 
-    return;
   }
-  
-  // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
-  // are filled with a Q and Qbar
   Int_t line1 = mu1->GetQuarkPythiaLine(2); //[2] ... very first quark
   Int_t line2 = mu2->GetQuarkPythiaLine(2); 
-  Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(2)); //[2] ... very first quark
-  Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(2)); 
-  if ((line1 == 6 || line1 == 7) && (line2 == 6 || line2 == 7)) { 
-    if((flavour1 == 4 || flavour1 == 5) && (flavour2 == 4 || flavour2 == 5)){
-      this->SetCorrelated(kTRUE);
-      fCreationProcess = 0; 
-      return;
+
+  Int_t line6or7[2] = {-1, -1}; //holds the index of quark in line 6 or 7
+  Int_t flavourLine6or7[2] = {-1, -1};
+  // 2.) pair creation: if pythia line 6 of one track *and* pythia line 7 of second track
+  // are filled with a Q and Qbar
+  for (Int_t imoth1 = 3; imoth1>=0; imoth1--) { 
+    Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
+    Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
+    if(lineMo1 == 6 || lineMo1 == 7){ //track 0 has a mother in line 6 or 7
+      line6or7[0] = imoth1;
+      flavourLine6or7[0] = flavour1;
+    }
+    for (Int_t imoth2 = 3; imoth2>=0; imoth2--) { 
+      Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
+      Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
+      if(lineMo2 == 6 || lineMo2 == 7){ //track 1 has a mother in line 6 or 7
+       line6or7[1] = imoth2;
+       flavourLine6or7[1] = flavour2;
+      }
+      if((line6or7[0] > 0 && line6or7[1] > 0) && //both tracks must have an entry in line 6 or 7
+        (flavourLine6or7[0] == 4 || flavourLine6or7[0] == 5) && //this entry must be a c or b quark
+        (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5) && // == " ==
+        (flavPar1 == flavPar2)){ //make sure that the first hadronised parents of the 2 tracks are of the same flavour 
+       this->SetCorrelated(kTRUE);
+       fCreationProcess = 0; 
+       return;
+      }
     }
   }
 
+
   // 3.)flavour excitation:
-  if((((line1 == 6 || line1 == 7) && (flavour1 == 4 || flavour1 == 5)) && //first track has Q in line 6 or 7
-      ((line2 == 2 || line2 == 3) && (flavour2 == 21 || flavour2 < 10))) || //second track has a g/q in line 2 or 3
-      (((line2 == 6 || line2 == 7) && (flavour2 == 4 || flavour2 == 5)) && //or the same,
-      ((line1 == 2 || line1 == 3) && (flavour1 == 21 || flavour1 < 10)))){ // swapping the track's indices
-
-//     printf("candidate for flavour excitation\n");
-
-    //Hermine: is it needed to check the lines 4-7???
-
-    //now, we have a candidate for flavour excitation
-    //we must also verify if in the Pythia listing
-    //the "incoming" lines 4 and 5 and "outgoing" lines 6 and 7
-    //are filled with g Q each: e.g. 4(g),5(Q),5(g),6(Q)
-//     Int_t pdg4 = TMath::Abs(stack()->Particle(4)->GetPdgCode());
-//     Int_t pdg5 = TMath::Abs(stack()->Particle(5)->GetPdgCode());
-//     Int_t pdg6 = TMath::Abs(stack()->Particle(6)->GetPdgCode());
-//     Int_t pdg7 = TMath::Abs(stack()->Particle(7)->GetPdgCode());
-//     if((pdg4 == 21 && pdg5 < 10 || pdg5 == 21 && pdg4 < 10 ) &&
-//        (pdg6 == 21 && pdg7 < 10 || pdg7 == 21 && pdg6 < 10 )){
-//       this->PrintInfo("H");
-      this->SetCorrelated(kTRUE);
-      fCreationProcess = 2;
-      return;
-//     }
+  Int_t line2or3[2] = {-1, -1}; //holds the index of g/q in line 2 or 3
+  Int_t flavourLine2or3[2] = {-1, -1};
+  for (Int_t imoth1 = 3; imoth1>=0; imoth1--) { 
+    Int_t lineMo1 = mu1->GetQuarkPythiaLine(imoth1);
+    Int_t flavour1 = TMath::Abs(mu1->GetQuarkPDGCode(imoth1));
+    if(lineMo1 == 2 || lineMo1 == 3){ //track 0 has a mother in line 2 or 3
+      line2or3[0] = imoth1;
+      flavourLine2or3[0] = flavour1;
+    }
+    for (Int_t imoth2 = 3; imoth2>=0; imoth2--) { 
+      Int_t lineMo2 = mu2->GetQuarkPythiaLine(imoth2);
+      Int_t flavour2 = TMath::Abs(mu2->GetQuarkPDGCode(imoth2));
+      if(lineMo2 == 2 || lineMo2 == 3){ //track 1 has a mother in line 2 or 3
+       line2or3[1] = imoth2;
+       flavourLine2or3[1] = flavour2;
+      }
+      if(((line6or7[0] > 0 && (flavourLine6or7[0] == 4  || flavourLine6or7[0] == 5)) && //first track has Q in line 6 or 7
+         (line2or3[1] > 0 && (flavourLine2or3[1] == 21 || flavourLine2or3[1] < 10))) || //second track has a g/q in line 2 or 3
+        ((line6or7[1] > 0 && (flavourLine6or7[1] == 4 || flavourLine6or7[1] == 5)) &&  //or the same,
+         (line2or3[0] > 0 && (flavourLine2or3[0] == 21 || flavourLine2or3[0] < 10)))){ // swapping the track's indices
+       //now, check also that the string fragmented into two hadrons
+       //of the same flavour (string usually splits into many hadrons
+       //among which there are mostly soft particles)
+       if(flavPar1 == flavPar2){
+         this->SetCorrelated(kTRUE);
+         fCreationProcess = 2;
+         return;
+       }
+      }
+    }
+  } 
+
+  //now flag (rare) processes in which only the incoming parton in line 2 or 3
+  //radiates a gluon which produces a QQbar pair:
+  //exclude the light quarks
+  if(line1 == line2 && (line1 == 2 || line1 == 3)){
+    if((TMath::Abs(mu1->GetQuarkPDGCode(1)) == 4 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 4) ||
+       (TMath::Abs(mu1->GetQuarkPDGCode(1)) == 5 && TMath::Abs(mu2->GetQuarkPDGCode(1)) == 5)){
+
+      //now, check also that the string fragmented into two hadrons
+      //of the same flavour (string usually splits into many hadrons
+      //among which there are mostly soft particles)
+      if(flavPar1 == flavPar2){
+
+       this->SetCorrelated(kTRUE);
+       fCreationProcess = 1;
+       if(GetCauseOfCorrelation() == -1){
+         this->SetCauseOfCorrelation(mu1->GetQuarkPDGCode(1));
+       }
+       return;
+      }
+    }
   }
 }
 
@@ -265,8 +323,8 @@ void AliMUONPairLight::PrintInfo(Option_t* opt){
          printf("\n");
       }
     }
-    IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(fIsFeedDown)) : printf("(*) it is not a resonance\n");
-    fIsFeedDown ? printf("(*) mother has feed-down: %d --> %d\n", this->GetMuonMotherPDG(1), this->GetMuonMotherPDG(0)) : printf("(*) no feed-down\n");
+    IsAResonance() ? printf("(*) it is a resonance: %d\n", this->GetMuonMotherPDG(0, fIsFeedDown)) : printf("(*) it is not a resonance\n");
+    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");
     printf("=====================================\n");
   }
   if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
@@ -286,8 +344,9 @@ void AliMUONPairLight::PrintInfo(Option_t* opt){
   }
 }
 
+//====================================
 Double_t AliMUONPairLight::GetOpeningAngle() { 
-  // opening angle between the two muons in the lab frame (in degrees)
+  /// opening angle between the two muons in the lab frame (in degrees)
   TLorentzVector pRecMu0 =  fMu0.GetPRec();
   TLorentzVector pRecMu1 =  fMu1.GetPRec();
   TVector3 pRecMu03 = pRecMu0.Vect();
index e8fec47..178375b 100644 (file)
@@ -37,36 +37,55 @@ public:
   AliMUONPairLight(AliMUONPairLight &dimuCopy);
   virtual ~AliMUONPairLight();
   virtual void SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1);
-  AliMUONTrackLight* GetMuon(Int_t index) ; 
+  AliMUONTrackLight* GetMuon(Int_t index) ;   
   Int_t GetMuonMotherPDG(Int_t imuon, Int_t mother=0) ; 
+  
+  /// \todo add comment
   Bool_t IsOpenCharm() {return (TMath::Abs(fMu0.GetParentFlavour(0))==4 && TMath::Abs(fMu1.GetParentFlavour(0))==4 && fIsCorrelated && !IsAResonance());}
+  /// \todo add comment
   Bool_t IsOpenBeauty() {return (TMath::Abs(fMu0.GetParentFlavour(0))==5 && TMath::Abs(fMu1.GetParentFlavour(0))==5 && fIsCorrelated  && !IsAResonance());}
   Bool_t IsAResonance(); 
+  /// \todo add comment
+  Bool_t IsOneMuonFromPionOrKaon(){return (fMu0.IsParentPionOrKaon(0) || fMu1.IsParentPionOrKaon(0));}
+  /// Return the info if the two muons are of correlated origin
   Bool_t IsCorrelated() const {return fIsCorrelated;}   
+  /// Return the pdg of common mother
   Int_t GetCauseOfCorrelation() const {return fCauseOfCorrelation;}
+  /// Return the info if the process is from feeddown 
   Bool_t IsFeedDown() const {return fIsFeedDown;}   
+  /// \todo add comment
   Bool_t IsOneTrackNotAMuon() {return (!( fMu0.IsAMuon() && fMu1.IsAMuon() )) ;}
+  /// \todo add comment
   Int_t GetCharge() {return fMu0.GetCharge() + fMu1.GetCharge();}
+  /// \brief Return the info ablout creation process
+  ///0: pair creation, 1: gluon splitting, 2: flavour excitation, 3: same fragmented mother, -1: resonance
   Int_t GetCreationProcess() const  {return fCreationProcess;}
+  /// Set the info ablout creation process
   void SetCorrelated(Bool_t answer) {fIsCorrelated = answer; }
+  /// Set the pdg of common mother
   void SetCauseOfCorrelation(Int_t pdg) {fCauseOfCorrelation = pdg; }
+  /// Set the info if the process is from feeddown
   void SetFeedDown(Int_t answer) {fIsFeedDown = answer;}
+  /// \todo add comment
   TLorentzVector GetPRec(){return fMu0.GetPRec()+fMu1.GetPRec();}
+  /// \todo add comment
   TLorentzVector GetPGen(){return fMu0.GetPGen()+fMu1.GetPGen();}
   Double_t GetOpeningAngle(); 
   virtual void PrintInfo(Option_t* opt);//"H" single muons' decay histories
                                            //"K" dimuon kinematics
                                            //"F" dimuon flags
                                            //"A" all variables
+
 protected: 
-  AliMUONTrackLight fMu0;   /// first muon  
-  AliMUONTrackLight fMu1;   /// second muon
-  Int_t fCreationProcess; ///0: pair creation, 1: gluon splitting, 2: flavour excitation, 3: same fragmented mother, -1: resonance
-  Bool_t fIsCorrelated; ///tells if the two muons are of correlated origin
-  Int_t fCauseOfCorrelation; ///pdg of common mother
-  Int_t fIsFeedDown;     /// tells if the process is from feeddown 
+  /// Checks if muons are correlated and assigns
+  void SetProcess();
 
-  void SetProcess();///checks if muons are correlated and assigns
+  AliMUONTrackLight fMu0;   ///< first muon  
+  AliMUONTrackLight fMu1;   ///< second muon
+  Int_t fCreationProcess; ///<0: pair creation, 1: gluon splitting, 2: flavour excitation, 3: same fragmented mother, -1: resonance
+  Bool_t fIsCorrelated; ///<tells if the two muons are of correlated origin
+  Int_t fCauseOfCorrelation; ///<pdg of common mother
+  Int_t fIsFeedDown;     ///< tells if the process is from feeddown 
 
   ClassDef(AliMUONPairLight,1)  /// Dimuon carrying info at generation and reconstruction level
 }; 
index 29102fe..16d551a 100644 (file)
@@ -1,6 +1,6 @@
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *      SigmaEffect_thetadegrees                                                                  *
+ *      SigmaEffect_thetadegrees                                          *
  * Author: The ALICE Off-line Project.                                    *
  * Contributors are mentioned in the code where appropriate.              *
  *                                                                        *
@@ -38,6 +38,7 @@
 #include "AliRunLoader.h"
 #include "AliStack.h"
 #include "AliHeader.h"
+#include "AliMUONTrackExtrap.h"
 
 #include "TDatabasePDG.h"
 #include "TClonesArray.h"
@@ -53,12 +54,14 @@ AliMUONTrackLight::AliMUONTrackLight()
     fPrec(), 
     fIsTriggered(kFALSE),
     fCharge(-999), 
+    fChi2(-1), 
     fCentr(-1),
     fPgen(), 
     fTrackPythiaLine(-999),
     fTrackPDGCode(-999),
     fOscillation(kFALSE), 
-    fNParents(0)
+    fNParents(0),
+    fWeight(1)    
 {
   /// default constructor
   fPgen.SetPxPyPzE(0.,0.,0.,0.); 
@@ -80,12 +83,14 @@ AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
     fPrec(muonCopy.fPrec), 
     fIsTriggered(muonCopy.fIsTriggered),
     fCharge(muonCopy.fCharge), 
+    fChi2(muonCopy.fChi2), 
     fCentr(muonCopy.fCentr),
     fPgen(muonCopy.fPgen), 
     fTrackPythiaLine(muonCopy.fTrackPythiaLine),
     fTrackPDGCode(muonCopy.fTrackPDGCode),
     fOscillation(muonCopy.fOscillation), 
-    fNParents(muonCopy.fNParents)
+    fNParents(muonCopy.fNParents),
+    fWeight(muonCopy.fWeight)
 {
   /// copy constructor
   for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i]; 
@@ -105,20 +110,45 @@ AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
     fPrec(), 
     fIsTriggered(kFALSE),
     fCharge(-999), 
+    fChi2(-1),
     fCentr(-1),
     fPgen(), 
     fTrackPythiaLine(-999),
     fTrackPDGCode(-999),
     fOscillation(kFALSE), 
-    fNParents(0)
+    fNParents(0),
+    fWeight(1)
 { 
   /// constructor
   //AliMUONTrackLight(); 
-  ComputePRecAndChargeFromESD(muonTrack); 
+  FillFromESD(muonTrack); 
 }
 
 //============================================
-void AliMUONTrackLight::ComputePRecAndChargeFromESD(AliESDMuonTrack* muonTrack){ 
+AliMUONTrackLight::~AliMUONTrackLight()
+{
+/// Destructor
+} 
+
+//============================================
+
+void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
+  /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
+  AliMUONTrackParam trPar(*((AliMUONTrackParam*) (trackReco->GetTrackParamAtHit()->First())));
+  //  AliMUONTrackParam *trPar  = trackReco->GetTrackParamAtVertex();
+  AliMUONTrackExtrap::ExtrapToVertex(&trPar,0.,0.,0.);
+  this->SetCharge(Int_t(TMath::Sign(1.,trPar.GetInverseBendingMomentum())));
+  this->SetPxPyPz(trPar.Px(),trPar.Py(), trPar.Pz()); 
+  this->SetTriggered(trackReco->GetMatchTrigger()); 
+  
+  Double_t xyz[3] = { trPar.GetNonBendingCoor(), 
+                     trPar.GetBendingCoor(),
+                     zvert};
+  this->SetVertex(xyz); 
+}
+
+//============================================
+void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
   /// computes prec and charge from ESD track
   Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass(); 
   Double_t thetaX = muonTrack->GetThetaX();
@@ -132,6 +162,14 @@ void AliMUONTrackLight::ComputePRecAndChargeFromESD(AliESDMuonTrack* muonTrack){
   fCharge   = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
   Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
   fPrec.SetPxPyPzE(px,py,pz,energy);
+  // get the position
+  fXYZ[0] = muonTrack->GetNonBendingCoor();
+  fXYZ[1] = muonTrack->GetBendingCoor();
+  if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
+  else fXYZ[2] = zvert;
+  // get the chi2 per d.o.f.
+  fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
+  fIsTriggered = muonTrack->GetMatchTrigger();
 }
 
 //============================================
@@ -194,18 +232,27 @@ Int_t AliMUONTrackLight::TrackCheck(Bool_t *compTrack){
 //============================================
 void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, TParticle *part){
   /// scans the muon history to determine parents pdg code and pythia line
+  // kept for backward compatibility
+  // see the overloaded method FillMuonHistory(AliStack *stack, TParticle *part)
+  AliStack *stack = runLoader->GetHeader()->Stack();
+  FillMuonHistory(stack,part);
+}
+
+//============================================
+void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
+  /// scans the muon history to determine parents pdg code and pythia line
   Int_t countP = -1;
   Int_t parents[10], parLine[10];
   Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
 
-  AliStack *stack = runLoader->GetHeader()->Stack();
   TParticle *mother;
   Int_t status=-1, pdg=-1;
   while(lineM >= 0){
-
+    
     mother = stack->Particle(lineM); //direct mother of rec. track
     pdg = mother->GetPdgCode();//store PDG code of first mother
-    if(pdg == 92) break;  // break if a string is found 
+    // break if a string, gluon, quark or diquark is found 
+    if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
     parents[++countP] = pdg;
     parLine[countP] = lineM;
     status = mother->GetStatusCode();//get its status code to check if oscillation occured
@@ -217,9 +264,9 @@ void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, TParticle *part
     this->SetParentPDGCode(i,parents[countP-i]);
     this->SetParentPythiaLine(i,parLine[countP-i]);
   }
-
   fNParents = countP+1;
   countP = -1;
+
   //and store the lines of the string and further quarks in another array:
   while(lineM >= 0){
     mother = stack->Particle(lineM);
@@ -229,14 +276,21 @@ void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, TParticle *part
     this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
     lineM = mother->GetFirstMother();
   }
-
+  
   //check if in case of HF production, the string points to the correct end
   //and correct it in case of need:
   countP = 1;
+  for(int par = 0; par < 4; par++){
+    if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
+      countP = par; //get the quark just before hadronisation
+      break;
+    }
+  }
   if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
     if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
 
-      Int_t pdg = this->GetQuarkPDGCode(1), line = this->GetQuarkPythiaLine(1);
+      printf("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n", this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)));
+      Int_t pdg = this->GetQuarkPDGCode(countP), line = this->GetQuarkPythiaLine(countP);
       this->ResetQuarkInfo();
       while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
                                                         //must coincide with the flavour of the last fragmented mother
@@ -252,9 +306,11 @@ void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, TParticle *part
        this->SetQuarkPDGCode(countP++, pdg);
        line = mother->GetFirstMother();
       }
-    }
+      this->PrintInfo("h");
+    }//mismatch
   }
 }
+
 //====================================
 void AliMUONTrackLight::ResetQuarkInfo(){
   /// resets parton information
@@ -346,5 +402,25 @@ void AliMUONTrackLight::PrintInfo(Option_t* opt){
           momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
   }
 }
-
-
+//====================================
+Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
+  /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
+  Int_t pdg = this->GetParentPDGCode(idparent); 
+  if (TMath::Abs(pdg)==211 || //pi+
+      TMath::Abs(pdg)==321 || //K+
+      TMath::Abs(pdg)==213 || //rho+
+      TMath::Abs(pdg)==311 || //K0
+      TMath::Abs(pdg)==313 || //K*0
+      TMath::Abs(pdg)==323    //K*+
+      ) { 
+    return kTRUE;
+  }
+  else return kFALSE;
+}
+//====================================
+Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg){
+  /// check if the provided pdg code corresponds to a diquark 
+  pdg = TMath::Abs(pdg);
+  if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
+  else return kFALSE;
+}
index ad89a01..d7afc3e 100644 (file)
 /// (authors: H.Woehri, A.de Falco)
 
 // ROOT classes
-
-#include <TClonesArray.h>
-#include <TLorentzVector.h>
+#include "TLorentzVector.h"
 
 class AliMUONTrack;
 class AliESDMuonTrack;
 class AliRunLoader;
-
+class AliStack; 
 class TParticle;
 
 class AliMUONTrackLight : public TObject { 
@@ -38,60 +36,107 @@ class AliMUONTrackLight : public TObject {
   AliMUONTrackLight(); 
   AliMUONTrackLight(AliESDMuonTrack* muonTrack); 
   AliMUONTrackLight(const AliMUONTrackLight &muonCopy);
-  virtual ~AliMUONTrackLight(){;} 
+  virtual ~AliMUONTrackLight(); 
   
+  /// Set 4-momentum of the generated particle
   void SetPGen(TLorentzVector pgen) {fPgen = pgen;}
+  /// Return 4-momentum of the generated particle  
   TLorentzVector GetPGen() const {return fPgen;}
+  /// Set reconstructed 4-momentum
   void SetPRec(TLorentzVector prec) {fPrec = prec;}
+  /// Return reconstructed 4-momentum 
   TLorentzVector GetPRec() const {return fPrec;}
+  /// Set primary vertex position from the ITS
   void SetVertex(Double_t *xyz) {for (Int_t i=0; i<3; i++) fXYZ[i]=xyz[i];}
+  /// Return primary vertex x position from the ITS 
+  Double_t GetX() { return fXYZ[0]; } 
+  /// Return primary vertex y position from the ITS 
+  Double_t GetY() { return fXYZ[1]; } 
+  /// Return primary vertex z position from the ITS 
+  Double_t GetZ() { return fXYZ[2]; } 
+  /// Return  primary vertex position from the ITS
   Double_t* GetVertex() { return fXYZ; } 
+  /// Set chi2 / ndf in the MUON track fit
+  void SetChi2(Double_t chi2) {fChi2=chi2;}
+  /// Return chi2 / ndf in the MUON track fit 
+  Double_t GetChi2() { return fChi2; } 
+  /// Set weight assigned to the muon
+  void SetWeight(Double_t w) {fWeight=w;}
+  /// Return weight assigned to the muon 
+  Double_t GetWeight() { return fWeight; } 
+  
+  /// Set muon charge 
   void SetCharge(Int_t charge) {fCharge = charge;}
+  /// Return muon charge 
   Int_t GetCharge() const {return fCharge;}
+  /// Return hadronised parents and grandparents 
   Int_t GetParentPDGCode(Int_t index = 0) const { return fParentPDGCode[index]; } 
+  /// Return line of Pythia output for hadronised parents & grandparents 
   Int_t GetParentPythiaLine(Int_t index = 0) const { return fParentPythiaLine[index]; } 
+  /// Return pdg of the string [0], quarks/gluons [1,2], sometimes proton [3] 
   Int_t GetQuarkPDGCode(Int_t index = 0) const { return fQuarkPDGCode[index]; } 
+  /// Return line of Pythia output for string [0] and quarks [1,2], sometimes proton [3] 
   Int_t GetQuarkPythiaLine(Int_t index = 0) const { return fQuarkPythiaLine[index]; } 
+  /// Return line of Pythia output for string [0] and quarks [1,2], sometimes proton [3] 
   Int_t GetTrackPythiaLine() const {return fTrackPythiaLine;}
+  /// Return pdg code of the rec. track (in general will be a muon) 
   Int_t GetTrackPDGCode() const {return fTrackPDGCode;}
+  /// Set line of kin. stack where rec. track (in general, the muon) is stored
   void SetTrackPythiaLine(Int_t trackLine) {fTrackPythiaLine = trackLine;}
+  /// Set pdg code of the rec. track (in general will be a muon)
   void SetTrackPDGCode(Int_t trackPdg) {fTrackPDGCode = trackPdg;}
-  void ComputePRecAndChargeFromESD(AliESDMuonTrack* muonTrack);
+  void FillFromESD(AliESDMuonTrack* muonTrack, Double_t zvert=-9999);
+  void FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert=-9999);  
+  /// Return info if is a muon 
   Bool_t IsAMuon() const { return (TMath::Abs(fTrackPDGCode)==13); }
+  Bool_t IsParentPionOrKaon(Int_t idParent = 0);
   void SetPxPyPz(Double_t px, Double_t py, Double_t pz); 
+  /// Set flag for trigger 
   void SetTriggered(Bool_t isTriggered) { fIsTriggered = isTriggered; } 
+  /// Return flag for trigger  
   Bool_t IsTriggered() const { return fIsTriggered; } 
   TParticle* FindRefTrack(AliMUONTrack* trackReco, TClonesArray* trackRefArray, AliRunLoader *runLoader); 
   Int_t TrackCheck(Bool_t *compTrack);
+  /// Return acually filled no. of *fragmented* parents 
   Int_t GetNParents() const {return fNParents;}
+  void FillMuonHistory(AliStack *stack, TParticle *part);
   void FillMuonHistory(AliRunLoader *runLoader, TParticle *part);
   Bool_t IsB0(Int_t intTest) const;//checks if the provided PDG code corresponds to a neutral B meson
   Bool_t IsMotherAResonance(Int_t index=0) const;
+  /// Return flag for oscillation 
   Bool_t GetOscillation() const {return fOscillation;}
   virtual void PrintInfo(Option_t* opt); //"H" muon's decay history
   //"K" muon kinematics
   //"A" all variables
   Int_t GetParentFlavour(Int_t idParent=0) const;
+  Bool_t IsDiquark(Int_t pdg);
 protected:
-  static const Int_t fgkNParentsMax = 5;   /// maximum number of parents
-  TLorentzVector fPrec; /// reconstructed 4-momentum
-  Double_t fXYZ[3];     /// primary vertex position from the ITS 
-  Bool_t fIsTriggered;  /// flag for trigger 
-  Int_t fCharge;        /// muon charge 
-  Float_t fCentr;       /// centrality 
-  TLorentzVector fPgen;   /// 4-momentum of the generated particle
-  Int_t fTrackPythiaLine; ///line of kin. stack where rec. track (in general, the muon) is stored
-  Int_t fTrackPDGCode; ///pdg code of the rec. track (in general will be a muon)
-  Int_t fParentPDGCode[fgkNParentsMax]; /// hadronised parents and grandparents 
-  Int_t fParentPythiaLine[fgkNParentsMax];///line of Pythia output for hadronised parents & grandparents
-  Int_t fQuarkPDGCode[4]; /// pdg of the string [0], quarks/gluons [1,2], sometimes proton [3] 
-  Int_t fQuarkPythiaLine[4]; ///line of Pythia output for string [0] and quarks [1,2], sometimes proton [3]
-  Bool_t fOscillation; /// flag for oscillation
-  Int_t fNParents; ///acually filled no. of *fragmented* parents
+  static const Int_t fgkNParentsMax = 5;   ///< maximum number of parents
+  TLorentzVector fPrec; ///< reconstructed 4-momentum
+  Double_t fXYZ[3];     ///< primary vertex position from the ITS 
+  Bool_t fIsTriggered;  ///< flag for trigger 
+  Int_t fCharge;        ///< muon charge 
+  Double_t fChi2;       ///< chi2 / ndf in the MUON track fit 
+  Float_t fCentr;       ///< centrality 
+  TLorentzVector fPgen;   ///< 4-momentum of the generated particle
+  Int_t fTrackPythiaLine; ///< line of kin. stack where rec. track (in general, the muon) is stored
+  Int_t fTrackPDGCode; ///< pdg code of the rec. track (in general will be a muon)
+  Int_t fParentPDGCode[fgkNParentsMax]; ///< hadronised parents and grandparents 
+  Int_t fParentPythiaLine[fgkNParentsMax];///< line of Pythia output for hadronised parents & grandparents
+  Int_t fQuarkPDGCode[4]; ///< pdg of the string [0], quarks/gluons [1,2], sometimes proton [3] 
+  Int_t fQuarkPythiaLine[4]; ///< line of Pythia output for string [0] and quarks [1,2], sometimes proton [3]
+  Bool_t fOscillation; ///< flag for oscillation
+  Int_t fNParents; ///< acually filled no. of *fragmented* parents
+  Double_t fWeight; ///< weight assigned to the muon
+  /// Set flag for oscillation
   void SetOscillation(Bool_t oscillation) { fOscillation = oscillation; }
+  /// Set hadronised parents and grandparents 
   void SetParentPDGCode(Int_t index, Int_t pdg) { fParentPDGCode[index] = pdg; } 
+  /// Set line of Pythia output for hadronised parents & grandparents
   void SetParentPythiaLine(Int_t index, Int_t line) { fParentPythiaLine[index] = line; } 
+  /// Set pdg of the string [0], quarks/gluons [1,2], sometimes proton [3] 
   void SetQuarkPDGCode(Int_t index, Int_t pdg){ fQuarkPDGCode[index] = pdg; }
+  /// Set line of Pythia output for string [0] and quarks [1,2], sometimes proton [3]
   void SetQuarkPythiaLine(Int_t index, Int_t line){ fQuarkPythiaLine[index] = line; }
   void ResetQuarkInfo();