rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2019_I1744604

$t$-channel single top-quark differential cross sections and charge ratios at 13 TeV
Experiment: CMS (LHC)
Inspire ID: 1744604
Status: VALIDATED
Authors:
  • Matthias Komm
References:
  • Eur.Phys.J. C80 (2020) 370, 2020
  • DOI:10.1140/epjc/s10052-020-7858-1
  • arXiv: 1907.08330
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • t-channel single top quark and antiquark events at $\sqrt{s} = 13$ TeV; charged lepton (including those from intermediate $\tau$ lepton decays) + two jets selected at particle level

Abstract: A measurement is presented of differential cross sections for $t$-channel single top quark and antiquark production in proton-proton collisions at a centre-of-mass energy of 13 TeV by the CMS experiment at the LHC. From a data set corresponding to an integrated luminosity of 35.9 $\textrm{fb}^{-1}$, events containing one muon or electron and two or three jets are analysed. The cross section is measured as a function of the top quark transverse momentum ($p_\textrm{T}$), rapidity, and polarisation angle, the charged lepton $p_\textrm{T}$ and rapidity, and the $p_\textrm{T}$ of the W boson from the top quark decay. In addition, the charge ratio is measured differentially as a function of the top quark, charged lepton, and W boson kinematic observables. The results are found to be in agreement with standard model predictions using various next-to-leading-order event generators and sets of parton distribution functions. Additionally, the spin asymmetry, sensitive to the top quark polarisation, is determined from the differential distribution of the polarisation angle at parton level to be $0.440 \pm 0.070$, in agreement with the standard model prediction. Rivet: This analysis is to be run separately on t-channel single top quark and top antiquark simulation and the resulting outputs have to be merged for caculating the top quark and antiquark cross section sums and ratios. The particle level selection requires: exactly one dressed lepton (electron or muon; including those from intermediate $\tau$ lepton decay) with $p_{T} > 26$ GeV and $|\eta| < 2.4$; exactly two jets with $p_{T} > 40$ GeV and $|\eta| < 4.7$ that are separated from the dressed lepton by $\Delta R>0.4$. A neutrino candidate reconstructed from the summed transverse momenta of all prompt neutrinos using a W boson mass constraint ($m_\textrm{W}=80.4$ GeV). The four momentum of the top quark is found by summing the dressed lepton, neutrino, and the jet momenta, where the latter is chosen such that a top quark mass closest to 172.5 GeV is obtained. The other jet is taken to be the one that recoils against the W boson in production and used for calculating the top quark polarization angle.

Source code: CMS_2019_I1744604.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/ChargedLeptons.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/MissingMomentum.hh"
  8#include "Rivet/Projections/PromptFinalState.hh"
  9#include "Rivet/Projections/VetoedFinalState.hh"
 10#include "Rivet/Projections/PartonicTops.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// Differential cross sections and charge ratios for 13 TeV t-channel single top-quark production
 16  class CMS_2019_I1744604 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2019_I1744604);
 21
 22
 23    /// @brief Book histograms and initialise projections before the run
 24    void init() override {
 25
 26      // final state of all stable particles
 27      Cut particle_cut = Cuts::abseta < 5.0;
 28      FinalState fs(particle_cut);
 29
 30      // select charged leptons
 31      ChargedLeptons charged_leptons(fs);
 32      // select final state photons for dressed lepton clustering
 33      IdentifiedFinalState photons(fs);
 34      photons.acceptIdPair(PID::PHOTON);
 35
 36      // select final state prompt charged leptons
 37      PromptFinalState prompt_leptons(charged_leptons);
 38      prompt_leptons.acceptMuonDecays(true);
 39      prompt_leptons.acceptTauDecays(true);
 40      // select final state prompt photons
 41      PromptFinalState prompt_photons(photons);
 42      prompt_photons.acceptMuonDecays(true);
 43      prompt_photons.acceptTauDecays(true);
 44
 45      // Dressed leptons from selected prompt charged leptons and photons
 46      Cut lepton_cut = Cuts::abseta < 2.4 and Cuts::pT > 26*GeV;
 47      LeptonFinder dressed_leptons(prompt_leptons, prompt_photons, 0.1, lepton_cut);
 48      declare(dressed_leptons, "LeptonFinder");
 49
 50      // Jets
 51      VetoedFinalState fsForJets(fs);
 52      fsForJets.addVetoOnThisFinalState(dressed_leptons);
 53      declare(FastJets(fsForJets, JetAlg::ANTIKT, 0.4), "Jets"); //< excludes all neutrinos by default
 54
 55      // Neutrinos
 56      IdentifiedFinalState neutrinos(fs);
 57      neutrinos.acceptNeutrinos();
 58      PromptFinalState prompt_neutrinos(neutrinos);
 59      prompt_neutrinos.acceptMuonDecays(true);
 60      prompt_neutrinos.acceptTauDecays(true);
 61      declare(prompt_neutrinos, "Neutrinos");
 62
 63      // Partonic top for differentiating between t and tbar events
 64      declare(PartonicTops(),"TopQuarks");
 65
 66
 67      // book top quark pt histograms
 68      book(_hist_abs_top_pt, "d13-x01-y01"); // absolute cross section
 69      book(_hist_norm_top_pt, "d37-x01-y01");  // normalized cross section
 70      book(_hist_ratio_top_pt, "d59-x01-y01"); // charge ratio
 71      // temporary histograms of absolute top quark and antiquark cross sections
 72      book(_hist_t_top_pt, "t_top_pt", refData(13, 1, 1));
 73      book(_hist_tbar_top_pt, "tbar_top_pt", refData(13, 1, 1));
 74
 75
 76      // book top quark rapidity histograms
 77      book(_hist_abs_top_y, "d15-x01-y01"); // absolute cross section
 78      book(_hist_norm_top_y, "d39-x01-y01"); // normalized cross section
 79      book(_hist_ratio_top_y, "d61-x01-y01"); // charge ratio
 80      // temporary histograms of absolute top quark and antiquark cross sections
 81      book(_hist_t_top_y, "t_top_y", refData(15, 1, 1));
 82      book(_hist_tbar_top_y, "tbar_top_y", refData(15, 1, 1));
 83
 84
 85      // book charged lepton pt histograms
 86      book(_hist_abs_lepton_pt,"d17-x01-y01"); // absolute cross section
 87      book(_hist_norm_lepton_pt,"d41-x01-y01"); // normalized cross section
 88      book(_hist_ratio_lepton_pt,"d63-x01-y01");  // charge ratio
 89      // temporary histograms of absolute top quark and antiquark cross sections
 90      book(_hist_t_lepton_pt,"t_lepton_pt",refData(17, 1, 1));
 91      book(_hist_tbar_lepton_pt,"tbar_lepton_pt",refData(17, 1, 1));
 92
 93
 94      // book charged lepton rapidity histograms
 95      book(_hist_abs_lepton_y, "d19-x01-y01"); // absolute cross section
 96      book(_hist_norm_lepton_y, "d43-x01-y01"); // normalized cross section
 97      book(_hist_ratio_lepton_y, "d65-x01-y01"); // charge ratio
 98      // temporary histograms of absolute top quark and antiquark cross sections
 99      book(_hist_t_lepton_y, "t_lepton_y", refData(19, 1, 1));
100      book(_hist_tbar_lepton_y, "tbar_lepton_y", refData(19, 1, 1));
101
102
103      // book W boson pt histograms
104      book(_hist_abs_w_pt, "d21-x01-y01"); // absolute cross section
105      book(_hist_norm_w_pt, "d45-x01-y01"); // normalized cross section
106      book(_hist_ratio_w_pt, "d67-x01-y01"); // charge ratio
107      // temporary histograms of absolute top quark and antiquark cross sections
108      book(_hist_t_w_pt, "t_w_pt", refData(21, 1, 1));
109      book(_hist_tbar_w_pt, "tbar_w_pt", refData(21, 1, 1));
110
111
112      // book top quark polarization angle histograms
113      book(_hist_abs_top_cos, "d23-x01-y01"); // absolute cross section
114      book(_hist_norm_top_cos, "d47-x01-y01"); // normalized cross section
115      // temporary histograms of absolute top quark and antiquark cross sections
116      book(_hist_t_top_cos, "t_top_cos", refData(23, 1, 1));
117      book(_hist_tbar_top_cos, "tbar_top_cos", refData(23, 1, 1));
118
119    }
120
121
122    /// @brief Perform the per-event analysis
123    void analyze(const Event& event) override {
124      vector<Particle> topQuarks = apply<PartonicTops>(
125        event,
126        "TopQuarks"
127      ).tops();
128
129
130      // skip events with no partonic top quark
131      if (topQuarks.size() != 1) {
132        return;
133      }
134
135      DressedLeptons dressedLeptons = apply<LeptonFinder>(
136        event,
137        "LeptonFinder"
138      ).dressedLeptons();
139
140      // only analyze events with one dressed lepton (muon or electron)
141      if (dressedLeptons.size()!=1) {
142        return;
143      }
144
145      Cut jet_cut((Cuts::abseta < 4.7) and (Cuts::pT > 40.*GeV));
146      vector<Jet> jets = apply<FastJets>(
147        event,
148        "Jets"
149      ).jets(jet_cut);
150
151      // ignore jets that overlap with dressed leptons within dR<0.4
152      Jets cleanedJets;
153      DeltaRLess dRFct(dressedLeptons[0], 0.4);
154      for (const Jet& jet: jets) {
155        if (not dRFct(jet)) cleanedJets.push_back(jet);
156      }
157
158      // select events with exactly two jets
159      if (cleanedJets.size() != 2) {
160        return;
161      }
162
163      Particles neutrinos = apply<PromptFinalState>(event, "Neutrinos").particles();
164      // construct missing transverse momentum by summing over all prompt neutrinos
165      FourMomentum met;
166      for (const Particle& neutrino: neutrinos) {
167        met += neutrino.momentum();
168      }
169
170      /* find unknown pz component of missing transverse momentum by imposing
171         a W boson mass constraint */
172      pair<FourMomentum,FourMomentum> nuMomentum = NuMomentum(
173        dressedLeptons[0].px(), dressedLeptons[0].py(), dressedLeptons[0].pz(),
174        dressedLeptons[0].E(), met.px(), met.py()
175      );
176
177      // define the W boson momentum as the sum of the dressed lepton + neutrino
178      FourMomentum wboson = nuMomentum.first + dressedLeptons[0].momentum();
179
180      /** construct the pseudo top quark momentum by summing the W boson and
181       *  the jet that yields a top quark mass closest to TOPMASS
182       */
183      FourMomentum topQuark(0, 0, 0, 0);
184      int bjetIndex = -1;
185      for (size_t i = 0; i < cleanedJets.size(); ++i) {
186        const auto& jet = cleanedJets[i];
187        FourMomentum topCandidate = jet.momentum() + wboson;
188        if (fabs(topQuark.mass() - TOPMASS) > fabs(topCandidate.mass() - TOPMASS)) {
189          bjetIndex = i;
190          topQuark = topCandidate;
191        }
192      }
193
194      if (bjetIndex < 0) {
195        return;
196      }
197
198      // define the jet used to construct the pseudo top quark as the b jet
199      Jet bjet = cleanedJets[bjetIndex];
200
201      // define the other jet as the spectator jet
202      Jet lightjet = cleanedJets[(bjetIndex + 1) % 2];
203
204      // calculate the cosine of the polarization angle that is defined as the
205      //   angle between the charged lepton and the spectator jet in the top
206      //   quark rest frame
207      LorentzTransform boostToTopFrame = LorentzTransform::mkFrameTransform(topQuark);
208      Vector3 ljetInTopFrame = boostToTopFrame.transform(lightjet.momentum()).vector3().unit();
209      Vector3 leptonInTopFrame = boostToTopFrame.transform(dressedLeptons[0].momentum()).vector3().unit();
210      double polarizationAngle = ljetInTopFrame.dot(leptonInTopFrame);
211
212      // fill the histograms depending on the partonic top quark charge
213      if (topQuarks[0].charge() > 0) {
214        _hist_t_top_pt->fill(topQuark.pt() / GeV);
215        _hist_t_top_y->fill(topQuark.absrapidity());
216        _hist_t_lepton_pt->fill(dressedLeptons[0].pt() / GeV);
217        _hist_t_lepton_y->fill(dressedLeptons[0].absrapidity());
218        _hist_t_w_pt->fill(wboson.pt() / GeV);
219        _hist_t_top_cos->fill(polarizationAngle);
220
221      } else {
222        _hist_tbar_top_pt->fill(topQuark.pt() / GeV);
223        _hist_tbar_top_y->fill(topQuark.absrapidity());
224        _hist_tbar_lepton_pt->fill(dressedLeptons[0].pt() / GeV);
225        _hist_tbar_lepton_y->fill(dressedLeptons[0].absrapidity());
226        _hist_tbar_w_pt->fill(wboson.pt() / GeV);
227        _hist_tbar_top_cos->fill(polarizationAngle);
228      }
229    }
230
231
232    /// @brief Normalise histograms etc., after the run
233    void finalize() override {
234
235      // multiply by 0.5 to average electron/muon decay channels
236      scale(_hist_t_top_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
237      scale(_hist_t_top_y, 0.5 * crossSection() / picobarn / sumOfWeights());
238      scale(_hist_t_lepton_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
239      scale(_hist_t_lepton_y, 0.5 * crossSection() / picobarn / sumOfWeights());
240      scale(_hist_t_w_pt,0.5 * crossSection() / picobarn / sumOfWeights());
241      scale(_hist_t_top_cos, 0.5 * crossSection() / picobarn / sumOfWeights());
242
243      scale(_hist_tbar_top_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
244      scale(_hist_tbar_top_y, 0.5 * crossSection() / picobarn / sumOfWeights());
245      scale(_hist_tbar_lepton_pt,0.5 * crossSection() / picobarn / sumOfWeights());
246      scale(_hist_tbar_lepton_y, 0.5 * crossSection() / picobarn / sumOfWeights());
247      scale(_hist_tbar_w_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
248      scale(_hist_tbar_top_cos, 0.5 * crossSection() /picobarn / sumOfWeights());
249
250      // populate absolute, normalized, and ratio histograms once top quark and
251      // antiquark histograms have been populated
252      if (_hist_t_top_pt->numEntries() > 0 and _hist_tbar_top_pt->numEntries() > 0) {
253        fillAbsHist(_hist_abs_top_pt, _hist_t_top_pt, _hist_tbar_top_pt);
254        fillNormHist(_hist_norm_top_pt, _hist_t_top_pt, _hist_tbar_top_pt);
255        divide(_hist_t_top_pt, _hist_abs_top_pt, _hist_ratio_top_pt);
256
257        fillAbsHist(_hist_abs_top_y, _hist_t_top_y, _hist_tbar_top_y);
258        fillNormHist(_hist_norm_top_y, _hist_t_top_y, _hist_tbar_top_y);
259        divide(_hist_t_top_y, _hist_abs_top_y, _hist_ratio_top_y);
260
261        fillAbsHist(_hist_abs_lepton_pt, _hist_t_lepton_pt, _hist_tbar_lepton_pt);
262        fillNormHist(_hist_norm_lepton_pt, _hist_t_lepton_pt, _hist_tbar_lepton_pt);
263        divide(_hist_t_lepton_pt, _hist_abs_lepton_pt, _hist_ratio_lepton_pt);
264
265        fillAbsHist(_hist_abs_lepton_y, _hist_t_lepton_y, _hist_tbar_lepton_y);
266        fillNormHist(_hist_norm_lepton_y, _hist_t_lepton_y, _hist_tbar_lepton_y);
267        divide(_hist_t_lepton_y, _hist_abs_lepton_y, _hist_ratio_lepton_y);
268
269        fillAbsHist(_hist_abs_w_pt, _hist_t_w_pt, _hist_tbar_w_pt);
270        fillNormHist(_hist_norm_w_pt, _hist_t_w_pt, _hist_tbar_w_pt);
271        divide(_hist_t_w_pt, _hist_abs_w_pt, _hist_ratio_w_pt);
272
273        fillAbsHist(_hist_abs_top_cos, _hist_t_top_cos, _hist_tbar_top_cos);
274        fillNormHist(_hist_norm_top_cos, _hist_t_top_cos, _hist_tbar_top_cos);
275      }
276    }
277
278
279  private:
280
281    // for reconstruction only
282    const double WMASS = 80.399;
283    const double TOPMASS = 172.5;
284
285    // Top quark pt histograms and ratio
286    Histo1DPtr _hist_abs_top_pt;
287    Histo1DPtr _hist_norm_top_pt;
288    Estimate1DPtr _hist_ratio_top_pt;
289    Histo1DPtr _hist_t_top_pt;
290    Histo1DPtr _hist_tbar_top_pt;
291
292    // Top quark rapidity histograms and ratio
293    Histo1DPtr _hist_abs_top_y;
294    Histo1DPtr _hist_norm_top_y;
295    Estimate1DPtr _hist_ratio_top_y;
296    Histo1DPtr _hist_t_top_y;
297    Histo1DPtr _hist_tbar_top_y;
298
299    // Charged lepton pt histograms and ratio
300    Histo1DPtr _hist_abs_lepton_pt;
301    Histo1DPtr _hist_norm_lepton_pt;
302    Estimate1DPtr _hist_ratio_lepton_pt;
303    Histo1DPtr _hist_t_lepton_pt;
304    Histo1DPtr _hist_tbar_lepton_pt;
305
306    // Charged lepton rapidity histograms and ratio
307    Histo1DPtr _hist_abs_lepton_y;
308    Histo1DPtr _hist_norm_lepton_y;
309    Estimate1DPtr _hist_ratio_lepton_y;
310    Histo1DPtr _hist_t_lepton_y;
311    Histo1DPtr _hist_tbar_lepton_y;
312
313    // W boson pt histograms and ratio
314    Histo1DPtr _hist_abs_w_pt;
315    Histo1DPtr _hist_norm_w_pt;
316    Estimate1DPtr _hist_ratio_w_pt;
317    Histo1DPtr _hist_t_w_pt;
318    Histo1DPtr _hist_tbar_w_pt;
319
320    // Top quark polarization angle histograms
321    Histo1DPtr _hist_abs_top_cos;
322    Histo1DPtr _hist_norm_top_cos;
323    Histo1DPtr _hist_t_top_cos;
324    Histo1DPtr _hist_tbar_top_cos;
325
326
327    /// @brief helper function to fill absolute cross section histograms
328    void fillAbsHist(
329      Histo1DPtr& hist_abs, const Histo1DPtr& hist_t,
330      const Histo1DPtr& hist_tbar
331    ) {
332      (*hist_abs) += (*hist_t);
333      (*hist_abs) += (*hist_tbar);
334    }
335
336    /// @brief helper function to fill normalized cross section histograms
337    void fillNormHist(
338      Histo1DPtr& hist_norm, const Histo1DPtr& hist_t,
339      const Histo1DPtr& hist_tbar
340    ) {
341      (*hist_norm) += (*hist_t);
342      (*hist_norm) += (*hist_tbar);
343      hist_norm->normalize();
344    }
345
346
347    /** @brief helper function to solve for the unknown neutrino pz momentum
348     *  using a W boson mass constraint
349     */
350    pair<FourMomentum,FourMomentum> NuMomentum(
351      double pxlep, double pylep, double pzlep,
352      double elep, double metpx, double metpy
353    ) {
354
355      FourMomentum result(0, 0, 0, 0);
356      FourMomentum result2(0, 0, 0, 0);
357
358      double misET2 = (metpx * metpx + metpy * metpy);
359      double mu = (WMASS * WMASS) / 2 + metpx * pxlep + metpy * pylep;
360      double a  = (mu * pzlep) / (elep * elep - pzlep * pzlep);
361      double a2 = pow(a, 2);
362
363      double b  = (pow(elep, 2.) * (misET2) - pow(mu, 2.))
364                  / (pow(elep, 2) - pow(pzlep, 2));
365
366      double pz1(0), pz2(0), pznu(0), pznu2(0);
367
368      FourMomentum p4W_rec;
369      FourMomentum p4b_rec;
370      FourMomentum p4Top_rec;
371      FourMomentum p4lep_rec;
372
373      p4lep_rec.setXYZE(pxlep, pylep, pzlep, elep);
374
375      FourMomentum p40_rec(0, 0, 0, 0);
376
377      // there are two real solutions
378      if (a2 - b > 0 ) {
379        double root = sqrt(a2 - b);
380        pz1 = a + root;
381        pz2 = a - root;
382
383        pznu = pz1;
384        pznu2 = pz2;
385
386        // first solution is the one with the smallest |pz|
387        if (fabs(pz1) > fabs(pz2)) {
388          pznu = pz2;
389          pznu2 = pz1;
390        }
391
392        double Enu = sqrt(misET2 + pznu * pznu);
393        double Enu2 = sqrt(misET2 + pznu2 * pznu2);
394
395        result.setXYZE(metpx, metpy, pznu, Enu);
396        result2.setXYZE(metpx, metpy, pznu2, Enu2);
397
398      } else {
399
400        // there are only complex solutions; set pz=0 and vary px/py such
401        // that mT=mW while keeping px^2+py^2 close to the original pT^2
402        double ptlep = sqrt(pxlep * pxlep + pylep * pylep);
403
404        double EquationA = 1;
405        double EquationB = -3 * pylep * WMASS / (ptlep);
406
407        double EquationC = WMASS * WMASS * (2 * pylep * pylep) / (ptlep * ptlep)
408                           + WMASS * WMASS
409                           - 4 * pxlep * pxlep * pxlep * metpx / (ptlep * ptlep)
410                           - 4 * pxlep * pxlep * pylep * metpy / (ptlep * ptlep);
411
412        double EquationD = 4 * pxlep * pxlep * WMASS * metpy / (ptlep)
413                           - pylep * WMASS * WMASS * WMASS / ptlep;
414
415        vector<double> solutions = EquationSolve(EquationA, EquationB, EquationC, EquationD);
416
417        vector<double> solutions2 = EquationSolve(EquationA, -EquationB, EquationC, -EquationD);
418
419        double deltaMin = 14000 * 14000;
420        double zeroValue = -WMASS * WMASS / (4 * pxlep);
421        double minPx = 0;
422        double minPy = 0;
423
424        for ( size_t i = 0; i < solutions.size(); ++i) {
425          if (solutions[i] < 0) continue;
426          double p_x = (solutions[i] * solutions[i] - WMASS * WMASS) / (4 * pxlep);
427          double p_y = (WMASS * WMASS * pylep
428                        + 2 * pxlep * pylep * p_x
429                        - WMASS * ptlep * solutions[i]
430                        ) / (2 * pxlep * pxlep);
431          double Delta2 = (p_x - metpx) * (p_x - metpx) + (p_y - metpy) * (p_y - metpy);
432
433          if (Delta2 < deltaMin && Delta2 > 0) {
434            deltaMin = Delta2;
435            minPx = p_x;
436            minPy = p_y;
437          }
438
439        }
440
441        for ( size_t i = 0; i < solutions2.size(); ++i) {
442          if (solutions2[i] < 0) continue;
443          double p_x = (solutions2[i] * solutions2[i] - WMASS * WMASS) / (4 * pxlep);
444          double p_y = (WMASS * WMASS * pylep
445                        + 2 * pxlep * pylep * p_x
446                        + WMASS * ptlep * solutions2[i]
447                       ) / (2 * pxlep * pxlep);
448          double Delta2 = (p_x - metpx) * (p_x - metpx) + (p_y - metpy) * (p_y - metpy);
449          if (Delta2 < deltaMin && Delta2 > 0) {
450            deltaMin = Delta2;
451            minPx = p_x;
452            minPy = p_y;
453          }
454        }
455
456        double pyZeroValue = (WMASS * WMASS * pxlep + 2 * pxlep * pylep * zeroValue);
457        double delta2ZeroValue = (zeroValue - metpx) * (zeroValue - metpx)
458                                 + (pyZeroValue - metpy) * (pyZeroValue - metpy);
459
460        if (deltaMin == 14000 * 14000) {
461          return make_pair(result, result2);
462        }
463
464        if (delta2ZeroValue < deltaMin) {
465          deltaMin = delta2ZeroValue;
466          minPx = zeroValue;
467          minPy = pyZeroValue;
468        }
469
470
471        double mu_Minimum = (WMASS * WMASS) / 2 + minPx * pxlep + minPy * pylep;
472        double a_Minimum  = (mu_Minimum * pzlep) /
473                            (elep * elep - pzlep * pzlep);
474        pznu = a_Minimum;
475
476        double Enu = sqrt(minPx * minPx + minPy * minPy + pznu * pznu);
477        result.setXYZE(minPx, minPy, pznu , Enu);
478      }
479      return make_pair(result, result2);
480    }
481
482
483    /// @brief helper function find root of the cubic equation a*x^3 + b*x^2 + c*x + d = 0
484    vector<double> EquationSolve(
485      double a, double b,
486      double c, double d
487    ) {
488      vector<double> result;
489
490      complex<double> x1;
491      complex<double> x2;
492      complex<double> x3;
493
494      double q = (3 * a * c - b * b) / (9 * a * a);
495      double r = (9 * a * b * c - 27 * a * a * d - 2 * b * b * b
496                 ) / (54 * a * a * a);
497      double Delta = q * q * q + r * r;
498
499      complex<double> s;
500      complex<double> t;
501
502      double rho = 0;
503      double theta = 0;
504
505      if (Delta <= 0) {
506        rho = sqrt(-(q * q * q));
507
508        theta = acos(r / rho);
509
510        s = polar<double>(sqrt(-q), theta / 3.0);
511        t = polar<double>(sqrt(-q), -theta / 3.0);
512      }
513
514      if (Delta > 0) {
515        s = complex<double>(cbrt(r + sqrt(Delta)), 0);
516        t = complex<double>(cbrt(r - sqrt(Delta)), 0);
517      }
518
519      complex<double> i(0, 1.0);
520
521
522      x1 = s + t + complex<double>(-b / (3.0 * a), 0);
523
524      x2 = (s + t) * complex<double>(-0.5, 0)
525           - complex<double>(b / (3.0 * a), 0)
526           + (s - t) * i * complex<double>(sqrt(3) / 2.0, 0);
527
528      x3 = (s + t) * complex<double>(-0.5, 0)
529           - complex<double>(b / (3.0 * a), 0)
530           - (s - t) * i * complex<double>(sqrt(3) / 2.0, 0);
531
532      if (fabs(x1.imag()) < 0.0001) result.push_back(x1.real());
533      if (fabs(x2.imag()) < 0.0001) result.push_back(x2.real());
534      if (fabs(x3.imag()) < 0.0001) result.push_back(x3.real());
535
536      return result;
537    }
538
539  };
540
541
542  RIVET_DECLARE_PLUGIN(CMS_2019_I1744604);
543
544}