Skip to content

Commit 004de26

Browse files
hanseoparkhanseopark
authored
[PWGJE] Add validation function of heavy-flavour definition and fix response matrix (#10385)
1 parent cae8811 commit 004de26

File tree

1 file changed

+94
-2
lines changed

1 file changed

+94
-2
lines changed

‎PWGJE/Tasks/jetTaggerHFQA.cxx‎

Lines changed: 94 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@ struct JetTaggerHFQA {
5959
Configurable<float> trackPtMax{"trackPtMax", 100.0, "maximum pT acceptance for tracks"};
6060
Configurable<float> trackDcaXYMax{"trackDcaXYMax", 1, "minimum DCA xy acceptance for tracks [cm]"};
6161
Configurable<float> trackDcaZMax{"trackDcaZMax", 2, "minimum DCA z acceptance for tracks [cm]"};
62+
Configurable<float> maxDeltaR{"maxDeltaR", 0.25, "maximum distance of jet axis from flavour initiating parton"};
6263
Configurable<float> jetEtaMin{"jetEtaMin", -99.0, "minimum jet pseudorapidity"};
6364
Configurable<float> jetEtaMax{"jetEtaMax", 99.0, "maximum jet pseudorapidity"};
6465
Configurable<float> prongChi2PCAMin{"prongChi2PCAMin", 1, "minimum Chi2 PCA of decay length of prongs"};
@@ -155,6 +156,18 @@ struct JetTaggerHFQA {
155156
registry.add("h_impact_parameter_xyz", "", {HistType::kTH1F, {{axisImpactParameterXYZ}}});
156157
registry.add("h_impact_parameter_xyz_significance", "", {HistType::kTH1F, {{axisImpactParameterXYZSignificance}}});
157158
}
159+
if (doprocessValFlavourDefMCD) {
160+
registry.add("h2_flavour_dist_quark_flavour_dist_hadron", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
161+
registry.add("h2_flavour_const_quark_flavour_const_hadron", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
162+
registry.add("h2_flavour_const_hadron_flavour_dist_hadron", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
163+
registry.add("h2_flavour_const_quark_flavour_dist_quark", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
164+
}
165+
if (doprocessValFlavourDefMCP) {
166+
registry.add("h2_part_flavour_dist_quark_part_flavour_dist_hadron", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
167+
registry.add("h2_part_flavour_const_quark_part_flavour_const_hadron", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
168+
registry.add("h2_part_flavour_const_hadron_part_flavour_dist_hadron", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
169+
registry.add("h2_part_flavour_const_quark_part_flavour_dist_quark", "", {HistType::kTH2F, {{axisJetFlavour}, {axisJetFlavour}}});
170+
}
158171
if (doprocessIPsData) {
159172
registry.add("h_jet_pt", "", {HistType::kTH1F, {{axisJetPt}}});
160173
registry.add("h_jet_eta", "", {HistType::kTH1F, {{axisEta}}});
@@ -443,6 +456,47 @@ struct JetTaggerHFQA {
443456
return true;
444457
}
445458

459+
template <typename U, typename T, typename V, typename W, typename X>
460+
void fillValidationFlavourDefMCD(T const& mcdjet, V const& tracks, W const& particles, X const& particlesPerColl, float eventWeight = 1.0)
461+
{
462+
float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
463+
if (mcdjet.pt() > pTHatMaxMCD * pTHat) {
464+
return;
465+
}
466+
typename V::iterator hftrack;
467+
int jetflavourConstQuark = jettaggingutilities::mcdJetFromHFShower(mcdjet, tracks, particles, maxDeltaR, true);
468+
int jetflavourConstHadron = jettaggingutilities::mcdJetFromHFShower(mcdjet, tracks, particles, maxDeltaR, false);
469+
int jetflavourDistQuark = -1;
470+
int jetflavourDistHadron = -1;
471+
for (auto const& mcpjet : mcdjet.template matchedJetGeo_as<U>()) {
472+
jetflavourDistQuark = jettaggingutilities::getJetFlavor(mcpjet, particlesPerColl);
473+
jetflavourDistHadron = jettaggingutilities::getJetFlavorHadron(mcpjet, particlesPerColl);
474+
}
475+
if (jetflavourDistQuark < 0 || jetflavourDistHadron < 0)
476+
return;
477+
registry.fill(HIST("h2_flavour_dist_quark_flavour_dist_hadron"), jetflavourDistQuark, jetflavourDistHadron, eventWeight);
478+
registry.fill(HIST("h2_flavour_const_quark_flavour_const_hadron"), jetflavourConstQuark, jetflavourConstHadron, eventWeight);
479+
registry.fill(HIST("h2_flavour_const_hadron_flavour_dist_hadron"), jetflavourConstHadron, jetflavourDistHadron, eventWeight);
480+
registry.fill(HIST("h2_flavour_const_quark_flavour_dist_quark"), jetflavourConstQuark, jetflavourDistQuark, eventWeight);
481+
}
482+
483+
template <typename T, typename U, typename V>
484+
void fillValidationFlavourDefMCP(T const& mcpjet, U const& particles, V const& particlesPerColl, float eventWeight = 1.0)
485+
{
486+
float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
487+
if (mcpjet.pt() > pTHatMaxMCP * pTHat) {
488+
return;
489+
}
490+
int jetflavourConstQuark = jettaggingutilities::mcpJetFromHFShower(mcpjet, particles, maxDeltaR, true);
491+
int jetflavourConstHadron = jettaggingutilities::mcpJetFromHFShower(mcpjet, particles, maxDeltaR, false);
492+
int jetflavourDistQuark = jettaggingutilities::getJetFlavor(mcpjet, particlesPerColl);
493+
int jetflavourDistHadron = jettaggingutilities::getJetFlavorHadron(mcpjet, particlesPerColl);
494+
registry.fill(HIST("h2_part_flavour_dist_quark_part_flavour_dist_hadron"), jetflavourDistQuark, jetflavourDistHadron, eventWeight);
495+
registry.fill(HIST("h2_part_flavour_const_quark_part_flavour_const_hadron"), jetflavourConstQuark, jetflavourConstHadron, eventWeight);
496+
registry.fill(HIST("h2_part_flavour_const_hadron_part_flavour_dist_hadron"), jetflavourConstHadron, jetflavourDistHadron, eventWeight);
497+
registry.fill(HIST("h2_part_flavour_const_quark_part_flavour_dist_quark"), jetflavourConstQuark, jetflavourDistQuark, eventWeight);
498+
}
499+
446500
template <typename T, typename U>
447501
void fillHistogramIPsData(T const& jet, U const& /*tracks*/)
448502
{
@@ -1013,6 +1067,44 @@ struct JetTaggerHFQA {
10131067
}
10141068
PROCESS_SWITCH(JetTaggerHFQA, processTracksDca, "Fill inclusive tracks' imformation for data", false);
10151069

1070+
void processValFlavourDefMCD(soa::Filtered<soa::Join<aod::JCollisions, aod::JCollisionPIs, aod::JMcCollisionLbs>>::iterator const& collision, soa::Join<JetTableMCD, TagTableMCD, JetTableMCDMCP, weightMCD> const& mcdjets, soa::Join<JetTableMCP, JetTableMCPMCD> const& /*mcpjets*/, JetTagTracksMCD const& tracks, aod::JetParticles const& particles)
1071+
{
1072+
if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
1073+
return;
1074+
}
1075+
for (auto const& mcdjet : mcdjets) {
1076+
auto const particlesPerColl = particles.sliceBy(particlesPerCollision, collision.mcCollisionId());
1077+
if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
1078+
continue;
1079+
}
1080+
if (!isAcceptedJet<aod::JetTracks>(mcdjet)) {
1081+
continue;
1082+
}
1083+
fillValidationFlavourDefMCD<soa::Join<JetTableMCP, JetTableMCPMCD>>(mcdjet, tracks, particles, particlesPerColl, mcdjet.eventWeight());
1084+
}
1085+
}
1086+
PROCESS_SWITCH(JetTaggerHFQA, processValFlavourDefMCD, "to check the validation of jet-flavour definition when compared to distance for mcd jets", false);
1087+
1088+
void processValFlavourDefMCP(soa::Join<JetTableMCP, weightMCP> const& mcpjets, aod::JetParticles const& particles, aod::JetMcCollisions const&)
1089+
{
1090+
for (auto const& mcpjet : mcpjets) {
1091+
auto const particlesPerColl = particles.sliceBy(particlesPerCollision, mcpjet.globalIndex());
1092+
if (!jetfindingutilities::isInEtaAcceptance(mcpjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
1093+
continue;
1094+
}
1095+
if (!isAcceptedJet<aod::JetParticles>(mcpjet)) {
1096+
continue;
1097+
}
1098+
int eventWeight = mcpjet.eventWeight();
1099+
float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
1100+
if (mcpjet.pt() > pTHatMaxMCD * pTHat) {
1101+
return;
1102+
}
1103+
fillValidationFlavourDefMCP(mcpjet, particles, particlesPerColl);
1104+
}
1105+
}
1106+
PROCESS_SWITCH(JetTaggerHFQA, processValFlavourDefMCP, "to check the validation of jet-flavour definition when compared to distance for mcp jets", false);
1107+
10161108
void processIPsData(soa::Filtered<aod::JetCollisions>::iterator const& collision, soa::Join<JetTableData, TagTableData> const& jets, JetTagTracksData const& tracks)
10171109
{
10181110
if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
@@ -1119,7 +1211,7 @@ struct JetTaggerHFQA {
11191211
if (!mcdjet.has_matchedJetGeo())
11201212
continue;
11211213
for (auto const& mcpjet : mcdjet.template matchedJetGeo_as<soa::Join<JetTableMCP, JetTableMCPMCD>>()) {
1122-
registry.fill(HIST("h3_jet_pt_jet_pt_part_matchedgeo_flavour"), mcpjet.pt(), mcdjet.pt(), mcdjet.origin());
1214+
registry.fill(HIST("h3_jet_pt_jet_pt_part_matchedgeo_flavour"), mcdjet.pt(), mcpjet.pt(), mcdjet.origin());
11231215
}
11241216
if (!doprocessIPsMCD)
11251217
fillHistogramIPsMCD(mcdjet, tracks);
@@ -1150,7 +1242,7 @@ struct JetTaggerHFQA {
11501242
if (mcpjet.pt() > pTHatMaxMCP * pTHat) {
11511243
continue;
11521244
}
1153-
registry.fill(HIST("h3_jet_pt_jet_pt_part_matchedgeo_flavour"), mcpjet.pt(), mcdjet.pt(), mcdjet.origin(), mcdjet.eventWeight());
1245+
registry.fill(HIST("h3_jet_pt_jet_pt_part_matchedgeo_flavour"), mcdjet.pt(), mcpjet.pt(), mcdjet.origin(), mcdjet.eventWeight());
11541246
}
11551247
if (!doprocessIPsMCDWeighted)
11561248
fillHistogramIPsMCD(mcdjet, tracks, mcdjet.eventWeight());

0 commit comments

Comments
 (0)