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>(event, "LeptonFinder").dressedLeptons();
136
137 // only analyze events with one dressed lepton (muon or electron)
138 if (dressedLeptons.size()!=1) {
139 return;
140 }
141
142 Cut jet_cut((Cuts::abseta < 4.7) and (Cuts::pT > 40.*GeV));
143 vector<Jet> jets = apply<FastJets>(
144 event,
145 "Jets"
146 ).jets(jet_cut);
147
148 // ignore jets that overlap with dressed leptons within dR<0.4
149 Jets cleanedJets;
150 DeltaRLess dRFct(dressedLeptons[0], 0.4);
151 for (const Jet& jet: jets) {
152 if (not dRFct(jet)) cleanedJets.push_back(jet);
153 }
154
155 // select events with exactly two jets
156 if (cleanedJets.size() != 2) {
157 return;
158 }
159
160 Particles neutrinos = apply<PromptFinalState>(event, "Neutrinos").particles();
161 // construct missing transverse momentum by summing over all prompt neutrinos
162 FourMomentum met;
163 for (const Particle& neutrino: neutrinos) {
164 met += neutrino.momentum();
165 }
166
167 /* find unknown pz component of missing transverse momentum by imposing
168 a W boson mass constraint */
169 pair<FourMomentum,FourMomentum> nuMomentum = NuMomentum(
170 dressedLeptons[0].px(), dressedLeptons[0].py(), dressedLeptons[0].pz(),
171 dressedLeptons[0].E(), met.px(), met.py()
172 );
173
174 // define the W boson momentum as the sum of the dressed lepton + neutrino
175 FourMomentum wboson = nuMomentum.first + dressedLeptons[0].momentum();
176
177 /** construct the pseudo top quark momentum by summing the W boson and
178 * the jet that yields a top quark mass closest to TOPMASS
179 */
180 FourMomentum topQuark(0, 0, 0, 0);
181 int bjetIndex = -1;
182 for (size_t i = 0; i < cleanedJets.size(); ++i) {
183 const auto& jet = cleanedJets[i];
184 FourMomentum topCandidate = jet.momentum() + wboson;
185 if (fabs(topQuark.mass() - TOPMASS) > fabs(topCandidate.mass() - TOPMASS)) {
186 bjetIndex = i;
187 topQuark = topCandidate;
188 }
189 }
190
191 if (bjetIndex < 0) {
192 return;
193 }
194
195 // define the jet used to construct the pseudo top quark as the b jet
196 Jet bjet = cleanedJets[bjetIndex];
197
198 // define the other jet as the spectator jet
199 Jet lightjet = cleanedJets[(bjetIndex + 1) % 2];
200
201 // calculate the cosine of the polarization angle that is defined as the
202 // angle between the charged lepton and the spectator jet in the top
203 // quark rest frame
204 LorentzTransform boostToTopFrame = LorentzTransform::mkFrameTransform(topQuark);
205 Vector3 ljetInTopFrame = boostToTopFrame.transform(lightjet.momentum()).vector3().unit();
206 Vector3 leptonInTopFrame = boostToTopFrame.transform(dressedLeptons[0].momentum()).vector3().unit();
207 double polarizationAngle = ljetInTopFrame.dot(leptonInTopFrame);
208
209 // fill the histograms depending on the partonic top quark charge
210 if (topQuarks[0].charge() > 0) {
211 _hist_t_top_pt->fill(topQuark.pt() / GeV);
212 _hist_t_top_y->fill(topQuark.absrapidity());
213 _hist_t_lepton_pt->fill(dressedLeptons[0].pt() / GeV);
214 _hist_t_lepton_y->fill(dressedLeptons[0].absrapidity());
215 _hist_t_w_pt->fill(wboson.pt() / GeV);
216 _hist_t_top_cos->fill(polarizationAngle);
217
218 } else {
219 _hist_tbar_top_pt->fill(topQuark.pt() / GeV);
220 _hist_tbar_top_y->fill(topQuark.absrapidity());
221 _hist_tbar_lepton_pt->fill(dressedLeptons[0].pt() / GeV);
222 _hist_tbar_lepton_y->fill(dressedLeptons[0].absrapidity());
223 _hist_tbar_w_pt->fill(wboson.pt() / GeV);
224 _hist_tbar_top_cos->fill(polarizationAngle);
225 }
226 }
227
228
229 /// @brief Normalise histograms etc., after the run
230 void finalize() override {
231
232 // multiply by 0.5 to average electron/muon decay channels
233 scale(_hist_t_top_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
234 scale(_hist_t_top_y, 0.5 * crossSection() / picobarn / sumOfWeights());
235 scale(_hist_t_lepton_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
236 scale(_hist_t_lepton_y, 0.5 * crossSection() / picobarn / sumOfWeights());
237 scale(_hist_t_w_pt,0.5 * crossSection() / picobarn / sumOfWeights());
238 scale(_hist_t_top_cos, 0.5 * crossSection() / picobarn / sumOfWeights());
239
240 scale(_hist_tbar_top_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
241 scale(_hist_tbar_top_y, 0.5 * crossSection() / picobarn / sumOfWeights());
242 scale(_hist_tbar_lepton_pt,0.5 * crossSection() / picobarn / sumOfWeights());
243 scale(_hist_tbar_lepton_y, 0.5 * crossSection() / picobarn / sumOfWeights());
244 scale(_hist_tbar_w_pt, 0.5 * crossSection() / picobarn / sumOfWeights());
245 scale(_hist_tbar_top_cos, 0.5 * crossSection() /picobarn / sumOfWeights());
246
247 // populate absolute, normalized, and ratio histograms once top quark and
248 // antiquark histograms have been populated
249 if (_hist_t_top_pt->numEntries() > 0 and _hist_tbar_top_pt->numEntries() > 0) {
250 fillAbsHist(_hist_abs_top_pt, _hist_t_top_pt, _hist_tbar_top_pt);
251 fillNormHist(_hist_norm_top_pt, _hist_t_top_pt, _hist_tbar_top_pt);
252 divide(_hist_t_top_pt, _hist_abs_top_pt, _hist_ratio_top_pt);
253
254 fillAbsHist(_hist_abs_top_y, _hist_t_top_y, _hist_tbar_top_y);
255 fillNormHist(_hist_norm_top_y, _hist_t_top_y, _hist_tbar_top_y);
256 divide(_hist_t_top_y, _hist_abs_top_y, _hist_ratio_top_y);
257
258 fillAbsHist(_hist_abs_lepton_pt, _hist_t_lepton_pt, _hist_tbar_lepton_pt);
259 fillNormHist(_hist_norm_lepton_pt, _hist_t_lepton_pt, _hist_tbar_lepton_pt);
260 divide(_hist_t_lepton_pt, _hist_abs_lepton_pt, _hist_ratio_lepton_pt);
261
262 fillAbsHist(_hist_abs_lepton_y, _hist_t_lepton_y, _hist_tbar_lepton_y);
263 fillNormHist(_hist_norm_lepton_y, _hist_t_lepton_y, _hist_tbar_lepton_y);
264 divide(_hist_t_lepton_y, _hist_abs_lepton_y, _hist_ratio_lepton_y);
265
266 fillAbsHist(_hist_abs_w_pt, _hist_t_w_pt, _hist_tbar_w_pt);
267 fillNormHist(_hist_norm_w_pt, _hist_t_w_pt, _hist_tbar_w_pt);
268 divide(_hist_t_w_pt, _hist_abs_w_pt, _hist_ratio_w_pt);
269
270 fillAbsHist(_hist_abs_top_cos, _hist_t_top_cos, _hist_tbar_top_cos);
271 fillNormHist(_hist_norm_top_cos, _hist_t_top_cos, _hist_tbar_top_cos);
272 }
273 }
274
275
276 private:
277
278 // for reconstruction only
279 const double WMASS = 80.399;
280 const double TOPMASS = 172.5;
281
282 // Top quark pt histograms and ratio
283 Histo1DPtr _hist_abs_top_pt;
284 Histo1DPtr _hist_norm_top_pt;
285 Estimate1DPtr _hist_ratio_top_pt;
286 Histo1DPtr _hist_t_top_pt;
287 Histo1DPtr _hist_tbar_top_pt;
288
289 // Top quark rapidity histograms and ratio
290 Histo1DPtr _hist_abs_top_y;
291 Histo1DPtr _hist_norm_top_y;
292 Estimate1DPtr _hist_ratio_top_y;
293 Histo1DPtr _hist_t_top_y;
294 Histo1DPtr _hist_tbar_top_y;
295
296 // Charged lepton pt histograms and ratio
297 Histo1DPtr _hist_abs_lepton_pt;
298 Histo1DPtr _hist_norm_lepton_pt;
299 Estimate1DPtr _hist_ratio_lepton_pt;
300 Histo1DPtr _hist_t_lepton_pt;
301 Histo1DPtr _hist_tbar_lepton_pt;
302
303 // Charged lepton rapidity histograms and ratio
304 Histo1DPtr _hist_abs_lepton_y;
305 Histo1DPtr _hist_norm_lepton_y;
306 Estimate1DPtr _hist_ratio_lepton_y;
307 Histo1DPtr _hist_t_lepton_y;
308 Histo1DPtr _hist_tbar_lepton_y;
309
310 // W boson pt histograms and ratio
311 Histo1DPtr _hist_abs_w_pt;
312 Histo1DPtr _hist_norm_w_pt;
313 Estimate1DPtr _hist_ratio_w_pt;
314 Histo1DPtr _hist_t_w_pt;
315 Histo1DPtr _hist_tbar_w_pt;
316
317 // Top quark polarization angle histograms
318 Histo1DPtr _hist_abs_top_cos;
319 Histo1DPtr _hist_norm_top_cos;
320 Histo1DPtr _hist_t_top_cos;
321 Histo1DPtr _hist_tbar_top_cos;
322
323
324 /// @brief helper function to fill absolute cross section histograms
325 void fillAbsHist(
326 Histo1DPtr& hist_abs, const Histo1DPtr& hist_t,
327 const Histo1DPtr& hist_tbar
328 ) {
329 (*hist_abs) += (*hist_t);
330 (*hist_abs) += (*hist_tbar);
331 }
332
333 /// @brief helper function to fill normalized cross section histograms
334 void fillNormHist(
335 Histo1DPtr& hist_norm, const Histo1DPtr& hist_t,
336 const Histo1DPtr& hist_tbar
337 ) {
338 (*hist_norm) += (*hist_t);
339 (*hist_norm) += (*hist_tbar);
340 hist_norm->normalize();
341 }
342
343
344 /** @brief helper function to solve for the unknown neutrino pz momentum
345 * using a W boson mass constraint
346 */
347 pair<FourMomentum,FourMomentum> NuMomentum(
348 double pxlep, double pylep, double pzlep,
349 double elep, double metpx, double metpy
350 ) {
351
352 FourMomentum result(0, 0, 0, 0);
353 FourMomentum result2(0, 0, 0, 0);
354
355 double misET2 = (metpx * metpx + metpy * metpy);
356 double mu = (WMASS * WMASS) / 2 + metpx * pxlep + metpy * pylep;
357 double a = (mu * pzlep) / (elep * elep - pzlep * pzlep);
358 double a2 = pow(a, 2);
359
360 double b = (pow(elep, 2.) * (misET2) - pow(mu, 2.))
361 / (pow(elep, 2) - pow(pzlep, 2));
362
363 double pz1(0), pz2(0), pznu(0), pznu2(0);
364
365 FourMomentum p4W_rec;
366 FourMomentum p4b_rec;
367 FourMomentum p4Top_rec;
368 FourMomentum p4lep_rec;
369
370 p4lep_rec.setXYZE(pxlep, pylep, pzlep, elep);
371
372 FourMomentum p40_rec(0, 0, 0, 0);
373
374 // there are two real solutions
375 if (a2 - b > 0 ) {
376 double root = sqrt(a2 - b);
377 pz1 = a + root;
378 pz2 = a - root;
379
380 pznu = pz1;
381 pznu2 = pz2;
382
383 // first solution is the one with the smallest |pz|
384 if (fabs(pz1) > fabs(pz2)) {
385 pznu = pz2;
386 pznu2 = pz1;
387 }
388
389 double Enu = sqrt(misET2 + pznu * pznu);
390 double Enu2 = sqrt(misET2 + pznu2 * pznu2);
391
392 result.setXYZE(metpx, metpy, pznu, Enu);
393 result2.setXYZE(metpx, metpy, pznu2, Enu2);
394
395 } else {
396
397 // there are only complex solutions; set pz=0 and vary px/py such
398 // that mT=mW while keeping px^2+py^2 close to the original pT^2
399 double ptlep = sqrt(pxlep * pxlep + pylep * pylep);
400
401 double EquationA = 1;
402 double EquationB = -3 * pylep * WMASS / (ptlep);
403
404 double EquationC = WMASS * WMASS * (2 * pylep * pylep) / (ptlep * ptlep)
405 + WMASS * WMASS
406 - 4 * pxlep * pxlep * pxlep * metpx / (ptlep * ptlep)
407 - 4 * pxlep * pxlep * pylep * metpy / (ptlep * ptlep);
408
409 double EquationD = 4 * pxlep * pxlep * WMASS * metpy / (ptlep)
410 - pylep * WMASS * WMASS * WMASS / ptlep;
411
412 vector<double> solutions = EquationSolve(EquationA, EquationB, EquationC, EquationD);
413
414 vector<double> solutions2 = EquationSolve(EquationA, -EquationB, EquationC, -EquationD);
415
416 double deltaMin = 14000 * 14000;
417 double zeroValue = -WMASS * WMASS / (4 * pxlep);
418 double minPx = 0;
419 double minPy = 0;
420
421 for ( size_t i = 0; i < solutions.size(); ++i) {
422 if (solutions[i] < 0) continue;
423 double p_x = (solutions[i] * solutions[i] - WMASS * WMASS) / (4 * pxlep);
424 double p_y = (WMASS * WMASS * pylep
425 + 2 * pxlep * pylep * p_x
426 - WMASS * ptlep * solutions[i]
427 ) / (2 * pxlep * pxlep);
428 double Delta2 = (p_x - metpx) * (p_x - metpx) + (p_y - metpy) * (p_y - metpy);
429
430 if (Delta2 < deltaMin && Delta2 > 0) {
431 deltaMin = Delta2;
432 minPx = p_x;
433 minPy = p_y;
434 }
435
436 }
437
438 for ( size_t i = 0; i < solutions2.size(); ++i) {
439 if (solutions2[i] < 0) continue;
440 double p_x = (solutions2[i] * solutions2[i] - WMASS * WMASS) / (4 * pxlep);
441 double p_y = (WMASS * WMASS * pylep
442 + 2 * pxlep * pylep * p_x
443 + WMASS * ptlep * solutions2[i]
444 ) / (2 * pxlep * pxlep);
445 double Delta2 = (p_x - metpx) * (p_x - metpx) + (p_y - metpy) * (p_y - metpy);
446 if (Delta2 < deltaMin && Delta2 > 0) {
447 deltaMin = Delta2;
448 minPx = p_x;
449 minPy = p_y;
450 }
451 }
452
453 double pyZeroValue = (WMASS * WMASS * pxlep + 2 * pxlep * pylep * zeroValue);
454 double delta2ZeroValue = (zeroValue - metpx) * (zeroValue - metpx)
455 + (pyZeroValue - metpy) * (pyZeroValue - metpy);
456
457 if (deltaMin == 14000 * 14000) {
458 return make_pair(result, result2);
459 }
460
461 if (delta2ZeroValue < deltaMin) {
462 deltaMin = delta2ZeroValue;
463 minPx = zeroValue;
464 minPy = pyZeroValue;
465 }
466
467
468 double mu_Minimum = (WMASS * WMASS) / 2 + minPx * pxlep + minPy * pylep;
469 double a_Minimum = (mu_Minimum * pzlep) /
470 (elep * elep - pzlep * pzlep);
471 pznu = a_Minimum;
472
473 double Enu = sqrt(minPx * minPx + minPy * minPy + pznu * pznu);
474 result.setXYZE(minPx, minPy, pznu , Enu);
475 }
476 return make_pair(result, result2);
477 }
478
479
480 /// @brief helper function find root of the cubic equation a*x^3 + b*x^2 + c*x + d = 0
481 vector<double> EquationSolve(
482 double a, double b,
483 double c, double d
484 ) {
485 vector<double> result;
486
487 complex<double> x1;
488 complex<double> x2;
489 complex<double> x3;
490
491 double q = (3 * a * c - b * b) / (9 * a * a);
492 double r = (9 * a * b * c - 27 * a * a * d - 2 * b * b * b
493 ) / (54 * a * a * a);
494 double Delta = q * q * q + r * r;
495
496 complex<double> s;
497 complex<double> t;
498
499 double rho = 0;
500 double theta = 0;
501
502 if (Delta <= 0) {
503 rho = sqrt(-(q * q * q));
504
505 theta = acos(r / rho);
506
507 s = polar<double>(sqrt(-q), theta / 3.0);
508 t = polar<double>(sqrt(-q), -theta / 3.0);
509 }
510
511 if (Delta > 0) {
512 s = complex<double>(cbrt(r + sqrt(Delta)), 0);
513 t = complex<double>(cbrt(r - sqrt(Delta)), 0);
514 }
515
516 complex<double> i(0, 1.0);
517
518
519 x1 = s + t + complex<double>(-b / (3.0 * a), 0);
520
521 x2 = (s + t) * complex<double>(-0.5, 0)
522 - complex<double>(b / (3.0 * a), 0)
523 + (s - t) * i * complex<double>(sqrt(3) / 2.0, 0);
524
525 x3 = (s + t) * complex<double>(-0.5, 0)
526 - complex<double>(b / (3.0 * a), 0)
527 - (s - t) * i * complex<double>(sqrt(3) / 2.0, 0);
528
529 if (fabs(x1.imag()) < 0.0001) result.push_back(x1.real());
530 if (fabs(x2.imag()) < 0.0001) result.push_back(x2.real());
531 if (fabs(x3.imag()) < 0.0001) result.push_back(x3.real());
532
533 return result;
534 }
535
536 };
537
538
539 RIVET_DECLARE_PLUGIN(CMS_2019_I1744604);
540
541}
|