Skip to content

Commit 2d3f7fd

Browse files
comrademarvincomrademarvin
committed
revamp of HF task to look at full decay chain
1 parent fc08591 commit 2d3f7fd

2 files changed

Lines changed: 104 additions & 47 deletions

File tree

‎MyTasks/HF/HFMuonFwdTracks.cxx‎

Lines changed: 100 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -9,88 +9,141 @@
99
using namespace o2;
1010
using namespace o2::framework;
1111

12-
// output textfile for muon tracks
13-
std::ofstream muonTracksOut("muontracks.txt");
14-
1512
struct HFMuonFwdTracks {
1613
// Histogram registry: an object to hold your histograms
1714
HistogramRegistry histos{"histos", {},
1815
OutputObjHandlingPolicy::AnalysisObject};
1916

20-
// store id of selected tracks to account for ambiguous tracks
21-
std::vector<int64_t> selectedTracksID;
17+
// output textfile for muon tracks
18+
std::ofstream muonTracksOut;
2219

2320
void init(InitContext const&)
2421
{
22+
// open output file and write header
23+
muonTracksOut.open("muonTracks_HF.csv");
24+
if (!muonTracksOut.is_open()) {
25+
LOGF(fatal, "Failed to open muonTracks.csv for writing");
26+
}
27+
muonTracksOut << "trackID,trackType,phi,tgl,signed1Pt,nClusters,pDCA,rAtAbsorberEnd,sign,chi2,chi2MatchMCHMID,chi2MatchMCHMFT,trackTime,eta,pt,p" << std::endl;
28+
2529
// define axes you want to use
2630
const AxisSpec axisCounter{1, 0, +1, ""};
2731
const AxisSpec axisEta{10, -4.0, -2.5, "#eta"};
2832
const AxisSpec axisPt{10, 0.0, +20.0, "p_{T} (GeV/c)"};
33+
const AxisSpec axisTrackType{5, 0, 5, "Track Type"};
2934

3035
// create histograms
3136
histos.add("eventCounterReco", "eventCounterReco", kTH1F, {axisCounter});
3237
histos.add("muPtHistReco", "muPtHistReco", kTH1F, {axisPt});
33-
histos.add("muPtHistRecoD", "muPtHistRecoD", kTH1F, {axisPt});
38+
histos.add("trackType", "trackType", kTH1D, {axisTrackType});
39+
}
3440

35-
// set textfile header
36-
muonTracksOut << "ID,eta,pt,p,phi,motherPDG,nClusters,pDca,chi2,chi2MatchMCHMID,chi2MatchMCHMFT,isPrompt\n";
41+
Int_t GetFlavour(Int_t pdgCode)
42+
{
43+
//
44+
// return the flavour of a particle
45+
// input: pdg code of the particle
46+
// output: Int_t
47+
// 3 in case of strange (open and hidden)
48+
// 4 in case of charm (")
49+
// 5 in case of beauty (")
50+
//
51+
Int_t pdg = TMath::Abs(pdgCode);
52+
// Resonance
53+
if (pdg > 100000)
54+
pdg %= 100000;
55+
if (pdg > 10000)
56+
pdg %= 10000;
57+
// meson ?
58+
if (pdg > 10)
59+
pdg /= 100;
60+
// baryon ?
61+
if (pdg > 10)
62+
pdg /= 10;
63+
return pdg;
3764
}
3865

3966
using muonTracks = soa::Join<aod::FwdTracks, aod::McFwdTrackLabels>;
4067

41-
void process(aod::Collision const& collision, muonTracks const& tracks, aod::McParticles const&)
68+
void process(muonTracks const& tracks, aod::McParticles const&)
4269
{
4370
histos.fill(HIST("eventCounterReco"), 0.5);
4471
for (auto& track : tracks) {
4572
if(track.has_mcParticle()){
4673
auto mcParticle = track.mcParticle();
74+
int64_t recoTrackID = track.globalIndex(); // Reconstructed track ID
4775
//auto statusCode = mcParticle.getGenStatusCode();
4876
if (abs(mcParticle.pdgCode())==13) {
49-
auto muID = mcParticle.globalIndex();
50-
auto muEta = mcParticle.eta();
51-
if (muEta >= -4.0 && muEta <= -2.5) histos.fill(HIST("muPtHistReco"), mcParticle.pt()); // forward eta region
52-
53-
// check if duplicate (ambiguous track)
54-
int occuranceCount = count(selectedTracksID.begin(), selectedTracksID.end(), muID);
55-
56-
// write to textfile
57-
if (occuranceCount < 1) {
58-
auto muRecoEta = track.eta();
59-
auto muRecoPt = track.pt();
60-
auto muRecoP = track.p();
61-
auto muRecoPhi = track.phi();
62-
auto muCluster = track.nClusters();
63-
auto muDCA = track.pDca();
64-
auto muChi2 = track.chi2();
65-
auto muChi2MCHMID = track.chi2MatchMCHMID();
66-
auto muChi2MCHMFT = track.chi2MatchMCHMFT();
67-
68-
// look for HF mother
69-
auto muMother = mcParticle.mothers_first_as<aod::McParticles>();
70-
auto muMotherPDG = abs(muMother.pdgCode());
71-
auto motherStatusCode = muMother.getGenStatusCode();
72-
if (abs(muMotherPDG) >= 411 && abs(muMotherPDG) <= 435) { // D meson parent - can be prompt/no-prompt
73-
histos.fill(HIST("muPtHistRecoD"), mcParticle.pt());
74-
}
77+
int64_t mcTrackID = mcParticle.globalIndex(); // Simulated track ID
78+
79+
// check if muon is from HF decay
80+
auto muMother = mcParticle.mothers_first_as<aod::McParticles>();
81+
auto muMotherPDG = abs(muMother.pdgCode());
82+
auto muMotherFlavour = GetFlavour(muMotherPDG);
7583

76-
// check whether is prompt HF hadron
84+
bool hasHFmother = false;
85+
std::cout << "==== Mu decay chain: 13 <- " << muMotherPDG;
86+
87+
if (muMotherFlavour != 4 && muMotherFlavour != 5) { // direct mother is not HF hadron
88+
// check if HF hadron mother in decay chain
7789
auto mcPart(muMother);
78-
auto prevMcPart(muMother);
7990
auto mcPartPDG = abs(mcPart.pdgCode());
80-
int isPrompt = 1;
81-
std::cout << "==== Forward muon decay chain: mu";
82-
while (mcPart.has_mothers() && (abs(mcPart.getGenStatusCode()) > 80 || mcPart.getGenStatusCode() == 0)) { // print out mother chain
83-
std::cout << " <- " << mcPartPDG;
84-
prevMcPart = *(mcPart);
91+
auto mcPartFlavour = GetFlavour(mcPartPDG);
92+
93+
while (mcPart.has_mothers()) {
8594
mcPart = *(mcPart.mothers_first_as<aod::McParticles>());
8695
mcPartPDG = abs(mcPart.pdgCode());
87-
}
88-
if (div(abs(prevMcPart.pdgCode()), 100).quot != div(muMotherPDG, 100).quot) isPrompt = 0;
89-
std::cout << "; isPrompt = " << isPrompt << std::endl;
96+
mcPartFlavour = GetFlavour(mcPartPDG);
97+
std::cout << " <- " << mcPartPDG;
9098

91-
muonTracksOut << muID << "," << muRecoEta << "," << muRecoPt << "," << muRecoP << "," << muRecoPhi << "," << muMotherPDG << "," << std::to_string(muCluster) << "," << muDCA << "," << muChi2 << "," << muChi2MCHMID << "," << muChi2MCHMFT << "," << isPrompt << std::endl;
92-
selectedTracksID.emplace(selectedTracksID.end(), muID);
99+
if (mcPartFlavour == 4 || mcPartFlavour == 5) {
100+
hasHFmother = true;
101+
break;
102+
}
103+
}
104+
} else {
105+
hasHFmother = true;
93106
}
107+
if (hasHFmother) std::cout << " (PICKED)";
108+
std::cout << std::endl;
109+
110+
if (hasHFmother) {
111+
auto muTrackType = static_cast<int64_t>(track.trackType());
112+
auto muChi2 = track.chi2();
113+
auto muChi2MatchMCHMID = track.chi2MatchMCHMID();
114+
auto muChi2MatchMCHMFT = track.chi2MatchMCHMFT();
115+
auto muMatchScoreMCHMFT = track.matchScoreMCHMFT();
116+
auto muMatchMFTTrackId = track.matchMFTTrackId();
117+
auto muDca = track.pDca();
118+
auto muPt = track.pt();
119+
auto muEta = track.eta();
120+
auto ptResolution = abs(mcParticle.pt() - track.pt());
121+
122+
histos.fill(HIST("trackType"), muTrackType);
123+
124+
if (muTrackType == 3) histos.fill(HIST("muPtHistReco"), muPt); // look at pT for standalone tracks
125+
126+
// save all muon tracks to the output file
127+
muonTracksOut << recoTrackID << "," << muTrackType << "," << track.phi() << "," << track.tgl() << "," << track.signed1Pt() << ","
128+
<< static_cast<int64_t>(track.nClusters()) << "," << muDca << "," << track.rAtAbsorberEnd() << ","
129+
<< static_cast<int64_t>(track.sign()) << "," << muChi2 << "," << muChi2MatchMCHMID << "," << muChi2MatchMCHMFT << ","
130+
<< track.trackTime() << "," << muEta << "," << muPt << "," << track.p() << std::endl;
131+
};
132+
133+
// check whether is prompt HF hadron
134+
// auto mcPart(muMother);
135+
// auto prevMcPart(muMother);
136+
// auto mcPartPDG = abs(mcPart.pdgCode());
137+
// int isPrompt = 1;
138+
// std::cout << "==== Forward muon decay chain: mu";
139+
// while (mcPart.has_mothers() && (abs(mcPart.getGenStatusCode()) > 80 || mcPart.getGenStatusCode() == 0)) { // print out mother chain
140+
// std::cout << " <- " << mcPartPDG;
141+
// prevMcPart = *(mcPart);
142+
// mcPart = *(mcPart.mothers_first_as<aod::McParticles>());
143+
// mcPartPDG = abs(mcPart.pdgCode());
144+
// }
145+
// if (div(abs(prevMcPart.pdgCode()), 100).quot != div(muMotherPDG, 100).quot) isPrompt = 0;
146+
// std::cout << "; isPrompt = " << isPrompt << std::endl;
94147
}
95148
}
96149
}

‎MyTasks/W/wMuonFwdEfficiency.cxx‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ struct wMuonFwdEfficiency {
5252
const AxisSpec axisChi2{50, 0.0, +10.0, "#chi^{2}"};
5353
const AxisSpec axisChi2Global{50, 0.0, +100.0, "#chi^{2}"};
5454
const AxisSpec axisDCA{100, 0.0, +1000.0, "pDCA"};
55+
const AxisSpec axisTrackType{5, 0, 5, "Track Type"};
5556

5657
// create histograms
5758
histos.add("eventCounterReco", "eventCounterReco", kTH1F, {axisCounter});
@@ -66,6 +67,7 @@ struct wMuonFwdEfficiency {
6667
histos.add("chi2MatchMCHMID", "chi2MatchMCHMID", kTH1F, {axisChi2});
6768
histos.add("chi2MatchMCHMFT", "chi2MatchMCHMFT", kTH1F, {axisChi2});
6869
histos.add("pDCA", "pDCA", kTH1F, {axisDCA});
70+
histos.add("trackType", "trackType", kTH1D, {axisTrackType});
6971
}
7072

7173
using muonTracks = soa::Join<aod::FwdTracks, aod::McFwdTrackLabels>;
@@ -134,6 +136,8 @@ struct wMuonFwdEfficiency {
134136
<< static_cast<int64_t>(track.sign()) << "," << muChi2 << "," << muChi2MatchMCHMID << "," << muChi2MatchMCHMFT << ","
135137
<< track.trackTime() << "," << muEta << "," << muPt << "," << track.p() << std::endl;
136138

139+
histos.fill(HIST("trackType"), muTrackType);
140+
137141
// basic cuts before plotting
138142
if (muTrackType == 3 && muChi2MatchMCHMID > 0) { // only look at standalone tracks for now
139143
// write to histograms

0 commit comments

Comments
 (0)