Skip to content

Commit 7a87e70

Browse files
comrademarvincomrademarvin
committed
updated filtering for duplicate tracks
1 parent 9aa8f01 commit 7a87e70

1 file changed

Lines changed: 93 additions & 24 deletions

File tree

‎MyTasks/W/wMuonFwdEfficiency.cxx‎

Lines changed: 93 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "Framework/AnalysisDataModel.h"
44

55
#include <iostream>
6+
#include <fstream>
67
#include <vector>
78

89
using namespace o2;
@@ -13,17 +14,37 @@ struct wMuonFwdEfficiency {
1314
HistogramRegistry histos{"histos", {},
1415
OutputObjHandlingPolicy::AnalysisObject};
1516

16-
// store id of selected tracks to account for ambiguous tracks
17-
std::vector<int64_t> selectedTracksID;
17+
// definition of reconstructed info to be saved
18+
struct RecoTrackInfo {
19+
int64_t trackID;
20+
float chi2;
21+
float chi2MatchMCHMID;
22+
float pDca;
23+
float pt;
24+
float eta;
25+
};
26+
// Map to group reconstructed tracks by their associated simulated track ID
27+
std::unordered_map<int64_t, std::vector<RecoTrackInfo>> trackGroups;
28+
29+
// output textfile for muon tracks
30+
std::ofstream muonTracksOut;
1831

1932
void init(InitContext const&)
2033
{
34+
// open output file and write header
35+
muonTracksOut.open("muonTracks.csv");
36+
if (!muonTracksOut.is_open()) {
37+
LOGF(fatal, "Failed to open muonTracks.csv for writing");
38+
}
39+
muonTracksOut << "trackID,chi2,chi2MatchMCHMID,pDCA,pt,eta" << std::endl;
40+
2141
// define axes you want to use
2242
const AxisSpec axisCounter{1, 0, +1, ""};
2343
const AxisSpec axisEta{20, -4.0, -2.5, "#eta"};
2444
const AxisSpec axisPt{20, 0.0, +80.0, "p_{T} (GeV/c)"};
2545
const AxisSpec axisDeltaPt{40, 0.0, +20.0, "|p_{T}^{true} - p_{T}^{reco}|(GeV/c)"};
26-
const AxisSpec axisChi2{100, 0.0, +100.0, "#chi^2"};
46+
const AxisSpec axisChi2{50, 0.0, +10.0, "#chi^{2}"};
47+
const AxisSpec axisDCA{100, 0.0, +1000.0, "pDCA"};
2748

2849
// create histograms
2950
histos.add("eventCounterReco", "eventCounterReco", kTH1F, {axisCounter});
@@ -32,36 +53,33 @@ struct wMuonFwdEfficiency {
3253
histos.add("yPtTruthHist", "yPtTruthHist", kTH2F, {axisPt, axisEta});
3354
histos.add("PtRecoHist", "PtRecoHist", kTH1F, {axisPt});
3455
histos.add("PtTruthHist", "PtTruthHist", kTH1F, {axisPt});
35-
histos.add("PtResolution", "PtResolution", kTH1F, {axisDeltaPt});
56+
//histos.add("PtResolution", "PtResolution", kTH1F, {axisDeltaPt});
3657
histos.add("chi2", "chi2", kTH1F, {axisChi2});
58+
histos.add("chi2MatchMCHMID", "chi2MatchMCHMID", kTH1F, {axisChi2});
59+
histos.add("pDCA", "pDCA", kTH1F, {axisDCA});
3760
}
3861

3962
using muonTracks = soa::Join<aod::FwdTracks, aod::McFwdTrackLabels>;
4063

41-
void processReco(aod::Collision const& collision, muonTracks const& tracks, aod::McParticles const&) // run with collisions
42-
//void processReco(muonTracks const& tracks, aod::McParticles const&) // run without collisions
64+
//void processReco(aod::Collision const& collision, muonTracks const& tracks, aod::McParticles const&) // run with collisions
65+
void processReco(muonTracks const& tracks, aod::McParticles const&) // run without collisions
4366
{
4467
histos.fill(HIST("eventCounterReco"), 0.5);
68+
4569
for (auto& track : tracks) {
46-
if(track.has_mcParticle()){
70+
if(track.has_mcParticle() && track.chi2MatchMCHMID() > 0) { // check if track has associated MC particle and has a valid chi2 match
4771
auto mcParticle = track.mcParticle();
48-
auto statusCode = abs(mcParticle.getGenStatusCode());
72+
int64_t recoTrackID = track.globalIndex(); // Reconstructed track ID
4973

5074
if (abs(mcParticle.pdgCode())==13) { // check if associated MC track is a muon
51-
// check if duplicate (ambiguous track)
52-
auto muID = mcParticle.globalIndex();
53-
int occuranceCount = count(selectedTracksID.begin(), selectedTracksID.end(), muID);
54-
55-
if (occuranceCount > 0) { // if already selected, move on
56-
continue;
57-
}
75+
int64_t mcTrackID = mcParticle.globalIndex(); // Simulated track ID
5876

5977
// check if muon is from W decay
6078
auto muMother = mcParticle.mothers_first_as<aod::McParticles>();
6179
auto muMotherPDG = abs(muMother.pdgCode());
6280

6381
bool hasWmother = false;
64-
std::cout << "==== Mu decay chain: 13 <- " << muMotherPDG;
82+
//std::cout << "==== Mu decay chain: 13 <- " << muMotherPDG;
6583

6684
if (muMotherPDG != 24) {
6785
// check if W mother in decay chain
@@ -71,7 +89,7 @@ struct wMuonFwdEfficiency {
7189
while (mcPart.has_mothers()) {
7290
mcPart = *(mcPart.mothers_first_as<aod::McParticles>());
7391
mcPartPDG = abs(mcPart.pdgCode());
74-
std::cout << " <- " << mcPartPDG;
92+
//std::cout << " <- " << mcPartPDG;
7593

7694
if (mcPartPDG == 24) {
7795
hasWmother = true;
@@ -81,17 +99,68 @@ struct wMuonFwdEfficiency {
8199
} else {
82100
hasWmother = true;
83101
}
84-
if (hasWmother) std::cout << " (PICKED)";
85-
std::cout << std::endl;
102+
// if (hasWmother) std::cout << " (PICKED)";
103+
// std::cout << std::endl;
86104

87105
if (hasWmother) {
88-
selectedTracksID.emplace(selectedTracksID.end(), muID); // add track ID to already processed selected tracks
89-
histos.fill(HIST("PtRecoHist"), mcParticle.pt());
90-
histos.fill(HIST("yPtRecoHist"), mcParticle.pt(), mcParticle.eta());
91-
histos.fill(HIST("PtResolution"), abs(mcParticle.pt()-track.pt()));
92-
histos.fill(HIST("chi2"), track.chi2());
106+
auto muChi2 = track.chi2();
107+
auto muChi2MatchMCHMID = track.chi2MatchMCHMID();
108+
auto pDca = track.pDca();
109+
auto pt = track.pt();
110+
auto eta = track.eta();
111+
112+
// Store the reconstructed track information in the map
113+
RecoTrackInfo trackInfo = {recoTrackID, muChi2, muChi2MatchMCHMID, pDca, pt, eta};
114+
trackGroups[mcTrackID].push_back(trackInfo);
115+
}
116+
}
117+
}
118+
}
119+
120+
for (auto& [mcTrackID, recoTracks] : trackGroups) {
121+
// prcess tracks to choose track with lowest chi2
122+
int trackIndex = 0;
123+
int chosenIndex = trackIndex;
124+
float minChi2 = 100000.0;
125+
for (const auto& trackInfo : recoTracks) {
126+
if (trackInfo.chi2 < minChi2) {
127+
minChi2 = trackInfo.chi2;
128+
chosenIndex = trackIndex; // remember the index of the track with the lowest chi2
129+
}
130+
trackIndex++;
131+
}
132+
133+
// write to histograms
134+
histos.fill(HIST("PtRecoHist"), recoTracks[chosenIndex].pt);
135+
histos.fill(HIST("yPtRecoHist"), recoTracks[chosenIndex].pt, recoTracks[chosenIndex].eta);
136+
//histos.fill(HIST("PtResolution"), abs(mcParticle.pt-track.pt));
137+
histos.fill(HIST("chi2"), recoTracks[chosenIndex].chi2);
138+
histos.fill(HIST("chi2MatchMCHMID"), recoTracks[chosenIndex].chi2MatchMCHMID);
139+
histos.fill(HIST("pDCA"), recoTracks[chosenIndex].pDca);
140+
141+
// output to textfile
142+
muonTracksOut << recoTracks[chosenIndex].trackID << ","
143+
<< recoTracks[chosenIndex].chi2 << ","
144+
<< recoTracks[chosenIndex].chi2MatchMCHMID << ","
145+
<< recoTracks[chosenIndex].pDca << ","
146+
<< recoTracks[chosenIndex].pt << ","
147+
<< recoTracks[chosenIndex].eta << std::endl;
148+
149+
// only print if the first track isn't the best track
150+
if (chosenIndex > 0) {
151+
// printing
152+
if (recoTracks.size() > 1) {
153+
// Output a heading for the simulated track
154+
LOGF(info, "Simulated track ID %ld has multiple associated reconstructed tracks:", mcTrackID);
155+
156+
// List the associated reconstructed tracks
157+
for (const auto& trackInfo : recoTracks) {
158+
// Option 1: Using RecoTrackInfo struct
159+
LOGF(info, "RecoTrack: ID = %ld, chi2 = %.2f, chi2MatchMCHMID = %.2f, pDCA = %.2f",
160+
trackInfo.trackID, trackInfo.chi2, trackInfo.chi2MatchMCHMID, trackInfo.pDca);
93161
}
94162
}
163+
LOGF(info, "Chosen track has index %d with track ID %ld", chosenIndex, recoTracks[chosenIndex].trackID);
95164
}
96165
}
97166
}

0 commit comments

Comments
 (0)