Skip to content

[PWGJE] include functionality for outlier rejection task #14127

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

Open
jaimenorman wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[PWGJE] include functionality for outlier rejection task #14127

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
36 changes: 29 additions & 7 deletions PWGJE/Tasks/jetHadronRecoil.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,7 @@
Configurable<bool> skipMBGapEvents{"skipMBGapEvents", false, "flag to choose to reject min. bias gap events; jet-level rejection applied at the jet finder level, here rejection is applied for collision and track process functions"};
Configurable<bool> outlierRejectEvent{"outlierRejectEvent", true, "where outliers are found, reject event (true) or just reject the single track/jet (false)"};
Configurable<bool> doSumw{"doSumw", false, "enable sumw2 for weighted histograms"};
Configurable<bool> rejectOutlierCollisions{"rejectOutlierCollisions", false, "flag to reject outlier collisions for weighed MC"};

TRandom3* rand = new TRandom3(0);

Expand All @@ -102,7 +103,7 @@
std::vector<double> dRBinning = {0.0, 1.0e-9, 0.003, 0.006, 0.009, 0.012, 0.015, 0.018, 0.021, 0.024,
0.027, 0.03, 0.033, 0.036, 0.039, 0.042, 0.045, 0.048, 0.051, 0.054,
0.057, 0.06, 0.063, 0.066, 0.069, 0.072, 0.075, 0.078, 0.081, 0.084,
0.087, 0.09, 0.093, 0.096, 0.099, 0.102, 0.105, 0.108, 0.111, 0.114,

Check failure on line 106 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-mass] 

Avoid hard-coded particle masses. Use o2::constants::physics::Mass... instead.
0.117, 0.12, 0.123, 0.126, 0.129, 0.132, 0.135, 0.138, 0.141, 0.144,
0.147, 0.15, 0.153, 0.156, 0.159, 0.162, 0.165, 0.168, 0.171, 0.174,
0.177, 0.18, 0.183, 0.186, 0.189, 0.192, 0.195, 0.198, 0.201, 0.204,
Expand Down Expand Up @@ -281,7 +282,7 @@
registry.fill(HIST("hNtrig"), 0.5, weight);
registry.fill(HIST("hRefEventTriggers"), nTT, weight);
registry.fill(HIST("hRhoReference"), rhoReference, weight);
for (double shift = 0.0; shift <= 2.0; shift += 0.1) {

Check failure on line 285 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hRhoReferenceShift"), rho + shift, shift, weight);
}
registry.fill(HIST("hReferenceTriggersPtHard"), ptTT / pTHat, weight);
Expand Down Expand Up @@ -316,31 +317,31 @@
float dphi = RecoDecay::constrainAngle(jet.phi() - phiTT);
double dR = getWTAaxisDifference(jet, tracks);
if (isSigCol) {
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {

Check failure on line 320 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hDeltaRpTSignal"), jet.pt() - (rho * jet.area()), dR, weight);
registry.fill(HIST("hDeltaRSignal"), dR, weight);
}
registry.fill(HIST("hDeltaRpTDPhiSignal"), jet.pt() - (rho * jet.area()), dphi, dR, weight);
registry.fill(HIST("hSignalPtDPhi"), dphi, jet.pt() - (rho * jet.area()), weight);
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {

Check failure on line 326 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hSignalPt"), jet.pt() - (rho * jet.area()), weight);
registry.fill(HIST("hSignalPtHard"), jet.pt() - (rho * jet.area()), ptTT / pTHat, weight);
}
}
if (!isSigCol) {
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {

Check failure on line 332 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hDeltaRpTReference"), jet.pt() - (rhoReference * jet.area()), dR, weight);
registry.fill(HIST("hDeltaRReference"), dR, weight);
}
registry.fill(HIST("hDeltaRpTDPhiReference"), jet.pt() - (rhoReference * jet.area()), dphi, dR, weight);
for (double shift = 0.0; shift <= 2.0; shift += 0.1) {

Check failure on line 337 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hDeltaRpTDPhiReferenceShifts"), jet.pt() - ((rho + shift) * jet.area()), dphi, dR, shift, weight);
}
registry.fill(HIST("hReferencePtDPhi"), dphi, jet.pt() - (rhoReference * jet.area()), weight);
for (double shift = 0.0; shift <= 2.0; shift += 0.1) {

Check failure on line 341 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hReferencePtDPhiShifts"), dphi, jet.pt() - ((rho + shift) * jet.area()), shift, weight);
}
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {

Check failure on line 344 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hReferencePt"), jet.pt() - (rhoReference * jet.area()), weight);
registry.fill(HIST("hReferencePtHard"), jet.pt() - (rhoReference * jet.area()), ptTT / pTHat, weight);
}
Expand Down Expand Up @@ -423,7 +424,7 @@
registry.fill(HIST("hNtrig"), 0.5, weight);
registry.fill(HIST("hRefEventTriggers"), nTT, weight);
registry.fill(HIST("hRhoReference"), rhoReference, weight);
for (double shift = 0.0; shift <= 2.0; shift += 0.1) {

Check failure on line 427 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
registry.fill(HIST("hRhoReferenceShift"), rho + shift, shift, weight);
}
registry.fill(HIST("hReferenceTriggersPtHard"), ptTT / pTHat, weight);
Expand Down Expand Up @@ -811,7 +812,7 @@
}
PROCESS_SWITCH(JetHadronRecoil, processMCDWithRhoSubtraction, "process MC detector level with rho subtraction", false);

void processMCDWeighted(soa::Filtered<soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>>::iterator const& collision,
void processMCDWeighted(soa::Filtered<soa::Join<aod::JetCollisions, aod::JMcCollisionLbs, aod::JCollisionOutliers>>::iterator const& collision,
aod::JMcCollisions const&,
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents>> const& jets,
soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras, aod::JMcTrackLbs>> const& tracks)
Expand All @@ -831,12 +832,15 @@
if (collision.mcCollision().ptHard() < pTHatMinEvent) {
return;
}
if (rejectOutlierCollisions && collision.isOutlier()) {
return;
}
registry.fill(HIST("hZvtxSelected"), collision.posZ(), collision.mcCollision().weight());
fillHistogramsMCD(jets, tracks, collision.mcCollision().weight(), 0.0, collision.mcCollision().ptHard(), collision.mcCollisionId());
}
PROCESS_SWITCH(JetHadronRecoil, processMCDWeighted, "process MC detector level with event weights", false);

void processMCDWeightedWithRhoSubtraction(soa::Filtered<soa::Join<aod::JetCollisions, aod::JMcCollisionLbs, aod::BkgChargedRhos>>::iterator const& collision,
void processMCDWeightedWithRhoSubtraction(soa::Filtered<soa::Join<aod::JetCollisions, aod::JMcCollisionLbs, aod::BkgChargedRhos, aod::JCollisionOutliers>>::iterator const& collision,
aod::JMcCollisions const&,
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents>> const& jets,
soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras, aod::JMcTrackLbs>> const& tracks)
Expand All @@ -856,6 +860,9 @@
if (collision.mcCollision().ptHard() < pTHatMinEvent) {
return;
}
if (rejectOutlierCollisions && collision.isOutlier()) {
return;
}
registry.fill(HIST("hZvtxSelected"), collision.posZ(), collision.mcCollision().weight());
fillHistogramsMCD(jets, tracks, collision.mcCollision().weight(), collision.rho(), collision.mcCollision().ptHard(), collision.mcCollisionId());
}
Expand Down Expand Up @@ -892,7 +899,7 @@
PROCESS_SWITCH(JetHadronRecoil, processMCP, "process MC particle level", false);

void processMCPWeighted(aod::JetMcCollision const& mccollision,
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
soa::SmallGroups<soa::Join<aod::JetCollisionsMCD, aod::JCollisionOutliers>> const& collisions,
soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents>> const& jets,
soa::Filtered<aod::JetParticles> const& particles)
{
Expand All @@ -915,6 +922,9 @@
if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) {
return;
}
if (rejectOutlierCollisions && collision.isOutlier()) {
return;
}
}
registry.fill(HIST("hZvtxSelected"), mccollision.posZ(), mccollision.weight());
fillMCPHistograms(jets, particles, mccollision.weight(), mccollision.ptHard());
Expand Down Expand Up @@ -969,7 +979,7 @@
}
PROCESS_SWITCH(JetHadronRecoil, processJetsMCPMCDMatchedWithRhoSubtraction, "process MC matched (inc jets) with rho subtraction", false);

void processJetsMCPMCDMatchedWeighted(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JMcCollisionLbs>>::iterator const& collision,
void processJetsMCPMCDMatchedWeighted(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JMcCollisionLbs, aod::JCollisionOutliers>>::iterator const& collision,
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const& mcdjets,
aod::JetTracks const& tracks,
aod::JetParticles const& particles,
Expand All @@ -988,12 +998,15 @@
if (collision.mcCollision().ptHard() < pTHatMinEvent) {
return;
}
if (rejectOutlierCollisions && collision.isOutlier()) {
return;
}
registry.fill(HIST("hZvtxSelected"), collision.posZ());
fillMatchedHistograms(mcdjets, mcpjets, tracks, particles, collision.mcCollision().weight(), 0.0, collision.mcCollision().ptHard());
}
PROCESS_SWITCH(JetHadronRecoil, processJetsMCPMCDMatchedWeighted, "process MC matched with event weights (inc jets)", false);

void processJetsMCPMCDMatchedWeightedWithRhoSubtraction(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JMcCollisionLbs, aod::BkgChargedRhos>>::iterator const& collision,
void processJetsMCPMCDMatchedWeightedWithRhoSubtraction(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JMcCollisionLbs, aod::BkgChargedRhos, aod::JCollisionOutliers>>::iterator const& collision,
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const& mcdjets,
aod::JetTracks const& tracks,
aod::JetParticles const& particles,
Expand All @@ -1012,6 +1025,9 @@
if (collision.mcCollision().ptHard() < pTHatMinEvent) {
return;
}
if (rejectOutlierCollisions && collision.isOutlier()) {
return;
}
registry.fill(HIST("hZvtxSelected"), collision.posZ());
fillMatchedHistograms(mcdjets, mcpjets, tracks, particles, collision.mcCollision().weight(), collision.rho(), collision.mcCollision().ptHard());
}
Expand Down Expand Up @@ -1050,7 +1066,7 @@
PROCESS_SWITCH(JetHadronRecoil, processRecoilJetsMCPMCDMatched, "process MC matched (recoil jets)", false);

void processRecoilJetsMCPMCDMatchedWeighted(aod::JetMcCollisions::iterator const& mccollision,
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
soa::SmallGroups<soa::Join<aod::JetCollisionsMCD, aod::JCollisionOutliers>> const& collisions,
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const& mcdjets,
soa::Filtered<aod::JetTracksMCD> const& tracks,
soa::Filtered<aod::JetParticles> const& particles,
Expand All @@ -1075,14 +1091,17 @@
if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) {
return;
}
if (rejectOutlierCollisions && collision.isOutlier()) {
return;
}
}
registry.fill(HIST("hZvtxSelected"), mccollision.posZ(), mccollision.weight());
fillRecoilJetMatchedHistograms(mcdjets, mcpjets, tracks, particles, mccollision.weight(), 0.0, mccollision.ptHard());
}
PROCESS_SWITCH(JetHadronRecoil, processRecoilJetsMCPMCDMatchedWeighted, "process MC matched with event weights (recoil jets)", false);

void processRecoilJetsMCPMCDMatchedWeightedWithRhoSubtraction(soa::Join<aod::JetMcCollisions, aod::BkgChargedRhos>::iterator const& mccollision,
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
soa::SmallGroups<soa::Join<aod::JetCollisionsMCD, aod::JCollisionOutliers>> const& collisions,
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const& mcdjets,
soa::Filtered<aod::JetTracksMCD> const& tracks,
soa::Filtered<aod::JetParticles> const& particles,
Expand All @@ -1107,6 +1126,9 @@
if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) {
return;
}
if (rejectOutlierCollisions && collision.isOutlier()) {
return;
}
}
registry.fill(HIST("hZvtxSelected"), mccollision.posZ(), mccollision.weight());
fillRecoilJetMatchedHistograms(mcdjets, mcpjets, tracks, particles, mccollision.weight(), mccollision.rho(), mccollision.ptHard());
Expand All @@ -1122,7 +1144,7 @@
jetConstituents.clear();
for (auto& jetConstituent : jet.template tracks_as<X>()) {
fastjetutilities::fillTracks(jetConstituent, jetConstituents, jetConstituent.globalIndex());
}

Check failure on line 1147 in PWGJE/Tasks/jetHadronRecoil.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
jetReclustered.clear();
fastjet::ClusterSequenceArea clusterSeq(jetReclusterer.findJets(jetConstituents, jetReclustered));
jetReclustered = sorted_by_pt(jetReclustered);
Expand Down
Loading