Skip to content

Commit b875773

Browse files
ariedel-cernariedel-cern
committed
Feat: add event shape table
1 parent cb509c0 commit b875773

3 files changed

Lines changed: 139 additions & 44 deletions

File tree

‎PWGCF/Femto/Core/collisionBuilder.h‎

Lines changed: 71 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525

2626
#include "Common/CCDB/EventSelectionParams.h"
2727
#include "Common/CCDB/RCTSelectionFlags.h"
28+
#include "Common/Core/RecoDecay.h"
2829
#include "Common/Core/Zorro.h"
2930

3031
#include <DataFormatsParameters/GRPMagField.h>
@@ -33,6 +34,10 @@
3334
#include <Framework/HistogramRegistry.h>
3435
#include <Framework/Logger.h>
3536

37+
#include <sys/stat.h>
38+
39+
#include <Rtypes.h>
40+
3641
#include <algorithm>
3742
#include <cmath>
3843
#include <cstddef>
@@ -80,6 +85,8 @@ struct ConfCollisionBits : o2::framework::ConfigurableGroup {
8085
o2::framework::Configurable<std::vector<float>> sphericityMin{"sphericityMin", {}, "Minimum sphericity"};
8186
o2::framework::Configurable<std::vector<float>> sphericityMax{"sphericityMax", {}, "Maximum sphericity"};
8287
o2::framework::Configurable<std::vector<std::string>> triggers{"triggers", {}, "List of all triggers to be used"};
88+
o2::framework::Configurable<datatypes::QvecDetectorType> qvecDetector{"qvecDetector", 0, "Detector used to estimate the Q-vector: 0 -> FT0C, 1 -> FT0A"};
89+
o2::framework::Configurable<datatypes::QvecHarmonicType> qvecHarmonic{"qvecHarmonic", 2, "Harmonic n of the Q-vector and event plane angle Psi_n: 2 -> elliptic, 3 -> triangular"};
8390
};
8491

8592
struct ConfCcdb : o2::framework::ConfigurableGroup {
@@ -240,6 +247,10 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
240247
mZorro.setBaseCCDBPath(confCcdb.triggerPath.value);
241248
}
242249

250+
// event shape
251+
mQvecDetector = static_cast<modes::QvecDetector>(config.qvecDetector.value);
252+
mQvecHarmonic = static_cast<modes::QvecHarmonic>(config.qvecHarmonic.value);
253+
243254
this->addSelection(kSel8, collisionSelectionNames.at(kSel8), config.sel8.value);
244255
this->addSelection(kNoSameBunchPileUp, collisionSelectionNames.at(kNoSameBunchPileUp), config.noSameBunchPileup.value);
245256
this->addSelection(kIsVertexItsTpc, collisionSelectionNames.at(kIsVertexItsTpc), config.isVertexItsTpc.value);
@@ -345,6 +356,35 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
345356
}
346357
[[nodiscard]] float getMultiplicity() const { return mMultiplicity; }
347358

359+
template <modes::System system, typename T>
360+
void setQvector(T const& col)
361+
{
362+
switch (mQvecDetector) {
363+
case modes::QvecDetector::kFT0C:
364+
mQvec = std::hypot(col.qvecFT0CReVec()[0], col.qvecFT0CImVec()[0]) * std::sqrt(col.sumAmplFT0C());
365+
break;
366+
case modes::QvecDetector::kFT0A:
367+
mQvec = std::hypot(col.qvecFT0AReVec()[0], col.qvecFT0AImVec()[0]) * std::sqrt(col.sumAmplFT0A());
368+
break;
369+
}
370+
}
371+
[[nodiscard]] float getQvector() const { return mQvec; }
372+
373+
template <modes::System system, typename T>
374+
void setEventPlane(T const& col)
375+
{
376+
float harmonic = static_cast<float>(mQvecHarmonic);
377+
switch (mQvecDetector) {
378+
case modes::QvecDetector::kFT0C:
379+
mEventPlane = RecoDecay::constrainAngle((std::atan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])) / harmonic, 0, 2); // constrain between 0 and pi
380+
break;
381+
case modes::QvecDetector::kFT0A:
382+
mEventPlane = RecoDecay::constrainAngle((std::atan2(col.qvecFT0AImVec()[0], col.qvecFT0AReVec()[0])) / harmonic, 0, 2); // constrain between 0 and pi
383+
break;
384+
}
385+
}
386+
[[nodiscard]] float getEventPlane() const { return mEventPlane; }
387+
348388
/// \brief Evaluate all pre-filters (kinematics, quality, RCT flags) for a collision candidate,
349389
/// filling one filter-histogram bin per bound plus the "All analyzed"/"All passed" summary bins.
350390
template <typename T>
@@ -415,6 +455,7 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
415455
this->evaluateObservable(kIsGoodZvtxFt0VsPv, static_cast<float>(col.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)));
416456
this->evaluateObservable(kNoCollInTimeRangeNarrow, static_cast<float>(col.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow)));
417457
this->evaluateObservable(kNoCollInTimeRangeStrict, static_cast<float>(col.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStrict)));
458+
this->evaluateObservable(kNoCollInTimeRangeStandard, static_cast<float>(col.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)));
418459
this->evaluateObservable(kNoCollInRofStrict, static_cast<float>(col.selection_bit(o2::aod::evsel::kNoCollInRofStrict)));
419460
this->evaluateObservable(kNoCollInRofStandard, static_cast<float>(col.selection_bit(o2::aod::evsel::kNoCollInRofStandard)));
420461
this->evaluateObservable(kNoHighMultCollInPrevRof, static_cast<float>(col.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)));
@@ -487,6 +528,11 @@ class CollisionSelection : public baseselection::BaseSelection<float, o2::analys
487528
float mSphericity = 0.f;
488529
float mCentrality = 0.f;
489530
float mMultiplicity = 0.f;
531+
float mQvec = 0.f;
532+
float mEventPlane = 0.f;
533+
534+
modes::QvecDetector mQvecDetector = modes::QvecDetector::kFT0C;
535+
modes::QvecHarmonic mQvecHarmonic = modes::QvecHarmonic::kN2;
490536

491537
// RCT flags
492538
mutable aod::rctsel::RCTFlagsChecker mRctFlagsChecker;
@@ -506,7 +552,7 @@ struct CollisionBuilderProducts : o2::framework::ProducesGroup {
506552
o2::framework::Produces<o2::aod::FColSphericities> producedSphericities;
507553
o2::framework::Produces<o2::aod::FColMults> producedMultiplicityEstimators;
508554
o2::framework::Produces<o2::aod::FColCents> producedCentralityEstimators;
509-
o2::framework::Produces<o2::aod::FColQns> producedQns;
555+
o2::framework::Produces<o2::aod::FColShapes> producedShapes;
510556
};
511557

512558
struct ConfCollisionTables : o2::framework::ConfigurableGroup {
@@ -518,7 +564,7 @@ struct ConfCollisionTables : o2::framework::ConfigurableGroup {
518564
o2::framework::Configurable<int> produceSphericities{"produceSphericities", -1, "Produce Sphericity (-1: auto; 0 off; 1 on)"};
519565
o2::framework::Configurable<int> produceMults{"produceMults", -1, "Produce Multiplicities (-1: auto; 0 off; 1 on)"};
520566
o2::framework::Configurable<int> produceCents{"produceCents", -1, "Produce Centralities (-1: auto; 0 off; 1 on)"};
521-
o2::framework::Configurable<int> produceQns{"produceQns", -1, "Produce Qn (-1: auto; 0 off; 1 on)"};
567+
o2::framework::Configurable<int> produceShapes{"produceShapes", -1, "Produce Event shape variables (-1: auto; 0 off; 1 on)"};
522568
};
523569

524570
template <auto& SelectionHistName, auto& FilterHistName>
@@ -544,15 +590,17 @@ class CollisionBuilder
544590
mProducedSphericities = utils::enableTable("FColSphericities_001", confTable.produceSphericities.value, initContext);
545591
mProducedMultiplicities = utils::enableTable("FColMults_001", confTable.produceMults.value, initContext);
546592
mProducedCentralities = utils::enableTable("FColCents_001", confTable.produceCents.value, initContext);
547-
mProduceQns = utils::enableTable("FColQnBins_001", confTable.produceQns.value, initContext);
593+
mProducedShapes = utils::enableTable("FColShapes_001", confTable.produceShapes.value, initContext);
548594

549595
if (mProducedCollisions && mProducedLiteCollisions) {
550596
LOG(fatal) << "FCols and FLiteCols are mutually exclusive -- enable only one. "
551597
<< "FLiteCols is meant to only replace FCols at the producer stage (for better compression in derived data); "
552598
<< "use the dedicated converter task to reconstruct FCols from FLiteCols downstream.";
553599
}
554600

555-
if (mProducedCollisions || mProducedLiteCollisions || mProducedCollisionMasks || mProducedPositions || mProducedSphericities || mProducedMultiplicities || mProducedCentralities) {
601+
if (mProducedCollisions || mProducedLiteCollisions || mProducedCollisionMasks ||
602+
mProducedPositions || mProducedSphericities || mProducedMultiplicities ||
603+
mProducedCentralities) {
556604
mFillAnyTable = true;
557605
} else {
558606
LOG(info) << "No tables configured, Selection object will not be configured...";
@@ -566,7 +614,7 @@ class CollisionBuilder
566614
}
567615

568616
template <modes::System system, typename T1, typename T2, typename T3, typename T4, typename T5>
569-
void initCollision(T1& bc, T2& col, T3& tracks, T4& ccdb, T5& histRegistry)
617+
void initCollision(T1 const& bc, T2 const& col, T3 const& tracks, T4& ccdb, T5& histRegistry)
570618
{
571619
if (mRunNumber != bc.runNumber()) {
572620
mRunNumber = bc.runNumber();
@@ -591,6 +639,11 @@ class CollisionBuilder
591639
mCollisionSelection.template setMultiplicity<system>(col);
592640
mCollisionSelection.template setCentrality<system>(col);
593641

642+
if constexpr (utils::HasQvectors<T2>) {
643+
mCollisionSelection.template setQvector<system>(col);
644+
mCollisionSelection.template setEventPlane<system>(col);
645+
}
646+
594647
std::vector<bool> triggerDecisions = mCollisionSelection.getTriggerDecisions(bc.globalBC());
595648

596649
mCollisionSelection.applySelections(col, triggerDecisions);
@@ -671,18 +724,18 @@ class CollisionBuilder
671724
col.ft0cOccupancyInTimeRange());
672725
}
673726

674-
// TODO: enable later for better QA
675-
// if (mProducedCentralities) {
676-
// collisionProducts.producedCentralityEstimators(
677-
// col.centFT0A(),
678-
// col.centFT0C());
679-
// }
680-
// PbPb specific columns
681-
// if constexpr (modes::isFlagSet(system, modes::System::kPbPb)) {
682-
// if (mProduceQns) {
683-
// collisionProducts.producedQns(utils::qn(col));
684-
// }
685-
// }
727+
if (mProducedCentralities) {
728+
collisionProducts.producedCentralityEstimators(
729+
col.centFT0A(),
730+
col.centFT0C(),
731+
col.centFT0M());
732+
}
733+
734+
if (mProducedShapes) {
735+
collisionProducts.producedShapes(
736+
mCollisionSelection.getQvector(),
737+
mCollisionSelection.getEventPlane());
738+
}
686739

687740
mCollisionAlreadyFilled = true;
688741
}
@@ -726,7 +779,7 @@ class CollisionBuilder
726779
bool mProducedSphericities = false;
727780
bool mProducedMultiplicities = false;
728781
bool mProducedCentralities = false;
729-
bool mProduceQns = false;
782+
bool mProducedShapes = false;
730783
};
731784

732785
struct CollisionBuilderDerivedToDerivedProducts : o2::framework::ProducesGroup {

‎PWGCF/Femto/Core/collisionHistManager.h‎

Lines changed: 36 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -16,14 +16,17 @@
1616
#ifndef PWGCF_FEMTO_CORE_COLLISIONHISTMANAGER_H_
1717
#define PWGCF_FEMTO_CORE_COLLISIONHISTMANAGER_H_
1818

19+
#include "PWGCF/Femto/Core/femtoUtils.h"
1920
#include "PWGCF/Femto/Core/histManager.h"
2021
#include "PWGCF/Femto/Core/modes.h"
2122

23+
#include <CommonConstants/MathConstants.h>
2224
#include <Framework/Configurable.h>
2325
#include <Framework/HistogramRegistry.h>
2426
#include <Framework/HistogramSpec.h>
2527

2628
#include <array>
29+
#include <cmath>
2730
#include <map>
2831
#include <string>
2932
#include <string_view>
@@ -48,6 +51,10 @@ enum ColHist {
4851
kCentVsMult,
4952
kCentVsSphericity,
5053
kMultVsSphericity,
54+
kFT0AvsFT0C,
55+
// event shape
56+
kQvector,
57+
kEventPlaneAngle,
5158
// mc
5259
kTruePosZ, // pure mc-truth, no reco collision (kMc without kReco)
5360
kTrueCent, // pure mc-truth, no reco collision (kMc without kReco)
@@ -80,6 +87,10 @@ constexpr std::array<histmanager::HistInfo<ColHist>, kColHistLast> HistTable = {
8087
{kCentVsMult, o2::framework::HistType::kTH2F, "hCentVsMult", "Centrality vs Multiplicity; Centrality (%); Multiplicity"},
8188
{kMultVsSphericity, o2::framework::HistType::kTH2F, "hMultVsSphericity", "Multiplicity vs Sphericity; Multiplicity; Sphericity"},
8289
{kCentVsSphericity, o2::framework::HistType::kTH2F, "hCentVsSphericity", "Centrality vs Sphericity; Centrality (%); Sphericity"},
90+
{kFT0AvsFT0C, o2::framework::HistType::kTH2F, "hFT0AvsFT0C", "FT0A centrality vs FT0C centrality; Centrality_{FT0A} (%); Centrality_{FT0C}"},
91+
// event shape
92+
{kQvector, o2::framework::HistType::kTH1F, "hQvector", "Q-vector; Q-vector; Entries"},
93+
{kEventPlaneAngle, o2::framework::HistType::kTH1F, "hEventPlaneAngle", "Event Plane angle; #Psi_{n}; Entries"},
8394
// mc
8495
{kTruePosZ, o2::framework::HistType::kTH1F, "hTruePosZ", "True vertex Z (mc-truth collision); V_{Z,True} (cm); Entries"},
8596
{kTrueCent, o2::framework::HistType::kTH1F, "hTrueCent", "True centrality (mc-truth collision); Centrality_{True} (%); Entries"},
@@ -94,7 +105,9 @@ constexpr std::array<histmanager::HistInfo<ColHist>, kColHistLast> HistTable = {
94105
{kPosZ, {(conf).vtxZ}}, \
95106
{kMult, {(conf).mult}}, \
96107
{kCent, {(conf).cent}}, \
97-
{kMagField, {(conf).magField}},
108+
{kMagField, {(conf).magField}}, \
109+
{kQvector, {(conf).qvector}}, \
110+
{kEventPlaneAngle, {(conf).eventPlaneAngle}},
98111

99112
// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
100113
#define COL_HIST_QA_MAP(confAnalysis, confQa) \
@@ -107,7 +120,8 @@ constexpr std::array<histmanager::HistInfo<ColHist>, kColHistLast> HistTable = {
107120
{kPoszVsCent, {(confAnalysis).vtxZ, (confAnalysis).cent}}, \
108121
{kCentVsMult, {(confAnalysis).cent, (confAnalysis).mult}}, \
109122
{kMultVsSphericity, {(confAnalysis).mult, (confQa).sphericity}}, \
110-
{kCentVsSphericity, {confBinningAnalysis.cent, (confQa).sphericity}},
123+
{kCentVsSphericity, {(confAnalysis).cent, (confQa).sphericity}}, \
124+
{kFT0AvsFT0C, {(confAnalysis).cent, (confAnalysis).cent}},
111125

112126
// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
113127
#define COL_HIST_MC_MAP(conf) \
@@ -164,6 +178,9 @@ struct ConfCollisionBinning : o2::framework::ConfigurableGroup {
164178
o2::framework::ConfigurableAxis mult{"mult", {200, 0, 200}, "Multiplicity binning"};
165179
o2::framework::ConfigurableAxis cent{"cent", {100, 0.0f, 100.0f}, "Centrality (multiplicity percentile) binning"};
166180
o2::framework::ConfigurableAxis magField{"magField", {11, -5.5, 5.5}, "Magnetic field binning"};
181+
o2::framework::Configurable<bool> plotEventShape{"plotEventShape", false, "Activate histograms for event shape (qvector, event plane angle)"};
182+
o2::framework::ConfigurableAxis qvector{"qvector", {100, 0.0f, 100.0f}, "Q-vector binning"};
183+
o2::framework::ConfigurableAxis eventPlaneAngle{"eventPlaneAngle", {720, 0, 1.f * o2::constants::math::TwoPI}, "Event plane angle binning"};
167184
};
168185

169186
struct ConfCollisionQaBinning : o2::framework::ConfigurableGroup {
@@ -184,9 +201,10 @@ class CollisionHistManager
184201
template <modes::Mode mode, typename T>
185202
void init(o2::framework::HistogramRegistry* registry,
186203
std::map<ColHist, std::vector<o2::framework::AxisSpec>> const& Specs,
187-
T const& /*ConfCollisionBinning*/)
204+
T const& ConfCollisionBinning)
188205
{
189206
mHistogramRegistry = registry;
207+
mPlotEventShape = ConfCollisionBinning.plotEventShape.value;
190208
if constexpr (isFlagSet(mode, modes::Mode::kReco)) {
191209
initAnalysis(Specs);
192210
}
@@ -257,6 +275,11 @@ class CollisionHistManager
257275
mHistogramRegistry->add(analysisDir + getHistNameV2(kMult, HistTable), getHistDesc(kMult, HistTable), getHistType(kMult, HistTable), {Specs.at(kMult)});
258276
mHistogramRegistry->add(analysisDir + getHistNameV2(kCent, HistTable), getHistDesc(kCent, HistTable), getHistType(kCent, HistTable), {Specs.at(kCent)});
259277
mHistogramRegistry->add(analysisDir + getHistNameV2(kMagField, HistTable), getHistDesc(kMagField, HistTable), getHistType(kMagField, HistTable), {Specs.at(kMagField)});
278+
279+
if (mPlotEventShape) {
280+
mHistogramRegistry->add(analysisDir + getHistNameV2(kQvector, HistTable), getHistDesc(kQvector, HistTable), getHistType(kQvector, HistTable), {Specs.at(kQvector)});
281+
mHistogramRegistry->add(analysisDir + getHistNameV2(kEventPlaneAngle, HistTable), getHistDesc(kEventPlaneAngle, HistTable), getHistType(kEventPlaneAngle, HistTable), {Specs.at(kEventPlaneAngle)});
282+
}
260283
}
261284

262285
void initQa(std::map<ColHist, std::vector<o2::framework::AxisSpec>> const& Specs)
@@ -273,6 +296,7 @@ class CollisionHistManager
273296
mHistogramRegistry->add(qaDir + getHistNameV2(kCentVsMult, HistTable), getHistDesc(kCentVsMult, HistTable), getHistType(kCentVsMult, HistTable), {Specs.at(kCentVsMult)});
274297
mHistogramRegistry->add(qaDir + getHistNameV2(kMultVsSphericity, HistTable), getHistDesc(kMultVsSphericity, HistTable), getHistType(kMultVsSphericity, HistTable), {Specs.at(kMultVsSphericity)});
275298
mHistogramRegistry->add(qaDir + getHistNameV2(kCentVsSphericity, HistTable), getHistDesc(kCentVsSphericity, HistTable), getHistType(kCentVsSphericity, HistTable), {Specs.at(kCentVsSphericity)});
299+
mHistogramRegistry->add(qaDir + getHistNameV2(kFT0AvsFT0C, HistTable), getHistDesc(kFT0AvsFT0C, HistTable), getHistType(kFT0AvsFT0C, HistTable), {Specs.at(kFT0AvsFT0C)});
276300
}
277301
}
278302

@@ -301,6 +325,13 @@ class CollisionHistManager
301325
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kMult, HistTable)), col.mult());
302326
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kCent, HistTable)), col.cent());
303327
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kMagField, HistTable)), col.magField());
328+
329+
if (mPlotEventShape) {
330+
if constexpr (utils::HasEventShape<T>) {
331+
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kQvector, HistTable)), col.qvec());
332+
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kEventPlaneAngle, HistTable)), col.eventPlaneAngle());
333+
}
334+
}
304335
}
305336

306337
template <typename T>
@@ -317,6 +348,7 @@ class CollisionHistManager
317348
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kCentVsMult, HistTable)), col.cent(), col.mult());
318349
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kMultVsSphericity, HistTable)), col.mult(), col.sphericity());
319350
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kCentVsSphericity, HistTable)), col.cent(), col.sphericity());
351+
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kFT0AvsFT0C, HistTable)), col.centFT0A(), col.centFT0C());
320352
}
321353
}
322354

@@ -343,6 +375,7 @@ class CollisionHistManager
343375
}
344376

345377
o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
378+
bool mPlotEventShape = false;
346379
bool mPlot2d = false;
347380
};
348381
} // namespace o2::analysis::femto::colhistmanager

0 commit comments

Comments
 (0)