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 protoin-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
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/UnstableParticles.hh"
#include <unordered_map>

namespace Rivet {


  /// @brief W+D production in pp at 13 TeV
  class ATLAS_2023_I2628732 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2023_I2628732);

    // Book histograms and initialise before the run
    void init() {
      const FinalState fs;

      // lepton kinematic cuts
      Cut cuts(Cuts::pT > 30*GeV && Cuts::abseta < 2.5);

      // Get photons to dress leptons
      FinalState photons(Cuts::abspid == PID::PHOTON);

      // Get dressed leptons
      IdentifiedFinalState lepids(fs, {{PID::ELECTRON, PID::POSITRON, PID::MUON, PID::ANTIMUON}});
      PromptFinalState leptons(lepids, false);
      DressedLeptons dressedleptons(photons, leptons, 0.1, cuts, true);
      declare(dressedleptons, "DressedLeptons");

      // unstable final-state for Ds
      declare(UnstableParticles(), "UFS");

      // Fiducial cross sections vs species:
      // D+ and D* production fractions can be reweighted to the world average values.
      // We can calculate them for each MC sample separately with the `CharmSpecies`
      // histogram from the Rivet routine, e.g.:
      //   - f(D+) = CharmSpecies->GetBinContent(2) / CharmSpecies->Integral(2,8)
      //   - f(D*) = CharmSpecies->GetBinContent(1) / CharmSpecies->Integral(2,8)
      book(_h["CharmSpecies"], "_CharmSpecies", 8, 0, 8);

      // Differential cross sections per bin
      bookPair("lep_minus", "Dplus", "D_pt",         3);
      bookPair("lep_plus",  "Dplus", "D_pt",         4);
      bookPair("lep_minus", "Dplus", "lep_abs_eta",  5);
      bookPair("lep_plus",  "Dplus", "lep_abs_eta",  6);
      bookPair("lep_minus", "Dstar", "D_pt",         7);
      bookPair("lep_plus",  "Dstar", "D_pt",         8);
      bookPair("lep_minus", "Dstar", "lep_abs_eta",  9);
      bookPair("lep_plus",  "Dstar", "lep_abs_eta", 10);
    }

    /// Perform the per-event analysis
    void analyze(const Event &event) {
      // Retrieve the dressed electrons
      const Particles &signal_leptons = apply<DressedLeptons>(event, "DressedLeptons").particlesByPt();
      if (signal_leptons.size() != 1)  vetoEvent;

      const Particle &lepton = signal_leptons[0];
      const std::string lepton_name = _lepton_names.at(lepton.pid());

      // Get the charm hadrons
      const UnstableParticles &ufs = apply<UnstableFinalState>(event, "UFS");
      std::unordered_map<unsigned int, Particles> particles;

      // Loop over particles
      for (const Particle &p : ufs.particles()) {
        const int id = p.abspid();
        const double pt = p.pT() / GeV;
        const double eta = p.abseta();
        if (_charm_hadron_names.count(id) && pt > 8.0 && eta < 2.2) {
          particles[id].push_back(p);
        }
      }

      // Fill histograms
      for (auto &kv : particles) {
        const unsigned int absPdgId = kv.first;
        const std::string hadron_name = _charm_hadron_names.at(absPdgId);

        for (auto &p : kv.second) {
          // Weight: +1 for OS and -1 for SS
          float charm_charge = (absPdgId == 421) ? p.pid() : p.charge();
          double weight = (charm_charge * lepton.charge() < 0) ? +1.0 : -1.0;

          // Fill charm species for production fraction reweighting
          _h["CharmSpecies"]->fill(_charm_species_map.at(absPdgId), weight);

          // Fill only D+ and D* histograms
          if (absPdgId != PID::DPLUS && absPdgId != PID::DSTARPLUS)  continue;

          // pT(D) overflow
          // Last pT(D) bin extends to infinity (150 only for practical purposes)
          double pt = p.pT() / GeV;
          if (pt >= _max_D_pt)  pt = _max_D_pt - 10;

          // Fill histograms
          _h[histo_name(lepton_name, hadron_name, "lep_abs_eta")]->fill(std::abs(lepton.eta()), weight);
          _h[histo_name(lepton_name, hadron_name, "D_pt")]->fill(pt, weight);
          _h[histo_name(lepton_name, hadron_name, "lep_abs_eta") + "_norm"]->fill(std::abs(lepton.eta()), weight);
          _h[histo_name(lepton_name, hadron_name, "D_pt") + "_norm"]->fill(pt, weight);
        }
      }
    }

    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    /// Finalize
    void finalize() {

      scale(_h, crossSectionPerEvent());

      // D+ and D* production fractions
      const double sum = _h["CharmSpecies"]->integral(false);
      const double fDplus = safediv(_h["CharmSpecies"]->bin(1).sumW(), sum);
      const double fDstar = safediv(_h["CharmSpecies"]->bin(0).sumW(), sum);

      // Reweight to values used in the paper:
      // f(D+) = 0.2404
      // f(D*) = 0.2429
      for (const string lepton_name : {"lep_minus", "lep_plus"}) {
        for (const string hadron_name : {"Dplus", "Dstar"}) {
          const double sf = hadron_name == "Dplus"? (0.2404/fDplus) : (0.2429/fDstar);
          scale(_h[histo_name(lepton_name, hadron_name, "lep_abs_eta")], sf);
          scale(_h[histo_name(lepton_name, hadron_name, "D_pt")], sf);
          normalize(_h[histo_name(lepton_name, hadron_name, "lep_abs_eta") + "_norm"], 1, true);
          normalize(_h[histo_name(lepton_name, hadron_name, "D_pt") + "_norm"], 1, true);
        }
      }

      // The cross-sections from this are analysis are not differential
      for (auto& item : _h) {
        if (item.first != "CharmSpecies")  barchart(item.second, _s[item.first]);
      }
    }


  private:

    string histo_name(const string& lepton, const string& hadron, const string& val) {
      return lepton + "_" + hadron + "_" + val;
    }


    void bookPair(const string& lepton, const string& hadron,
                  const string& val, unsigned int d) {
      // absolute
      string label = histo_name(lepton, hadron, val);
      book(_h[label], "_"+label, refData(d, 1, 1));
      book(_s[label], d, 1, 1);
      // normalised
      label += "_norm";
      book(_h[label], "_"+label, refData(d, 1, 2));
      book(_s[label], d, 1, 2);
    }

    // Mappting for lepton names
    const std::unordered_map<int, std::string> _lepton_names = {
      {PID::ELECTRON, "lep_minus"},
      {PID::POSITRON, "lep_plus"},
      {PID::MUON, "lep_minus"},
      {PID::ANTIMUON, "lep_plus"},
    };

    // Mapping between pdg id an charm hadron names
    const std::unordered_map<unsigned int, std::string> _charm_hadron_names = {
      {PID::DPLUS, "Dplus"},
      {PID::DSTARPLUS, "Dstar"},
      {PID::D0, "Dzero"},
      {PID::DSPLUS, "Ds"},
      {PID::LAMBDACPLUS, "LambdaC"},
      {PID::XI0C, "XiCzero"},
      {PID::XICPLUS, "XiCplus"},
      {PID::OMEGA0C, "OmegaC"},
    };

    // Needed to fill the CharmSpecies histograms
    const std::unordered_map<unsigned int, float> _charm_species_map = {
      {PID::DPLUS, 1.5},
      {PID::DSTARPLUS, 0.5},
      {PID::D0, 2.5},
      {PID::DSPLUS, 3.5},
      {PID::LAMBDACPLUS, 4.5},
      {PID::XI0C, 5.5},
      {PID::XICPLUS, 6.5},
      {PID::OMEGA0C, 7.5},
    };

    // Histogram map
    map<string, Histo1DPtr> _h;
    map<string, Scatter2DPtr> _s;

    // Last pT(D) bin extends to infinity (150 only for practical purposes)
    double _max_D_pt = 150;
  };


  RIVET_DECLARE_PLUGIN(ATLAS_2023_I2628732);

}