Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

corrected the MCP outlier rejection for every event + identification of fake jet matches based on EMCAL fiducial cuts #7616

Merged
merged 13 commits into from
Sep 8, 2024
Merged

corrected the MCP outlier rejection for every event + identification of fake jet matches based on EMCAL fiducial cuts #7616

Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
7d6ee6f
proper jet matching+ more QA histos
Archita-Dash Aug 30, 2024
0f122bb
Merge branch 'master' of https://github.com/AliceO2Group/O2Physics
Archita-Dash Aug 31, 2024
acb2e6e
Merge branch 'master' of https://github.com/AliceO2Group/O2Physics
Archita-Dash Aug 31, 2024
7874102
fixed bugs
Archita-Dash Sep 2, 2024
ec60563
Merge branch 'master' of https://github.com/AliceO2Group/O2Physics
Archita-Dash Sep 2, 2024
c19f746
added more QA histos for matched jets
Archita-Dash Sep 2, 2024
489fbfb
fixed linter
Archita-Dash Sep 3, 2024
50f078b
Merge branch 'master' of https://github.com/AliceO2Group/O2Physics
Archita-Dash Sep 3, 2024
ef5a291
corrected the MCP outlier rejection for every event + identification …
Archita-Dash Sep 8, 2024
c692941
corrected the MCP outlier rejection for every event + identification …
Archita-Dash Sep 8, 2024
8d69ff1
corrected the MCP outlier rejection for every event + identification …
Archita-Dash Sep 8, 2024
99f056f
Please consider the following formatting changes
alibuild Sep 8, 2024
0275431
Merge pull request #10 from alibuild/alibot-cleanup-7616
Archita-Dash Sep 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 41 additions & 44 deletions PWGJE/Tasks/fulljetspectrapp.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,8 @@ struct FullJetSpectrapp {
h_collisions_unweighted->GetXaxis()->SetBinLabel(7, "JetsMCD w/o kTVXinEMC");
h_collisions_unweighted->GetXaxis()->SetBinLabel(8, "Tracks w/o kTVXinEMC");
h_collisions_unweighted->GetXaxis()->SetBinLabel(9, "JetsMCPMCDMatched w/o kTVXinEMC");
h_collisions_unweighted->GetXaxis()->SetBinLabel(10, "Fake Matched MCD Jets");
h_collisions_unweighted->GetXaxis()->SetBinLabel(11, "Fake Matched MCP Jets");
}

if (doprocessTracksWeighted) {
Expand Down Expand Up @@ -163,7 +165,7 @@ struct FullJetSpectrapp {

// Track QA histograms
if (doprocessTracks || doprocessTracksWeighted) {
registry.add("h_collisions_unweighted", "event status; event status;entries", {HistType::kTH1F, {{11, 0., 11.0}}});
registry.add("h_collisions_unweighted", "event status; event status;entries", {HistType::kTH1F, {{12, 0., 12.0}}});

registry.add("h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}});
registry.add("h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1., 1.}}});
Expand Down Expand Up @@ -230,7 +232,7 @@ struct FullJetSpectrapp {
registry.add("h_full_jet_pt_part", "jet pT;#it{p}_{T_jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}});
registry.add("h_full_jet_eta_part", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}});
registry.add("h_full_jet_phi_part", "jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}});
registry.add("h2_full_jet_NEF_part", "#it{p}_{T,jet} vs NEF at Part Level;#it{p}_{T,jet} (GeV/#it{c});NEF", {HistType::kTH2F, {{350, 0., 350.}, {100, 0.0, 1.5}}});
registry.add("h2_full_jet_NEF_part", "#it{p}_{T,jet} vs NEF at Part Level;#it{p}_{T,jet} (GeV/#it{c});NEF", {HistType::kTH2F, {{350, 0., 350.}, {105, 0., 1.05}}});

registry.add("h_Partjet_ntracks", "#it{p}_{T,constituent};#it{p}_{T_constituent} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}});
registry.add("h2_full_jet_chargedconstituents_part", "Number of charged constituents at Part Level;#it{p}_{T,jet} (GeV/#it{c});N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}});
Expand Down Expand Up @@ -305,8 +307,8 @@ struct FullJetSpectrapp {
using JetTableMCDMatchedJoined = soa::Join<aod::FullMCDetectorLevelJets, aod::FullMCDetectorLevelJetConstituents, aod::FullMCDetectorLevelJetsMatchedToFullMCParticleLevelJets>;
using JetTableMCPMatchedJoined = soa::Join<aod::FullMCParticleLevelJets, aod::FullMCParticleLevelJetConstituents, aod::FullMCParticleLevelJetsMatchedToFullMCDetectorLevelJets>;

using JetTableMCDMatchedWeightedJoined = soa::Join<aod::FullMCDetectorLevelJets, aod::FullMCDetectorLevelJetConstituents, aod::FullMCDetectorLevelJetsMatchedToFullMCParticleLevelJets, aod::FullMCDetectorLevelJetEventWeights>;
using JetTableMCPMatchedWeightedJoined = soa::Join<aod::FullMCParticleLevelJets, aod::FullMCParticleLevelJetConstituents, aod::FullMCParticleLevelJetsMatchedToFullMCDetectorLevelJets, aod::FullMCParticleLevelJetEventWeights>;
using JetTableMCDMatchedWeightedJoined = soa::Join<aod::FullMCDetectorLevelJets, aod::FullMCDetectorLevelJetConstituents, aod::FullMCDetectorLevelJetsMatchedToFullMCParticleLevelJets>;
using JetTableMCPMatchedWeightedJoined = soa::Join<aod::FullMCParticleLevelJets, aod::FullMCParticleLevelJetConstituents, aod::FullMCParticleLevelJetsMatchedToFullMCDetectorLevelJets>;

// Applying some cuts(filters) on collisions, tracks, clusters

Expand All @@ -315,6 +317,7 @@ struct FullJetSpectrapp {
Filter trackCuts = (aod::jtrack::pt >= trackpTMin && aod::jtrack::pt < trackpTMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax && aod::jtrack::phi >= trackPhiMin && aod::jtrack::phi <= trackPhiMax);
aod::EMCALClusterDefinition clusterDefinition = aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionS.value);
Filter clusterFilter = (aod::jcluster::definition == static_cast<int>(clusterDefinition) && aod::jcluster::eta > clusterEtaMin && aod::jcluster::eta < clusterEtaMax && aod::jcluster::phi >= clusterPhiMin && aod::jcluster::phi <= clusterPhiMax && aod::jcluster::energy >= clusterEnergyMin && aod::jcluster::time > clusterTimeMin && aod::jcluster::time < clusterTimeMax && (clusterRejectExotics && aod::jcluster::isExotic != true));
Preslice<JetTableMCPMatchedWeightedJoined> JetMCPPerMcCollision = aod::jet::mcCollisionId;

template <typename T, typename U>
bool isAcceptedJet(U const& jet)
Expand Down Expand Up @@ -344,7 +347,6 @@ struct FullJetSpectrapp {
void fillJetHistograms(T const& jet, float weight = 1.0)
{
float neutralEnergy = 0.0;
// std::cout << "jet r is " << jet.r() << " its rounded value is " << round(selectedJetsRadius * 100.0f) << std::endl;
if (jet.r() == round(selectedJetsRadius * 100.0f)) {
registry.fill(HIST("h_full_jet_pt"), jet.pt(), weight);
registry.fill(HIST("h_full_jet_eta"), jet.eta(), weight);
Expand Down Expand Up @@ -628,7 +630,7 @@ struct FullJetSpectrapp {
registry.fill(HIST("h_collisions_weighted"), 5.0); // JetsMCDWeighted w/o kTVXinEMC
for (auto const& jet : jets) {
if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedJet<JetTracks>(jet)) {
fillRejectedJetHistograms(jet, 1.0);
fillRejectedJetHistograms(jet, jet.eventWeight());
}
}
return;
Expand All @@ -644,10 +646,8 @@ struct FullJetSpectrapp {
if (!isAcceptedJet<JetTracks>(jet)) {
continue;
}
// std::cout << "jet pT" << jet.pt() << "jet eta"<< jet.eta()<< "jet phi"<< jet.phi()<< std::endl;
// std::cout << "Event weight: " << jet.eventWeight() << std::endl;

fillJetHistograms(jet, jet.eventWeight());
// std::cout << "jet pT" << jet.pt() << "jet eta"<< jet.eta()<< "jet phi"<< jet.phi()<< std::endl;
}
}
PROCESS_SWITCH(FullJetSpectrapp, processJetsMCDWeighted, "Full Jets at Detector Level on weighted events", false);
Expand All @@ -672,8 +672,7 @@ struct FullJetSpectrapp {
if (!isAcceptedJet<JetParticles>(jet)) {
return;
}
// std::cout << "jet pT" << jet.pt() << "jet eta"<< jet.eta()<< "jet phi"<< jet.phi()<< std::endl;
// std::cout << "Event weight: " << jet.eventWeight() << std::endl;

fillMCPHistograms(jet, jet.eventWeight());
}
PROCESS_SWITCH(FullJetSpectrapp, processJetsMCPWeighted, "Full Jets at Particle Level on weighted events", false);
Expand Down Expand Up @@ -718,6 +717,8 @@ struct FullJetSpectrapp {
{
registry.fill(HIST("h_collisions_unweighted"), 1.0); // total events
bool eventAccepted = false;
int fakemcdjet = 0;
int fakemcpjet = 0;

if (fabs(collision.posZ()) > VertexZCut) { // making double sure this condition is satisfied
return;
Expand Down Expand Up @@ -746,27 +747,28 @@ struct FullJetSpectrapp {
//**end of event selection**

for (const auto& mcdjet : mcdjets) {
if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
continue;
}
if (!isAcceptedJet<JetTracks>(mcdjet)) {
// Check if MCD jet is within the EMCAL fiducial region
if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax || mcdjet.eta() < jetEtaMin || mcdjet.eta() > jetEtaMax) {
fakemcdjet++;
registry.fill(HIST("h_collisions_unweighted"), 10.0); // Fake Matched MCD Jets
registry.fill(HIST("h2_full_fakemcdjets"), mcdjet.pt(), fakemcdjet, 1.0);
continue;
}
if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax) {
if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
continue;
}
if (!mcdjet.has_matchedJetGeo()) {
if (!isAcceptedJet<JetTracks>(mcdjet)) {
continue;
}
for (auto& mcpjet : mcdjet.template matchedJetGeo_as<JetTableMCPMatchedJoined>()) {

if (!mcpjet.has_matchedJetGeo()) {
continue;
}
// apply emcal fiducial cuts to the matched particle level jets
if (mcpjet.eta() > jetEtaMax || mcpjet.phi() > jetPhiMax || !jetfindingutilities::isInEtaAcceptance(mcpjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
if (mcpjet.eta() > jetEtaMax || mcpjet.eta() < jetEtaMin || mcpjet.phi() > jetPhiMax || mcpjet.phi() < jetPhiMin) {
fakemcpjet++;
registry.fill(HIST("h_collisions_unweighted"), 11.0); // Fake Matched MCP Jets
registry.fill(HIST("h2_full_fakemcpjets"), mcpjet.pt(), fakemcpjet, 1.0);
continue;
}
// Fill MCD jet histograms if a valid MCP jet match was found within the EMCAL region
registry.fill(HIST("h_full_matchedmcpjet_eta"), mcpjet.eta(), 1.0);
registry.fill(HIST("h_full_matchedmcpjet_phi"), mcpjet.phi(), 1.0);
fillMatchedHistograms<typename JetTableMCDMatchedJoined::iterator, JetTableMCPMatchedJoined>(mcdjet);
Expand All @@ -777,15 +779,21 @@ struct FullJetSpectrapp {
}
PROCESS_SWITCH(FullJetSpectrapp, processJetsMCPMCDMatched, "Full Jet finder MCP matched to MCD", false);

void processJetsMCPMCDMatchedWeighted(soa::Filtered<soa::Join<EMCCollisions, aod::JMcCollisionLbs>>::iterator const& collision, JetTableMCDMatchedWeightedJoined const& mcdjets, JetTableMCPMatchedWeightedJoined const&, aod::JMcCollisions const&, JetTracks const&, JetClusters const&, JetParticles const&)
void processJetsMCPMCDMatchedWeighted(soa::Filtered<soa::Join<EMCCollisions, aod::JMcCollisionLbs>>::iterator const& collision, JetTableMCDMatchedWeightedJoined const& mcdjets, JetTableMCPMatchedWeightedJoined const& mcpjets, aod::JMcCollisions const&, JetTracks const&, JetClusters const&, JetParticles const&)
{
float eventWeight = collision.mcCollision().weight();
registry.fill(HIST("h_collisions_weighted"), 1.0, eventWeight); // total events
bool eventAccepted = false;
int fakemcdjet = 0;
int fakemcpjet = 0;
float weight = 1.0;
float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent));
const auto mcpJetsPerMcCollision = mcpjets.sliceBy(JetMCPPerMcCollision, collision.mcCollisionId());

for (auto mcpjet : mcpJetsPerMcCollision) {
if (mcpjet.pt() > pTHatMaxMCP * pTHat) { // outlier rejection for MCP
return;
}
}

if (doEMCALEventWorkaround) {
if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content
Expand All @@ -807,44 +815,33 @@ struct FullJetSpectrapp {
}

for (const auto& mcdjet : mcdjets) {
if (mcdjet.pt() > pTHatMaxMCD * pTHat) {
eventAccepted = false; ////reject the whole event for outlier jets
// Check if MCD jet is within the EMCAL fiducial region
if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax || mcdjet.eta() < jetEtaMin || mcdjet.eta() > jetEtaMax) {
fakemcdjet++;
registry.fill(HIST("h_collisions_weighted"), 8.0); // Fake Matched Weighted MCD Jets
registry.fill(HIST("h2_full_fakemcdjets"), mcdjet.pt(), fakemcdjet, eventWeight);
break;
continue;
}
if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
continue;
}
if (!isAcceptedJet<JetTracks>(mcdjet)) {
continue;
}
if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax) {
continue;
}
if (!mcdjet.has_matchedJetGeo()) {
continue;
}
for (auto& mcpjet : mcdjet.template matchedJetGeo_as<JetTableMCPMatchedWeightedJoined>()) {
if (mcpjet.pt() > pTHatMaxMCP * pTHat) {
eventAccepted = false; // reject the whole event for outlier jets
// apply emcal fiducial cuts to the matched particle level jets
if (mcpjet.eta() > jetEtaMax || mcpjet.eta() < jetEtaMin || mcpjet.phi() > jetPhiMax || mcpjet.phi() < jetPhiMin) {
fakemcpjet++;
registry.fill(HIST("h_collisions_weighted"), 9.0); // Fake Matched Weighted MCP Jets
registry.fill(HIST("h2_full_fakemcpjets"), mcpjet.pt(), fakemcpjet, eventWeight);
break;
}
if (!mcpjet.has_matchedJetGeo()) {
continue;
}
// apply emcal fiducial cuts to the matched particle level jets
if (mcpjet.eta() > jetEtaMax || mcpjet.phi() > jetPhiMax || !jetfindingutilities::isInEtaAcceptance(mcpjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
continue;
}
// If both MCD-MCP matched jet pairs are within the EMCAL fiducial region, fill these histos
registry.fill(HIST("h_full_matchedmcpjet_eta"), mcpjet.eta(), eventWeight);
registry.fill(HIST("h_full_matchedmcpjet_phi"), mcpjet.phi(), eventWeight);
fillMatchedHistograms<typename JetTableMCDMatchedWeightedJoined::iterator, JetTableMCPMatchedWeightedJoined>(mcdjet, mcdjet.eventWeight());
fillMatchedHistograms<typename JetTableMCDMatchedWeightedJoined::iterator, JetTableMCPMatchedWeightedJoined>(mcdjet, eventWeight);
} // mcpjet
// Fill MCD jet histograms if a valid MCP jet match was found within the EMCAL region
registry.fill(HIST("h_full_matchedmcdjet_eta"), mcdjet.eta(), eventWeight);
registry.fill(HIST("h_full_matchedmcdjet_phi"), mcdjet.phi(), eventWeight);
} // mcdjet
Expand Down
Loading