Rivet analyses referenceATLAS_2023_I2628732W+D production in pp at 13 TeVExperiment: ATLAS (LHC) Inspire ID: 2628732 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
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// -*- 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}
|