rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2021_I1978840

Measurement of $\mathrm{W}^{\pm}\gamma$ differential cross sections in proton-proton collisions at $\sqrt{s}=13\,\mathrm{TeV}$ and effective field theory constraints
Experiment: CMS (LHC)
Inspire ID: 1978840
Status: VALIDATED
Authors:
  • cms-pag-conveners-smp@cern.ch
  • Andrew Gilbert
References:
  • Public page: SMP-20-005
  • arxiv:2111.13948
  • Phys. Rev. D 105 (2022) 052003
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp to $\mathrm{W}^{\pm}\gamma$ at $\sqrt{s}=13$ TeV.

Differential cross section measurements of $\mathrm{W}^{\pm}\gamma$ production in proton-proton collisions at $\sqrt{s} = 13\,\mathrm{TeV}$ are presented. The data set used in this study was collected with the CMS detector at the CERN LHC in 2016--2018 with an integrated luminosity of $138\,\mathrm{fb}^{-1}$. Candidate events containing an electron or muon, a photon, and missing transverse momentum are selected. The measurements are compared with standard model predictions computed at next-to-leading and next-to-next-to-leading orders in perturbative quantum chromodynamics. Constraints on the presence of TeV-scale new physics affecting the $\mathrm{WW}\gamma$ vertex are determined within an effective field theory framework, focusing on the $\mathcal{O}_{3W}$ operator. A simultaneous measurement of the photon transverse momentum and the azimuthal angle of the charged lepton in a special reference frame is performed.

Source code: CMS_2021_I1978840.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Event.hh"
  4#include "Rivet/Math/LorentzTrans.hh"
  5#include "Rivet/Particle.hh"
  6#include "Rivet/Projections/ChargedLeptons.hh"
  7#include "Rivet/Projections/DressedLeptons.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Projections/FinalState.hh"
 10#include "Rivet/Projections/IdentifiedFinalState.hh"
 11#include "Rivet/Projections/MissingMomentum.hh"
 12#include "Rivet/Projections/PromptFinalState.hh"
 13#include "Rivet/Projections/VetoedFinalState.hh"
 14#include "Rivet/Projections/InvisibleFinalState.hh"
 15
 16namespace Rivet {
 17
 18
 19  /// @brief Measurement of W/gamma differential cross sections at sqrt(s) = 13 TeV
 20  class CMS_2021_I1978840 : public Analysis {
 21  public:
 22
 23    struct WGammaRivetVariables {
 24      bool is_wg_gen;
 25
 26      double l0_pt;
 27      double l0_eta;
 28      double l0_phi;
 29      double l0_M;
 30      int l0_q;
 31      unsigned l0_abs_pdgid;
 32
 33      double p0_pt;
 34      double p0_eta;
 35      double p0_phi;
 36      double p0_M;
 37      bool p0_frixione;
 38      double p0_frixione_sum;
 39
 40      double l0p0_dr;
 41      double mt_cluster;
 42
 43      double n0_pt;
 44      double n0_eta;
 45      double n0_phi;
 46      double n0_M;
 47
 48      double met_pt;
 49      double met_phi;
 50
 51      double true_phi;
 52      double true_phi_f;
 53
 54      int n_jets;
 55
 56      WGammaRivetVariables() { resetVars(); }
 57
 58      void resetVars() {
 59        is_wg_gen = false;
 60        l0_pt = 0.;
 61        l0_eta = 0.;
 62        l0_phi = 0.;
 63        l0_M = 0.;
 64        l0_q = 0;
 65        l0_abs_pdgid = 0;
 66        p0_pt = 0.;
 67        p0_eta = 0.;
 68        p0_phi = 0.;
 69        p0_M = 0.;
 70        p0_frixione = false;
 71        n0_pt = 0.;
 72        n0_eta = 0.;
 73        n0_phi = 0.;
 74        n0_M = 0.;
 75        met_pt = 0.;
 76        met_phi = 0.;
 77        true_phi = 0.;
 78        true_phi_f = 0.;
 79        l0p0_dr = 0.;
 80        mt_cluster = 0.;
 81        n_jets = 0;
 82      }
 83    };
 84
 85    struct WGSystem {
 86      int lepton_charge;
 87
 88      FourMomentum wg_system;
 89
 90      FourMomentum c_w_boson;
 91      FourMomentum c_charged_lepton;
 92      FourMomentum c_neutrino;
 93      FourMomentum c_photon;
 94
 95      FourMomentum r_w_boson;
 96      FourMomentum r_charged_lepton;
 97      FourMomentum r_neutrino;
 98      FourMomentum r_photon;
 99
100      WGSystem(Particle const& lep, FourMomentum const& neu, Particle const& pho, bool verbose);
101
102      double phi();
103      double symphi();
104    };
105
106    double photon_iso_dr_ = 0.4;
107    double lepton_pt_cut_ = 30.;
108    double lepton_abs_eta_cut_ = 2.5;
109    double photon_pt_cut_ = 30.;
110    double photon_abs_eta_cut_ = 2.5;
111    double missing_pt_cut_ = 40.;
112    double lepton_photon_dr_cut_ = 0.7;
113    double dressed_lepton_cone_ = 0.1;
114
115    double eft_lepton_pt_cut_ = 80.;
116    double eft_photon_pt_cut_ = 150.;
117
118    double jet_pt_cut_ = 30.;
119    double jet_abs_eta_cut_ = 2.5;
120    double jet_dr_cut_ = 0.4;
121
122    WGammaRivetVariables vars_;
123    map<string, Histo1DPtr> _h;
124
125
126    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2021_I1978840);
127
128
129    void init() {
130      vars_.resetVars();
131
132      FinalState fs;
133      declare(fs, "FinalState");
134
135      // Jets - all final state particles excluding neutrinos
136      VetoedFinalState vfs;
137      vfs.vetoNeutrinos();
138      FastJets fastjets(vfs, FastJets::ANTIKT, 0.4);
139      declare(fastjets, "Jets");
140
141      // Dressed leptons
142      ChargedLeptons charged_leptons(fs);
143      PromptFinalState prompt_leptons(charged_leptons);
144      declare(prompt_leptons, "PromptLeptons");
145
146      IdentifiedFinalState photons(fs);
147      photons.acceptIdPair(PID::PHOTON);
148      PromptFinalState prompt_photons(photons);
149      prompt_photons.acceptMuonDecays(true);
150      prompt_photons.acceptTauDecays(true);
151
152      DressedLeptons dressed_leptons(prompt_photons, prompt_leptons, dressed_lepton_cone_,
153                                     Cuts::open(),
154                                     /*useDecayPhotons*/ false);
155      declare(dressed_leptons, "DressedLeptons");
156
157      // Photons
158      VetoedFinalState vetoed_prompt_photons(prompt_photons);
159      vetoed_prompt_photons.addVetoOnThisFinalState(dressed_leptons);
160      declare(vetoed_prompt_photons, "Photons");
161
162      // Invisibles
163      InvisibleFinalState invisibles(true, true, true);
164      declare(invisibles, "Invisibles");
165
166      // MET
167      declare(MissingMomentum(fs), "MET");
168
169      // Booking of histograms
170      book(_h["baseline_photon_pt"], 1, 1, 1);
171      book(_h["baseline_photon_eta"], 8, 1, 1);
172      book(_h["baseline_leppho_dr"], 15, 1, 1);
173      book(_h["baseline_leppho_deta"], 22, 1, 1);
174      book(_h["baseline_mt_cluster"], 29, 1, 1);
175      book(_h["baseline_njet"], 36, 1, 1);
176      book(_h["raz_leppho_deta"], 40, 1, 1);
177      book(_h["eft_photon_pt_phi_0"], 54, 1, 1);
178      book(_h["eft_photon_pt_phi_1"], 55, 1, 1);
179      book(_h["eft_photon_pt_phi_2"], 56, 1, 1);
180    }
181
182    /// Perform the per-event analysis
183    void analyze(const Event& event) {
184      vars_.resetVars();
185
186      const Particles leptons = apply<FinalState>(event, "DressedLeptons").particlesByPt();
187
188      if (leptons.size() == 0) {
189        vetoEvent;
190      }
191      auto l0 = leptons.at(0);
192
193      const Particles photons = apply<FinalState>(event, "Photons")
194        .particlesByPt(DeltaRGtr(l0, lepton_photon_dr_cut_));
195      if (photons.size() == 0) {
196        vetoEvent;
197      }
198      auto p0 = photons.at(0);
199
200
201      const Particles invisibles = apply<FinalState>(event, "Invisibles").particlesByPt();
202      FourMomentum n0;
203      for (auto const& inv : invisibles) {
204        n0 += inv.momentum();
205      }
206
207      if (invisibles.size() == 0) {
208        vetoEvent;
209      }
210
211      FourMomentum met = apply<MissingMomentum>(event, "MET").missingMomentum();
212      // Redefine the MET
213      met = FourMomentum(met.pt(), met.px(), met.py(), 0.);
214
215      // Filter jets on pT, eta and DR with lepton and photon
216      const Jets jets = apply<FastJets>(event, "Jets").jetsByPt([&](Jet const& j) {
217        return j.pt() > jet_pt_cut_ && j.abseta() < jet_abs_eta_cut_ && deltaR(j, l0) > jet_dr_cut_ && deltaR(j, p0) > jet_dr_cut_;
218      });
219
220      if (leptons.size() >= 1 && photons.size() >= 1 && invisibles.size() >= 1) {
221        // Populate variables
222        vars_.is_wg_gen = true;
223        vars_.l0_pt = l0.pt()/GeV;
224        vars_.l0_eta = l0.eta();
225        vars_.l0_phi = l0.phi(PhiMapping::MINUSPI_PLUSPI);
226        vars_.l0_M = l0.mass()/GeV;
227        vars_.l0_q = l0.charge();
228        vars_.l0_abs_pdgid = l0.abspid();
229        vars_.p0_pt = p0.pt()/GeV;
230        vars_.p0_eta = p0.eta();
231        vars_.p0_phi = p0.phi(PhiMapping::MINUSPI_PLUSPI);
232        vars_.p0_M = p0.mass()/GeV;
233        vars_.n_jets = jets.size();
234        vars_.l0p0_dr = deltaR(l0, p0);
235        vars_.n0_pt = n0.pt()/GeV;
236        vars_.n0_eta = n0.eta();
237        vars_.n0_phi = n0.phi(PhiMapping::MINUSPI_PLUSPI);
238        vars_.n0_M = n0.mass()/GeV;
239        vars_.met_pt = met.pt()/GeV;
240        vars_.met_phi = met.phi(PhiMapping::MINUSPI_PLUSPI);
241
242        // Here we build a list of particles to cluster jets, to be used in the photon isolation.
243        // The selection of particles that we want to veto from this isolation sum is non-trivial:
244        // the leading pT lepton, the leading pT photon that has DeltaR > 0.7 from the leading pT
245        // lepton, the invisibles  and any tau decay products. Therefore, the selection is done
246        // here instead of in the initialise method.
247        Particles finalparts_iso = apply<FinalState>(event, "FinalState").particles();
248        Particles filtered_iso;
249        for (Particle const& p : finalparts_iso) {
250          bool veto_particle = false;
251          if (p.genParticle() == l0.genParticle() || p.genParticle() == p0.genParticle() || p.fromTau()) {
252            veto_particle = true;
253          }
254          for (auto const& inv : invisibles) {
255            if (p.genParticle() == inv.genParticle()) {
256              veto_particle = true;
257            }
258          }
259          if (!veto_particle) {
260            filtered_iso.push_back(p);
261          }
262        }
263        auto proj = getProjection<FastJets>("Jets");
264        proj.reset();
265        proj.calc(filtered_iso);
266        auto jets_iso = proj.jets();
267
268        vars_.p0_frixione = true;
269        double frixione_sum = 0.;
270
271        // Apply Frixione isolation to the photon:
272        auto jparts = sortBy(jets_iso, [&](Jet const& part1, Jet const& part2) {
273          return deltaR(part1, p0) < deltaR(part2, p0);
274        });
275
276        for (auto const& ip : jparts) {
277          double dr = deltaR(ip, p0);
278          if (dr >= photon_iso_dr_) {
279            break;
280          }
281          frixione_sum += ip.pt();
282          if (frixione_sum > (p0.pt() * ((1. - cos(dr)) / (1. - cos(photon_iso_dr_))))) {
283            vars_.p0_frixione = false;
284          }
285        }
286
287        // Now calculate EFT phi observables
288        auto wg_system = WGSystem(l0, n0, p0, false);
289
290        vars_.true_phi = wg_system.phi();
291        vars_.true_phi_f = wg_system.symphi();
292
293        // Calculate mTcluster
294        auto cand1 = l0.momentum() + p0.momentum();
295        auto full_system = cand1 + met;
296        double mTcluster2 = sqr(sqrt(cand1.mass2() + cand1.pt2()) + met.pt()) - full_system.pt2();
297        if (mTcluster2 > 0) {
298          vars_.mt_cluster = sqrt(mTcluster2);
299        } else {
300          vars_.mt_cluster = 0.;
301        }
302
303        if (vars_.l0_pt > lepton_pt_cut_ && fabs(vars_.l0_eta) < lepton_abs_eta_cut_ &&
304            vars_.p0_pt > photon_pt_cut_ && fabs(vars_.p0_eta) < photon_abs_eta_cut_ &&
305            vars_.p0_frixione && vars_.l0p0_dr > lepton_photon_dr_cut_ &&
306            vars_.met_pt > missing_pt_cut_) {
307          _h["baseline_photon_pt"]->fill(vars_.p0_pt / GeV);
308          _h["baseline_photon_eta"]->fill(vars_.p0_eta);
309          _h["baseline_leppho_dr"]->fill(vars_.l0p0_dr);
310          _h["baseline_leppho_deta"]->fill(vars_.l0_eta - vars_.p0_eta);
311          _h["baseline_mt_cluster"]->fill(vars_.mt_cluster / GeV);
312          double fill_n_jets = vars_.n_jets >= 2 ? 2. : double(vars_.n_jets);
313          _h["baseline_njet"]->fill(fill_n_jets);
314          if (vars_.n_jets == 0 && vars_.mt_cluster > 150.) {
315            _h["raz_leppho_deta"]->fill(vars_.l0_eta - vars_.p0_eta);
316          }
317        }
318
319        if (vars_.l0_pt > eft_lepton_pt_cut_ && fabs(vars_.l0_eta) < lepton_abs_eta_cut_ &&
320            vars_.p0_pt > eft_photon_pt_cut_ && fabs(vars_.p0_eta) < photon_abs_eta_cut_ &&
321            vars_.p0_frixione && vars_.l0p0_dr > lepton_photon_dr_cut_ &&
322            vars_.met_pt > missing_pt_cut_ && vars_.n_jets == 0) {
323          double absphi = fabs(vars_.true_phi_f);
324          if (absphi > 0. && absphi <= (PI / 6.)) {
325            _h["eft_photon_pt_phi_0"]->fill(vars_.p0_pt / GeV);
326          } else if (absphi > (PI / 6.) && absphi <= (PI / 3.)) {
327            _h["eft_photon_pt_phi_1"]->fill(vars_.p0_pt / GeV);
328          } else if (absphi > (PI / 3.) && absphi <= (PI / 2.)) {
329            _h["eft_photon_pt_phi_2"]->fill(vars_.p0_pt / GeV);
330          }
331        }
332      }
333    }
334
335
336    void finalize() {
337      double flavor_factor = 3. / 2.;  // account for the fact that tau events are vetoed
338      // Scale according to cross section
339      for (string x :
340             {"baseline_photon_pt", "baseline_photon_eta", "baseline_leppho_dr", "baseline_leppho_deta",
341              "baseline_mt_cluster", "baseline_njet", "raz_leppho_deta", "eft_photon_pt_phi_0",
342              "eft_photon_pt_phi_1", "eft_photon_pt_phi_2"}) {
343        scale(_h[x], flavor_factor * crossSection() / femtobarn / sumOfWeights());
344      }
345
346      // Since these are really 2D, we need to divide by the y bin width:
347      for (string x : {"eft_photon_pt_phi_0", "eft_photon_pt_phi_1", "eft_photon_pt_phi_2"}) {
348        scale(_h[x], 1. / (PI / 6.));
349      }
350    }
351
352  };
353
354
355  CMS_2021_I1978840::WGSystem::WGSystem(Particle const& lep, FourMomentum const& neu,
356                                        Particle const& pho, bool verbose) {
357    lepton_charge = lep.charge3() / 3;
358    wg_system += lep.momentum();
359    wg_system += neu;
360    wg_system += pho.momentum();
361
362    if (verbose) {
363      cout << "> charge:    " << lepton_charge << "\n";
364      cout << "> wg_system: " << wg_system << "\n";
365      cout << "> lepton   : " << lep.momentum() << "\n";
366      cout << "> neutrino : " << neu << "\n";
367      cout << "> photon   : " << pho.momentum() << "\n";
368    }
369
370    auto boost = LorentzTransform::mkFrameTransformFromBeta(wg_system.betaVec());
371
372    c_charged_lepton = boost.transform(lep);
373    c_neutrino = boost.transform(neu);
374    c_photon = boost.transform(pho);
375
376    if (verbose) {
377      cout << "> c_lepton   : " << c_charged_lepton << "\n";
378      cout << "> c_neutrino : " << c_neutrino << "\n";
379      cout << "> c_photon   : " << c_photon << "\n";
380    }
381
382    FourMomentum c_w_boson;
383    c_w_boson += c_charged_lepton;
384    c_w_boson += c_neutrino;
385
386    Vector3 r_uvec = wg_system.vector3().unit();
387    if (verbose) {
388      cout << "> c_w_boson : " << c_w_boson << endl;
389      cout << "> r_uvec   : " << r_uvec << endl;
390    }
391
392    Vector3 z_uvec = c_w_boson.vector3().unit();
393    Vector3 y_uvec = z_uvec.cross(r_uvec).unit();
394    Vector3 x_uvec = y_uvec.cross(z_uvec).unit();
395    if (verbose) {
396      cout << "> x_uvec   : " << x_uvec << endl;
397      cout << "> y_uvec   : " << y_uvec << endl;
398      cout << "> z_uvec   : " << z_uvec << endl;
399    }
400
401    Matrix3 rot_matrix;
402    rot_matrix.setRow(0, x_uvec).setRow(1, y_uvec).setRow(2, z_uvec);
403    auto rotator = LorentzTransform();
404    rotator = rotator.postMult(rot_matrix);
405    if (verbose) {
406      cout << "> rotator   : " << rotator << endl;
407    }
408    r_w_boson = rotator.transform(c_w_boson);
409    r_charged_lepton = rotator.transform(c_charged_lepton);
410    r_neutrino = rotator.transform(c_neutrino);
411    r_photon = rotator.transform(c_photon);
412
413    if (verbose) {
414      cout << "> r_lepton   : " << r_charged_lepton << endl;
415      cout << "> r_neutrino : " << r_neutrino << endl;
416      cout << "> r_photon   : " << r_photon << endl;
417      cout << "> r_w_boson  : " << r_w_boson << endl;
418    }
419  }
420
421  double CMS_2021_I1978840::WGSystem::phi() {
422    double lep_phi = r_charged_lepton.phi(PhiMapping::MINUSPI_PLUSPI);
423    return mapAngleMPiToPi(lepton_charge > 0 ? (lep_phi) : (lep_phi + PI));
424  }
425
426  double CMS_2021_I1978840::WGSystem::symphi() {
427    double lep_phi = phi();
428    if (lep_phi > PI / 2.) {
429      return PI - lep_phi;
430    } else if (lep_phi < -1. * (PI / 2.)) {
431      return -1. * (PI + lep_phi);
432    } else {
433      return lep_phi;
434    }
435  }
436
437
438  RIVET_DECLARE_PLUGIN(CMS_2021_I1978840);
439
440}