Rivet analyses
W+D production in pp at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 2628732
Status: VALIDATED
Authors: - Miha Muskinja
References: - arXiv: 2302.00336 - Expt page: ATLAS-STDM-2019-22
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± → K∓π±π± and D*± → D0π± → (K∓π±)π±, where pT(e, μ) > 30GeV, |η(e, μ)| < 2.5, pT(D) > 8GeV, and |η(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 σfidOS − SS(W− + D+) = 50.2 ± 0.2 (stat.) −2.3+2.4 (syst.) pb, σfidOS − SS(W+ + D−) = 48.2 ± 0.2 (stat.) −2.2+2.3 (syst.) pb, σfidOS − SS(W− + D*+) = 51.1 ± 0.4 (stat.) −1.8+1.9 (syst.) pb, and σfidOS − SS(W+ + D*−) = 50.0 ± 0.4 (stat.) −1.8+1.9 (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 − s̄ quark asymmetry and is found to be Rc± = 0.971 ± 0.006 (stat.) ±0.011 (syst.).
Source
code:ATLAS_2023_I2628732.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/LeptonFinder.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, TauDecaysAs::NONPROMPT);
LeptonFinder dressedleptons(leptons, photons, 0.1, cuts);
declare(dressedleptons, "LeptonFinder");
// 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<LeptonFinder>(event, "LeptonFinder").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(lepton.abseta(), 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(lepton.abseta(), 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(2).sumW(), sum);
const double fDstar = safediv(_h["CharmSpecies"]->bin(1).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 analysis are not di_lep_minus_Dplus_D_ptfferential
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, Estimate1DPtr> _s;
// Last pT(D) bin extends to infinity (150 only for practical purposes)
double _max_D_pt = 150;
};
RIVET_DECLARE_PLUGIN(ATLAS_2023_I2628732);
}