Skip to content

Commit bee76a3

Browse files
dsekihatdsekihat
authored
[PWGEM/Dilepton] add a task to associate mc col to rec col (#9204)
1 parent 3e4aa5b commit bee76a3

File tree

6 files changed

+149
-63
lines changed

6 files changed

+149
-63
lines changed

‎PWGEM/Dilepton/Core/DileptonMC.h‎

Lines changed: 26 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,9 @@ using namespace o2::aod::pwgem::dilepton::utils::pairutil;
6565
using MyCollisions = soa::Join<aod::EMEvents, aod::EMEventsMult, aod::EMEventsCent, aod::EMMCEventLabels>;
6666
using MyCollision = MyCollisions::iterator;
6767

68+
using MyMCCollisions = soa::Join<aod::EMMCEvents, aod::MostProbableEMEventIdsInMC>;
69+
using MyMCCollision = MyMCCollisions::iterator;
70+
6871
using MyMCElectrons = soa::Join<aod::EMPrimaryElectrons, aod::EMPrimaryElectronsCov, aod::EMPrimaryElectronEMEventIds, aod::EMAmbiguousElectronSelfIds, aod::EMPrimaryElectronsPrefilterBit, aod::EMPrimaryElectronsPrefilterBitDerived, aod::EMPrimaryElectronMCLabels>;
6972
using MyMCElectron = MyMCElectrons::iterator;
7073
using FilteredMyMCElectrons = soa::Filtered<MyMCElectrons>;
@@ -1103,22 +1106,26 @@ struct DileptonMC {
11031106
}
11041107
fRegistry.fill(HIST("MCEvent/before/hZvtx"), mccollision.posZ());
11051108

1106-
auto rec_colls_per_mccoll = collisions.sliceBy(recColperMcCollision, mccollision.globalIndex());
1107-
// LOGF(info, "rec_colls_per_mccoll.size() = %d", rec_colls_per_mccoll.size());
1108-
if (rec_colls_per_mccoll.size() < 1) {
1109-
continue;
1110-
}
1109+
// auto rec_colls_per_mccoll = collisions.sliceBy(recColperMcCollision, mccollision.globalIndex());
1110+
// // LOGF(info, "rec_colls_per_mccoll.size() = %d", rec_colls_per_mccoll.size());
1111+
// if (rec_colls_per_mccoll.size() < 1) {
1112+
// continue;
1113+
// }
11111114

1112-
uint32_t maxNumContrib = 0;
1113-
int rec_col_globalIndex = -999;
1114-
for (auto& rec_col : rec_colls_per_mccoll) {
1115-
if (rec_col.numContrib() > maxNumContrib) {
1116-
rec_col_globalIndex = rec_col.globalIndex();
1117-
maxNumContrib = rec_col.numContrib(); // assign mc collision to collision where the number of contibutor is lager. LF/MM recommendation
1118-
}
1119-
}
1115+
// uint32_t maxNumContrib = 0;
1116+
// int rec_col_globalIndex = -999;
1117+
// for (auto& rec_col : rec_colls_per_mccoll) {
1118+
// if (rec_col.numContrib() > maxNumContrib) {
1119+
// rec_col_globalIndex = rec_col.globalIndex();
1120+
// maxNumContrib = rec_col.numContrib(); // assign mc collision to collision where the number of contibutor is lager. LF/MM recommendation
1121+
// }
1122+
// }
1123+
// auto collision = collisions.rawIteratorAt(rec_col_globalIndex);
11201124

1121-
auto collision = collisions.rawIteratorAt(rec_col_globalIndex);
1125+
if (mccollision.mpemeventId() < 0) {
1126+
continue;
1127+
}
1128+
auto collision = collisions.rawIteratorAt(mccollision.mpemeventId());
11221129

11231130
float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()};
11241131
if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) {
@@ -2101,9 +2108,9 @@ struct DileptonMC {
21012108
Partition<aod::EMMCParticles> negative_muonsMC = o2::aod::mcparticle::pdgCode == 13; // mu-
21022109
PresliceUnsorted<aod::EMMCParticles> perMcCollision = aod::emmcparticle::emmceventId;
21032110
PresliceUnsorted<aod::EMMCGenVectorMesons> perMcCollision_vm = aod::emmcgenvectormeson::emmceventId;
2104-
PresliceUnsorted<MyCollisions> recColperMcCollision = aod::emmceventlabel::emmceventId;
2111+
// PresliceUnsorted<MyCollisions> recColperMcCollision = aod::emmceventlabel::emmceventId;
21052112

2106-
void processAnalysis(FilteredMyCollisions const& collisions, aod::EMMCEvents const& mccollisions, aod::EMMCParticles const& mcparticles, TLeptons const& leptons)
2113+
void processAnalysis(FilteredMyCollisions const& collisions, MyMCCollisions const& mccollisions, aod::EMMCParticles const& mcparticles, TLeptons const& leptons)
21072114
{
21082115
// LOGF(info, "collisions.size() = %d, mccollisions.size() = %d, mcparticles.size() = %d", collisions.size(), mccollisions.size(), mcparticles.size());
21092116
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
@@ -2134,7 +2141,7 @@ struct DileptonMC {
21342141
Partition<MySmearedMuons> positive_muonsMC_smeared = o2::aod::mcparticle::pdgCode == -13; // mu+
21352142
Partition<MySmearedMuons> negative_muonsMC_smeared = o2::aod::mcparticle::pdgCode == 13; // mu-
21362143

2137-
void processAnalysis_Smeared(FilteredMyCollisions const& collisions, aod::EMMCEvents const& mccollisions, TLeptons const& leptons, TSmeardMCParitlces const& mcparticles_smeared)
2144+
void processAnalysis_Smeared(FilteredMyCollisions const& collisions, MyMCCollisions const& mccollisions, TLeptons const& leptons, TSmeardMCParitlces const& mcparticles_smeared)
21382145
{
21392146
// LOGF(info, "collisions.size() = %d, mccollisions.size() = %d, mcparticles.size() = %d", collisions.size(), mccollisions.size(), mcparticles.size());
21402147
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
@@ -2160,7 +2167,7 @@ struct DileptonMC {
21602167
}
21612168
PROCESS_SWITCH(DileptonMC, processAnalysis_Smeared, "run dilepton mc analysis with smearing", false);
21622169

2163-
void processGen_VM(FilteredMyCollisions const& collisions, aod::EMMCEvents const&, aod::EMMCGenVectorMesons const& mcparticles)
2170+
void processGen_VM(FilteredMyCollisions const& collisions, MyMCCollisions const&, aod::EMMCGenVectorMesons const& mcparticles)
21642171
{
21652172
// for oemga, phi efficiency
21662173
for (auto& collision : collisions) {
@@ -2172,7 +2179,7 @@ struct DileptonMC {
21722179
if (!fEMEventCut.IsSelected(collision)) {
21732180
continue;
21742181
}
2175-
auto mccollision = collision.template emmcevent_as<aod::EMMCEvents>();
2182+
auto mccollision = collision.template emmcevent_as<MyMCCollisions>();
21762183
if (cfgEventGeneratorType >= 0 && mccollision.getSubGeneratorId() != cfgEventGeneratorType) {
21772184
continue;
21782185
}

‎PWGEM/Dilepton/Core/SingleTrackQCMC.h‎

Lines changed: 22 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,9 @@ using namespace o2::aod::pwgem::dilepton::utils::emtrackutil;
5555
using MyCollisions = soa::Join<aod::EMEvents, aod::EMEventsMult, aod::EMEventsCent, aod::EMMCEventLabels>;
5656
using MyCollision = MyCollisions::iterator;
5757

58+
using MyMCCollisions = soa::Join<aod::EMMCEvents, aod::MostProbableEMEventIdsInMC>;
59+
using MyMCCollision = MyMCCollisions::iterator;
60+
5861
using MyMCElectrons = soa::Join<aod::EMPrimaryElectrons, aod::EMPrimaryElectronsCov, aod::EMPrimaryElectronEMEventIds, aod::EMAmbiguousElectronSelfIds, aod::EMPrimaryElectronsPrefilterBit, aod::EMPrimaryElectronMCLabels>;
5962
using MyMCElectron = MyMCElectrons::iterator;
6063
using FilteredMyMCElectrons = soa::Filtered<MyMCElectrons>;
@@ -763,20 +766,24 @@ struct SingleTrackQCMC {
763766
}
764767
fRegistry.fill(HIST("MCEvent/before/hZvtx"), mccollision.posZ());
765768

766-
auto rec_colls_per_mccoll = collisions.sliceBy(recColperMcCollision, mccollision.globalIndex());
767-
if (rec_colls_per_mccoll.size() < 1) {
768-
continue;
769-
}
769+
// auto rec_colls_per_mccoll = collisions.sliceBy(recColperMcCollision, mccollision.globalIndex());
770+
// if (rec_colls_per_mccoll.size() < 1) {
771+
// continue;
772+
// }
773+
// uint32_t maxNumContrib = 0;
774+
// int rec_col_globalIndex = -999;
775+
// for (auto& rec_col : rec_colls_per_mccoll) {
776+
// if (rec_col.numContrib() > maxNumContrib) {
777+
// rec_col_globalIndex = rec_col.globalIndex();
778+
// maxNumContrib = rec_col.numContrib(); // assign mc collision to collision where the number of contibutor is lager. LF/MM recommendation
779+
// }
780+
// }
781+
// auto collision = collisions.rawIteratorAt(rec_col_globalIndex);
770782

771-
uint32_t maxNumContrib = 0;
772-
int rec_col_globalIndex = -999;
773-
for (auto& rec_col : rec_colls_per_mccoll) {
774-
if (rec_col.numContrib() > maxNumContrib) {
775-
rec_col_globalIndex = rec_col.globalIndex();
776-
maxNumContrib = rec_col.numContrib(); // assign mc collision to collision where the number of contibutor is lager. LF/MM recommendation
777-
}
783+
if (mccollision.mpemeventId() < 0) {
784+
continue;
778785
}
779-
auto collision = collisions.rawIteratorAt(rec_col_globalIndex);
786+
auto collision = collisions.rawIteratorAt(mccollision.mpemeventId());
780787

781788
float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()};
782789
if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) {
@@ -967,9 +974,9 @@ struct SingleTrackQCMC {
967974
Partition<aod::EMMCParticles> muonsMC = nabs(o2::aod::mcparticle::pdgCode) == 13; // mu+, mu-
968975
PresliceUnsorted<aod::EMMCParticles> perMcCollision = aod::emmcparticle::emmceventId;
969976

970-
PresliceUnsorted<MyCollisions> recColperMcCollision = aod::emmceventlabel::emmceventId;
977+
// PresliceUnsorted<MyCollisions> recColperMcCollision = aod::emmceventlabel::emmceventId;
971978

972-
void processQCMC(FilteredMyCollisions const& collisions, aod::EMMCEvents const& mccollisions, aod::EMMCParticles const& mcparticles, TLeptons const& tracks)
979+
void processQCMC(FilteredMyCollisions const& collisions, MyMCCollisions const& mccollisions, aod::EMMCParticles const& mcparticles, TLeptons const& tracks)
973980
{
974981
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
975982
if (cfgApplyWeightTTCA) {
@@ -991,7 +998,7 @@ struct SingleTrackQCMC {
991998
Partition<MySmearedElectrons> electronsMC_smeared = nabs(o2::aod::mcparticle::pdgCode) == 11; // e+, e-
992999
Partition<MySmearedMuons> muonsMC_smeared = nabs(o2::aod::mcparticle::pdgCode) == 13; // mu+, mu-
9931000

994-
void processQCMC_Smeared(FilteredMyCollisions const& collisions, aod::EMMCEvents const& mccollisions, TLeptons const& tracks, TSmearedMCParticles const& mcparticles_smeared)
1001+
void processQCMC_Smeared(FilteredMyCollisions const& collisions, MyMCCollisions const& mccollisions, TLeptons const& tracks, TSmearedMCParticles const& mcparticles_smeared)
9951002
{
9961003
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
9971004
if (cfgApplyWeightTTCA) {

‎PWGEM/Dilepton/DataModel/dileptonTables.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -217,6 +217,7 @@ using EMEventNormInfo = EMEventNormInfos::iterator;
217217

218218
namespace emmcevent
219219
{
220+
DECLARE_SOA_INDEX_COLUMN(EMEvent, mpemevent); //! most propable emeventId
220221
DECLARE_SOA_COLUMN(McCollisionId, mcCollisionId, int);
221222
} // namespace emmcevent
222223

@@ -229,9 +230,11 @@ DECLARE_SOA_TABLE(EMMCEvents, "AOD", "EMMCEVENT", //! MC event information tab
229230
mccollision::GetGeneratorId<mccollision::GeneratorsID>,
230231
mccollision::GetSubGeneratorId<mccollision::GeneratorsID>,
231232
mccollision::GetSourceId<mccollision::GeneratorsID>);
232-
233233
using EMMCEvent = EMMCEvents::iterator;
234234

235+
DECLARE_SOA_TABLE(MostProbableEMEventIdsInMC, "AOD", "MPEMEVENTIDINMC", emmcevent::EMEventId); // To be joined with EMMCEvents table at analysis level.
236+
using MostProbableEMEventIdInMC = MostProbableEMEventIdsInMC::iterator;
237+
235238
namespace emmceventlabel
236239
{
237240
DECLARE_SOA_INDEX_COLUMN(EMMCEvent, emmcevent); //! MC collision

‎PWGEM/Dilepton/Tasks/CMakeLists.txt‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ add_subdirectory(Converters)
1313

1414
o2physics_add_dpl_workflow(efficiency-ee
1515
SOURCES emEfficiencyEE.cxx
16-
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore
16+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::PWGDQCore
1717
COMPONENT_NAME Analysis)
1818

1919
o2physics_add_dpl_workflow(lmee-lf-cocktail
@@ -28,12 +28,12 @@ o2physics_add_dpl_workflow(lmee-hf-cocktail
2828

2929
o2physics_add_dpl_workflow(mc-templates
3030
SOURCES MCtemplates.cxx
31-
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore
31+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::PWGDQCore
3232
COMPONENT_NAME Analysis)
3333

3434
o2physics_add_dpl_workflow(smearing
3535
SOURCES smearing.cxx
36-
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore
36+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
3737
COMPONENT_NAME Analysis)
3838

3939
o2physics_add_dpl_workflow(create-resolution-map
@@ -43,7 +43,7 @@ o2physics_add_dpl_workflow(create-resolution-map
4343

4444
o2physics_add_dpl_workflow(table-reader-barrel
4545
SOURCES tableReaderBarrel.cxx
46-
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore
46+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::PWGDQCore
4747
COMPONENT_NAME Analysis)
4848

4949
o2physics_add_dpl_workflow(event-qc
@@ -121,3 +121,8 @@ o2physics_add_dpl_workflow(prefilter-dielectron
121121
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::MLCore O2Physics::PWGEMDileptonCore
122122
COMPONENT_NAME Analysis)
123123

124+
o2physics_add_dpl_workflow(associate-mccollision-to-collision
125+
SOURCES associateMCcollision.cxx
126+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
127+
COMPONENT_NAME Analysis)
128+

‎PWGEM/Dilepton/Tasks/associateMCcollision.cxx‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
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+
// ========================
13+
//
14+
// This code produces a table with an index between mc collision and rec. collision.
15+
// Please write to: daiki.sekihata@cern.ch
16+
17+
#include "Framework/runDataProcessing.h"
18+
#include "Framework/AnalysisTask.h"
19+
#include "Framework/ASoAHelpers.h"
20+
21+
#include "PWGEM/Dilepton/DataModel/dileptonTables.h"
22+
23+
using namespace o2;
24+
using namespace o2::aod;
25+
using namespace o2::framework;
26+
using namespace o2::framework::expressions;
27+
using namespace o2::soa;
28+
29+
struct associateMCcollision {
30+
Produces<aod::MostProbableEMEventIdsInMC> mpemeventIds;
31+
32+
HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
33+
34+
void init(InitContext&)
35+
{
36+
addhistograms();
37+
}
38+
39+
~associateMCcollision() {}
40+
41+
void addhistograms()
42+
{
43+
fRegistry.add("hReccollsPerMCcoll", "Rec. colls per MC coll;Rec. colls per MC coll;Number of MC collisions", kTH1F, {{21, -0.5, 20.5}}, false);
44+
}
45+
46+
using MyCollisions = soa::Join<aod::EMEvents, aod::EMMCEventLabels>;
47+
using MyCollision = MyCollisions::iterator;
48+
PresliceUnsorted<MyCollisions> recColperMcCollision = aod::emmceventlabel::emmceventId;
49+
50+
void processNcontrib(aod::EMMCEvents const& mccollisions, MyCollisions const& collisions)
51+
{
52+
for (auto& mccollision : mccollisions) {
53+
auto rec_colls_per_mccoll = collisions.sliceBy(recColperMcCollision, mccollision.globalIndex());
54+
fRegistry.fill(HIST("hReccollsPerMCcoll"), rec_colls_per_mccoll.size());
55+
uint32_t maxNumContrib = 0;
56+
int rec_col_globalIndex = -999;
57+
for (auto& rec_col : rec_colls_per_mccoll) {
58+
if (rec_col.numContrib() > maxNumContrib) {
59+
rec_col_globalIndex = rec_col.globalIndex();
60+
maxNumContrib = rec_col.numContrib(); // assign mc collision to collision where the number of contibutors is the lagest. LF/MM recommendation
61+
}
62+
}
63+
// LOGF(info, "rec_col_globalIndex = %d", rec_col_globalIndex);
64+
mpemeventIds(rec_col_globalIndex);
65+
} // end of mc collision
66+
} // end of process
67+
PROCESS_SWITCH(associateMCcollision, processNcontrib, "produce most probable emeventId based on Ncontrib to PV", true);
68+
};
69+
70+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
71+
{
72+
return WorkflowSpec{adaptAnalysisTask<associateMCcollision>(cfgc, TaskName{"associate-mccollision-to-collision"})};
73+
}

0 commit comments

Comments
 (0)