--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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());
+ }
+}
+
+
--- /dev/null
+#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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff