|
| 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 emcalGammaGammaBcWise.cxx |
| 13 | +/// |
| 14 | +/// \brief This code loops over BCs to pair photons using only EMCal + TVX (no central barrel, no collisions) |
| 15 | +/// |
| 16 | +/// \author Nicolas Strangmann ([email protected]) - Goethe University Frankfurt |
| 17 | +/// |
| 18 | + |
| 19 | +#include <CommonConstants/MathConstants.h> |
| 20 | +#include <vector> |
| 21 | +#include "TLorentzVector.h" |
| 22 | +#include "TVector3.h" |
| 23 | + |
| 24 | +#include "Framework/runDataProcessing.h" |
| 25 | +#include "Framework/AnalysisTask.h" |
| 26 | +#include "Framework/HistogramRegistry.h" |
| 27 | +#include "Common/DataModel/EventSelection.h" |
| 28 | +#include "PWGJE/DataModel/EMCALClusters.h" |
| 29 | + |
| 30 | +using namespace o2; |
| 31 | +using namespace o2::framework; |
| 32 | +using namespace o2::framework::expressions; |
| 33 | + |
| 34 | +using MyBCs = soa::Join<aod::BCs, aod::BcSels, aod::Timestamps, aod::Run3MatchedToBCSparse>; |
| 35 | +using MyCollisions = soa::Join<aod::Collisions, aod::EvSels>; |
| 36 | + |
| 37 | +using SelectedUniqueClusters = soa::Filtered<aod::EMCALClusters>; // Clusters from collisions with only one collision in the BC |
| 38 | +using SelectedAmbiguousClusters = soa::Filtered<aod::EMCALAmbiguousClusters>; // Clusters from BCs with multiple collisions (no vertex assignment possible) |
| 39 | + |
| 40 | +struct Photon { // Struct to store photons (unique and ambiguous clusters that passed the cuts) |
| 41 | + Photon(float eta, float phi, float energy) : eta(eta), phi(phi), e(energy), theta(2 * std::atan2(std::exp(-eta), 1)), px(e * std::sin(theta) * std::cos(phi)), py(e * std::sin(theta) * std::sin(phi)), pz(e * std::cos(theta)), pt(std::sqrt(px * px + py * py)) |
| 42 | + { |
| 43 | + photon.SetPxPyPzE(px, py, pz, e); |
| 44 | + } |
| 45 | + |
| 46 | + TLorentzVector photon; |
| 47 | + float eta, phi, e; |
| 48 | + float theta; |
| 49 | + float px, py, pz, pt; |
| 50 | +}; |
| 51 | + |
| 52 | +struct Meson { |
| 53 | + Meson(Photon p1, Photon p2) : p1(p1), p2(p2) |
| 54 | + { |
| 55 | + pMeson = p1.photon + p2.photon; |
| 56 | + } |
| 57 | + Photon p1, p2; |
| 58 | + TLorentzVector pMeson; |
| 59 | + |
| 60 | + float m() const { return pMeson.M(); } |
| 61 | + float pT() const { return pMeson.Pt(); } |
| 62 | + float openingAngle() const { return p1.photon.Angle(p2.photon.Vect()); } |
| 63 | +}; |
| 64 | + |
| 65 | +struct EmcalGammaGammaBcWise { |
| 66 | + HistogramRegistry mHistManager{"EmcalGammaGammaBcWiseHistograms"}; |
| 67 | + |
| 68 | + Configurable<float> cfgMinClusterEnergy{"cfgMinClusterEnergy", 0.7, "Minimum energy of selected clusters (GeV)"}; |
| 69 | + Configurable<float> cfgMinM02{"cfgMinM02", 0.1, "Minimum M02 of selected clusters"}; |
| 70 | + Configurable<float> cfgMaxM02{"cfgMaxM02", 0.7, "Maximum M02 of selected clusters"}; |
| 71 | + Configurable<float> cfgMinTime{"cfgMinTime", -15, "Minimum time of selected clusters (ns)"}; |
| 72 | + Configurable<float> cfgMaxTime{"cfgMaxTime", 15, "Maximum time of selected clusters (ns)"}; |
| 73 | + |
| 74 | + Filter energyFilter = aod::emcalcluster::energy > cfgMinClusterEnergy; |
| 75 | + Filter m02Filter = (aod::emcalcluster::nCells == 1 || (aod::emcalcluster::m02 > cfgMinM02 && aod::emcalcluster::m02 < cfgMaxM02)); |
| 76 | + Filter timeFilter = (aod::emcalcluster::time > cfgMinTime && aod::emcalcluster::time < cfgMaxTime); |
| 77 | + |
| 78 | + std::vector<Photon> mPhotons; |
| 79 | + |
| 80 | + void init(InitContext const&) |
| 81 | + { |
| 82 | + mHistManager.add("nBCs", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{3, -0.5, 2.5}}); |
| 83 | + |
| 84 | + mHistManager.add("nCollisionsVsClusters", "Number of collisions vs number of clusters;N_{Collisions};N_{Clusters}", HistType::kTH2F, {{4, -0.5, 3.5}, {21, -0.5, 20.5}}); |
| 85 | + |
| 86 | + mHistManager.add("clusterE", "Energy of cluster;#bf{#it{E} (GeV)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 20}}); |
| 87 | + mHistManager.add("clusterM02", "Shape of cluster;#bf{#it{M}_{02}};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 2}}); |
| 88 | + mHistManager.add("clusterTime", "Time of cluster;#bf{#it{t} (ns)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, -100, 100}}); |
| 89 | + |
| 90 | + mHistManager.add("invMassVsPt", "Invariant mass and pT of meson candidates", HistType::kTH2F, {{400, 0., 0.8}, {200, 0., 20.}}); |
| 91 | + mHistManager.add("invMassVsPtBackground", "Invariant mass and pT of background meson candidates", HistType::kTH2F, {{400, 0., 0.8}, {200, 0., 20.}}); |
| 92 | + } |
| 93 | + |
| 94 | + PresliceUnsorted<MyCollisions> perFoundBC = aod::evsel::foundBCId; |
| 95 | + Preslice<aod::EMCALClusters> perCol = aod::emcalcluster::collisionId; |
| 96 | + Preslice<aod::EMCALAmbiguousClusters> perBC = aod::emcalcluster::bcId; |
| 97 | + |
| 98 | + void process(MyBCs const& bcs, MyCollisions const& collisions, SelectedUniqueClusters const& uClusters, SelectedAmbiguousClusters const& aClusters) |
| 99 | + { |
| 100 | + for (const auto& bc : bcs) { |
| 101 | + mHistManager.fill(HIST("nBCs"), 0.); |
| 102 | + if (!bc.selection_bit(aod::evsel::kIsTriggerTVX)) |
| 103 | + continue; |
| 104 | + mHistManager.fill(HIST("nBCs"), 1.); |
| 105 | + if (!bc.alias_bit(kTVXinEMC)) |
| 106 | + continue; |
| 107 | + mHistManager.fill(HIST("nBCs"), 2.); |
| 108 | + |
| 109 | + auto collisionsInFoundBC = collisions.sliceBy(perFoundBC, bc.globalIndex()); |
| 110 | + |
| 111 | + if (collisionsInFoundBC.size() == 1) { // Unique |
| 112 | + auto clustersInCollision = uClusters.sliceBy(perCol, collisionsInFoundBC.begin().globalIndex()); |
| 113 | + processClusters(clustersInCollision); |
| 114 | + } else { // Ambiguous |
| 115 | + auto clustersInBC = aClusters.sliceBy(perBC, bc.globalIndex()); |
| 116 | + processClusters(clustersInBC); |
| 117 | + } |
| 118 | + |
| 119 | + mHistManager.fill(HIST("nCollisionsVsClusters"), collisionsInFoundBC.size(), mPhotons.size()); |
| 120 | + |
| 121 | + processMesons(); |
| 122 | + } |
| 123 | + } |
| 124 | + |
| 125 | + /// \brief Process EMCAL clusters (either ambigous or unique) |
| 126 | + template <typename Clusters> |
| 127 | + void processClusters(Clusters const& clusters) |
| 128 | + { |
| 129 | + mPhotons.clear(); |
| 130 | + |
| 131 | + // loop over all clusters from accepted collision |
| 132 | + for (const auto& cluster : clusters) { |
| 133 | + |
| 134 | + mHistManager.fill(HIST("clusterE"), cluster.energy()); |
| 135 | + mHistManager.fill(HIST("clusterM02"), cluster.m02()); |
| 136 | + mHistManager.fill(HIST("clusterTime"), cluster.time()); |
| 137 | + |
| 138 | + mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy())); |
| 139 | + } |
| 140 | + } |
| 141 | + |
| 142 | + /// \brief Process meson candidates, calculate invariant mass and pT and fill histograms |
| 143 | + void processMesons() |
| 144 | + { |
| 145 | + if (mPhotons.size() < 2) // if less then 2 clusters are found, skip event |
| 146 | + return; |
| 147 | + |
| 148 | + // loop over all photon combinations and build meson candidates |
| 149 | + for (unsigned int ig1 = 0; ig1 < mPhotons.size(); ++ig1) { |
| 150 | + for (unsigned int ig2 = ig1 + 1; ig2 < mPhotons.size(); ++ig2) { |
| 151 | + // build meson from photons |
| 152 | + Meson meson(mPhotons[ig1], mPhotons[ig2]); |
| 153 | + mHistManager.fill(HIST("invMassVsPt"), meson.m(), meson.pT()); |
| 154 | + |
| 155 | + calculateBackground(meson, ig1, ig2); // calculate background candidates (rotation background) |
| 156 | + } |
| 157 | + } |
| 158 | + } |
| 159 | + |
| 160 | + /// \brief Calculate background (using rotation background method) |
| 161 | + void calculateBackground(const Meson& meson, const unsigned int ig1, const unsigned int ig2) |
| 162 | + { |
| 163 | + if (mPhotons.size() < 3) // if less than 3 clusters are present, skip event |
| 164 | + return; |
| 165 | + |
| 166 | + TVector3 lvRotationPion = (meson.pMeson).Vect(); // calculate rotation axis |
| 167 | + for (unsigned int ig3 = 0; ig3 < mPhotons.size(); ++ig3) { |
| 168 | + for (const unsigned int ig : {ig1, ig2}) { |
| 169 | + if (ig == ig3) |
| 170 | + continue; |
| 171 | + |
| 172 | + TLorentzVector lvRotationPhoton(mPhotons[ig].px, mPhotons[ig].py, mPhotons[ig].pz, mPhotons[ig].e); |
| 173 | + lvRotationPhoton.Rotate(constants::math::PIHalf, lvRotationPion); |
| 174 | + Photon rotPhoton(lvRotationPhoton.Eta(), lvRotationPhoton.Phi(), lvRotationPhoton.E()); |
| 175 | + Meson mesonRotated(rotPhoton, mPhotons[ig3]); |
| 176 | + mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated.m(), mesonRotated.pT()); |
| 177 | + } |
| 178 | + } |
| 179 | + } |
| 180 | +}; |
| 181 | + |
| 182 | +WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc) |
| 183 | +{ |
| 184 | + return WorkflowSpec{adaptAnalysisTask<EmcalGammaGammaBcWise>(cfgc)}; |
| 185 | +} |
0 commit comments