diff --git a/PWGLF/DataModel/LFEventTopologyTables.h b/PWGLF/DataModel/LFEventTopologyTables.h new file mode 100644 index 00000000000..a42baa8c94d --- /dev/null +++ b/PWGLF/DataModel/LFEventTopologyTables.h @@ -0,0 +1,120 @@ +// 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. + +/// \file LFEventTopologyTables.h +/// \brief Per-collision forward/backward sub-event tables, reconstructed and generated level. +/// \author Cristian Andrei + +#ifndef PWGLF_DATAMODEL_LFEVENTTOPOLOGYTABLES_H_ +#define PWGLF_DATAMODEL_LFEVENTTOPOLOGYTABLES_H_ + +#include + +#include + +namespace o2::aod +{ + +namespace evshapecoex +{ +DECLARE_SOA_COLUMN(CollisionId, collisionId, int); //! source collision global index (DF-local; not an index column) +DECLARE_SOA_COLUMN(RunNumber, runNumber, int); //! run number +DECLARE_SOA_COLUMN(GlobalBC, globalBC, uint64_t); //! global bunch crossing +DECLARE_SOA_COLUMN(PosZ, posZ, float); //! primary-vertex z (cm) +DECLARE_SOA_COLUMN(MultFT0A, multFT0A, float); //! FT0-A amplitude +DECLARE_SOA_COLUMN(MultFT0C, multFT0C, float); //! FT0-C amplitude +DECLARE_SOA_COLUMN(NF, nF, uint16_t); //! track count, forward sub-event +DECLARE_SOA_COLUMN(NB, nB, uint16_t); //! track count, backward sub-event +DECLARE_SOA_COLUMN(NGap, nGap, uint16_t); //! track count, central gap +DECLARE_SOA_COLUMN(AF, aF, float); //! mean pT, forward (GeV/c); NaN if empty +DECLARE_SOA_COLUMN(AB, aB, float); //! mean pT, backward (GeV/c); NaN if empty +DECLARE_SOA_COLUMN(SumPt2F, sumPt2F, float); //! sum of pT^2, forward (GeV^2/c^2) +DECLARE_SOA_COLUMN(SumPt2B, sumPt2B, float); //! sum of pT^2, backward (GeV^2/c^2) +DECLARE_SOA_COLUMN(QxF, qxF, float); //! sum of cos(2 phi), forward +DECLARE_SOA_COLUMN(QyF, qyF, float); //! sum of sin(2 phi), forward +DECLARE_SOA_COLUMN(QxB, qxB, float); //! sum of cos(2 phi), backward +DECLARE_SOA_COLUMN(QyB, qyB, float); //! sum of sin(2 phi), backward +DECLARE_SOA_COLUMN(Qx4F, qx4F, float); //! sum of cos(4 phi), forward +DECLARE_SOA_COLUMN(Qy4F, qy4F, float); //! sum of sin(4 phi), forward +DECLARE_SOA_COLUMN(Qx4B, qx4B, float); //! sum of cos(4 phi), backward +DECLARE_SOA_COLUMN(Qy4B, qy4B, float); //! sum of sin(4 phi), backward +DECLARE_SOA_COLUMN(NPlus, nPlus, uint16_t); //! positive tracks, forward + backward +DECLARE_SOA_COLUMN(NMinus, nMinus, uint16_t); //! negative tracks, forward + backward +DECLARE_SOA_COLUMN(QaBits, qaBits, uint16_t); //! event-selection bits, order as QaBitOrder in eventShapeCoex.cxx +DECLARE_SOA_COLUMN(OccTracks, occTracks, int); //! track occupancy in time range +DECLARE_SOA_COLUMN(OccFT0C, occFT0C, float); //! FT0C occupancy in time range +DECLARE_SOA_COLUMN(NumContrib, numContrib, uint16_t); //! number of PV contributors +DECLARE_SOA_COLUMN(CollTimeRes, collTimeRes, float); //! collision time resolution (ns) +DECLARE_SOA_COLUMN(NPVC, nPVC, uint16_t); //! PV contributors among the selected tracks +} // namespace evshapecoex + +DECLARE_SOA_TABLE(EvShapeCoex, "AOD", "EVSHAPECOEX", //! per-collision forward/backward sub-events, reconstructed level + o2::soa::Index<>, + evshapecoex::CollisionId, evshapecoex::RunNumber, evshapecoex::GlobalBC, evshapecoex::PosZ, + evshapecoex::MultFT0A, evshapecoex::MultFT0C, + evshapecoex::NF, evshapecoex::NB, evshapecoex::NGap, + evshapecoex::AF, evshapecoex::AB, evshapecoex::SumPt2F, evshapecoex::SumPt2B, + evshapecoex::QxF, evshapecoex::QyF, evshapecoex::QxB, evshapecoex::QyB, + evshapecoex::Qx4F, evshapecoex::Qy4F, evshapecoex::Qx4B, evshapecoex::Qy4B, + evshapecoex::NPlus, evshapecoex::NMinus, + evshapecoex::QaBits, evshapecoex::OccTracks, evshapecoex::OccFT0C, + evshapecoex::NumContrib, evshapecoex::CollTimeRes, evshapecoex::NPVC); +using EvShapeCoexRow = EvShapeCoex::iterator; + +// Generator level: charged physical primaries, one row per McCollision (no event selection). +namespace evshapecoexgen +{ +DECLARE_SOA_COLUMN(McCollisionId, mcCollisionId, int); //! source McCollision global index (DF-local; not an index column) +DECLARE_SOA_COLUMN(PosZ, posZ, float); //! generated primary-vertex z (cm) +DECLARE_SOA_COLUMN(NFwdA, nFwdA, uint16_t); //! charged primaries in the FT0-A acceptance +DECLARE_SOA_COLUMN(NFwdC, nFwdC, uint16_t); //! charged primaries in the FT0-C acceptance +DECLARE_SOA_COLUMN(NF, nF, uint16_t); //! particle count, forward sub-event +DECLARE_SOA_COLUMN(NB, nB, uint16_t); //! particle count, backward sub-event +DECLARE_SOA_COLUMN(NGap, nGap, uint16_t); //! particle count, central gap +DECLARE_SOA_COLUMN(AF, aF, float); //! mean pT, forward (GeV/c); NaN if empty +DECLARE_SOA_COLUMN(AB, aB, float); //! mean pT, backward (GeV/c); NaN if empty +DECLARE_SOA_COLUMN(SumPt2F, sumPt2F, float); //! sum of pT^2, forward (GeV^2/c^2) +DECLARE_SOA_COLUMN(SumPt2B, sumPt2B, float); //! sum of pT^2, backward (GeV^2/c^2) +DECLARE_SOA_COLUMN(QxF, qxF, float); //! sum of cos(2 phi), forward +DECLARE_SOA_COLUMN(QyF, qyF, float); //! sum of sin(2 phi), forward +DECLARE_SOA_COLUMN(QxB, qxB, float); //! sum of cos(2 phi), backward +DECLARE_SOA_COLUMN(QyB, qyB, float); //! sum of sin(2 phi), backward +DECLARE_SOA_COLUMN(Qx4F, qx4F, float); //! sum of cos(4 phi), forward +DECLARE_SOA_COLUMN(Qy4F, qy4F, float); //! sum of sin(4 phi), forward +DECLARE_SOA_COLUMN(Qx4B, qx4B, float); //! sum of cos(4 phi), backward +DECLARE_SOA_COLUMN(Qy4B, qy4B, float); //! sum of sin(4 phi), backward +DECLARE_SOA_COLUMN(NPlus, nPlus, uint16_t); //! positive particles, forward + backward +DECLARE_SOA_COLUMN(NMinus, nMinus, uint16_t); //! negative particles, forward + backward +} // namespace evshapecoexgen + +DECLARE_SOA_TABLE(EvShapeCoexGen, "AOD", "EVSHAPECOEXGEN", //! per-McCollision forward/backward sub-events, generator level + o2::soa::Index<>, + evshapecoexgen::McCollisionId, evshapecoexgen::PosZ, + evshapecoexgen::NFwdA, evshapecoexgen::NFwdC, + evshapecoexgen::NF, evshapecoexgen::NB, evshapecoexgen::NGap, + evshapecoexgen::AF, evshapecoexgen::AB, evshapecoexgen::SumPt2F, evshapecoexgen::SumPt2B, + evshapecoexgen::QxF, evshapecoexgen::QyF, evshapecoexgen::QxB, evshapecoexgen::QyB, + evshapecoexgen::Qx4F, evshapecoexgen::Qy4F, evshapecoexgen::Qx4B, evshapecoexgen::Qy4B, + evshapecoexgen::NPlus, evshapecoexgen::NMinus); +using EvShapeCoexGenRow = EvShapeCoexGen::iterator; + +namespace evshapecoexmclabel +{ +DECLARE_SOA_INDEX_COLUMN_FULL(EvShapeCoexGen, evShapeCoexGen, int, EvShapeCoexGen, ""); //! matching EvShapeCoexGen row; negative if unlabelled +} // namespace evshapecoexmclabel + +DECLARE_SOA_TABLE(EvShapeCoexMcLabels, "AOD", "EVSHAPECOEXLBL", //! MC label, joinable with EvShapeCoex + o2::soa::Index<>, evshapecoexmclabel::EvShapeCoexGenId); +using EvShapeCoexMcLabel = EvShapeCoexMcLabels::iterator; + +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_LFEVENTTOPOLOGYTABLES_H_ diff --git a/PWGLF/TableProducer/CMakeLists.txt b/PWGLF/TableProducer/CMakeLists.txt index ec21b553563..745d984c81d 100644 --- a/PWGLF/TableProducer/CMakeLists.txt +++ b/PWGLF/TableProducer/CMakeLists.txt @@ -12,6 +12,7 @@ add_subdirectory(QC) add_subdirectory(Common) +add_subdirectory(GlobalEventProperties) add_subdirectory(Nuspex) add_subdirectory(Strangeness) add_subdirectory(Resonances) diff --git a/PWGLF/TableProducer/GlobalEventProperties/CMakeLists.txt b/PWGLF/TableProducer/GlobalEventProperties/CMakeLists.txt new file mode 100644 index 00000000000..6161ead0a4b --- /dev/null +++ b/PWGLF/TableProducer/GlobalEventProperties/CMakeLists.txt @@ -0,0 +1,15 @@ +# 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. + +o2physics_add_dpl_workflow(event-shape-coex + SOURCES eventShapeCoex.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/GlobalEventProperties/eventShapeCoex.cxx b/PWGLF/TableProducer/GlobalEventProperties/eventShapeCoex.cxx new file mode 100644 index 00000000000..1dfd790903e --- /dev/null +++ b/PWGLF/TableProducer/GlobalEventProperties/eventShapeCoex.cxx @@ -0,0 +1,247 @@ +// 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. + +/// \file eventShapeCoex.cxx +/// \brief Per-collision forward/backward sub-event quantities of the mid-rapidity tracks. +/// \author Cristian Andrei +/// +/// First O2Physics task of the AliPhysics PWGLF/SPECTRA/MultEvShape analyses. + +#include "PWGLF/DataModel/LFEventTopologyTables.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace +{ +/// Sums over one sub-event, shared by the reconstructed and generated levels. +struct SubEventSums { + uint32_t n = 0; + double sumPt = 0.; + double sumPt2 = 0.; + double qx2 = 0., qy2 = 0.; + double qx4 = 0., qy4 = 0.; + + void add(double pt, float phi) + { + // phi stays float, as stored in the AOD + const double c = std::cos(phi); + const double s = std::sin(phi); + const double cos2 = c * c - s * s; + const double sin2 = 2. * s * c; + ++n; + sumPt += pt; + sumPt2 += pt * pt; + qx2 += cos2; + qy2 += sin2; + qx4 += cos2 * cos2 - sin2 * sin2; + qy4 += 2. * sin2 * cos2; + } + + [[nodiscard]] float meanPt() const { return (n > 0) ? static_cast(sumPt / n) : std::nanf(""); } +}; +} // namespace + +struct EventShapeCoex { + Produces coex; + Produces coexGen; + Produces coexMcLabels; + + Service pdg{}; + + Configurable vtxZCut{"vtxZCut", 10.0f, "Max |PV z| (cm)"}; + Configurable ptMin{"ptMin", 0.15f, "Minimum track pT (GeV/c)"}; + Configurable ptMax{"ptMax", 2.0f, "Maximum track pT (GeV/c)"}; + Configurable etaMax{"etaMax", 0.8f, "Max |eta| (mid-rapidity measurement region)"}; + Configurable gapHalf{"gapHalf", 0.2f, "Central-gap half-width: |eta| <= gapHalf excluded from F/B"}; + Configurable etaFwdAMin{"etaFwdAMin", 3.5f, "FT0-A proxy: lower eta edge (generator level)"}; + Configurable etaFwdAMax{"etaFwdAMax", 4.9f, "FT0-A proxy: upper eta edge (generator level)"}; + Configurable etaFwdCMin{"etaFwdCMin", -3.3f, "FT0-C proxy: lower eta edge (generator level)"}; + Configurable etaFwdCMax{"etaFwdCMax", -2.1f, "FT0-C proxy: upper eta edge (generator level)"}; + + Filter trackFilter = (aod::track::pt > ptMin) && (aod::track::pt < ptMax) && (nabs(aod::track::eta) < etaMax) && requireGlobalTrackInFilter(); + + using MyCollisions = soa::Join; + using MyCollisionsMc = soa::Join; + using MyTracks = soa::Filtered>; + + void init(InitContext const&) + { + if (doprocessReco && doprocessRecoMC) { + LOGF(fatal, "processReco and processRecoMC both fill EvShapeCoex; enable only one of them."); + } + if (doprocessRecoMC && !doprocessMC) { + LOGF(fatal, "processRecoMC links to EvShapeCoexGen rows and requires processMC in the same job."); + } + } + + /// Fills one EvShapeCoex row; returns false if the collision is rejected. + template + bool fillReco(TCollision const& coll, MyTracks const& tracks) + { + if (!coll.sel8() || std::abs(coll.posZ()) > vtxZCut) { + return false; + } + + const float g = gapHalf; + SubEventSums fwd, bwd; + uint32_t nGap = 0, nPlus = 0, nMinus = 0, nPVC = 0; + for (const auto& track : tracks) { + if (track.isPVContributor()) { + ++nPVC; + } + const float eta = track.eta(); + if (eta > g) { + fwd.add(track.pt(), track.phi()); + } else if (eta < -g) { + bwd.add(track.pt(), track.phi()); + } else { + ++nGap; + continue; + } + if (track.sign() > 0) { + ++nPlus; + } else if (track.sign() < 0) { + ++nMinus; + } + } + + // bit i of qaBits = QaBitOrder[i] + static constexpr int NQaBits = 12; + static constexpr std::array QaBitOrder = { + o2::aod::evsel::kNoSameBunchPileup, o2::aod::evsel::kIsGoodZvtxFT0vsPV, + o2::aod::evsel::kIsVertexITSTPC, o2::aod::evsel::kIsVertexTOFmatched, + o2::aod::evsel::kNoCollInTimeRangeNarrow, o2::aod::evsel::kNoCollInTimeRangeStrict, + o2::aod::evsel::kNoCollInTimeRangeStandard, o2::aod::evsel::kNoCollInRofStrict, + o2::aod::evsel::kNoCollInRofStandard, o2::aod::evsel::kNoHighMultCollInPrevRof, + o2::aod::evsel::kNoITSROFrameBorder, o2::aod::evsel::kNoTimeFrameBorder}; + uint16_t qaBits = 0; + for (int i = 0; i < NQaBits; ++i) { + if (coll.selection_bit(QaBitOrder[i])) { + qaBits |= static_cast(1u << i); + } + } + + const auto& bc = coll.template bc_as(); + + coex(static_cast(coll.globalIndex()), bc.runNumber(), bc.globalBC(), coll.posZ(), + coll.multFT0A(), coll.multFT0C(), + static_cast(fwd.n), static_cast(bwd.n), static_cast(nGap), + fwd.meanPt(), bwd.meanPt(), static_cast(fwd.sumPt2), static_cast(bwd.sumPt2), + static_cast(fwd.qx2), static_cast(fwd.qy2), + static_cast(bwd.qx2), static_cast(bwd.qy2), + static_cast(fwd.qx4), static_cast(fwd.qy4), + static_cast(bwd.qx4), static_cast(bwd.qy4), + static_cast(nPlus), static_cast(nMinus), + qaBits, coll.trackOccupancyInTimeRange(), coll.ft0cOccupancyInTimeRange(), + coll.numContrib(), coll.collisionTimeRes(), static_cast(nPVC)); + return true; + } + + void processReco(MyCollisions::iterator const& coll, MyTracks const& tracks, aod::BCs const&) + { + fillReco(coll, tracks); + } + PROCESS_SWITCH(EventShapeCoex, processReco, "Reconstructed-level reduction (EvShapeCoex)", true); + + void processRecoMC(MyCollisionsMc::iterator const& coll, MyTracks const& tracks, aod::BCs const&) + { + if (fillReco(coll, tracks)) { + // EvShapeCoexGen has one row per McCollision, in order: the McCollision index is its row index + coexMcLabels(coll.mcCollisionId()); + } + } + PROCESS_SWITCH(EventShapeCoex, processRecoMC, "Reconstructed-level reduction with MC labels (EvShapeCoex + EvShapeCoexMcLabels)", false); + + void processMC(aod::McCollision const& mcCollision, aod::McParticles const& particles) + { + static constexpr double ChargeTolerance = 1.e-3; + + const float g = gapHalf; + const float em = etaMax; + SubEventSums fwd, bwd; + uint32_t nFwdA = 0, nFwdC = 0, nGap = 0, nPlus = 0, nMinus = 0; + for (const auto& particle : particles) { + if (!particle.isPhysicalPrimary()) { + continue; + } + const auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle == nullptr) { + continue; // unknown PDG code, charge undetermined + } + const double charge = pdgParticle->Charge(); + if (std::abs(charge) < ChargeTolerance) { + continue; + } + const float eta = particle.eta(); + // forward counts: no pT window + if (eta > etaFwdAMin && eta < etaFwdAMax) { + ++nFwdA; + } + if (eta > etaFwdCMin && eta < etaFwdCMax) { + ++nFwdC; + } + const float pt = particle.pt(); + if (pt <= ptMin || pt >= ptMax || std::abs(eta) >= em) { + continue; + } + if (eta > g) { + fwd.add(pt, particle.phi()); + } else if (eta < -g) { + bwd.add(pt, particle.phi()); + } else { + ++nGap; + continue; + } + if (charge > 0.) { + ++nPlus; + } else { + ++nMinus; + } + } + + coexGen(static_cast(mcCollision.globalIndex()), mcCollision.posZ(), + static_cast(nFwdA), static_cast(nFwdC), + static_cast(fwd.n), static_cast(bwd.n), static_cast(nGap), + fwd.meanPt(), bwd.meanPt(), static_cast(fwd.sumPt2), static_cast(bwd.sumPt2), + static_cast(fwd.qx2), static_cast(fwd.qy2), + static_cast(bwd.qx2), static_cast(bwd.qy2), + static_cast(fwd.qx4), static_cast(fwd.qy4), + static_cast(bwd.qx4), static_cast(bwd.qy4), + static_cast(nPlus), static_cast(nMinus)); + } + PROCESS_SWITCH(EventShapeCoex, processMC, "Generator-level reduction (EvShapeCoexGen)", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}