Skip to content

Commit ec98923

Browse files
committed
Adding Outlier rejection task for MC
1 parent f10c7ea commit ec98923

File tree

5 files changed

+151
-4
lines changed

5 files changed

+151
-4
lines changed

PWGJE/DataModel/JetReducedData.h

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -99,6 +99,7 @@ DECLARE_SOA_COLUMN(ReadCountsWithTVXAndZVertexAndSel7KINT7, readCountsWithTVXAnd
9999
DECLARE_SOA_COLUMN(ReadCountsWithCustom, readCountsWithCustom, std::vector<int>);
100100
DECLARE_SOA_COLUMN(IsAmbiguous, isAmbiguous, bool);
101101
DECLARE_SOA_COLUMN(IsEMCALReadout, isEmcalReadout, bool);
102+
DECLARE_SOA_COLUMN(IsOutlier, isOutlier, bool);
102103
} // namespace jcollision
103104

104105
DECLARE_SOA_TABLE_STAGED(JCollisions, "JCOLLISION",
@@ -122,6 +123,9 @@ DECLARE_SOA_TABLE_STAGED(JCollisionMcInfos, "JCOLLISIONMCINFO",
122123
jcollision::Weight,
123124
jcollision::SubGeneratorId);
124125

126+
DECLARE_SOA_TABLE_STAGED(JCollisionOutliers, "JCOLLISIONOUTLR",
127+
jcollision::IsOutlier);
128+
125129
DECLARE_SOA_TABLE_STAGED(JEMCCollisionLbs, "JEMCCOLLISIONLB",
126130
jcollision::IsAmbiguous,
127131
jcollision::IsEMCALReadout);
@@ -170,6 +174,8 @@ DECLARE_SOA_COLUMN(Accepted, accepted, uint64_t);
170174
DECLARE_SOA_COLUMN(Attempted, attempted, uint64_t);
171175
DECLARE_SOA_COLUMN(XsectGen, xsectGen, float);
172176
DECLARE_SOA_COLUMN(XsectErr, xsectErr, float);
177+
DECLARE_SOA_COLUMN(PtHard, ptHard, float);
178+
DECLARE_SOA_COLUMN(IsOutlier, isOutlier, bool);
173179
} // namespace jmccollision
174180
DECLARE_SOA_TABLE_STAGED(JMcCollisions, "JMCCOLLISION",
175181
o2::soa::Index<>,
@@ -181,7 +187,8 @@ DECLARE_SOA_TABLE_STAGED(JMcCollisions, "JMCCOLLISION",
181187
jmccollision::Accepted,
182188
jmccollision::Attempted,
183189
jmccollision::XsectGen,
184-
jmccollision::XsectErr);
190+
jmccollision::XsectErr,
191+
jmccollision::PtHard);
185192

186193
using JMcCollision = JMcCollisions::iterator;
187194
using StoredJMcCollision = StoredJMcCollisions::iterator;
@@ -197,6 +204,9 @@ DECLARE_SOA_INDEX_COLUMN(JMcCollision, mcCollision);
197204
DECLARE_SOA_TABLE_STAGED(JMcCollisionLbs, "JMCCOLLISIONLB",
198205
jmccollisionlb::JMcCollisionId);
199206

207+
DECLARE_SOA_TABLE_STAGED(JMcCollisionOutliers, "JMCCOLLISIONOUTLR",
208+
jmccollision::IsOutlier);
209+
200210
namespace jtrack
201211
{
202212
DECLARE_SOA_INDEX_COLUMN(JCollision, collision);

PWGJE/TableProducer/CMakeLists.txt

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,11 @@ o2physics_add_dpl_workflow(jet-eventweight-mcp
5353
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
5454
COMPONENT_NAME Analysis)
5555

56+
o2physics_add_dpl_workflow(mc-outlier-rejector
57+
SOURCES mcOutlierRejector.cxx
58+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
59+
COMPONENT_NAME Analysis)
60+
5661
o2physics_add_dpl_workflow(jet-track-derived
5762
SOURCES jetTrackDerived.cxx
5863
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore

PWGJE/TableProducer/derivedDataProducer.cxx

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -256,14 +256,14 @@ struct JetDerivedDataProducerTask {
256256

257257
void processMcCollisions(aod::McCollision const& mcCollision)
258258
{
259-
products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), 1, 1, 1.0, 1.0);
259+
products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), 1, 1, 1.0, 1.0, 999.0);
260260
products.jMcCollisionsParentIndexTable(mcCollision.globalIndex());
261261
}
262262
PROCESS_SWITCH(JetDerivedDataProducerTask, processMcCollisions, "produces derived MC collision table", false);
263263

264264
void processMcCollisionsWithXsection(soa::Join<aod::McCollisions, aod::HepMCXSections>::iterator const& mcCollision)
265265
{
266-
products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr());
266+
products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr(), mcCollision.ptHard());
267267
products.jMcCollisionsParentIndexTable(mcCollision.globalIndex());
268268
}
269269
PROCESS_SWITCH(JetDerivedDataProducerTask, processMcCollisionsWithXsection, "produces derived MC collision table with cross section information", false);

PWGJE/TableProducer/derivedDataWriter.cxx

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -432,7 +432,7 @@ struct JetDerivedDataWriter {
432432
mcCollisionMapping.resize(mcCollisions.size(), -1);
433433
for (auto const& mcCollision : mcCollisions) {
434434
if (mcCollision.isMcCollisionSelected()) {
435-
products.storedJMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.subGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr());
435+
products.storedJMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.subGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr(), mcCollision.ptHard());
436436
products.storedJMcCollisionsParentIndexTable(mcCollision.mcCollisionId());
437437
mcCollisionMapping[mcCollision.globalIndex()] = products.storedJMcCollisionsTable.lastIndex();
438438
}
Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,132 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
// Task to produce a table joinable to the jcollision table which contains an outlier rejector
13+
//
14+
/// \author Nima Zardoshti <nima.zardoshti@cern.ch>
15+
16+
#include "PWGJE/DataModel/Jet.h"
17+
#include "PWGJE/DataModel/JetReducedData.h"
18+
19+
#include "Framework/ASoA.h"
20+
#include "Framework/AnalysisTask.h"
21+
#include "Framework/O2DatabasePDGPlugin.h"
22+
#include <Framework/AnalysisHelpers.h>
23+
#include <Framework/Configurable.h>
24+
#include <Framework/InitContext.h>
25+
#include <Framework/runDataProcessing.h>
26+
27+
#include <vector>
28+
29+
using namespace o2;
30+
using namespace o2::framework;
31+
using namespace o2::framework::expressions;
32+
33+
struct McOutlierRejectorTask {
34+
Produces<aod::JCollisionOutliers> collisionOutliers;
35+
Produces<aod::JMcCollisionOutliers> mcCollisionOutliers;
36+
37+
Configurable<bool> checkmcCollisionForCollision{"checkmcCollisionForCollision", true, "additionally reject collision based on mcCollision"};
38+
Configurable<float> ptHatMax{"ptHatMax", 4.0, "maximum factor of pt hat the leading jet in the event is allowed"};
39+
40+
std::vector<bool> collisionFlag;
41+
std::vector<bool> mcCollisionFlag;
42+
43+
void processSetupCollisionSelection(aod::JCollisions const& collisions)
44+
{
45+
collisionFlag.clear();
46+
collisionFlag.resize(collisions.size(), false);
47+
}
48+
PROCESS_SWITCH(McOutlierRejectorTask, processSetupCollisionSelection, "Setup collision processing", true);
49+
50+
void processSetupMcCollisionSelection(aod::JMcCollisions const& mcCollisions)
51+
{
52+
mcCollisionFlag.clear();
53+
mcCollisionFlag.resize(mcCollisions.size(), false);
54+
}
55+
PROCESS_SWITCH(McOutlierRejectorTask, processSetupMcCollisionSelection, "Setup MC Collision processing", true);
56+
57+
template <typename T>
58+
void collisionSelection(int32_t collisionIndex, T const& selectionObjects, float ptHard, std::vector<bool>& flagArray)
59+
{
60+
61+
if (selectionObjects.size() != 0) {
62+
float maxSelectionObjectPt = 0.0;
63+
if constexpr (std::is_same_v<std::decay_t<T>, aod::JetTracks> || std::is_same_v<std::decay_t<T>, aod::JetParticles>) {
64+
for (auto selectionObject : selectionObjects) {
65+
if (selectionObject.pt() > maxSelectionObjectPt) {
66+
maxSelectionObjectPt = selectionObject.pt();
67+
}
68+
}
69+
} else {
70+
maxSelectionObjectPt = selectionObjects.iteratorAt(0).pt();
71+
}
72+
73+
if (maxSelectionObjectPt > ptHatMax * ptHard) {
74+
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
75+
}
76+
}
77+
}
78+
79+
template <typename T>
80+
void processSelectionObjects(aod::JetCollisionMCD const& collision, T const& selectionObjects, aod::JetMcCollisions const&)
81+
{
82+
auto mcCollision = collision.mcCollision_as<aod::JetMcCollisions>();
83+
collisionSelection(collision.globalIndex(), selectionObjects, mcCollision.ptHard(), collisionFlag);
84+
}
85+
86+
template <typename T>
87+
void processSelectionMcObjects(aod::JetMcCollision const& mcCollision, T const& selectionMcObjects)
88+
{
89+
collisionSelection(mcCollision.globalIndex(), selectionMcObjects, mcCollision.ptHard(), mcCollisionFlag);
90+
}
91+
92+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::ChargedMCDetectorLevelJets>, processSelectingChargedMCDetectorLevelJets, "process mc detector level charged jets", true);
93+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::NeutralMCDetectorLevelJets>, processSelectingNeutralMCDetectorLevelJets, "process mc detector level neutral jets", false);
94+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::FullMCDetectorLevelJets>, processSelectingFullMCDetectorLevelJets, "process mc detector level full jets", false);
95+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::D0ChargedMCDetectorLevelJets>, processSelectingD0ChargedMCDetectorLevelJets, "process mc detector level D0 charged jets", false);
96+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::DplusChargedMCDetectorLevelJets>, processSelectingDplusChargedMCDetectorLevelJets, "process mc detector level Dplus charged jets", false);
97+
// PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::DstarChargedMCDetectorLevelJets>, processSelectingDstarChargedMCDetectorLevelJets, "process mc detector level Dstar charged jets", false);
98+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::LcChargedMCDetectorLevelJets>, processSelectingLcChargedMCDetectorLevelJets, "process mc detector level Lc charged jets", false);
99+
// PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::B0ChargedMCDetectorLevelJets>, processSelectingB0ChargedMCDetectorLevelJets, "process mc detector level B0 charged jets", false);
100+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::BplusChargedMCDetectorLevelJets>, processSelectingBplusChargedMCDetectorLevelJets, "process mc detector level Bplus charged jets", false);
101+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::DielectronChargedMCDetectorLevelJets>, processSelectingDielectronChargedMCDetectorLevelJets, "process mc detector level Dielectron charged jets", false);
102+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::JetTracks>, processSelectingTracks, "process tracks", false);
103+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::ChargedMCParticleLevelJets>, processSelectingChargedMCParticleLevelJets, "process mc particle level charged jets", true);
104+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::NeutralMCParticleLevelJets>, processSelectingNeutralMCParticleLevelJets, "process mc particle level neutral jets", false);
105+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::FullMCParticleLevelJets>, processSelectingFullMCParticleLevelJets, "process mc particle level full jets", false);
106+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::D0ChargedMCParticleLevelJets>, processSelectingD0ChargedMCParticleLevelJets, "process mc particle level D0 charged jets", false);
107+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::DplusChargedMCParticleLevelJets>, processSelectingDplusChargedMCParticleLevelJets, "process mc particle level Dplus charged jets", false);
108+
// PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::DstarChargedMCParticleLevelJets>, processSelectingDstarChargedMCParticleLevelJets, "process mc particle level Dstar charged jets", false);
109+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::LcChargedMCParticleLevelJets>, processSelectingLcChargedMCParticleLevelJets, "process mc particle level Lc charged jets", false);
110+
// PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::B0ChargedMCParticleLevelJets>, processSelectingB0ChargedMCParticleLevelJets, "process mc particle level B0 charged jets", false);
111+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::BplusChargedMCParticleLevelJets>, processSelectingBplusChargedMCParticleLevelJets, "process mc particle level Bplus charged jets", false);
112+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::DielectronChargedMCParticleLevelJets>, processSelectingDielectronChargedMCParticleLevelJets, "process mc particle level Dielectron charged jets", false);
113+
PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::JetParticles>, processSelectingParticles, "process mc particles", false);
114+
115+
void processStoreCollisionDecision(aod::JetCollisionMCD const& collision)
116+
{
117+
bool rejectCollision = collisionFlag[collision.globalIndex()];
118+
if (!rejectCollision && checkmcCollisionForCollision && collision.has_mcCollision()) {
119+
rejectCollision = mcCollisionFlag[collision.mcCollisionId()];
120+
}
121+
collisionOutliers(rejectCollision);
122+
}
123+
PROCESS_SWITCH(McOutlierRejectorTask, processStoreCollisionDecision, "write out decision of rejecting collision", true);
124+
125+
void processStoreMcCollisionDecision(aod::JetMcCollision const& mcCollision)
126+
{
127+
mcCollisionOutliers(mcCollisionFlag[mcCollision.globalIndex()]);
128+
}
129+
PROCESS_SWITCH(McOutlierRejectorTask, processStoreMcCollisionDecision, "write out decision of rejecting mcCollision", true);
130+
};
131+
132+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<McOutlierRejectorTask>(cfgc, TaskName{"mc-outlier-rejector"})}; }

0 commit comments

Comments
 (0)