|
| 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 | +/// \file eventShapeCoex.cxx |
| 13 | +/// \brief Per-collision forward/backward sub-event quantities of the mid-rapidity tracks. |
| 14 | +/// \author Cristian Andrei <Cristian.Andrei@cern.ch> |
| 15 | +/// |
| 16 | +/// First O2Physics task of the AliPhysics PWGLF/SPECTRA/MultEvShape analyses. |
| 17 | + |
| 18 | +#include "PWGLF/DataModel/LFEventTopologyTables.h" |
| 19 | + |
| 20 | +#include "Common/CCDB/EventSelectionParams.h" |
| 21 | +#include "Common/DataModel/EventSelection.h" |
| 22 | +#include "Common/DataModel/Multiplicity.h" |
| 23 | +#include "Common/DataModel/TrackSelectionTables.h" |
| 24 | + |
| 25 | +#include <Framework/ASoA.h> |
| 26 | +#include <Framework/AnalysisDataModel.h> |
| 27 | +#include <Framework/AnalysisTask.h> |
| 28 | +#include <Framework/Configurable.h> |
| 29 | +#include <Framework/InitContext.h> |
| 30 | +#include <Framework/Logger.h> |
| 31 | +#include <Framework/O2DatabasePDGPlugin.h> |
| 32 | +#include <Framework/runDataProcessing.h> |
| 33 | + |
| 34 | +#include <cmath> |
| 35 | +#include <cstdint> |
| 36 | + |
| 37 | +using namespace o2; |
| 38 | +using namespace o2::framework; |
| 39 | +using namespace o2::framework::expressions; |
| 40 | + |
| 41 | +namespace |
| 42 | +{ |
| 43 | +/// Sums over one sub-event, shared by the reconstructed and generated levels. |
| 44 | +struct SubEventSums { |
| 45 | + uint32_t n = 0; |
| 46 | + double sumPt = 0.; |
| 47 | + double sumPt2 = 0.; |
| 48 | + double qx2 = 0., qy2 = 0.; |
| 49 | + double qx4 = 0., qy4 = 0.; |
| 50 | + |
| 51 | + void add(double pt, float phi) |
| 52 | + { |
| 53 | + // phi stays float, as stored in the AOD |
| 54 | + const double c = std::cos(phi); |
| 55 | + const double s = std::sin(phi); |
| 56 | + const double cos2 = c * c - s * s; |
| 57 | + const double sin2 = 2. * s * c; |
| 58 | + ++n; |
| 59 | + sumPt += pt; |
| 60 | + sumPt2 += pt * pt; |
| 61 | + qx2 += cos2; |
| 62 | + qy2 += sin2; |
| 63 | + qx4 += cos2 * cos2 - sin2 * sin2; |
| 64 | + qy4 += 2. * sin2 * cos2; |
| 65 | + } |
| 66 | + |
| 67 | + float meanPt() const { return (n > 0) ? static_cast<float>(sumPt / n) : std::nanf(""); } |
| 68 | +}; |
| 69 | +} // namespace |
| 70 | + |
| 71 | +struct EventShapeCoex { |
| 72 | + Produces<aod::EvShapeCoex> coex; |
| 73 | + Produces<aod::EvShapeCoexGen> coexGen; |
| 74 | + Produces<aod::EvShapeCoexMcLabels> coexMcLabels; |
| 75 | + |
| 76 | + Service<o2::framework::O2DatabasePDG> pdg; |
| 77 | + |
| 78 | + Configurable<float> vtxZCut{"vtxZCut", 10.0f, "Max |PV z| (cm)"}; |
| 79 | + Configurable<float> ptMin{"ptMin", 0.15f, "Minimum track pT (GeV/c)"}; |
| 80 | + Configurable<float> ptMax{"ptMax", 2.0f, "Maximum track pT (GeV/c)"}; |
| 81 | + Configurable<float> etaMax{"etaMax", 0.8f, "Max |eta| (mid-rapidity measurement region)"}; |
| 82 | + Configurable<float> gapHalf{"gapHalf", 0.2f, "Central-gap half-width: |eta| <= gapHalf excluded from F/B"}; |
| 83 | + Configurable<float> etaFwdAMin{"etaFwdAMin", 3.5f, "FT0-A proxy: lower eta edge (generator level)"}; |
| 84 | + Configurable<float> etaFwdAMax{"etaFwdAMax", 4.9f, "FT0-A proxy: upper eta edge (generator level)"}; |
| 85 | + Configurable<float> etaFwdCMin{"etaFwdCMin", -3.3f, "FT0-C proxy: lower eta edge (generator level)"}; |
| 86 | + Configurable<float> etaFwdCMax{"etaFwdCMax", -2.1f, "FT0-C proxy: upper eta edge (generator level)"}; |
| 87 | + |
| 88 | + Filter trackFilter = (aod::track::pt > ptMin) && (aod::track::pt < ptMax) && (nabs(aod::track::eta) < etaMax) && requireGlobalTrackInFilter(); |
| 89 | + |
| 90 | + using MyCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults>; |
| 91 | + using MyCollisionsMc = soa::Join<MyCollisions, aod::McCollisionLabels>; |
| 92 | + using MyTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection>>; |
| 93 | + |
| 94 | + void init(InitContext const&) |
| 95 | + { |
| 96 | + if (doprocessReco && doprocessRecoMC) { |
| 97 | + LOGF(fatal, "processReco and processRecoMC both fill EvShapeCoex; enable only one of them."); |
| 98 | + } |
| 99 | + if (doprocessRecoMC && !doprocessMC) { |
| 100 | + LOGF(fatal, "processRecoMC links to EvShapeCoexGen rows and requires processMC in the same job."); |
| 101 | + } |
| 102 | + } |
| 103 | + |
| 104 | + /// Fills one EvShapeCoex row; returns false if the collision is rejected. |
| 105 | + template <typename TCollision> |
| 106 | + bool fillReco(TCollision const& coll, MyTracks const& tracks) |
| 107 | + { |
| 108 | + if (!coll.sel8() || std::abs(coll.posZ()) > vtxZCut) { |
| 109 | + return false; |
| 110 | + } |
| 111 | + |
| 112 | + const float g = gapHalf; |
| 113 | + SubEventSums fwd, bwd; |
| 114 | + uint32_t nGap = 0, nPlus = 0, nMinus = 0, nPVC = 0; |
| 115 | + for (const auto& track : tracks) { |
| 116 | + if (track.isPVContributor()) { |
| 117 | + ++nPVC; |
| 118 | + } |
| 119 | + const float eta = track.eta(); |
| 120 | + if (eta > g) { |
| 121 | + fwd.add(track.pt(), track.phi()); |
| 122 | + } else if (eta < -g) { |
| 123 | + bwd.add(track.pt(), track.phi()); |
| 124 | + } else { |
| 125 | + ++nGap; |
| 126 | + continue; |
| 127 | + } |
| 128 | + if (track.sign() > 0) { |
| 129 | + ++nPlus; |
| 130 | + } else if (track.sign() < 0) { |
| 131 | + ++nMinus; |
| 132 | + } |
| 133 | + } |
| 134 | + |
| 135 | + // bit i of qaBits = QaBitOrder[i] |
| 136 | + static constexpr int NQaBits = 12; |
| 137 | + static constexpr o2::aod::evsel::EventSelectionFlags QaBitOrder[NQaBits] = { |
| 138 | + o2::aod::evsel::kNoSameBunchPileup, o2::aod::evsel::kIsGoodZvtxFT0vsPV, |
| 139 | + o2::aod::evsel::kIsVertexITSTPC, o2::aod::evsel::kIsVertexTOFmatched, |
| 140 | + o2::aod::evsel::kNoCollInTimeRangeNarrow, o2::aod::evsel::kNoCollInTimeRangeStrict, |
| 141 | + o2::aod::evsel::kNoCollInTimeRangeStandard, o2::aod::evsel::kNoCollInRofStrict, |
| 142 | + o2::aod::evsel::kNoCollInRofStandard, o2::aod::evsel::kNoHighMultCollInPrevRof, |
| 143 | + o2::aod::evsel::kNoITSROFrameBorder, o2::aod::evsel::kNoTimeFrameBorder}; |
| 144 | + uint16_t qaBits = 0; |
| 145 | + for (int i = 0; i < NQaBits; ++i) { |
| 146 | + if (coll.selection_bit(QaBitOrder[i])) { |
| 147 | + qaBits |= static_cast<uint16_t>(1u << i); |
| 148 | + } |
| 149 | + } |
| 150 | + |
| 151 | + const auto& bc = coll.template bc_as<aod::BCs>(); |
| 152 | + |
| 153 | + coex(static_cast<int>(coll.globalIndex()), bc.runNumber(), bc.globalBC(), coll.posZ(), |
| 154 | + coll.multFT0A(), coll.multFT0C(), |
| 155 | + static_cast<uint16_t>(fwd.n), static_cast<uint16_t>(bwd.n), static_cast<uint16_t>(nGap), |
| 156 | + fwd.meanPt(), bwd.meanPt(), static_cast<float>(fwd.sumPt2), static_cast<float>(bwd.sumPt2), |
| 157 | + static_cast<float>(fwd.qx2), static_cast<float>(fwd.qy2), |
| 158 | + static_cast<float>(bwd.qx2), static_cast<float>(bwd.qy2), |
| 159 | + static_cast<float>(fwd.qx4), static_cast<float>(fwd.qy4), |
| 160 | + static_cast<float>(bwd.qx4), static_cast<float>(bwd.qy4), |
| 161 | + static_cast<uint16_t>(nPlus), static_cast<uint16_t>(nMinus), |
| 162 | + qaBits, coll.trackOccupancyInTimeRange(), coll.ft0cOccupancyInTimeRange(), |
| 163 | + coll.numContrib(), coll.collisionTimeRes(), static_cast<uint16_t>(nPVC)); |
| 164 | + return true; |
| 165 | + } |
| 166 | + |
| 167 | + void processReco(MyCollisions::iterator const& coll, MyTracks const& tracks, aod::BCs const&) |
| 168 | + { |
| 169 | + fillReco(coll, tracks); |
| 170 | + } |
| 171 | + PROCESS_SWITCH(EventShapeCoex, processReco, "Reconstructed-level reduction (EvShapeCoex)", true); |
| 172 | + |
| 173 | + void processRecoMC(MyCollisionsMc::iterator const& coll, MyTracks const& tracks, aod::BCs const&) |
| 174 | + { |
| 175 | + if (fillReco(coll, tracks)) { |
| 176 | + // EvShapeCoexGen has one row per McCollision, in order: the McCollision index is its row index |
| 177 | + coexMcLabels(coll.mcCollisionId()); |
| 178 | + } |
| 179 | + } |
| 180 | + PROCESS_SWITCH(EventShapeCoex, processRecoMC, "Reconstructed-level reduction with MC labels (EvShapeCoex + EvShapeCoexMcLabels)", false); |
| 181 | + |
| 182 | + void processMC(aod::McCollision const& mcCollision, aod::McParticles const& particles) |
| 183 | + { |
| 184 | + static constexpr double ChargeTolerance = 1.e-3; |
| 185 | + |
| 186 | + const float g = gapHalf; |
| 187 | + const float em = etaMax; |
| 188 | + SubEventSums fwd, bwd; |
| 189 | + uint32_t nFwdA = 0, nFwdC = 0, nGap = 0, nPlus = 0, nMinus = 0; |
| 190 | + for (const auto& particle : particles) { |
| 191 | + if (!particle.isPhysicalPrimary()) { |
| 192 | + continue; |
| 193 | + } |
| 194 | + const auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); |
| 195 | + if (pdgParticle == nullptr) { |
| 196 | + continue; // unknown PDG code, charge undetermined |
| 197 | + } |
| 198 | + const double charge = pdgParticle->Charge(); |
| 199 | + if (std::abs(charge) < ChargeTolerance) { |
| 200 | + continue; |
| 201 | + } |
| 202 | + const float eta = particle.eta(); |
| 203 | + // forward counts: no pT window |
| 204 | + if (eta > etaFwdAMin && eta < etaFwdAMax) { |
| 205 | + ++nFwdA; |
| 206 | + } |
| 207 | + if (eta > etaFwdCMin && eta < etaFwdCMax) { |
| 208 | + ++nFwdC; |
| 209 | + } |
| 210 | + const float pt = particle.pt(); |
| 211 | + if (pt <= ptMin || pt >= ptMax || std::abs(eta) >= em) { |
| 212 | + continue; |
| 213 | + } |
| 214 | + if (eta > g) { |
| 215 | + fwd.add(pt, particle.phi()); |
| 216 | + } else if (eta < -g) { |
| 217 | + bwd.add(pt, particle.phi()); |
| 218 | + } else { |
| 219 | + ++nGap; |
| 220 | + continue; |
| 221 | + } |
| 222 | + if (charge > 0.) { |
| 223 | + ++nPlus; |
| 224 | + } else { |
| 225 | + ++nMinus; |
| 226 | + } |
| 227 | + } |
| 228 | + |
| 229 | + coexGen(static_cast<int>(mcCollision.globalIndex()), mcCollision.posZ(), |
| 230 | + static_cast<uint16_t>(nFwdA), static_cast<uint16_t>(nFwdC), |
| 231 | + static_cast<uint16_t>(fwd.n), static_cast<uint16_t>(bwd.n), static_cast<uint16_t>(nGap), |
| 232 | + fwd.meanPt(), bwd.meanPt(), static_cast<float>(fwd.sumPt2), static_cast<float>(bwd.sumPt2), |
| 233 | + static_cast<float>(fwd.qx2), static_cast<float>(fwd.qy2), |
| 234 | + static_cast<float>(bwd.qx2), static_cast<float>(bwd.qy2), |
| 235 | + static_cast<float>(fwd.qx4), static_cast<float>(fwd.qy4), |
| 236 | + static_cast<float>(bwd.qx4), static_cast<float>(bwd.qy4), |
| 237 | + static_cast<uint16_t>(nPlus), static_cast<uint16_t>(nMinus)); |
| 238 | + } |
| 239 | + PROCESS_SWITCH(EventShapeCoex, processMC, "Generator-level reduction (EvShapeCoexGen)", false); |
| 240 | +}; |
| 241 | + |
| 242 | +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
| 243 | +{ |
| 244 | + return WorkflowSpec{adaptAnalysisTask<EventShapeCoex>(cfgc)}; |
| 245 | +} |
0 commit comments