rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2023_I2628732

W+D production in pp at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 2628732
Status: VALIDATED
Authors:
  • Miha Muskinja
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • W+charm production

The production of a $W$ boson in association with a single charm quark is studied using 140 ifb$^{-1}$ of i$\sqrt{s} = 13$TeV proton-proton collision data collected with the ATLAS detector at the Large Hadron Collider. The charm quark is tagged by a charmed hadron, reconstructed with a secondary-vertex fit. The $W$ boson is reconstructed from an electron/muon decay and the missing transverse momentum. The mesons reconstructed are $D^\pm\to K^\mp\pi^\pm\pi^\pm$ and $D^{\ast\pm} \to D^0\pi^\pm \to (K^\mp\pi^\pm)\pi^\pm$, where $p_\text{T}(e,\mu)>30$GeV, $|\eta(e,\mu)|<2.5$, $p_\text{T}(D)>8$GeV, and $|\eta(D)|<2.2$. The integrated and normalized differential cross-sections as a function of the pseudorapidity of the lepton from the $W$ boson decay, and of the transverse momentum of the meson, are extracted from the data using a profile likelihood fit. The measured fiducial cross-sections are $\sigma_\text{fid}^{OS-SS}(W^- + D^+) = 50.2\pm 0.2$ (stat.) $^{+2.4}_{-2.3}$ (syst.) pb, $\sigma_\text{fid}^{OS-SS}(W^+ + D^-) = 48.2\pm 0.2$ (stat.) $^{+2.3}_{-2.2}$ (syst.) pb, $\sigma_\text{fid}^{OS-SS}(W^- + D^{\ast+}) = 51.1\pm 0.4$ (stat.) $^{+1.9}_{-1.8}$ (syst.) pb, and $\sigma_\text{fid}^{OS-SS}(W^+ + D^{\ast-}) = 50.0\pm 0.4$ (stat.) $^{+1.9}_{-1.8}$ (syst.) pb. Results are compared with the predictions of next-to-leading-order quantum chromodynamics calculations performed using state-of-the-art parton distribution functions. The ratio of charm to anti-charm production cross-sections is studied to probe the $s-\bar{s}$ quark asymmetry and is found to be $R^\pm_c = 0.971 \pm 0.006$ (stat.) $\pm0.011$ (syst.).

Source code: ATLAS_2023_I2628732.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/LeptonFinder.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/UnstableParticles.hh"
  8#include <unordered_map>
  9
 10namespace Rivet {
 11
 12
 13  /// @brief W+D production in pp at 13 TeV
 14  class ATLAS_2023_I2628732 : public Analysis {
 15  public:
 16
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2023_I2628732);
 18
 19    // Book histograms and initialise before the run
 20    void init() {
 21      const FinalState fs;
 22
 23      // lepton kinematic cuts
 24      Cut cuts(Cuts::pT > 30*GeV && Cuts::abseta < 2.5);
 25
 26      // Get photons to dress leptons
 27      FinalState photons(Cuts::abspid == PID::PHOTON);
 28
 29      // Get dressed leptons
 30      IdentifiedFinalState lepids(fs, {{PID::ELECTRON, PID::POSITRON, PID::MUON, PID::ANTIMUON}});
 31      PromptFinalState leptons(lepids, TauDecaysAs::NONPROMPT);
 32      LeptonFinder dressedleptons(leptons, photons, 0.1, cuts);
 33      declare(dressedleptons, "LeptonFinder");
 34
 35      // unstable final-state for Ds
 36      declare(UnstableParticles(), "UFS");
 37
 38      // Fiducial cross sections vs species:
 39      // D+ and D* production fractions can be reweighted to the world average values.
 40      // We can calculate them for each MC sample separately with the `CharmSpecies`
 41      // histogram from the Rivet routine, e.g.:
 42      //   - f(D+) = CharmSpecies->GetBinContent(2) / CharmSpecies->Integral(2,8)
 43      //   - f(D*) = CharmSpecies->GetBinContent(1) / CharmSpecies->Integral(2,8)
 44      book(_h["CharmSpecies"], "_CharmSpecies", 8, 0, 8);
 45
 46      // Differential cross sections per bin
 47      bookPair("lep_minus", "Dplus", "D_pt",         3);
 48      bookPair("lep_plus",  "Dplus", "D_pt",         4);
 49      bookPair("lep_minus", "Dplus", "lep_abs_eta",  5);
 50      bookPair("lep_plus",  "Dplus", "lep_abs_eta",  6);
 51      bookPair("lep_minus", "Dstar", "D_pt",         7);
 52      bookPair("lep_plus",  "Dstar", "D_pt",         8);
 53      bookPair("lep_minus", "Dstar", "lep_abs_eta",  9);
 54      bookPair("lep_plus",  "Dstar", "lep_abs_eta", 10);
 55    }
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event &event) {
 59      // Retrieve the dressed electrons
 60      const Particles &signal_leptons = apply<LeptonFinder>(event, "LeptonFinder").particlesByPt();
 61      if (signal_leptons.size() != 1)  vetoEvent;
 62
 63      const Particle &lepton = signal_leptons[0];
 64      const std::string lepton_name = _lepton_names.at(lepton.pid());
 65
 66      // Get the charm hadrons
 67      const UnstableParticles &ufs = apply<UnstableFinalState>(event, "UFS");
 68      std::unordered_map<unsigned int, Particles> particles;
 69
 70      // Loop over particles
 71      for (const Particle &p : ufs.particles()) {
 72        const int id = p.abspid();
 73        const double pt = p.pT() / GeV;
 74        const double eta = p.abseta();
 75        if (_charm_hadron_names.count(id) && pt > 8.0 && eta < 2.2) {
 76          particles[id].push_back(p);
 77        }
 78      }
 79
 80      // Fill histograms
 81      for (auto &kv : particles) {
 82        const unsigned int absPdgId = kv.first;
 83        const std::string hadron_name = _charm_hadron_names.at(absPdgId);
 84
 85        for (auto &p : kv.second) {
 86          // Weight: +1 for OS and -1 for SS
 87          float charm_charge = (absPdgId == 421) ? p.pid() : p.charge();
 88          double weight = (charm_charge * lepton.charge() < 0) ? +1.0 : -1.0;
 89
 90          // Fill charm species for production fraction reweighting
 91          _h["CharmSpecies"]->fill(_charm_species_map.at(absPdgId), weight);
 92
 93          // Fill only D+ and D* histograms
 94          if (absPdgId != PID::DPLUS && absPdgId != PID::DSTARPLUS)  continue;
 95
 96          // pT(D) overflow
 97          // Last pT(D) bin extends to infinity (150 only for practical purposes)
 98          double pt = p.pT() / GeV;
 99          if (pt >= _max_D_pt)  pt = _max_D_pt - 10;
100
101          // Fill histograms
102          _h[histo_name(lepton_name, hadron_name, "lep_abs_eta")]->fill(lepton.abseta(), weight);
103          _h[histo_name(lepton_name, hadron_name, "D_pt")]->fill(pt, weight);
104          _h[histo_name(lepton_name, hadron_name, "lep_abs_eta") + "_norm"]->fill(lepton.abseta(), weight);
105          _h[histo_name(lepton_name, hadron_name, "D_pt") + "_norm"]->fill(pt, weight);
106        }
107      }
108    }
109
110    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
111    /// Finalize
112    void finalize() {
113
114      scale(_h, crossSectionPerEvent());
115
116      // D+ and D* production fractions
117      const double sum = _h["CharmSpecies"]->integral(false);
118      const double fDplus = safediv(_h["CharmSpecies"]->bin(2).sumW(), sum);
119      const double fDstar = safediv(_h["CharmSpecies"]->bin(1).sumW(), sum);
120
121      // Reweight to values used in the paper:
122      // f(D+) = 0.2404
123      // f(D*) = 0.2429
124      for (const string lepton_name : {"lep_minus", "lep_plus"}) {
125        for (const string hadron_name : {"Dplus", "Dstar"}) {
126          const double sf = hadron_name == "Dplus"? (0.2404/fDplus) : (0.2429/fDstar);
127          scale(_h[histo_name(lepton_name, hadron_name, "lep_abs_eta")], sf);
128          scale(_h[histo_name(lepton_name, hadron_name, "D_pt")], sf);
129          normalize(_h[histo_name(lepton_name, hadron_name, "lep_abs_eta") + "_norm"], 1, true);
130          normalize(_h[histo_name(lepton_name, hadron_name, "D_pt") + "_norm"], 1, true);
131        }
132      }
133
134      // The cross-sections from this analysis are not di_lep_minus_Dplus_D_ptfferential
135      for (auto& item : _h) {
136        if (item.first != "CharmSpecies")  barchart(item.second, _s[item.first]);
137      }
138    }
139
140
141  private:
142
143    string histo_name(const string& lepton, const string& hadron, const string& val) {
144      return lepton + "_" + hadron + "_" + val;
145    }
146
147
148    void bookPair(const string& lepton, const string& hadron,
149                  const string& val, unsigned int d) {
150      // absolute
151      string label = histo_name(lepton, hadron, val);
152      book(_h[label], "_"+label, refData(d, 1, 1));
153      book(_s[label], d, 1, 1);
154      // normalised
155      label += "_norm";
156      book(_h[label], "_"+label, refData(d, 1, 2));
157      book(_s[label], d, 1, 2);
158    }
159
160    // Mappting for lepton names
161    const std::unordered_map<int, std::string> _lepton_names = {
162      {PID::ELECTRON, "lep_minus"},
163      {PID::POSITRON, "lep_plus"},
164      {PID::MUON, "lep_minus"},
165      {PID::ANTIMUON, "lep_plus"},
166    };
167
168    // Mapping between pdg id an charm hadron names
169    const std::unordered_map<unsigned int, std::string> _charm_hadron_names = {
170      {PID::DPLUS, "Dplus"},
171      {PID::DSTARPLUS, "Dstar"},
172      {PID::D0, "Dzero"},
173      {PID::DSPLUS, "Ds"},
174      {PID::LAMBDACPLUS, "LambdaC"},
175      {PID::XI0C, "XiCzero"},
176      {PID::XICPLUS, "XiCplus"},
177      {PID::OMEGA0C, "OmegaC"},
178    };
179
180    // Needed to fill the CharmSpecies histograms
181    const std::unordered_map<unsigned int, float> _charm_species_map = {
182      {PID::DPLUS, 1.5},
183      {PID::DSTARPLUS, 0.5},
184      {PID::D0, 2.5},
185      {PID::DSPLUS, 3.5},
186      {PID::LAMBDACPLUS, 4.5},
187      {PID::XI0C, 5.5},
188      {PID::XICPLUS, 6.5},
189      {PID::OMEGA0C, 7.5},
190    };
191
192    // Histogram map
193    map<string, Histo1DPtr> _h;
194    map<string, Estimate1DPtr> _s;
195
196    // Last pT(D) bin extends to infinity (150 only for practical purposes)
197    double _max_D_pt = 150;
198  };
199
200
201  RIVET_DECLARE_PLUGIN(ATLAS_2023_I2628732);
202
203}