Skip to content

Commit e4228e0

Browse files
nzardoshnzardosh
committed
adding possibility to select multiple cluster definitions for hadronic correction
1 parent a3dcedd commit e4228e0

File tree

1 file changed

+137
-123
lines changed

1 file changed

+137
-123
lines changed

‎PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx‎

Lines changed: 137 additions & 123 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ struct EmcalClusterHadronicCorrectionTask {
4545
PresliceUnsorted<aod::JEMCTracks> perTrackMatchedTrack = aod::jemctrack::trackId;
4646

4747
// define configurables here
48-
Configurable<std::string> clusterDefinitionS{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"};
48+
Configurable<std::string> clusterDefinitions{"clusterDefinitions", "kV3Default", "cluster definitions to be selected, e.g. V3Default"};
4949

5050
Configurable<float> minTime{"minTime", -25., "Minimum cluster time for time cut"};
5151
Configurable<float> maxTime{"maxTime", 20., "Maximum cluster time for time cut"};
@@ -75,6 +75,8 @@ struct EmcalClusterHadronicCorrectionTask {
7575
Configurable<bool> useFraction1{"useFraction1", false, "Fractional momentum subtraction for clusterE1 and clusterEAll1"};
7676
Configurable<bool> useFraction2{"useFraction2", false, "Fractional momentum subtraction for clusterE2 and clusterEAll2"};
7777

78+
std::vector<int> clusterDefinitionsVec;
79+
7880
void init(o2::framework::InitContext&)
7981
{
8082
// Event histograms
@@ -97,15 +99,20 @@ struct EmcalClusterHadronicCorrectionTask {
9799

98100
// Matched-Track histograms
99101
registry.add("h_matchedtracks", "Total matched tracks; track status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}});
100-
}
101102

102-
aod::EMCALClusterDefinition clusterDefinition = aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionS.value);
103-
Filter clusterDefinitionSelection = (aod::jcluster::definition == static_cast<int>(clusterDefinition));
103+
std::string clusterDefinitionsString = clusterDefinitions.value;
104+
size_t start = 0;
105+
size_t end;
106+
while ((end = clusterDefinitionsString.find(',', start)) != std::string::npos) {
107+
clusterDefinitionsVec.push_back(static_cast<int>(aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionsString.substr(start, end - start))));
108+
start = end + 1;
109+
}
110+
}
104111

105112
// The matching of clusters and tracks is already centralised in the EMCAL framework.
106113
// One only needs to apply a filter on matched clusters
107114
// Here looping over all collisions matched to EMCAL clusters
108-
void processMatchedCollisions(aod::JetCollision const&, soa::Filtered<soa::Join<aod::JClusters, aod::JClusterTracks>> const& clusters, aod::JEMCTracks const& emcTracks, aod::JetTracks const&)
115+
void processMatchedCollisions(aod::JetCollision const&, soa::Join<aod::JClusters, aod::JClusterTracks> const& clusters, aod::JEMCTracks const& emcTracks, aod::JetTracks const&)
109116
{
110117
registry.fill(HIST("h_allcollisions"), 1);
111118

@@ -114,35 +121,132 @@ struct EmcalClusterHadronicCorrectionTask {
114121
return;
115122
}
116123

117-
// Looping over all clusters matched to the collision
118-
for (const auto& cluster : clusters) {
119-
registry.fill(HIST("h_matchedclusters"), 1);
120-
121-
double clusterE1;
122-
double clusterE2;
123-
double clusterEAll1;
124-
double clusterEAll2;
125-
clusterE1 = clusterE2 = clusterEAll1 = clusterEAll2 = cluster.energy();
124+
// Looping over all cluster definitions
125+
for (const auto& clusterDefinition : clusterDefinitionsVec) {
126+
// Looping over all clusters matched to the collision
127+
for (const auto& cluster : clusters) {
128+
if (cluster.definition() != clusterDefinition) {
129+
continue; // Skip clusters that do not match the current cluster definition
130+
}
126131

127-
registry.fill(HIST("h_ClsE"), cluster.energy());
128-
registry.fill(HIST("h_ClsM02"), cluster.m02());
129-
registry.fill(HIST("h_ClsTime"), cluster.time());
132+
registry.fill(HIST("h_matchedclusters"), 1);
133+
134+
double clusterE1;
135+
double clusterE2;
136+
double clusterEAll1;
137+
double clusterEAll2;
138+
clusterE1 = clusterE2 = clusterEAll1 = clusterEAll2 = cluster.energy();
139+
140+
registry.fill(HIST("h_ClsE"), cluster.energy());
141+
registry.fill(HIST("h_ClsM02"), cluster.m02());
142+
registry.fill(HIST("h_ClsTime"), cluster.time());
143+
144+
int nMatches = 0; // counter for closest matched track
145+
double closestTrkP = 0.0; // closest track momentum
146+
double totalTrkP = 0.0; // total track momentum
147+
148+
// pT-dependent track-matching instead of PID based track-matching to be adapted from Run 2 - suggested by Markus Fasel
149+
150+
TF1 funcPtDepEta("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
151+
funcPtDepEta.SetParameters(eta0, eta1, eta2);
152+
TF1 funcPtDepPhi("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
153+
funcPtDepEta.SetParameters(phi0, phi1, phi2);
154+
155+
// No matched tracks (trackless case)
156+
if (cluster.matchedTracks().size() == 0) {
157+
// Use original cluster energy values, no subtraction needed.
158+
registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), 0);
159+
registry.fill(HIST("h_Ecluster1"), clusterE1);
160+
registry.fill(HIST("h_Ecluster2"), clusterE2);
161+
registry.fill(HIST("h_EclusterAll1"), clusterEAll1);
162+
registry.fill(HIST("h_EclusterAll2"), clusterEAll2);
163+
registry.fill(HIST("h2_ClsEvsEcluster1"), cluster.energy(), clusterE1);
164+
registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), clusterE2);
165+
registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), clusterEAll1);
166+
registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), clusterEAll2);
167+
clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2);
168+
continue;
169+
}
130170

131-
int nMatches = 0; // counter for closest matched track
132-
double closestTrkP = 0.0; // closest track momentum
133-
double totalTrkP = 0.0; // total track momentum
171+
// Looping over all matched tracks for the cluster
172+
// Total number of matched tracks = 20 (hard-coded)
173+
for (const auto& matchedTrack : cluster.matchedTracks_as<aod::JetTracks>()) {
174+
if (matchedTrack.pt() < minTrackPt) {
175+
continue;
176+
}
177+
double mom = std::abs(matchedTrack.p());
178+
registry.fill(HIST("h_matchedtracks"), 1);
179+
180+
// CASE 1: skip tracks with a very low pT
181+
constexpr double kMinMom = 1e-6;
182+
if (mom < kMinMom) {
183+
continue;
184+
} // end CASE 1
185+
186+
// CASE 2:
187+
// a) If one matched track -> 100% energy subtraction
188+
// b) If more than one matched track -> 100% energy subtraction
189+
// c) If you want to do systematic studies -> perform the above two checks a) and b), and then subtract 70% energy instead of 100%
190+
191+
// Perform dEta/dPhi matching
192+
auto emcTrack = (emcTracks.sliceBy(perTrackMatchedTrack, matchedTrack.globalIndex())).iteratorAt(0);
193+
double dEta = emcTrack.etaDiff();
194+
double dPhi = emcTrack.phiDiff();
195+
196+
// Apply the eta and phi matching thresholds
197+
// dEta and dPhi cut : ensures that the matched track is within the desired eta/phi window
198+
199+
// Do pT-dependent track matching
200+
if (doMomDepMatching) {
201+
auto trackEtaMax = funcPtDepEta.Eval(mom);
202+
auto trackPhiHigh = +funcPtDepPhi.Eval(mom);
203+
auto trackPhiLow = -funcPtDepPhi.Eval(mom);
204+
205+
if ((dPhi < trackPhiHigh && dPhi > trackPhiLow) && std::fabs(dEta) < trackEtaMax) {
206+
if (nMatches == 0) {
207+
closestTrkP = mom;
208+
}
209+
totalTrkP += mom;
210+
nMatches++;
211+
}
212+
} else {
213+
// Do fixed dEta/dPhi matching (non-pT dependent)
214+
if (std::fabs(dEta) >= minDEta || std::fabs(dPhi) >= minDPhi) {
215+
continue; // Skip this track if outside the fixed cut region
216+
}
134217

135-
// pT-dependent track-matching instead of PID based track-matching to be adapted from Run 2 - suggested by Markus Fasel
218+
// If track passes fixed dEta/dPhi cuts, process it
219+
if (nMatches == 0) {
220+
closestTrkP = mom; // Closest track match
221+
}
222+
totalTrkP += mom; // Accumulate momentum
223+
nMatches++; // Count this match
224+
}
136225

137-
TF1 funcPtDepEta("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
138-
funcPtDepEta.SetParameters(eta0, eta1, eta2);
139-
TF1 funcPtDepPhi("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
140-
funcPtDepEta.SetParameters(phi0, phi1, phi2);
226+
} // End of track loop
227+
registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), nMatches);
228+
229+
if (nMatches >= 1) {
230+
if (useM02SubtractionScheme1) {
231+
// Do M02-based correction if enabled
232+
clusterE1 = subtractM02ClusterEnergy(cluster.m02(), clusterE1, nMatches, closestTrkP, hadCorr1, useFraction1);
233+
clusterEAll1 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll1, nMatches, totalTrkP, hadCorralltrks1, useFraction1);
234+
} else {
235+
// Default energy subtraction (100% and 70%)
236+
clusterE1 = subtractClusterEnergy(clusterE1, closestTrkP, hadCorr1, nMatches, useFraction1);
237+
clusterEAll1 = subtractClusterEnergy(clusterEAll1, totalTrkP, hadCorralltrks1, nMatches, useFraction1);
238+
}
141239

142-
// No matched tracks (trackless case)
143-
if (cluster.matchedTracks().size() == 0) {
144-
// Use original cluster energy values, no subtraction needed.
145-
registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), 0);
240+
if (useM02SubtractionScheme2) {
241+
// Do M02-based correction if enabled
242+
clusterE2 = subtractM02ClusterEnergy(cluster.m02(), clusterE2, nMatches, closestTrkP, hadCorr2, useFraction2);
243+
clusterEAll2 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll2, nMatches, totalTrkP, hadCorralltrks2, useFraction2);
244+
} else {
245+
// Default energy subtraction (100% and 70%)
246+
clusterE2 = subtractClusterEnergy(clusterE2, closestTrkP, hadCorr2, nMatches, useFraction2);
247+
clusterEAll2 = subtractClusterEnergy(clusterEAll2, totalTrkP, hadCorralltrks2, nMatches, useFraction2);
248+
}
249+
}
146250
registry.fill(HIST("h_Ecluster1"), clusterE1);
147251
registry.fill(HIST("h_Ecluster2"), clusterE2);
148252
registry.fill(HIST("h_EclusterAll1"), clusterEAll1);
@@ -151,102 +255,12 @@ struct EmcalClusterHadronicCorrectionTask {
151255
registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), clusterE2);
152256
registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), clusterEAll1);
153257
registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), clusterEAll2);
154-
clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2);
155-
continue;
156-
}
157-
158-
// Looping over all matched tracks for the cluster
159-
// Total number of matched tracks = 20 (hard-coded)
160-
for (const auto& matchedTrack : cluster.matchedTracks_as<aod::JetTracks>()) {
161-
if (matchedTrack.pt() < minTrackPt) {
162-
continue;
163-
}
164-
double mom = std::abs(matchedTrack.p());
165-
registry.fill(HIST("h_matchedtracks"), 1);
166-
167-
// CASE 1: skip tracks with a very low pT
168-
constexpr double kMinMom = 1e-6;
169-
if (mom < kMinMom) {
170-
continue;
171-
} // end CASE 1
172-
173-
// CASE 2:
174-
// a) If one matched track -> 100% energy subtraction
175-
// b) If more than one matched track -> 100% energy subtraction
176-
// c) If you want to do systematic studies -> perform the above two checks a) and b), and then subtract 70% energy instead of 100%
177-
178-
// Perform dEta/dPhi matching
179-
auto emcTrack = (emcTracks.sliceBy(perTrackMatchedTrack, matchedTrack.globalIndex())).iteratorAt(0);
180-
double dEta = emcTrack.etaDiff();
181-
double dPhi = emcTrack.phiDiff();
182-
183-
// Apply the eta and phi matching thresholds
184-
// dEta and dPhi cut : ensures that the matched track is within the desired eta/phi window
185-
186-
// Do pT-dependent track matching
187-
if (doMomDepMatching) {
188-
auto trackEtaMax = funcPtDepEta.Eval(mom);
189-
auto trackPhiHigh = +funcPtDepPhi.Eval(mom);
190-
auto trackPhiLow = -funcPtDepPhi.Eval(mom);
191-
192-
if ((dPhi < trackPhiHigh && dPhi > trackPhiLow) && std::fabs(dEta) < trackEtaMax) {
193-
if (nMatches == 0) {
194-
closestTrkP = mom;
195-
}
196-
totalTrkP += mom;
197-
nMatches++;
198-
}
199-
} else {
200-
// Do fixed dEta/dPhi matching (non-pT dependent)
201-
if (std::fabs(dEta) >= minDEta || std::fabs(dPhi) >= minDPhi) {
202-
continue; // Skip this track if outside the fixed cut region
203-
}
204-
205-
// If track passes fixed dEta/dPhi cuts, process it
206-
if (nMatches == 0) {
207-
closestTrkP = mom; // Closest track match
208-
}
209-
totalTrkP += mom; // Accumulate momentum
210-
nMatches++; // Count this match
211-
}
212258

213-
} // End of track loop
214-
registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), nMatches);
215-
216-
if (nMatches >= 1) {
217-
if (useM02SubtractionScheme1) {
218-
// Do M02-based correction if enabled
219-
clusterE1 = subtractM02ClusterEnergy(cluster.m02(), clusterE1, nMatches, closestTrkP, hadCorr1, useFraction1);
220-
clusterEAll1 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll1, nMatches, totalTrkP, hadCorralltrks1, useFraction1);
221-
} else {
222-
// Default energy subtraction (100% and 70%)
223-
clusterE1 = subtractClusterEnergy(clusterE1, closestTrkP, hadCorr1, nMatches, useFraction1);
224-
clusterEAll1 = subtractClusterEnergy(clusterEAll1, totalTrkP, hadCorralltrks1, nMatches, useFraction1);
225-
}
259+
// Fill the table with all four corrected energies
260+
clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2);
226261

227-
if (useM02SubtractionScheme2) {
228-
// Do M02-based correction if enabled
229-
clusterE2 = subtractM02ClusterEnergy(cluster.m02(), clusterE2, nMatches, closestTrkP, hadCorr2, useFraction2);
230-
clusterEAll2 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll2, nMatches, totalTrkP, hadCorralltrks2, useFraction2);
231-
} else {
232-
// Default energy subtraction (100% and 70%)
233-
clusterE2 = subtractClusterEnergy(clusterE2, closestTrkP, hadCorr2, nMatches, useFraction2);
234-
clusterEAll2 = subtractClusterEnergy(clusterEAll2, totalTrkP, hadCorralltrks2, nMatches, useFraction2);
235-
}
236-
}
237-
registry.fill(HIST("h_Ecluster1"), clusterE1);
238-
registry.fill(HIST("h_Ecluster2"), clusterE2);
239-
registry.fill(HIST("h_EclusterAll1"), clusterEAll1);
240-
registry.fill(HIST("h_EclusterAll2"), clusterEAll2);
241-
registry.fill(HIST("h2_ClsEvsEcluster1"), cluster.energy(), clusterE1);
242-
registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), clusterE2);
243-
registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), clusterEAll1);
244-
registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), clusterEAll2);
245-
246-
// Fill the table with all four corrected energies
247-
clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2);
248-
249-
} // End of cluster loop
262+
} // End of cluster loop
263+
} // End of cluster definition loop
250264
} // process function ends
251265
PROCESS_SWITCH(EmcalClusterHadronicCorrectionTask, processMatchedCollisions, "hadronic correction", true);
252266

0 commit comments

Comments
 (0)