|
| 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 | +// O2 includes |
| 12 | + |
| 13 | +#include "ReconstructionDataFormats/Track.h" |
| 14 | +#include "Framework/AnalysisTask.h" |
| 15 | +#include "Framework/AnalysisDataModel.h" |
| 16 | +#include "Framework/ASoAHelpers.h" |
| 17 | +#include "AnalysisCore/MC.h" |
| 18 | +#include "AnalysisDataModel/PID/PIDResponse.h" |
| 19 | +#include "AnalysisDataModel/TrackSelectionTables.h" |
| 20 | + |
| 21 | +#include "AnalysisDataModel/EventSelection.h" |
| 22 | +#include "AnalysisDataModel/TrackSelectionTables.h" |
| 23 | +#include "AnalysisDataModel/Centrality.h" |
| 24 | + |
| 25 | +#include "Framework/HistogramRegistry.h" |
| 26 | + |
| 27 | +#include <TLorentzVector.h> |
| 28 | +#include <TMath.h> |
| 29 | +#include <TObjArray.h> |
| 30 | + |
| 31 | +#include <cmath> |
| 32 | + |
| 33 | +using namespace o2; |
| 34 | +using namespace o2::framework; |
| 35 | +using namespace o2::framework::expressions; |
| 36 | + |
| 37 | +void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions) |
| 38 | +{ |
| 39 | + std::vector<ConfigParamSpec> options{ |
| 40 | + {"add-vertex", VariantType::Int, 1, {"Vertex plots"}}, |
| 41 | + {"add-gen", VariantType::Int, 1, {"Generated plots"}}, |
| 42 | + {"add-rec", VariantType::Int, 1, {"Reconstructed plots"}}}; |
| 43 | + std::swap(workflowOptions, options); |
| 44 | +} |
| 45 | + |
| 46 | +#include "Framework/runDataProcessing.h" // important to declare after the options |
| 47 | + |
| 48 | +struct NucleiSpectraEfficienctyVtx { |
| 49 | + OutputObj<TH1F> histVertexTrueZ{TH1F("histVertexTrueZ", "MC true z position of z-vertex; vertex z (cm)", 100, -20., 20.)}; |
| 50 | + |
| 51 | + void process(aod::McCollision const& mcCollision) |
| 52 | + { |
| 53 | + histVertexTrueZ->Fill(mcCollision.posZ()); |
| 54 | + } |
| 55 | +}; |
| 56 | + |
| 57 | +struct NucleiSpectraEfficiencyGen { |
| 58 | + |
| 59 | + HistogramRegistry spectra{"spectraGen", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; |
| 60 | + |
| 61 | + void init(o2::framework::InitContext&) |
| 62 | + { |
| 63 | + std::vector<double> ptBinning = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, |
| 64 | + 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4., 5., 6., 8., 10., 12., 14.}; |
| 65 | + std::vector<double> centBinning = {0., 1., 5., 10., 20., 30., 40., 50., 70., 100.}; |
| 66 | + // |
| 67 | + AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"}; |
| 68 | + AxisSpec centAxis = {centBinning, "V0M (%)"}; |
| 69 | + // |
| 70 | + spectra.add("histGenPt", "generated particles", HistType::kTH1F, {ptAxis}); |
| 71 | + } |
| 72 | + |
| 73 | + void process(aod::McCollision const& mcCollision, aod::McParticles& mcParticles) |
| 74 | + { |
| 75 | + // |
| 76 | + // loop over generated particles and fill generated particles |
| 77 | + // |
| 78 | + for (auto& mcParticleGen : mcParticles) { |
| 79 | + if (mcParticleGen.pdgCode() != -1000020030) { |
| 80 | + continue; |
| 81 | + } |
| 82 | + if (!MC::isPhysicalPrimary(mcParticles, mcParticleGen)) { |
| 83 | + continue; |
| 84 | + } |
| 85 | + if (abs(mcParticleGen.y()) > 0.5) { |
| 86 | + continue; |
| 87 | + } |
| 88 | + spectra.fill(HIST("histGenPt"), mcParticleGen.pt()); |
| 89 | + } |
| 90 | + } |
| 91 | +}; |
| 92 | + |
| 93 | +struct NucleiSpectraEfficiencyRec { |
| 94 | + |
| 95 | + HistogramRegistry spectra{"spectraRec", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; |
| 96 | + |
| 97 | + void init(o2::framework::InitContext&) |
| 98 | + { |
| 99 | + std::vector<double> ptBinning = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, |
| 100 | + 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4., 5., 6., 8., 10., 12., 14.}; |
| 101 | + std::vector<double> centBinning = {0., 1., 5., 10., 20., 30., 40., 50., 70., 100.}; |
| 102 | + // |
| 103 | + AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"}; |
| 104 | + AxisSpec centAxis = {centBinning, "V0M (%)"}; |
| 105 | + // |
| 106 | + spectra.add("histRecVtxZ", "collision z position", HistType::kTH1F, {{600, -20., +20., "z position (cm)"}}); |
| 107 | + spectra.add("histRecPt", "reconstructed particles", HistType::kTH1F, {ptAxis}); |
| 108 | + spectra.add("histTpcSignal", "Specific energy loss", HistType::kTH2F, {{600, -6., 6., "#it{p} (GeV/#it{c})"}, {1400, 0, 1400, "d#it{E} / d#it{X} (a. u.)"}}); |
| 109 | + spectra.add("histTpcNsigma", "n-sigma TPC", HistType::kTH2F, {ptAxis, {200, -100., +100., "n#sigma_{He} (a. u.)"}}); |
| 110 | + } |
| 111 | + |
| 112 | + Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; |
| 113 | + Configurable<float> cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; |
| 114 | + Configurable<float> nsigmacutLow{"nsigmacutLow", -10.0, "Value of the Nsigma cut"}; |
| 115 | + Configurable<float> nsigmacutHigh{"nsigmacutHigh", +10.0, "Value of the Nsigma cut"}; |
| 116 | + |
| 117 | + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; |
| 118 | + Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::isGlobalTrack == (uint8_t) true); |
| 119 | + |
| 120 | + using TrackCandidates = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksExtended, aod::McTrackLabels, aod::pidTPCFullHe, aod::pidTOFFullHe, aod::TrackSelection>>; |
| 121 | + |
| 122 | + void process(soa::Filtered<soa::Join<aod::Collisions, aod::McCollisionLabels>>::iterator const& collision, |
| 123 | + TrackCandidates const& tracks, aod::McParticles& mcParticles, aod::McCollisions const& mcCollisions) |
| 124 | + { |
| 125 | + // |
| 126 | + // check the vertex-z distribution |
| 127 | + // |
| 128 | + spectra.fill(HIST("histRecVtxZ"), collision.posZ()); |
| 129 | + // |
| 130 | + // loop over reconstructed particles and fill reconstructed tracks |
| 131 | + // |
| 132 | + for (auto track : tracks) { |
| 133 | + TLorentzVector lorentzVector{}; |
| 134 | + lorentzVector.SetPtEtaPhiM(track.pt() * 2.0, track.eta(), track.phi(), constants::physics::MassHelium3); |
| 135 | + if (lorentzVector.Rapidity() < -0.5 || lorentzVector.Rapidity() > 0.5) { |
| 136 | + continue; |
| 137 | + } |
| 138 | + // |
| 139 | + // fill QA histograms |
| 140 | + // |
| 141 | + float nSigmaHe3 = track.tpcNSigmaHe(); |
| 142 | + nSigmaHe3 += 94.222101 * TMath::Exp(-0.905203 * track.tpcInnerParam()); |
| 143 | + // |
| 144 | + spectra.fill(HIST("histTpcSignal"), track.tpcInnerParam() * track.sign(), track.tpcSignal()); |
| 145 | + spectra.fill(HIST("histTpcNsigma"), track.tpcInnerParam(), nSigmaHe3); |
| 146 | + // |
| 147 | + // fill histograms |
| 148 | + // |
| 149 | + if (nSigmaHe3 > nsigmacutLow && nSigmaHe3 < nsigmacutHigh) { |
| 150 | + // check on perfect PID |
| 151 | + if (track.mcParticle().pdgCode() != -1000020030) { |
| 152 | + continue; |
| 153 | + } |
| 154 | + // fill reconstructed histogram |
| 155 | + spectra.fill(HIST("histRecPt"), track.pt() * 2.0); |
| 156 | + } |
| 157 | + } |
| 158 | + } |
| 159 | +}; |
| 160 | + |
| 161 | +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
| 162 | +{ |
| 163 | + const bool vertex = cfgc.options().get<int>("add-vertex"); |
| 164 | + const bool gen = cfgc.options().get<int>("add-gen"); |
| 165 | + const bool rec = cfgc.options().get<int>("add-rec"); |
| 166 | + // |
| 167 | + WorkflowSpec workflow{}; |
| 168 | + // |
| 169 | + if (vertex) { |
| 170 | + workflow.push_back(adaptAnalysisTask<NucleiSpectraEfficienctyVtx>(cfgc, TaskName{"nuclei-efficiency-vtx"})); |
| 171 | + } |
| 172 | + if (gen) { |
| 173 | + workflow.push_back(adaptAnalysisTask<NucleiSpectraEfficiencyGen>(cfgc, TaskName{"nuclei-efficiency-gen"})); |
| 174 | + } |
| 175 | + if (rec) { |
| 176 | + workflow.push_back(adaptAnalysisTask<NucleiSpectraEfficiencyRec>(cfgc, TaskName{"nuclei-efficiency-rec"})); |
| 177 | + } |
| 178 | + // |
| 179 | + return workflow; |
| 180 | +} |
0 commit comments