|
| 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 | +#include <cmath> |
| 13 | +#include <vector> |
| 14 | + |
| 15 | +#include "Framework/runDataProcessing.h" |
| 16 | +#include "Framework/AnalysisTask.h" |
| 17 | +#include "Framework/AnalysisDataModel.h" |
| 18 | +#include "Framework/ASoAHelpers.h" |
| 19 | +#include "Framework/HistogramRegistry.h" |
| 20 | +#include "Framework/Configurable.h" |
| 21 | +#include "Framework/InitContext.h" |
| 22 | +#include "Framework/O2DatabasePDGPlugin.h" |
| 23 | + |
| 24 | +#include "ReconstructionDataFormats/Track.h" |
| 25 | +#include "Common/DataModel/Multiplicity.h" |
| 26 | +#include "Common/DataModel/EventSelection.h" |
| 27 | +#include "Common/DataModel/TrackSelectionTables.h" |
| 28 | +#include "Common/Core/TrackSelection.h" |
| 29 | +#include "Common/Core/TrackSelectionDefaults.h" |
| 30 | + |
| 31 | +using namespace o2; |
| 32 | +using namespace o2::framework; |
| 33 | +using namespace o2::framework::expressions; |
| 34 | + |
| 35 | +struct flattenicityTask { |
| 36 | + |
| 37 | + // --- Flattenicity constants --- |
| 38 | + static constexpr int N_CH_A = 96; |
| 39 | + static constexpr int N_CH_C = 112; |
| 40 | + static constexpr int N_CELL = N_CH_A + N_CH_C; |
| 41 | + static constexpr int N_PHI_SECTORS = 8; |
| 42 | + static constexpr int N_ETA_A = N_CH_A / N_PHI_SECTORS; |
| 43 | + static constexpr int N_ETA_C = N_CH_C / N_PHI_SECTORS; |
| 44 | + |
| 45 | + // FT0 acceptance |
| 46 | + static constexpr float FT0A_ETA_MIN = 3.5; |
| 47 | + static constexpr float FT0A_ETA_MAX = 4.9; |
| 48 | + static constexpr float FT0C_ETA_MIN = -3.3; |
| 49 | + static constexpr float FT0C_ETA_MAX = -2.1; |
| 50 | + |
| 51 | + // --- Configurables --- |
| 52 | + Configurable<float> cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, "Eta range for tracks"}; |
| 53 | + Configurable<float> cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.15f, "Minimum pT"}; |
| 54 | + |
| 55 | + Configurable<bool> isRun3{"isRun3", true, "is Run3 dataset"}; |
| 56 | + Configurable<bool> timeEvsel{"timeEvsel", true, "TPC Time frame boundary cut"}; |
| 57 | + Configurable<bool> piluprejection{"piluprejection", true, "Pileup rejection"}; |
| 58 | + Configurable<bool> goodzvertex{"goodzvertex", true, "Good Z vertex"}; |
| 59 | + |
| 60 | + // --- Track selection --- |
| 61 | + TrackSelection mySelectionPrim; |
| 62 | + |
| 63 | + // --- Histograms --- |
| 64 | + HistogramRegistry registry; |
| 65 | + |
| 66 | + // --- Init --- |
| 67 | + void init(InitContext const&) |
| 68 | + { |
| 69 | + // Initialize track selection |
| 70 | + mySelectionPrim = myTrackSelectionPrim(); |
| 71 | + |
| 72 | + // Define histograms |
| 73 | + AxisSpec flatBins = {40, 0.0, 1.0, "#rho"}; |
| 74 | + AxisSpec nchBins = {100, -0.5, 99.5, "N_{ch}"}; |
| 75 | + |
| 76 | + registry.add("hFlattenicityTruth", "Truth flattenicity; 1-#rho; Events", |
| 77 | + HistType::kTH1D, {flatBins}); |
| 78 | + registry.add("hFlattenicityReco", "Reco flattenicity; 1-#rho; Events", |
| 79 | + HistType::kTH1D, {flatBins}); |
| 80 | + registry.add("hFlattenicityCorrelation", "Truth vs Reco; 1-#rho_{truth}; 1-#rho_{reco}", |
| 81 | + HistType::kTH2D, {flatBins, flatBins}); |
| 82 | + registry.add("hNch", "Nch distribution; N_{ch}; Events", |
| 83 | + HistType::kTH1D, {nchBins}); |
| 84 | + } |
| 85 | + |
| 86 | + // --- Track selection function --- |
| 87 | + TrackSelection myTrackSelectionPrim() |
| 88 | + { |
| 89 | + TrackSelection selectedTracks; |
| 90 | + selectedTracks.SetPtRange(0.1f, 1e10f); |
| 91 | + selectedTracks.SetEtaRange(-0.8f, 0.8f); |
| 92 | + selectedTracks.SetRequireITSRefit(true); |
| 93 | + selectedTracks.SetRequireTPCRefit(true); |
| 94 | + selectedTracks.SetMinNCrossedRowsTPC(70); |
| 95 | + selectedTracks.SetMinNCrossedRowsOverFindableClustersTPC(0.4f); |
| 96 | + selectedTracks.SetMaxChi2PerClusterTPC(4.f); |
| 97 | + selectedTracks.SetRequireHitsInITSLayers(1, {0, 1}); |
| 98 | + selectedTracks.SetMaxChi2PerClusterITS(36.f); |
| 99 | + selectedTracks.SetMaxDcaXYPtDep([](float pt) { return 0.0105f + 0.0350f / pow(pt, 1.1f); }); |
| 100 | + selectedTracks.SetMaxDcaZ(2.f); |
| 101 | + return selectedTracks; |
| 102 | + } |
| 103 | + |
| 104 | + // --- Flattenicity calculation --- |
| 105 | + float calculateFlattenicity(const std::vector<float>& counts) |
| 106 | + { |
| 107 | + if (counts.size() != N_CELL) return -1; |
| 108 | + |
| 109 | + float total = 0; |
| 110 | + for (auto c : counts) total += c; |
| 111 | + if (total <= 0) return -1; |
| 112 | + |
| 113 | + float mean = total / N_CELL; |
| 114 | + if (mean <= 0) return -1; |
| 115 | + |
| 116 | + float sumSq = 0; |
| 117 | + for (auto c : counts) { |
| 118 | + sumSq += (c - mean) * (c - mean); |
| 119 | + } |
| 120 | + |
| 121 | + float rho = std::sqrt(sumSq / (N_CELL * N_CELL)) / mean; |
| 122 | + return 1.0f - rho; |
| 123 | + } |
| 124 | + |
| 125 | + // --- Process Data --- |
| 126 | + void processData(aod::Collision const& collision, |
| 127 | + soa::Filtered<aod::Tracks> const& tracks, |
| 128 | + aod::FT0s const& ft0s) |
| 129 | + { |
| 130 | + // Event selection (Paola/Jesus) |
| 131 | + if (!collision.sel8()) return; |
| 132 | + if (std::abs(collision.posZ()) >= 10.0f) return; |
| 133 | + |
| 134 | + // Track loop for Nch |
| 135 | + int nch = 0; |
| 136 | + for (auto& track : tracks) { |
| 137 | + if (!mySelectionPrim.IsSelected(track)) continue; |
| 138 | + nch++; |
| 139 | + } |
| 140 | + registry.fill(HIST("hNch"), nch); |
| 141 | + |
| 142 | + // FT0 flattenicity |
| 143 | + auto ft0 = collision.ft0(); |
| 144 | + if (ft0.hasAmplitudeA() && ft0.hasAmplitudeC()) { |
| 145 | + auto ampA = ft0.amplitudeA(); |
| 146 | + auto ampC = ft0.amplitudeC(); |
| 147 | + |
| 148 | + std::vector<float> counts(N_CELL, 0.0f); |
| 149 | + for (int i = 0; i < ampA.size() && i < N_CH_A; ++i) { |
| 150 | + counts[i] = ampA[i]; |
| 151 | + } |
| 152 | + for (int i = 0; i < ampC.size() && i < N_CH_C; ++i) { |
| 153 | + counts[N_CH_A + i] = ampC[i]; |
| 154 | + } |
| 155 | + |
| 156 | + float flat = calculateFlattenicity(counts); |
| 157 | + if (flat >= 0) { |
| 158 | + registry.fill(HIST("hFlattenicityReco"), flat); |
| 159 | + } |
| 160 | + } |
| 161 | + } |
| 162 | + PROCESS_SWITCH(flattenicityTask, processData, "Process data", true); |
| 163 | + |
| 164 | + // --- Process MC --- |
| 165 | + void processMC(aod::McCollision const& mcCollision, |
| 166 | + aod::McParticles const& particles, |
| 167 | + soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions>> const& collisions, |
| 168 | + aod::FT0s const& ft0s) |
| 169 | + { |
| 170 | + // MC processing - to be implemented |
| 171 | + } |
| 172 | + PROCESS_SWITCH(flattenicityTask, processMC, "Process MC", false); |
| 173 | +}; |
| 174 | + |
| 175 | +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
| 176 | +{ |
| 177 | + WorkflowSpec workflow{}; |
| 178 | + workflow.push_back(adaptAnalysisTask<flattenicityTask>(cfgc)); |
| 179 | + return workflow; |
| 180 | +} |
0 commit comments