Rivet analyses referenceCMS_2019_I1744604$t$-channel single top-quark differential cross sections and charge ratios at 13 TeVExperiment: CMS (LHC) Inspire ID: 1744604 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
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}
|