New classes for analysis
authorivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 8 Sep 2006 15:33:56 +0000 (15:33 +0000)
committerivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 8 Sep 2006 15:33:56 +0000 (15:33 +0000)
(Hermine, Alessandro)

MUON/AliMUONPairLight.cxx [new file with mode: 0644]
MUON/AliMUONPairLight.h [new file with mode: 0644]
MUON/AliMUONTrackLight.cxx [new file with mode: 0644]
MUON/AliMUONTrackLight.h [new file with mode: 0644]

diff --git a/MUON/AliMUONPairLight.cxx b/MUON/AliMUONPairLight.cxx
new file mode 100644 (file)
index 0000000..d5bc186
--- /dev/null
@@ -0,0 +1,300 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *      SigmaEffect_thetadegrees                                                                  *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpeateose. It is      *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//===================================================================
+//This class was prepared by INFN Cagliari, July 2006
+//(authors: H.Woehri, A.de Falco)
+// 
+// Compact information for the generated muon pairs in the MUON arm 
+// useful at the last stage of the analysis chain
+// Pairs are built with two AliMUONTrackLight objects 
+// Using the class AliMUONTrackLight this class combines the decay
+// information ("history") of the reconstructed tracks and fills
+// a series of flags for the formed reconstructed dimuon:
+// fIsCorrelated, fCreationProcess, fIsFeedDown, ...
+// for information about the dimuon, use PrintInfo with the appropriate
+// printflag
+// To be used together with AliMUONTrackLight
+//===================================================================
+
+
+//MUON classes
+#include "AliMUONPairLight.h"
+//Root classes
+#include "TString.h"
+
+ClassImp(AliMUONPairLight) 
+
+//====================================
+AliMUONPairLight::AliMUONPairLight() : 
+  TObject(), 
+  fMu0(),
+  fMu1(), 
+  fCreationProcess(-999),
+  fIsCorrelated(kFALSE), 
+  fCauseOfCorrelation (-1),
+  fIsFeedDown(kFALSE)
+{
+  /// default constructor
+  ; 
+}
+
+//====================================
+
+AliMUONPairLight::AliMUONPairLight(AliMUONPairLight &dimuCopy) 
+  : TObject(dimuCopy),
+    fMu0(dimuCopy.fMu0),
+    fMu1(dimuCopy.fMu0), 
+    fCreationProcess(dimuCopy.fCreationProcess),
+    fIsCorrelated(dimuCopy.fIsCorrelated), 
+    fCauseOfCorrelation (dimuCopy.fCauseOfCorrelation),
+    fIsFeedDown(dimuCopy.fIsFeedDown)
+{ 
+/// copy constructor
+///   fMu0 = AliMUONTrackLight(dimuCopy.fMu0); 
+///   fMu1 = AliMUONTrackLight(dimuCopy.fMu1); 
+///   fIsCorrelated = dimuCopy.fIsCorrelated;
+///   fCauseOfCorrelation = dimuCopy.fCauseOfCorrelation;
+///   fCreationProcess = dimuCopy.fCreationProcess;
+///   fIsFeedDown = dimuCopy.fIsFeedDown;
+  ;
+}
+
+//====================================
+
+AliMUONPairLight::~AliMUONPairLight(){
+  /// destructor
+}
+
+//====================================
+
+Bool_t AliMUONPairLight::IsAResonance(){
+  /// checks if muon pair comes from a resonance decay  
+  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");
+
+}
+
+//====================================
+
+AliMUONTrackLight* AliMUONPairLight::GetMuon(Int_t index)  { 
+  /// return muon 0 or 1
+   if (index==0) return &fMu0;
+   else if (index==1) return &fMu1; 
+   else{ printf ("Index can be either 0 or 1\n"); return 0;}
+   //   else return &fMu1; 
+}
+
+//====================================
+
+Int_t AliMUONPairLight::GetMuonMotherPDG(Int_t imuon, Int_t mother) { 
+  /// return muon mother pdg code
+  if (imuon==0) return fMu0.GetParentPDGCode(mother); 
+  else if (imuon==1) return fMu1.GetParentPDGCode(mother); 
+  else { printf ("Index must be only 0 or 1\n"); return -999; } 
+}
+
+//====================================
+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) { 
+       this->SetCorrelated(kTRUE); 
+       this->SetCauseOfCorrelation(mu1->GetParentPDGCode(imoth1));
+       if(!IsAResonance()) fCreationProcess = 3; 
+       else fCreationProcess = -1;
+       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));
+    }
+    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;
+    }
+  }
+
+  // 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;
+//     }
+  }
+}
+
+//====================================
+void AliMUONPairLight::SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1){
+  /// set the two muons 
+  fMu0 = mu0; 
+  fMu1 = mu1; 
+  this->SetProcess();
+} 
+
+//====================================
+void AliMUONPairLight::PrintInfo(Option_t* opt){
+  /// print information about muon pairs
+  /// Options: 
+  /// - "H" single muons' decay histories
+  /// - "K" dimuon kinematics
+  /// - "F" dimuon flags
+  /// - "A" all variables
+  TString options(opt);
+  options.ToUpper();
+
+  if(options.Contains("H") || options.Contains("A")){//muon decay histories
+
+    AliMUONTrackLight *mu1 = &fMu0;
+    AliMUONTrackLight *mu2 = &fMu1;
+
+    printf("========= History =======================\n");
+    printf("first muon");
+    mu1->PrintInfo("H");
+    printf("second muon");
+    mu2->PrintInfo("H");
+    printf("=========================================\n");
+  }
+  if(options.Contains("F") || options.Contains("A")){//flags
+    printf("the flags set for this muon pair are:\n");
+    printf("=====================================\n");
+    if(this->IsOneTrackNotAMuon()) printf("(*) one rec. track is not a muon\n");
+    fIsCorrelated ? printf("(*) it is a correlated pair\n") : printf("(*) it is not a correlated pair\n");
+    if(IsOpenCharm()) printf("(*) correlated open charm: ");
+    if(IsOpenBeauty()) printf("(*) correlated open beauty: ");
+    if(IsOpenCharm() || IsOpenBeauty()){
+      switch(fCreationProcess){
+      case 0:
+       printf("pair creation");
+       break;
+      case 1:
+       printf("gluon splitting");
+       break;
+      case 2:
+       printf("flavour excitation");
+       break;
+      case 3:
+       printf("both muons come from same fragmented mother");
+       break;
+      }
+      if(this->GetMuon(0)->GetOscillation() || this->GetMuon(1)->GetOscillation()) 
+       printf("... where oscillation occured\n");
+      else{
+       if(IsOpenBeauty())
+         printf(" (no oscillation)\n");
+       else
+         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");
+    printf("=====================================\n");
+  }
+  if(options.Contains("K") || options.Contains("A")){//dimuon kinematics
+    Double_t *vtx = this->GetMuon(0)->GetVertex();
+    TLorentzVector momRec = this->GetPRec();
+    TLorentzVector momGen = this->GetPGen();
+    printf("the dimuon charge is %d\n", this->GetCharge());
+    printf("primary Vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
+    printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
+    printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
+    //rapidity, pT, angles, ...
+    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", 
+          momRec.Pt(), momRec.Eta(), 
+          TMath::Pi()/180.*this->GetOpeningAngle(), this->GetOpeningAngle(), 
+          momRec.Theta(), 180./TMath::Pi() * momRec.Theta(), 
+          momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
+  }
+}
+
+Double_t AliMUONPairLight::GetOpeningAngle() { 
+  // opening angle between the two muons in the lab frame (in degrees)
+  TLorentzVector pRecMu0 =  fMu0.GetPRec();
+  TLorentzVector pRecMu1 =  fMu1.GetPRec();
+  TVector3 pRecMu03 = pRecMu0.Vect();
+  TVector3 pRecMu13 = pRecMu1.Vect();
+  Double_t scalar = pRecMu03.Dot(pRecMu13);
+  Double_t modMu0 = pRecMu03.Mag();
+  Double_t modMu1 = pRecMu13.Mag();
+  Double_t theta = (TMath::ACos(scalar/(modMu0*modMu1)))*(180./TMath::Pi());
+  return theta; 
+}
diff --git a/MUON/AliMUONPairLight.h b/MUON/AliMUONPairLight.h
new file mode 100644 (file)
index 0000000..681576f
--- /dev/null
@@ -0,0 +1,73 @@
+#ifndef ALIMUONPAIRLIGHT_H
+#define ALIMUONPAIRLIGHT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+// Revision of includes 06/09/2006
+
+/// \ingroup analysis
+/// \class AliMUONPairLight
+/// \brief Compact information for the generated muon pairs 
+/// 
+/// Compact information for the generated muon pairs in the MUON arm 
+/// useful at the last stage of the analysis chain
+/// Pairs are built with two AliMUONTrackLight objects 
+/// Using the class AliMUONTrackLight this class combines the decay
+/// information ("history") of the reconstructed tracks and fills
+/// a series of flags for the formed reconstructed dimuon:
+/// fIsCorrelated, fCreationProcess, fIsFeedDown, ...
+/// for information about the dimuon, use PrintInfo with the appropriate
+/// printflag
+/// To be used together with AliMUONTrackLight
+///
+/// \author This class was prepared by INFN Cagliari, July 2006
+/// (authors: H.Woehri, A.de Falco)
+
+// MUON classes
+#include "AliMUONTrackLight.h"
+
+// ROOT classes
+//#include "TLorentzVector.h"
+class TLorentzVector;
+
+class AliMUONPairLight : public TObject { 
+public:
+  AliMUONPairLight();
+  AliMUONPairLight(AliMUONPairLight &dimuCopy);
+  virtual ~AliMUONPairLight();
+  virtual void SetMuons(AliMUONTrackLight mu0, AliMUONTrackLight mu1);
+  AliMUONTrackLight* GetMuon(Int_t index) ; 
+  Int_t GetMuonMotherPDG(Int_t imuon, Int_t mother=0) ; 
+  Bool_t IsOpenCharm() {return (TMath::Abs(fMu0.GetParentFlavour(0))==4 && TMath::Abs(fMu1.GetParentFlavour(0))==4 && fIsCorrelated && !IsAResonance());}
+  Bool_t IsOpenBeauty() {return (TMath::Abs(fMu0.GetParentFlavour(0))==5 && TMath::Abs(fMu1.GetParentFlavour(0))==5 && fIsCorrelated  && !IsAResonance());}
+  Bool_t IsAResonance(); 
+  Bool_t IsCorrelated() const {return fIsCorrelated;}   
+  Int_t GetCauseOfCorrelation() const {return fCauseOfCorrelation;}
+  Bool_t IsFeedDown() const {return fIsFeedDown;}   
+  Bool_t IsOneTrackNotAMuon() {return (!( fMu0.IsAMuon() && fMu1.IsAMuon() )) ;}
+  Int_t GetCharge() {return fMu0.GetCharge() + fMu1.GetCharge();}
+  Int_t GetCreationProcess() const  {return fCreationProcess;}
+  void SetCorrelated(Bool_t answer) {fIsCorrelated = answer; }
+  void SetCauseOfCorrelation(Int_t pdg) {fCauseOfCorrelation = pdg; }
+  void SetFeedDown(Int_t answer) {fIsFeedDown = answer;}
+  TLorentzVector GetPRec(){return fMu0.GetPRec()+fMu1.GetPRec();}
+  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 
+
+  void SetProcess();///checks if muons are correlated and assigns
+
+  ClassDef(AliMUONPairLight,1)  /// Dimuon carrying info at generation and reconstruction level
+}; 
+#endif
diff --git a/MUON/AliMUONTrackLight.cxx b/MUON/AliMUONTrackLight.cxx
new file mode 100644 (file)
index 0000000..29102fe
--- /dev/null
@@ -0,0 +1,350 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *      SigmaEffect_thetadegrees                                                                  *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpeateose. It is      *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//===================================================================
+//This class was prepared by INFN Cagliari, July 2006
+//(authors: H.Woehri, A.de Falco)
+// 
+// 
+// Compact information for the muon generated tracks in the MUON arm 
+// useful at the last stage of the analysis chain
+// provides a link between the reconstructed track and the generated particle 
+// stores kinematical information at gen. and rec. level and 
+// the decay history of the muon, allowing the identification of the 
+// mother process 
+// 
+// To be used together with AliMUONPairLight
+//===================================================================
+
+#include "AliMUONTrackLight.h"
+#include "AliMUONTrack.h"
+#include "AliMUONConstants.h"
+
+#include "AliESDMuonTrack.h"
+#include "AliRunLoader.h"
+#include "AliStack.h"
+#include "AliHeader.h"
+
+#include "TDatabasePDG.h"
+#include "TClonesArray.h"
+#include "TParticle.h"
+#include "TString.h"
+
+ClassImp(AliMUONTrackLight) 
+
+//===================================================================
+
+AliMUONTrackLight::AliMUONTrackLight() 
+  : TObject(), 
+    fPrec(), 
+    fIsTriggered(kFALSE),
+    fCharge(-999), 
+    fCentr(-1),
+    fPgen(), 
+    fTrackPythiaLine(-999),
+    fTrackPDGCode(-999),
+    fOscillation(kFALSE), 
+    fNParents(0)
+{
+  /// default constructor
+  fPgen.SetPxPyPzE(0.,0.,0.,0.); 
+  fPrec.SetPxPyPzE(0.,0.,0.,0.); 
+  for (Int_t i=0; i<3; i++) fXYZ[i]=-999; 
+  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
+    fParentPDGCode[npar] = -1; 
+    fParentPythiaLine[npar] = -1;
+  }
+  for (Int_t i = 0; i < 4; i++){
+    fQuarkPDGCode[i] = -1; 
+    fQuarkPythiaLine[i] = -1; 
+  }
+}
+
+//============================================
+AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy) 
+  : TObject(muonCopy), 
+    fPrec(muonCopy.fPrec), 
+    fIsTriggered(muonCopy.fIsTriggered),
+    fCharge(muonCopy.fCharge), 
+    fCentr(muonCopy.fCentr),
+    fPgen(muonCopy.fPgen), 
+    fTrackPythiaLine(muonCopy.fTrackPythiaLine),
+    fTrackPDGCode(muonCopy.fTrackPDGCode),
+    fOscillation(muonCopy.fOscillation), 
+    fNParents(muonCopy.fNParents)
+{
+  /// copy constructor
+  for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i]; 
+  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
+    fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar]; 
+    fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
+  }
+  for (Int_t i = 0; i < 4; i++){
+    fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i]; 
+    fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i]; 
+  }
+}
+
+//============================================
+AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
+  : TObject(), 
+    fPrec(), 
+    fIsTriggered(kFALSE),
+    fCharge(-999), 
+    fCentr(-1),
+    fPgen(), 
+    fTrackPythiaLine(-999),
+    fTrackPDGCode(-999),
+    fOscillation(kFALSE), 
+    fNParents(0)
+{ 
+  /// constructor
+  //AliMUONTrackLight(); 
+  ComputePRecAndChargeFromESD(muonTrack); 
+}
+
+//============================================
+void AliMUONTrackLight::ComputePRecAndChargeFromESD(AliESDMuonTrack* muonTrack){ 
+  /// computes prec and charge from ESD track
+  Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass(); 
+  Double_t thetaX = muonTrack->GetThetaX();
+  Double_t thetaY = muonTrack->GetThetaY();
+  Double_t tanthx = TMath::Tan(thetaX);
+  Double_t tanthy = TMath::Tan(thetaY);
+  Double_t pYZ    =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
+  Double_t pz     = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
+  Double_t px     = pz * tanthx;
+  Double_t py     = pz * tanthy;
+  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);
+}
+
+//============================================
+void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){ 
+  /// set the reconstructed 4-momentum, assuming the particle is a muon
+  Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass(); 
+  Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
+  fPrec.SetPxPyPzE(px,py,pz,energy);
+}
+
+//============================================
+TParticle* AliMUONTrackLight::FindRefTrack(AliMUONTrack* trackReco, TClonesArray* trackRefArray, AliRunLoader *runLoader){ 
+  /// find the MC particle that corresponds to a given rec track
+  TParticle *part = 0; 
+  const Double_t kSigma2Cut = 16;  // 4 sigmas cut, kSigma2Cut = 4*4
+  Int_t nTrackRef = trackRefArray->GetEntriesFast();
+  Int_t compPart = 0; 
+  for (Int_t iref = 0; iref < nTrackRef; iref++) {
+    AliMUONTrack *trackRef = (AliMUONTrack *)trackRefArray->At(iref);
+    // check if trackRef is compatible with trackReco:
+    //routine returns for each chamber a yes/no information if the
+    //hit of rec. track and hit of referenced track are compatible
+    Bool_t *compTrack = trackRef->CompatibleTrack(trackReco,kSigma2Cut);
+    Int_t iTrack = this->TrackCheck(compTrack); //returns number of validated conditions 
+    if (iTrack==4) { 
+      compPart++; 
+      Int_t trackID = trackRef->GetTrackID();
+      this->SetTrackPythiaLine(trackID);
+      part = ((AliStack *)(((AliHeader *) runLoader->GetHeader())->Stack()))->Particle(trackID);
+      fTrackPDGCode = part->GetPdgCode();
+    }
+  }
+  if (compPart>1) { 
+    printf ("<AliMUONTrackLight::FindRefTrack> ERROR: more than one particle compatible to the reconstructed track.\n"); 
+    Int_t i=0, j=1/i; 
+    printf ("j=%d \n",j); 
+  } 
+  return part; 
+}
+
+//============================================
+Int_t AliMUONTrackLight::TrackCheck(Bool_t *compTrack){
+  /// Apply reconstruction requirements
+  /// Return number of validated conditions 
+  /// If all the tests are verified then TrackCheck = 4 (good track)
+  Int_t iTrack = 0;
+  Int_t hitsInLastStations = 0;
+  
+  // apply reconstruction requirements
+  if (compTrack[0] || compTrack[1]) iTrack++; // at least one hit in st. 0
+  if (compTrack[2] || compTrack[3]) iTrack++; // at least one hit in st. 1
+  if (compTrack[4] || compTrack[5]) iTrack++; // at least one hit in st. 2
+  for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) {
+    if (compTrack[ch]) hitsInLastStations++; 
+  }
+  if (hitsInLastStations > 2) iTrack++; // at least 3 hits in st. 3 & 4
+  return iTrack;
+}
+
+//============================================
+void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, 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 
+    parents[++countP] = pdg;
+    parLine[countP] = lineM;
+    status = mother->GetStatusCode();//get its status code to check if oscillation occured
+    if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
+    lineM = mother->GetFirstMother();
+  }
+  //store all the fragmented parents in an array:
+  for(int i = 0; i <= countP; i++){
+    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);
+    pdg = mother->GetPdgCode();
+    //now, get information before the fragmentation
+    this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
+    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;
+  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);
+      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
+
+       pdg = stack->Particle(++line)->GetPdgCode();    
+      }
+      //now, we have the correct string end and the correct line
+      //continue to fill again all parents and corresponding lines
+      while(line >= 0){
+       mother = stack->Particle(line);//get again the mother
+       pdg = mother->GetPdgCode();
+       this->SetQuarkPythiaLine(countP, line);
+       this->SetQuarkPDGCode(countP++, pdg);
+       line = mother->GetFirstMother();
+      }
+    }
+  }
+}
+//====================================
+void AliMUONTrackLight::ResetQuarkInfo(){
+  /// resets parton information
+  for(int pos = 1; pos < 4; pos++){//[0] is the string
+    this->SetQuarkPDGCode(pos,-1);
+    this->SetQuarkPythiaLine(pos,-1);
+  }
+}
+
+//====================================
+Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
+  /// checks if the particle is a B0 
+  Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
+  Bool_t answer = kFALSE;
+  for(int i = 0; i < 2; i++){
+    if(TMath::Abs(intTest) == bMes0[i]){
+      answer = kTRUE;
+      break;
+    }
+  }
+  return answer;
+}
+//====================================
+Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
+  /// checks if mother is a resonance
+  Int_t intTest = GetParentPDGCode(index); 
+  // the resonance pdg code is built this way
+  // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour 
+  Int_t id=intTest%100000; 
+  return (!((id-id%10)%110));
+}
+//====================================
+Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
+  /// returns the flavour of parent idParent (idParent=0 is the oldest 
+  /// hadronized parent)
+  Int_t pdg = GetParentPDGCode(idParent); 
+  Int_t quark = TMath::Abs(pdg/100);
+  if(quark > 9) quark = quark/10;
+  return quark;
+}
+
+//====================================
+void AliMUONTrackLight::PrintInfo(Option_t* opt){
+  /// prints information about the track: 
+  /// - "H" muon's decay history
+  /// - "K" muon kinematics
+  /// - "A" all variables
+  TString options(opt);
+  options.ToUpper();
+
+  if(options.Contains("H") || options.Contains("A")){ //muon decay history
+    char *name= new char[100];
+    TString pdg = "", line = "";
+    for(int i = 3; i >= 0; i--){
+      if(this->GetQuarkPythiaLine(i)>= 0){
+       sprintf(name, "%4d --> ", this->GetQuarkPythiaLine(i));
+       line += name;
+       sprintf(name, "%4d --> ", this->GetQuarkPDGCode(i));
+       pdg += name;
+      }
+    }
+    for(int i = 0; i < fNParents; i++){ 
+      if(this->GetParentPythiaLine(i)>= 0){
+       sprintf(name, "%7d --> ", this->GetParentPythiaLine(i));
+       line += name;
+       sprintf(name, "%7d --> ", this->GetParentPDGCode(i));
+       pdg += name;
+      }
+    }
+    sprintf(name, "%4d", this->GetTrackPythiaLine()); line += name;
+    sprintf(name, "%4d", this->GetTrackPDGCode()); pdg += name;
+
+    printf("\nmuon's decay history:\n");
+    printf(" PDG: %s\n", pdg.Data());
+    printf("line: %s\n", line.Data());
+  }
+  if(options.Contains("K") || options.Contains("A")){ //muon kinematic
+
+    Int_t charge = this->GetCharge();
+    Double_t *vtx = this->GetVertex();
+    TLorentzVector momRec = this->GetPRec();
+    TLorentzVector momGen = this->GetPGen();
+    printf("the track's charge is %d\n", charge);
+    printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
+    printf("Generated:     Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
+    printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
+    printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n", 
+          momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(), 
+          momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
+  }
+}
+
+
diff --git a/MUON/AliMUONTrackLight.h b/MUON/AliMUONTrackLight.h
new file mode 100644 (file)
index 0000000..3471cbb
--- /dev/null
@@ -0,0 +1,99 @@
+#ifndef ALIMUONTRACKLIGHT_H
+#define ALIMUONTRACKLIGHT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+// Revision of includes 06/09/2006
+
+/// \ingroup analysis
+/// \class AliMUONTrackLight
+/// \brief Compact information for the muon generated tracks
+/// 
+/// Compact information for the muon generated tracks in the MUON arm 
+/// useful at the last stage of the analysis chain
+/// provides a link between the reconstructed track and the generated particle 
+/// stores kinematical information at gen. and rec. level and 
+/// the decay history of the muon, allowing the identification of the 
+/// mother process 
+/// 
+/// To be used together with AliMUONPairLight
+///
+/// \author This class was prepared by INFN Cagliari, July 2006
+/// (authors: H.Woehri, A.de Falco)
+
+// ROOT classes
+#include "TLorentzVector.h"
+
+class AliMUONTrack;
+class AliESDMuonTrack;
+class AliRunLoader;
+
+class TParticle;
+
+class AliMUONTrackLight : public TObject { 
+ public: 
+  AliMUONTrackLight(); 
+  AliMUONTrackLight(AliESDMuonTrack* muonTrack); 
+  AliMUONTrackLight(const AliMUONTrackLight &muonCopy);
+  virtual ~AliMUONTrackLight(){;} 
+  
+  void SetPGen(TLorentzVector pgen) {fPgen = pgen;}
+  TLorentzVector GetPGen() const {return fPgen;}
+  void SetPRec(TLorentzVector prec) {fPrec = prec;}
+  TLorentzVector GetPRec() const {return fPrec;}
+  void SetVertex(Double_t *xyz) {for (Int_t i=0; i<3; i++) fXYZ[i]=xyz[i];}
+  Double_t* GetVertex() { return fXYZ; } 
+  void SetCharge(Int_t charge) {fCharge = charge;}
+  Int_t GetCharge() const {return fCharge;}
+  Int_t GetParentPDGCode(Int_t index = 0) const { return fParentPDGCode[index]; } 
+  Int_t GetParentPythiaLine(Int_t index = 0) const { return fParentPythiaLine[index]; } 
+  Int_t GetQuarkPDGCode(Int_t index = 0) const { return fQuarkPDGCode[index]; } 
+  Int_t GetQuarkPythiaLine(Int_t index = 0) const { return fQuarkPythiaLine[index]; } 
+  Int_t GetTrackPythiaLine() const {return fTrackPythiaLine;}
+  Int_t GetTrackPDGCode() const {return fTrackPDGCode;}
+  void SetTrackPythiaLine(Int_t trackLine) {fTrackPythiaLine = trackLine;}
+  void SetTrackPDGCode(Int_t trackPdg) {fTrackPDGCode = trackPdg;}
+  void ComputePRecAndChargeFromESD(AliESDMuonTrack* muonTrack);
+  Bool_t IsAMuon() const { return (TMath::Abs(fTrackPDGCode)==13); }
+  void SetPxPyPz(Double_t px, Double_t py, Double_t pz); 
+  void SetTriggered(Bool_t isTriggered) { fIsTriggered = isTriggered; } 
+  Bool_t IsTriggered() const { return fIsTriggered; } 
+  TParticle* FindRefTrack(AliMUONTrack* trackReco, TClonesArray* trackRefArray, AliRunLoader *runLoader); 
+  Int_t TrackCheck(Bool_t *compTrack);
+  Int_t GetNParents() const {return fNParents;}
+  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;
+  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;
+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
+  void SetOscillation(Bool_t oscillation) { fOscillation = oscillation; }
+  void SetParentPDGCode(Int_t index, Int_t pdg) { fParentPDGCode[index] = pdg; } 
+  void SetParentPythiaLine(Int_t index, Int_t line) { fParentPythiaLine[index] = line; } 
+  void SetQuarkPDGCode(Int_t index, Int_t pdg){ fQuarkPDGCode[index] = pdg; }
+  void SetQuarkPythiaLine(Int_t index, Int_t line){ fQuarkPythiaLine[index] = line; }
+  void ResetQuarkInfo();
+
+  ClassDef(AliMUONTrackLight,1)  /// Muon Track for analysis
+}; 
+
+#endif