From ec98923ab7f5f0d8ae3bd822f23fa51ae283d4ab Mon Sep 17 00:00:00 2001 From: nzardosh Date: Sat, 21 Jun 2025 19:44:33 +0200 Subject: [PATCH] Adding Outlier rejection task for MC --- PWGJE/DataModel/JetReducedData.h | 12 +- PWGJE/TableProducer/CMakeLists.txt | 5 + PWGJE/TableProducer/derivedDataProducer.cxx | 4 +- PWGJE/TableProducer/derivedDataWriter.cxx | 2 +- PWGJE/TableProducer/mcOutlierRejector.cxx | 132 ++++++++++++++++++++ 5 files changed, 151 insertions(+), 4 deletions(-) create mode 100644 PWGJE/TableProducer/mcOutlierRejector.cxx diff --git a/PWGJE/DataModel/JetReducedData.h b/PWGJE/DataModel/JetReducedData.h index 4d5e401fe7d..08fba650bfe 100644 --- a/PWGJE/DataModel/JetReducedData.h +++ b/PWGJE/DataModel/JetReducedData.h @@ -99,6 +99,7 @@ DECLARE_SOA_COLUMN(ReadCountsWithTVXAndZVertexAndSel7KINT7, readCountsWithTVXAnd DECLARE_SOA_COLUMN(ReadCountsWithCustom, readCountsWithCustom, std::vector); DECLARE_SOA_COLUMN(IsAmbiguous, isAmbiguous, bool); DECLARE_SOA_COLUMN(IsEMCALReadout, isEmcalReadout, bool); +DECLARE_SOA_COLUMN(IsOutlier, isOutlier, bool); } // namespace jcollision DECLARE_SOA_TABLE_STAGED(JCollisions, "JCOLLISION", @@ -122,6 +123,9 @@ DECLARE_SOA_TABLE_STAGED(JCollisionMcInfos, "JCOLLISIONMCINFO", jcollision::Weight, jcollision::SubGeneratorId); +DECLARE_SOA_TABLE_STAGED(JCollisionOutliers, "JCOLLISIONOUTLR", + jcollision::IsOutlier); + DECLARE_SOA_TABLE_STAGED(JEMCCollisionLbs, "JEMCCOLLISIONLB", jcollision::IsAmbiguous, jcollision::IsEMCALReadout); @@ -170,6 +174,8 @@ DECLARE_SOA_COLUMN(Accepted, accepted, uint64_t); DECLARE_SOA_COLUMN(Attempted, attempted, uint64_t); DECLARE_SOA_COLUMN(XsectGen, xsectGen, float); DECLARE_SOA_COLUMN(XsectErr, xsectErr, float); +DECLARE_SOA_COLUMN(PtHard, ptHard, float); +DECLARE_SOA_COLUMN(IsOutlier, isOutlier, bool); } // namespace jmccollision DECLARE_SOA_TABLE_STAGED(JMcCollisions, "JMCCOLLISION", o2::soa::Index<>, @@ -181,7 +187,8 @@ DECLARE_SOA_TABLE_STAGED(JMcCollisions, "JMCCOLLISION", jmccollision::Accepted, jmccollision::Attempted, jmccollision::XsectGen, - jmccollision::XsectErr); + jmccollision::XsectErr, + jmccollision::PtHard); using JMcCollision = JMcCollisions::iterator; using StoredJMcCollision = StoredJMcCollisions::iterator; @@ -197,6 +204,9 @@ DECLARE_SOA_INDEX_COLUMN(JMcCollision, mcCollision); DECLARE_SOA_TABLE_STAGED(JMcCollisionLbs, "JMCCOLLISIONLB", jmccollisionlb::JMcCollisionId); +DECLARE_SOA_TABLE_STAGED(JMcCollisionOutliers, "JMCCOLLISIONOUTLR", + jmccollision::IsOutlier); + namespace jtrack { DECLARE_SOA_INDEX_COLUMN(JCollision, collision); diff --git a/PWGJE/TableProducer/CMakeLists.txt b/PWGJE/TableProducer/CMakeLists.txt index 2a82a7f0466..5413eef70db 100644 --- a/PWGJE/TableProducer/CMakeLists.txt +++ b/PWGJE/TableProducer/CMakeLists.txt @@ -53,6 +53,11 @@ o2physics_add_dpl_workflow(jet-eventweight-mcp PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(mc-outlier-rejector + SOURCES mcOutlierRejector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(jet-track-derived SOURCES jetTrackDerived.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore diff --git a/PWGJE/TableProducer/derivedDataProducer.cxx b/PWGJE/TableProducer/derivedDataProducer.cxx index a759b76a2b1..d7684932bad 100644 --- a/PWGJE/TableProducer/derivedDataProducer.cxx +++ b/PWGJE/TableProducer/derivedDataProducer.cxx @@ -256,14 +256,14 @@ struct JetDerivedDataProducerTask { void processMcCollisions(aod::McCollision const& mcCollision) { - products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), 1, 1, 1.0, 1.0); + products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), 1, 1, 1.0, 1.0, 999.0); products.jMcCollisionsParentIndexTable(mcCollision.globalIndex()); } PROCESS_SWITCH(JetDerivedDataProducerTask, processMcCollisions, "produces derived MC collision table", false); void processMcCollisionsWithXsection(soa::Join::iterator const& mcCollision) { - products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr()); + products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr(), mcCollision.ptHard()); products.jMcCollisionsParentIndexTable(mcCollision.globalIndex()); } PROCESS_SWITCH(JetDerivedDataProducerTask, processMcCollisionsWithXsection, "produces derived MC collision table with cross section information", false); diff --git a/PWGJE/TableProducer/derivedDataWriter.cxx b/PWGJE/TableProducer/derivedDataWriter.cxx index 2742c848f83..a1d2346fd7d 100644 --- a/PWGJE/TableProducer/derivedDataWriter.cxx +++ b/PWGJE/TableProducer/derivedDataWriter.cxx @@ -432,7 +432,7 @@ struct JetDerivedDataWriter { mcCollisionMapping.resize(mcCollisions.size(), -1); for (auto const& mcCollision : mcCollisions) { if (mcCollision.isMcCollisionSelected()) { - products.storedJMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.subGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr()); + products.storedJMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.subGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr(), mcCollision.ptHard()); products.storedJMcCollisionsParentIndexTable(mcCollision.mcCollisionId()); mcCollisionMapping[mcCollision.globalIndex()] = products.storedJMcCollisionsTable.lastIndex(); } diff --git a/PWGJE/TableProducer/mcOutlierRejector.cxx b/PWGJE/TableProducer/mcOutlierRejector.cxx new file mode 100644 index 00000000000..427a5afdfad --- /dev/null +++ b/PWGJE/TableProducer/mcOutlierRejector.cxx @@ -0,0 +1,132 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// Task to produce a table joinable to the jcollision table which contains an outlier rejector +// +/// \author Nima Zardoshti + +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetReducedData.h" + +#include "Framework/ASoA.h" +#include "Framework/AnalysisTask.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include +#include +#include +#include + +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct McOutlierRejectorTask { + Produces collisionOutliers; + Produces mcCollisionOutliers; + + Configurable checkmcCollisionForCollision{"checkmcCollisionForCollision", true, "additionally reject collision based on mcCollision"}; + Configurable ptHatMax{"ptHatMax", 4.0, "maximum factor of pt hat the leading jet in the event is allowed"}; + + std::vector collisionFlag; + std::vector mcCollisionFlag; + + void processSetupCollisionSelection(aod::JCollisions const& collisions) + { + collisionFlag.clear(); + collisionFlag.resize(collisions.size(), false); + } + PROCESS_SWITCH(McOutlierRejectorTask, processSetupCollisionSelection, "Setup collision processing", true); + + void processSetupMcCollisionSelection(aod::JMcCollisions const& mcCollisions) + { + mcCollisionFlag.clear(); + mcCollisionFlag.resize(mcCollisions.size(), false); + } + PROCESS_SWITCH(McOutlierRejectorTask, processSetupMcCollisionSelection, "Setup MC Collision processing", true); + + template + void collisionSelection(int32_t collisionIndex, T const& selectionObjects, float ptHard, std::vector& flagArray) + { + + if (selectionObjects.size() != 0) { + float maxSelectionObjectPt = 0.0; + if constexpr (std::is_same_v, aod::JetTracks> || std::is_same_v, aod::JetParticles>) { + for (auto selectionObject : selectionObjects) { + if (selectionObject.pt() > maxSelectionObjectPt) { + maxSelectionObjectPt = selectionObject.pt(); + } + } + } else { + maxSelectionObjectPt = selectionObjects.iteratorAt(0).pt(); + } + + if (maxSelectionObjectPt > ptHatMax * ptHard) { + flagArray[collisionIndex] = true; // Currently if running multiple different jet finders, then a single type of jet can veto an event for others. Decide if this is the best way + } + } + } + + template + void processSelectionObjects(aod::JetCollisionMCD const& collision, T const& selectionObjects, aod::JetMcCollisions const&) + { + auto mcCollision = collision.mcCollision_as(); + collisionSelection(collision.globalIndex(), selectionObjects, mcCollision.ptHard(), collisionFlag); + } + + template + void processSelectionMcObjects(aod::JetMcCollision const& mcCollision, T const& selectionMcObjects) + { + collisionSelection(mcCollision.globalIndex(), selectionMcObjects, mcCollision.ptHard(), mcCollisionFlag); + } + + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingChargedMCDetectorLevelJets, "process mc detector level charged jets", true); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingNeutralMCDetectorLevelJets, "process mc detector level neutral jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingFullMCDetectorLevelJets, "process mc detector level full jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingD0ChargedMCDetectorLevelJets, "process mc detector level D0 charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingDplusChargedMCDetectorLevelJets, "process mc detector level Dplus charged jets", false); + // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingDstarChargedMCDetectorLevelJets, "process mc detector level Dstar charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingLcChargedMCDetectorLevelJets, "process mc detector level Lc charged jets", false); + // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingB0ChargedMCDetectorLevelJets, "process mc detector level B0 charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingBplusChargedMCDetectorLevelJets, "process mc detector level Bplus charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingDielectronChargedMCDetectorLevelJets, "process mc detector level Dielectron charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects, processSelectingTracks, "process tracks", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingChargedMCParticleLevelJets, "process mc particle level charged jets", true); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingNeutralMCParticleLevelJets, "process mc particle level neutral jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingFullMCParticleLevelJets, "process mc particle level full jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingD0ChargedMCParticleLevelJets, "process mc particle level D0 charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingDplusChargedMCParticleLevelJets, "process mc particle level Dplus charged jets", false); + // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingDstarChargedMCParticleLevelJets, "process mc particle level Dstar charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingLcChargedMCParticleLevelJets, "process mc particle level Lc charged jets", false); + // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingB0ChargedMCParticleLevelJets, "process mc particle level B0 charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingBplusChargedMCParticleLevelJets, "process mc particle level Bplus charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingDielectronChargedMCParticleLevelJets, "process mc particle level Dielectron charged jets", false); + PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects, processSelectingParticles, "process mc particles", false); + + void processStoreCollisionDecision(aod::JetCollisionMCD const& collision) + { + bool rejectCollision = collisionFlag[collision.globalIndex()]; + if (!rejectCollision && checkmcCollisionForCollision && collision.has_mcCollision()) { + rejectCollision = mcCollisionFlag[collision.mcCollisionId()]; + } + collisionOutliers(rejectCollision); + } + PROCESS_SWITCH(McOutlierRejectorTask, processStoreCollisionDecision, "write out decision of rejecting collision", true); + + void processStoreMcCollisionDecision(aod::JetMcCollision const& mcCollision) + { + mcCollisionOutliers(mcCollisionFlag[mcCollision.globalIndex()]); + } + PROCESS_SWITCH(McOutlierRejectorTask, processStoreMcCollisionDecision, "write out decision of rejecting mcCollision", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"mc-outlier-rejector"})}; }