|
| 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 | +/// |
| 13 | +/// \file alice3-dilepton.cxx |
| 14 | +/// \author s.scheid@cern.ch, daiki.sekihata@cern.ch |
| 15 | +/// |
| 16 | + |
| 17 | +#include "ALICE3/DataModel/OTFCollision.h" |
| 18 | +#include "ALICE3/DataModel/OTFRICH.h" |
| 19 | +#include "ALICE3/DataModel/OTFTOF.h" |
| 20 | +#include "ALICE3/DataModel/prefilterDilepton.h" |
| 21 | +#include "ALICE3/DataModel/tracksAlice3.h" |
| 22 | +#include "Common/DataModel/TrackSelectionTables.h" |
| 23 | + |
| 24 | +#include <CommonConstants/PhysicsConstants.h> |
| 25 | +#include <Framework/ASoAHelpers.h> |
| 26 | +#include <Framework/AnalysisDataModel.h> |
| 27 | +#include <Framework/AnalysisTask.h> |
| 28 | +#include <Framework/HistogramRegistry.h> |
| 29 | +#include <Framework/O2DatabasePDGPlugin.h> |
| 30 | +#include <Framework/runDataProcessing.h> |
| 31 | +#include <MathUtils/Utils.h> |
| 32 | + |
| 33 | +#include <Math/GenVector/VectorUtil.h> |
| 34 | +#include <Math/Vector4D.h> |
| 35 | + |
| 36 | +#include <unordered_map> |
| 37 | +#include <vector> |
| 38 | + |
| 39 | +using namespace o2; |
| 40 | +using namespace o2::aod; |
| 41 | +using namespace o2::soa; |
| 42 | +using namespace o2::framework; |
| 43 | +using namespace o2::framework::expressions; |
| 44 | + |
| 45 | +struct Alice3DileptonEventCentSelection { |
| 46 | + |
| 47 | + Service<o2::framework::O2DatabasePDG> inspdg; |
| 48 | + |
| 49 | + Produces<aod::DiEventCentCuts> eventCentSel; |
| 50 | + Configurable<float> minFITPart{"minFITPart", 0.f, "Minimum number of charged particles in the FIT acceptance"}; |
| 51 | + Configurable<float> maxFITPart{"maxFITPart", 0.f, "Maximum number of charged particles in the FIT acceptance"}; |
| 52 | + |
| 53 | + HistogramRegistry registry{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; |
| 54 | + |
| 55 | + void init(InitContext&) |
| 56 | + { |
| 57 | + registry.add("Generated/Before/ParticlesFit", "Charged Particles in Fit acceptance per event", kTH1F, {{15000, 0, 15000}}); |
| 58 | + registry.add("Generated/After/ParticlesFit", "Charged Particles in Fit acceptance per event", kTH1F, {{15000, 0, 15000}}); |
| 59 | + registry.add("Generated/Before/ParticlesEta", "Charged Particles per event", kTH1F, {{1200, -6, 6}}); |
| 60 | + registry.add("Generated/After/ParticlesEta", "Charged Particles per event", kTH1F, {{1200, -6, 6}}); |
| 61 | + } |
| 62 | + |
| 63 | + using MyEvents = soa::Join<aod::Collisions, aod::McCollisionLabels>; |
| 64 | + Preslice<aod::McParticles> perMCCollision = o2::aod::mcparticle::mcCollisionId; |
| 65 | + |
| 66 | + void processGen(MyEvents::iterator const& event, o2::aod::McCollisions const&, aod::McParticles const& mcParticles) |
| 67 | + { |
| 68 | + if (!event.has_mcCollision()) { |
| 69 | + eventCentSel(0); |
| 70 | + } else { |
| 71 | + auto mccollision = event.mcCollision(); |
| 72 | + auto mcParticles_per_coll = mcParticles.sliceBy(perMCCollision, mccollision.globalIndex()); |
| 73 | + int nParticlesFIT = 0; |
| 74 | + for (const auto& mcParticle : mcParticles_per_coll) { |
| 75 | + if (mcParticle.isPhysicalPrimary()) { |
| 76 | + auto pdgParticle = inspdg->GetParticle(mcParticle.pdgCode()); |
| 77 | + if (pdgParticle) { |
| 78 | + float charge = pdgParticle->Charge() / 3.f; // Charge in units of |e| |
| 79 | + if (std::abs(charge) >= 1.) { |
| 80 | + registry.fill(HIST("Generated/Before/ParticlesEta"), mcParticle.eta()); |
| 81 | + if ((2.2 < mcParticle.eta() && mcParticle.eta() < 5.0) || (-3.4 < mcParticle.eta() && mcParticle.eta() < -2.3)) { |
| 82 | + nParticlesFIT++; |
| 83 | + } |
| 84 | + } |
| 85 | + } |
| 86 | + } |
| 87 | + } // end of mc particle loop |
| 88 | + registry.fill(HIST("Generated/Before/ParticlesFit"), nParticlesFIT); |
| 89 | + if (nParticlesFIT > minFITPart && nParticlesFIT < maxFITPart) { |
| 90 | + registry.fill(HIST("Generated/After/ParticlesFit"), nParticlesFIT); |
| 91 | + for (const auto& mcParticle : mcParticles_per_coll) { |
| 92 | + if (mcParticle.isPhysicalPrimary()) { |
| 93 | + auto pdgParticle = inspdg->GetParticle(mcParticle.pdgCode()); |
| 94 | + if (pdgParticle) { |
| 95 | + float charge = pdgParticle->Charge() / 3.f; // Charge in units of |e| |
| 96 | + if (std::abs(charge) >= 1.) { |
| 97 | + registry.fill(HIST("Generated/After/ParticlesEta"), mcParticle.eta()); |
| 98 | + } |
| 99 | + } |
| 100 | + } |
| 101 | + } // end of mc particle loop |
| 102 | + eventCentSel(1); |
| 103 | + } else { |
| 104 | + eventCentSel(0); |
| 105 | + } |
| 106 | + } // if mc collision loop |
| 107 | + } // end of process |
| 108 | + |
| 109 | + void processDummy(MyEvents::iterator const&) |
| 110 | + { |
| 111 | + eventCentSel(1); |
| 112 | + } |
| 113 | + |
| 114 | + PROCESS_SWITCH(Alice3DileptonEventCentSelection, processGen, "Select centrality", false); |
| 115 | + PROCESS_SWITCH(Alice3DileptonEventCentSelection, processDummy, "Dummy", true); |
| 116 | +}; |
| 117 | + |
| 118 | +struct Alice3DileptonPrefilter { |
| 119 | + |
| 120 | + Produces<aod::DiTrackPrefilter> pfb_derived; |
| 121 | + |
| 122 | + SliceCache cache_mc; |
| 123 | + SliceCache cache_rec; |
| 124 | + |
| 125 | + Service<o2::framework::O2DatabasePDG> inspdg; |
| 126 | + |
| 127 | + Configurable<float> ptMin{"ptMin", 0.f, "Lower limit in pT"}; |
| 128 | + Configurable<float> ptMax{"ptMax", 5.f, "Upper limit in pT"}; |
| 129 | + Configurable<float> etaMin{"etaMin", -5.f, "Lower limit in eta"}; |
| 130 | + Configurable<float> etaMax{"etaMax", 5.f, "Upper limit in eta"}; |
| 131 | + Configurable<float> maxMass{"maxMass", 5.f, "Upper limit in mass"}; |
| 132 | + Configurable<float> maxOp{"maxOp", 5.f, "Upper limit in opening angle"}; |
| 133 | + Configurable<float> nSigmaEleCutOuterTOF{"nSigmaEleCutOuterTOF", 3., "Electron inclusion in outer TOF"}; |
| 134 | + Configurable<float> nSigmaEleCutInnerTOF{"nSigmaEleCutInnerTOF", 3., "Electron inclusion in inner TOF"}; |
| 135 | + Configurable<float> nSigmaPionCutOuterTOF{"nSigmaPionCutOuterTOF", 3., "Pion exclusion in outer TOF"}; |
| 136 | + Configurable<float> nSigmaPionCutInnerTOF{"nSigmaPionCutInnerTOF", 3., "Pion exclusion in inner TOF"}; |
| 137 | + Configurable<float> nSigmaElectronRich{"nSigmaElectronRich", 3., "Electron inclusion RICH"}; |
| 138 | + Configurable<float> nSigmaPionRich{"nSigmaPionRich", 3., "Pion exclusion RICH"}; |
| 139 | + |
| 140 | + std::unordered_map<int, uint16_t> map_pfb; // map track.globalIndex -> prefilter bit |
| 141 | + |
| 142 | + HistogramRegistry registry{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; |
| 143 | + |
| 144 | + using MyTracksMCs = soa::Join<aod::Tracks, aod::UpgradeTofs, aod::UpgradeTofMCs, aod::UpgradeRichs, aod::TracksAlice3>; |
| 145 | + using MyTracksMC = MyTracksMCs::iterator; |
| 146 | + using DileptonCollisions = soa::Join<aod::Collisions, aod::DiEventCentCuts>; |
| 147 | + using DileptonCollision = DileptonCollisions::iterator; |
| 148 | + |
| 149 | + Preslice<MyTracksMCs> perCollision = aod::track::collisionId; |
| 150 | + Partition<MyTracksMCs> posTracks = o2::aod::track::signed1Pt > 0.f; |
| 151 | + Partition<MyTracksMCs> negTracks = o2::aod::track::signed1Pt < 0.f; |
| 152 | + |
| 153 | + void init(InitContext&) |
| 154 | + { |
| 155 | + const AxisSpec axisM{500, 0, 5, "#it{m}_{ll} (GeV/#it{c}^{2})"}; |
| 156 | + const AxisSpec axisPt{1000, 0, 10, "#it{p}_{T} (GeV/#it{c})"}; |
| 157 | + const AxisSpec axisSigmaEl{200, -10, 10, "n#sigma_{El}"}; |
| 158 | + const AxisSpec axisTrackLengthOuterTOF{300, 0., 300., "Track length (cm)"}; |
| 159 | + const AxisSpec axisEta{1000, -5, 5, "#it{#eta}"}; |
| 160 | + const AxisSpec axisPhi{360, 0, TMath::TwoPi(), "#it{#varphi} (rad.)"}; |
| 161 | + const AxisSpec axisProdx{2000, -100, 100, "Prod. Vertex X (cm)"}; |
| 162 | + const AxisSpec axisPrody{2000, -100, 100, "Prod. Vertex Y (cm)"}; |
| 163 | + const AxisSpec axisProdz{2000, -100, 100, "Prod. Vertex Z (cm)"}; |
| 164 | + |
| 165 | + registry.add("Reconstructed/Pair/ULS/Mass_Pt", "Pair Mass vs. Pt", kTH2F, {axisM, axisPt}, true); |
| 166 | + registry.add("ReconstructedFiltered/Pair/ULS/Mass_Pt", "Pair Mass vs. Pt", kTH2F, {axisM, axisPt}, true); |
| 167 | + registry.add("Reconstructed/Track/Pt", "Track Pt", kTH1F, {axisPt}); |
| 168 | + registry.add("Reconstructed/Track/Eta", "Track Eta", kTH1F, {axisEta}); |
| 169 | + registry.add("Reconstructed/Track/Phi", "Track Phi", kTH1F, {axisPhi}); |
| 170 | + registry.add("Reconstructed/Track/Eta_Pt", "Eta vs. Pt", kTH2F, {axisPt, axisEta}, true); |
| 171 | + registry.add("Reconstructed/Track/SigmaOTofvspt", "Track #sigma oTOF", kTH2F, {axisPt, axisSigmaEl}); |
| 172 | + registry.add("Reconstructed/Track/SigmaITofvspt", "Track #sigma iTOF", kTH2F, {axisPt, axisSigmaEl}); |
| 173 | + registry.add("Reconstructed/Track/SigmaRichvspt", "Track #sigma RICH", kTH2F, {axisPt, axisSigmaEl}); |
| 174 | + registry.add("Reconstructed/Track/outerTOFTrackLength", "Track length outer TOF", kTH1F, {axisTrackLengthOuterTOF}); |
| 175 | + } |
| 176 | + |
| 177 | + void processPreFilter(DileptonCollisions const& collisions, MyTracksMCs const& tracks) |
| 178 | + { |
| 179 | + for (const auto& track : tracks) { |
| 180 | + map_pfb[track.globalIndex()] = 0; |
| 181 | + } // end of track loop |
| 182 | + |
| 183 | + for (const auto& collision : collisions) { |
| 184 | + auto negTracks_coll = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache_rec); |
| 185 | + auto posTracks_coll = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache_rec); |
| 186 | + |
| 187 | + if (collision.isEventCentSelected() == 0) { |
| 188 | + for (const auto& pos : posTracks_coll) { |
| 189 | + map_pfb[pos.globalIndex()] = 0; |
| 190 | + } |
| 191 | + for (const auto& neg : negTracks_coll) { |
| 192 | + map_pfb[neg.globalIndex()] = 0; |
| 193 | + } |
| 194 | + continue; |
| 195 | + } |
| 196 | + |
| 197 | + for (const auto& pos : posTracks_coll) { |
| 198 | + if (!pos.isReconstructed()) { |
| 199 | + continue; |
| 200 | + } |
| 201 | + if (pos.eta() < etaMin || etaMax < pos.eta()) { |
| 202 | + continue; |
| 203 | + } |
| 204 | + if (pos.pt() < ptMin || ptMax < pos.pt()) { |
| 205 | + continue; |
| 206 | + } |
| 207 | + if ((std::abs(pos.nSigmaElectronRich()) < nSigmaElectronRich && nSigmaPionRich < std::abs(pos.nSigmaPionRich())) || (std::abs(pos.nSigmaElectronOuterTOF()) < nSigmaEleCutOuterTOF && nSigmaPionCutOuterTOF < std::abs(pos.nSigmaPionOuterTOF())) || (std::abs(pos.nSigmaElectronInnerTOF()) < nSigmaEleCutInnerTOF && nSigmaPionCutInnerTOF < std::abs(pos.nSigmaPionInnerTOF()))) { |
| 208 | + registry.fill(HIST("Reconstructed/Track/SigmaOTofvspt"), pos.pt(), pos.nSigmaElectronOuterTOF()); |
| 209 | + registry.fill(HIST("Reconstructed/Track/SigmaITofvspt"), pos.pt(), pos.nSigmaElectronInnerTOF()); |
| 210 | + registry.fill(HIST("Reconstructed/Track/SigmaRichvspt"), pos.pt(), pos.nSigmaElectronRich()); |
| 211 | + registry.fill(HIST("Reconstructed/Track/outerTOFTrackLength"), pos.outerTOFTrackLength()); |
| 212 | + registry.fill(HIST("Reconstructed/Track/Pt"), pos.pt()); |
| 213 | + registry.fill(HIST("Reconstructed/Track/Eta"), pos.eta()); |
| 214 | + registry.fill(HIST("Reconstructed/Track/Phi"), pos.phi()); |
| 215 | + registry.fill(HIST("Reconstructed/Track/Eta_Pt"), pos.pt(), pos.eta()); |
| 216 | + } |
| 217 | + } |
| 218 | + for (const auto& ele : negTracks_coll) { |
| 219 | + if (!ele.isReconstructed()) { |
| 220 | + continue; |
| 221 | + } |
| 222 | + if (ele.eta() < etaMin || etaMax < ele.eta()) { |
| 223 | + continue; |
| 224 | + } |
| 225 | + if (ele.pt() < ptMin || ptMax < ele.pt()) { |
| 226 | + continue; |
| 227 | + } |
| 228 | + if ((std::abs(ele.nSigmaElectronRich()) < nSigmaElectronRich && nSigmaPionRich < std::abs(ele.nSigmaPionRich())) || (std::abs(ele.nSigmaElectronOuterTOF()) < nSigmaEleCutOuterTOF && nSigmaPionCutOuterTOF < std::abs(ele.nSigmaPionOuterTOF())) || (std::abs(ele.nSigmaElectronInnerTOF()) < nSigmaEleCutInnerTOF && nSigmaPionCutInnerTOF < std::abs(ele.nSigmaPionInnerTOF()))) { |
| 229 | + registry.fill(HIST("Reconstructed/Track/SigmaOTofvspt"), ele.pt(), ele.nSigmaElectronOuterTOF()); |
| 230 | + registry.fill(HIST("Reconstructed/Track/SigmaITofvspt"), ele.pt(), ele.nSigmaElectronInnerTOF()); |
| 231 | + registry.fill(HIST("Reconstructed/Track/SigmaRichvspt"), ele.pt(), ele.nSigmaElectronRich()); |
| 232 | + registry.fill(HIST("Reconstructed/Track/outerTOFTrackLength"), ele.outerTOFTrackLength()); |
| 233 | + registry.fill(HIST("Reconstructed/Track/Pt"), ele.pt()); |
| 234 | + registry.fill(HIST("Reconstructed/Track/Eta"), ele.eta()); |
| 235 | + registry.fill(HIST("Reconstructed/Track/Phi"), ele.phi()); |
| 236 | + registry.fill(HIST("Reconstructed/Track/Eta_Pt"), ele.pt(), ele.eta()); |
| 237 | + } |
| 238 | + } |
| 239 | + |
| 240 | + for (const auto& [pos, ele] : combinations(CombinationsFullIndexPolicy(posTracks_coll, negTracks_coll))) { // ULS |
| 241 | + if (!pos.isReconstructed()) { |
| 242 | + continue; |
| 243 | + } |
| 244 | + if (pos.eta() < etaMin || etaMax < pos.eta()) { |
| 245 | + continue; |
| 246 | + } |
| 247 | + if (pos.pt() < ptMin || ptMax < pos.pt()) { |
| 248 | + continue; |
| 249 | + } |
| 250 | + if (!ele.isReconstructed()) { |
| 251 | + continue; |
| 252 | + } |
| 253 | + if (ele.eta() < etaMin || etaMax < ele.eta()) { |
| 254 | + continue; |
| 255 | + } |
| 256 | + if (ele.pt() < ptMin || ptMax < ele.pt()) { |
| 257 | + continue; |
| 258 | + } |
| 259 | + if ((std::abs(pos.nSigmaElectronRich()) < nSigmaElectronRich && nSigmaPionRich < std::abs(pos.nSigmaPionRich())) || (std::abs(pos.nSigmaElectronOuterTOF()) < nSigmaEleCutOuterTOF && nSigmaPionCutOuterTOF < std::abs(pos.nSigmaPionOuterTOF())) || (std::abs(pos.nSigmaElectronInnerTOF()) < nSigmaEleCutInnerTOF && nSigmaPionCutInnerTOF < std::abs(pos.nSigmaPionInnerTOF()))) { |
| 260 | + if ((std::abs(ele.nSigmaElectronRich()) < nSigmaElectronRich && nSigmaPionRich < std::abs(ele.nSigmaPionRich())) || (std::abs(ele.nSigmaElectronOuterTOF()) < nSigmaEleCutOuterTOF && nSigmaPionCutOuterTOF < std::abs(ele.nSigmaPionOuterTOF())) || (std::abs(ele.nSigmaElectronInnerTOF()) < nSigmaEleCutInnerTOF && nSigmaPionCutInnerTOF < std::abs(ele.nSigmaPionInnerTOF()))) { |
| 261 | + ROOT::Math::PtEtaPhiMVector v1(pos.pt(), pos.eta(), pos.phi(), o2::constants::physics::MassElectron); |
| 262 | + ROOT::Math::PtEtaPhiMVector v2(ele.pt(), ele.eta(), ele.phi(), o2::constants::physics::MassElectron); |
| 263 | + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; |
| 264 | + float angle = ROOT::Math::VectorUtil::Angle(v1, v2); |
| 265 | + o2::math_utils::bringToPMPi(angle); |
| 266 | + |
| 267 | + registry.fill(HIST("Reconstructed/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt()); |
| 268 | + |
| 269 | + if (v12.M() < maxMass && angle < maxOp) { |
| 270 | + map_pfb[pos.globalIndex()] = 1; |
| 271 | + map_pfb[ele.globalIndex()] = 1; |
| 272 | + registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt()); |
| 273 | + } |
| 274 | + } |
| 275 | + } |
| 276 | + } |
| 277 | + } // end of collision |
| 278 | + |
| 279 | + for (const auto& track : tracks) { |
| 280 | + pfb_derived(map_pfb[track.globalIndex()]); |
| 281 | + } // end of track loop |
| 282 | + map_pfb.clear(); |
| 283 | + } |
| 284 | + |
| 285 | + void processDummy(MyTracksMCs const& tracks) |
| 286 | + { |
| 287 | + for (int i = 0; i < tracks.size(); i++) { |
| 288 | + pfb_derived(0); |
| 289 | + } |
| 290 | + } |
| 291 | + |
| 292 | + PROCESS_SWITCH(Alice3DileptonPrefilter, processPreFilter, "Run prefilter", false); |
| 293 | + PROCESS_SWITCH(Alice3DileptonPrefilter, processDummy, "Dummy", true); |
| 294 | +}; |
| 295 | + |
| 296 | +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
| 297 | +{ |
| 298 | + return WorkflowSpec{ |
| 299 | + adaptAnalysisTask<Alice3DileptonEventCentSelection>(cfgc), |
| 300 | + adaptAnalysisTask<Alice3DileptonPrefilter>(cfgc)}; |
| 301 | +} |
0 commit comments