Skip to content

Commit 8ecc833

Browse files
mhemmer-cernmhemmer-cern
committed
[PWGEM] Add option to separate EMCal and DCal cluster pairing for flow task
1 parent 78c1280 commit 8ecc833

6 files changed

Lines changed: 51 additions & 16 deletions

File tree

‎PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ enum Observable {
4747
};
4848

4949
// Values in tables are stored in downscaled format to save disk space
50-
const std::array<float, nObservables> downscalingFactors{
50+
constexpr std::array<float, nObservables> downscalingFactors{
5151
1E0, // Cluster definition
5252
1E3, // Cluster energy
5353
1E4, // Cluster eta

‎PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx‎

Lines changed: 13 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -27,12 +27,14 @@
2727

2828
#include <CCDB/BasicCCDBManager.h>
2929
#include <CommonConstants/MathConstants.h>
30+
#include <CommonConstants/PhysicsConstants.h>
3031
#include <DataFormatsParameters/GRPMagField.h>
3132
#include <DataFormatsParameters/GRPObject.h>
3233
#include <Framework/ASoA.h>
3334
#include <Framework/ASoAHelpers.h>
3435
#include <Framework/AnalysisHelpers.h>
3536
#include <Framework/AnalysisTask.h>
37+
#include <Framework/Concepts.h>
3638
#include <Framework/Configurable.h>
3739
#include <Framework/HistogramRegistry.h>
3840
#include <Framework/HistogramSpec.h>
@@ -51,7 +53,12 @@
5153
#include <TPDGCode.h>
5254
#include <TTree.h>
5355

56+
#include <sys/types.h>
57+
58+
#include <algorithm>
5459
#include <cmath>
60+
#include <cstddef>
61+
#include <cstdint>
5562
#include <random>
5663
#include <string>
5764
#include <utility>
@@ -148,14 +155,16 @@ ClusterMcInfo classifyCluster(const TGroup& g, TIter& mcCluster, TIter& mcCluste
148155
info.isPhoton = true;
149156
}
150157

151-
info.photonId = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta});
158+
info.photonId = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime});
152159
info.isFromPi0 = info.photonId >= 0;
153160

154161
return info;
155162
}
163+
156164
struct EmcalPhotonMcTask {
157165
static constexpr float EMCALRadius = 440.f;
158166
static constexpr float PhiVUndefined = -999.f;
167+
static constexpr float Epsilon = 1.e-6f;
159168

160169
Produces<aod::ConvTagCandidates> convTagCandidates;
161170

@@ -241,7 +250,7 @@ struct EmcalPhotonMcTask {
241250
using McColls = o2::soa::Join<o2::aod::EMMCEvents, o2::aod::BinnedGenPts>;
242251
using McParticles = EMMCParticles;
243252

244-
PresliceOptional<EMCalPhotons> perCollisionEMC = o2::aod::emccluster::pmeventId;
253+
PresliceOptional<aod::EMCEMEventIds> perCollisionEMC = o2::aod::emccluster::pmeventId;
245254
PresliceOptional<MinMTracks> perEMCClusterMT = o2::aod::mintm::minClusterId;
246255
PresliceOptional<MinMSTracks> perEMCClusterMS = o2::aod::mintm::minClusterId;
247256

@@ -515,7 +524,7 @@ struct EmcalPhotonMcTask {
515524
ROOT::Math::XYZVector nRef = u.Cross(zAxis); // normal to the plane containing u and z
516525

517526
float phiV = PhiVUndefined; // sentinel for degenerate geometry
518-
if (nDecay.R() > 1e-6f && nRef.R() > 1e-6f) {
527+
if (nDecay.R() > Epsilon && nRef.R() > Epsilon) {
519528
float cosPhiV = static_cast<float>(nDecay.Unit().Dot(nRef.Unit()));
520529
cosPhiV = std::clamp(cosPhiV, -1.f, 1.f);
521530
phiV = std::acos(cosPhiV);
@@ -627,7 +636,7 @@ struct EmcalPhotonMcTask {
627636
}
628637

629638
mcCluster1.setCursor(cluster.emmcparticleIds()[0]);
630-
int photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster1, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta});
639+
int photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster1, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime});
631640
int motherId = -1;
632641
if (photonid1 >= 0) {
633642
mcPhoton1.setCursor(photonid1);

‎PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ using namespace o2::aod::pwgem::photon;
5050

5151
enum QvecEstimator {
5252
FT0M = 0,
53-
FT0A = 1,
53+
FT0A,
5454
FT0C,
5555
TPCPos,
5656
TPCNeg,
@@ -60,7 +60,7 @@ enum QvecEstimator {
6060

6161
enum CentralityEstimator {
6262
None = 0,
63-
CFT0A = 1,
63+
CFT0A,
6464
CFT0C,
6565
CFT0M,
6666
NCentralityEstimators
@@ -382,7 +382,7 @@ struct TaskFlowReso {
382382
// no selection on the centrality is applied on purpose to allow for the resolution study in post-processing
383383
return;
384384
}
385-
if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) {
385+
if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange()) || !(collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) {
386386
return;
387387
}
388388
float cent = getCentrality(collision);
@@ -530,7 +530,7 @@ struct TaskFlowReso {
530530

531531
}; // End struct TaskFlowReso
532532

533-
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
533+
WorkflowSpec defineDataProcessing(ConfigContext const& context)
534534
{
535-
return WorkflowSpec{adaptAnalysisTask<TaskFlowReso>(cfgc)};
535+
return WorkflowSpec{adaptAnalysisTask<TaskFlowReso>(context)};
536536
}

‎PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx‎

Lines changed: 29 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@
3838
#include <Framework/AnalysisHelpers.h>
3939
#include <Framework/AnalysisTask.h>
4040
#include <Framework/BinningPolicy.h>
41+
#include <Framework/Concepts.h>
4142
#include <Framework/Configurable.h>
4243
#include <Framework/Expressions.h>
4344
#include <Framework/GroupedCombinations.h>
@@ -172,6 +173,7 @@ struct TaskPi0FlowEMC {
172173
Configurable<bool> cfgEMCUseTM{"cfgEMCUseTM", false, "flag to use EMCal track matching cut or not"};
173174
Configurable<bool> emcUseSecondaryTM{"emcUseSecondaryTM", false, "flag to use EMCal secondary track matching cut or not"};
174175
Configurable<bool> cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"};
176+
Configurable<bool> separateEMCalDCal{"separateEMCalDCal", false, "flag to only pair EMCal with EMCal and DCal with DCal clusters"};
175177
} emccuts;
176178

177179
V0PhotonCut fV0PhotonCut;
@@ -251,6 +253,8 @@ struct TaskPi0FlowEMC {
251253
int runNow = 0;
252254
int runBefore = -1;
253255

256+
static constexpr float MaxPhiEMCal = 3.5f;
257+
254258
// Filter clusterFilter = aod::skimmedcluster::time >= emccuts.cfgEMCminTime && aod::skimmedcluster::time <= emccuts.cfgEMCmaxTime && aod::skimmedcluster::m02 >= emccuts.cfgEMCminM02 && aod::skimmedcluster::m02 <= emccuts.cfgEMCmaxM02 && aod::skimmedcluster::e >= emccuts.cfgEMCminE;
255259
Filter collisionFilter = (nabs(aod::collision::posZ) <= eventcuts.cfgZvtxMax) && (aod::evsel::ft0cOccupancyInTimeRange <= eventcuts.cfgFT0COccupancyMax) && (aod::evsel::ft0cOccupancyInTimeRange >= eventcuts.cfgFT0COccupancyMin);
256260
// using FilteredEMCalPhotons = soa::Filtered<soa::Join<aod::EMCEMEventIds, aod::MinClusters>>;
@@ -260,10 +264,13 @@ struct TaskPi0FlowEMC {
260264
using CollsWithQvecs = soa::Join<aod::PMEvents, aod::EMEventsAlias, aod::EMEventsMult_000, aod::EMEventsCent_000, aod::EMEventsQvec_001>;
261265
using Colls = soa::Join<aod::PMEvents, aod::EMEventsAlias, aod::EMEventsMult_000, aod::EMEventsCent_000>;
262266

267+
Partition<EMCalPhotons> emcalPhotons = aod::mincluster::storedPhi < std::lround(MaxPhiEMCal * emcdownscaling::downscalingFactors[emcdownscaling::kPhi]);
268+
Partition<EMCalPhotons> dcalPhotons = aod::mincluster::storedPhi >= std::lround(MaxPhiEMCal * emcdownscaling::downscalingFactors[emcdownscaling::kPhi]);
269+
263270
static constexpr std::size_t NQVecEntries = 6;
264271

265-
PresliceOptional<EMCalPhotons> perCollisionEMC = o2::aod::emccluster::pmeventId;
266-
PresliceOptional<PCMPhotons> perCollisionPCM = aod::v0photonkf::pmeventId;
272+
PresliceOptional<aod::EMCEMEventIds> perCollisionEMC = o2::aod::emccluster::pmeventId;
273+
PresliceOptional<aod::V0KFEMEventIds> perCollisionPCM = aod::v0photonkf::pmeventId;
267274
PresliceOptional<MinMTracks> perEMCClusterMT = o2::aod::mintm::minClusterId;
268275
PresliceOptional<MinMSTracks> perEMCClusterMS = o2::aod::mintm::minClusterId;
269276

@@ -481,6 +488,13 @@ struct TaskPi0FlowEMC {
481488

482489
}; // end init
483490

491+
/// \brief Check whether a photon (by its phi) falls in the EMCal or DCal acceptance
492+
/// \param phi azimuthal angle of the photon
493+
static bool isEMCalRegion(float phi)
494+
{
495+
return phi < MaxPhiEMCal;
496+
}
497+
484498
/// Change radians to degree
485499
/// \param angle in radians
486500
/// \return angle in degree
@@ -1047,7 +1061,7 @@ struct TaskPi0FlowEMC {
10471061
continue;
10481062
}
10491063
if (rotationConfig.cfgDoRotation.value && nColl % rotationConfig.cfgDownsampling == 0) {
1050-
rotationBackground<EMCalPhotons>(vMeson, v1, v2, photons1, g1.globalIndex(), g2.globalIndex(), collision);
1064+
rotationBackground<TPhotons1>(vMeson, v1, v2, photons1, g1.globalIndex(), g2.globalIndex(), collision);
10511065
}
10521066
if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) {
10531067
registry.fill(HIST("hMesonCuts"), 3);
@@ -1093,6 +1107,8 @@ struct TaskPi0FlowEMC {
10931107
}
10941108

10951109
auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex());
1110+
auto emcalPhotonsPerCollision = emcalPhotons.sliceBy(perCollisionEMC, collision.globalIndex());
1111+
auto dcalPhotonsPerCollision = dcalPhotons.sliceBy(perCollisionEMC, collision.globalIndex());
10961112

10971113
if (emccuts.cfgEnableQA.value) {
10981114
for (const auto& photon : photonsPerCollision) {
@@ -1108,7 +1124,12 @@ struct TaskPi0FlowEMC {
11081124
registry.fill(HIST("clusterQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // after cuts
11091125
}
11101126
}
1111-
runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags);
1127+
if (emccuts.separateEMCalDCal.value) {
1128+
runPairingLoop(collision, emcalPhotonsPerCollision, emcalPhotonsPerCollision, flags, flags);
1129+
runPairingLoop(collision, dcalPhotonsPerCollision, dcalPhotonsPerCollision, flags, flags);
1130+
} else {
1131+
runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags);
1132+
}
11121133
if (rotationConfig.cfgDoRotation.value) {
11131134
if (nColl % rotationConfig.cfgDownsampling == 0) {
11141135
nColl = 1; // reset counter
@@ -1163,6 +1184,10 @@ struct TaskPi0FlowEMC {
11631184
if (!(flags.test(g1.globalIndex())) || !(flags.test(g2.globalIndex()))) {
11641185
continue;
11651186
}
1187+
if (emccuts.separateEMCalDCal.value && isEMCalRegion(g1.phi()) != isEMCalRegion(g2.phi())) {
1188+
continue; // only pair EMCal-EMCal or DCal-DCal
1189+
}
1190+
11661191
// Cut edge clusters away, similar to rotation method to ensure same acceptance is used
11671192
if (cfgDistanceToEdge.value > 0) {
11681193
if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= cfgEMCalMapLevelBackground.value) {

‎PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -172,7 +172,7 @@ struct TestTaskEmc {
172172
PROCESS_SWITCH(TestTaskEmc, processEMCalCalib, "Process EMCal calibration same event", true);
173173
}; // End struct TestTaskEmc
174174

175-
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
175+
WorkflowSpec defineDataProcessing(ConfigContext const& context)
176176
{
177-
return WorkflowSpec{adaptAnalysisTask<TestTaskEmc>(cfgc)};
177+
return WorkflowSpec{adaptAnalysisTask<TestTaskEmc>(context)};
178178
}

‎PWGEM/PhotonMeson/Utils/MCUtilities.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#define PWGEM_PHOTONMESON_UTILS_MCUTILITIES_H_
1818

1919
#include <Framework/ASoA.h>
20+
#include <Framework/Concepts.h>
2021

2122
#include <TPDGCode.h>
2223

0 commit comments

Comments
 (0)