]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8175/src/SigmaProcess.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
CID 21256: Uninitialized pointer field (UNINIT_CTOR)
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8175 / src / SigmaProcess.cxx
CommitLineData
c6b60c38 1// SigmaProcess.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2013 Torbjorn Sjostrand.
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6// Function definitions (not found in the header) for the
7// SigmaProcess class, and classes derived from it.
8
9#include "SigmaProcess.h"
10
11namespace Pythia8 {
12
13//==========================================================================
14
15// The SigmaProcess class.
16// Base class for cross sections.
17
18//--------------------------------------------------------------------------
19
20// Constants: could be changed here if desired, but normally should not.
21// These are of technical nature, as described for each.
22
23// Conversion of GeV^{-2} to mb for cross section.
24const double SigmaProcess::CONVERT2MB = 0.389380;
25
26// The sum of outgoing masses must not be too close to the cm energy.
27const double SigmaProcess::MASSMARGIN = 0.1;
28
29// Parameters of momentum rescaling procedure: maximally allowed
30// relative energy error and number of iterations.
31const double SigmaProcess::COMPRELERR = 1e-10;
32const int SigmaProcess::NCOMPSTEP = 10;
33
34//--------------------------------------------------------------------------
35
36// Perform simple initialization and store pointers.
37
38void SigmaProcess::init(Info* infoPtrIn, Settings* settingsPtrIn,
39 ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, BeamParticle* beamAPtrIn,
40 BeamParticle* beamBPtrIn, Couplings* couplingsPtrIn,
41 SigmaTotal* sigmaTotPtrIn, SusyLesHouches* slhaPtrIn) {
42
43 // Store pointers.
44 infoPtr = infoPtrIn;
45 settingsPtr = settingsPtrIn;
46 particleDataPtr = particleDataPtrIn;
47 rndmPtr = rndmPtrIn;
48 beamAPtr = beamAPtrIn;
49 beamBPtr = beamBPtrIn;
50 couplingsPtr = couplingsPtrIn;
51 sigmaTotPtr = sigmaTotPtrIn;
52 slhaPtr = slhaPtrIn;
53
54 // Read out some properties of beams to allow shorthand.
55 idA = (beamAPtr != 0) ? beamAPtr->id() : 0;
56 idB = (beamBPtr != 0) ? beamBPtr->id() : 0;
57 mA = (beamAPtr != 0) ? beamAPtr->m() : 0.;
58 mB = (beamBPtr != 0) ? beamBPtr->m() : 0.;
59 isLeptonA = (beamAPtr != 0) ? beamAPtr->isLepton() : false;
60 isLeptonB = (beamBPtr != 0) ? beamBPtr->isLepton() : false;
61 hasLeptonBeams = isLeptonA || isLeptonB;
62
63 // K factor, multiplying resolved processes. (But not here for MPI.)
64 Kfactor = settingsPtr->parm("SigmaProcess:Kfactor");
65
66 // Maximum incoming quark flavour.
67 nQuarkIn = settingsPtr->mode("PDFinProcess:nQuarkIn");
68
69 // Medium heavy fermion masses set massless or not in ME expressions.
70 mcME = (settingsPtr->flag("SigmaProcess:cMassiveME"))
71 ? particleDataPtr->m0(4) : 0.;
72 mbME = (settingsPtr->flag("SigmaProcess:bMassiveME"))
73 ? particleDataPtr->m0(5) : 0.;
74 mmuME = (settingsPtr->flag("SigmaProcess:muMassiveME"))
75 ? particleDataPtr->m0(13) : 0.;
76 mtauME = (settingsPtr->flag("SigmaProcess:tauMassiveME"))
77 ? particleDataPtr->m0(15) : 0.;
78
79 // Renormalization scale choice.
80 renormScale1 = settingsPtr->mode("SigmaProcess:renormScale1");
81 renormScale2 = settingsPtr->mode("SigmaProcess:renormScale2");
82 renormScale3 = settingsPtr->mode("SigmaProcess:renormScale3");
83 renormScale3VV = settingsPtr->mode("SigmaProcess:renormScale3VV");
84 renormMultFac = settingsPtr->parm("SigmaProcess:renormMultFac");
85 renormFixScale = settingsPtr->parm("SigmaProcess:renormFixScale");
86
87 // Factorization scale choice.
88 factorScale1 = settingsPtr->mode("SigmaProcess:factorScale1");
89 factorScale2 = settingsPtr->mode("SigmaProcess:factorScale2");
90 factorScale3 = settingsPtr->mode("SigmaProcess:factorScale3");
91 factorScale3VV = settingsPtr->mode("SigmaProcess:factorScale3VV");
92 factorMultFac = settingsPtr->parm("SigmaProcess:factorMultFac");
93 factorFixScale = settingsPtr->parm("SigmaProcess:factorFixScale");
94
95 // CP violation parameters for the BSM Higgs sector.
96 higgsH1parity = settingsPtr->mode("HiggsH1:parity");
97 higgsH1eta = settingsPtr->parm("HiggsH1:etaParity");
98 higgsH2parity = settingsPtr->mode("HiggsH2:parity");
99 higgsH2eta = settingsPtr->parm("HiggsH2:etaParity");
100 higgsA3parity = settingsPtr->mode("HiggsA3:parity");
101 higgsA3eta = settingsPtr->parm("HiggsA3:etaParity");
102
103 // If BSM not switched on then H1 should have SM properties.
104 if (!settingsPtr->flag("Higgs:useBSM")){
105 higgsH1parity = 1;
106 higgsH1eta = 0.;
107 }
108
109}
110
111//--------------------------------------------------------------------------
112
113// Set up allowed flux of incoming partons.
114// addBeam: set up PDF's that need to be evaluated for the two beams.
115// addPair: set up pairs of incoming partons from the two beams.
116
117bool SigmaProcess::initFlux() {
118
119 // Reset arrays (in case of several init's in same run).
120 inBeamA.clear();
121 inBeamB.clear();
122 inPair.clear();
123
124 // Read in process-specific channel information.
125 string fluxType = inFlux();
126
127 // Case with g g incoming state.
128 if (fluxType == "gg") {
129 addBeamA(21);
130 addBeamB(21);
131 addPair(21, 21);
132 }
133
134 // Case with q g incoming state.
135 else if (fluxType == "qg") {
136 for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
137 int idNow = (i == 0) ? 21 : i;
138 addBeamA(idNow);
139 addBeamB(idNow);
140 }
141 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
142 if (idNow != 0) {
143 addPair(idNow, 21);
144 addPair(21, idNow);
145 }
146 }
147
148 // Case with q q', q qbar' or qbar qbar' incoming state.
149 else if (fluxType == "qq") {
150 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
151 if (idNow != 0) {
152 addBeamA(idNow);
153 addBeamB(idNow);
154 }
155 for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
156 if (id1Now != 0)
157 for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
158 if (id2Now != 0)
159 addPair(id1Now, id2Now);
160 }
161
162 // Case with q qbar incoming state.
163 else if (fluxType == "qqbarSame") {
164 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
165 if (idNow != 0) {
166 addBeamA(idNow);
167 addBeamB(idNow);
168 }
169 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
170 if (idNow != 0)
171 addPair(idNow, -idNow);
172 }
173
174 // Case with f f', f fbar', fbar fbar' incoming state.
175 else if (fluxType == "ff") {
176 // If beams are leptons then they are also the colliding partons.
177 if ( isLeptonA && isLeptonB ) {
178 addBeamA(idA);
179 addBeamB(idB);
180 addPair(idA, idB);
181 // First beam is lepton and second is hadron.
182 } else if ( isLeptonA ) {
183 addBeamA(idA);
184 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
185 if (idNow != 0) {
186 addBeamB(idNow);
187 addPair(idA, idNow);
188 }
189 // First beam is hadron and second is lepton.
190 } else if ( isLeptonB ) {
191 addBeamB(idB);
192 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
193 if (idNow != 0) {
194 addBeamA(idNow);
195 addPair(idNow, idB);
196 }
197 // Hadron beams gives quarks.
198 } else {
199 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
200 if (idNow != 0) {
201 addBeamA(idNow);
202 addBeamB(idNow);
203 }
204 for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
205 if (id1Now != 0)
206 for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
207 if (id2Now != 0)
208 addPair(id1Now, id2Now);
209 }
210 }
211
212 // Case with f fbar incoming state.
213 else if (fluxType == "ffbarSame") {
214 // If beams are antiparticle pair and leptons then also colliding partons.
215 if ( idA + idB == 0 && isLeptonA ) {
216 addBeamA(idA);
217 addBeamB(idB);
218 addPair(idA, idB);
219 // Else assume both to be hadrons, for better or worse.
220 } else {
221 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
222 if (idNow != 0) {
223 addBeamA(idNow);
224 addBeamB(idNow);
225 }
226 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
227 if (idNow != 0)
228 addPair(idNow, -idNow);
229 }
230 }
231
232 // Case with f fbar' charged(+-1) incoming state.
233 else if (fluxType == "ffbarChg") {
234 // If beams are leptons then also colliding partons.
235 if ( isLeptonA && isLeptonB && abs( particleDataPtr->chargeType(idA)
236 + particleDataPtr->chargeType(idB) ) == 3 ) {
237 addBeamA(idA);
238 addBeamB(idB);
239 addPair(idA, idB);
240 // Hadron beams gives quarks.
241 } else {
242 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
243 if (idNow != 0) {
244 addBeamA(idNow);
245 addBeamB(idNow);
246 }
247 for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
248 if (id1Now != 0)
249 for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
250 if (id2Now != 0 && id1Now * id2Now < 0
251 && (abs(id1Now) + abs(id2Now))%2 == 1) addPair(id1Now, id2Now);
252 }
253 }
254
255 // Case with f fbar' generic incoming state.
256 else if (fluxType == "ffbar") {
257 // If beams are leptons then also colliding partons.
258 if (isLeptonA && isLeptonB && idA * idB < 0) {
259 addBeamA(idA);
260 addBeamB(idB);
261 addPair(idA, idB);
262 // Hadron beams gives quarks.
263 } else {
264 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
265 if (idNow != 0) {
266 addBeamA(idNow);
267 addBeamB(idNow);
268 }
269 for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now)
270 if (id1Now != 0)
271 for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now)
272 if (id2Now != 0 && id1Now * id2Now < 0)
273 addPair(id1Now, id2Now);
274 }
275 }
276
277 // Case with f gamma incoming state.
278 else if (fluxType == "fgm") {
279 // Fermion from incoming side A.
280 if ( isLeptonA ) {
281 addBeamA(idA);
282 addPair(idA, 22);
283 } else {
284 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
285 if (idNow != 0) {
286 addBeamA(idNow);
287 addPair(idNow, 22);
288 }
289 }
290 // Fermion from incoming side B.
291 if ( isLeptonB ) {
292 addBeamB( idB);
293 addPair(22, idB);
294 } else {
295 for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow)
296 if (idNow != 0) {
297 addBeamB(idNow);
298 addPair(22, idNow);
299 }
300 }
301 // Photons in the beams.
302 addBeamA(22);
303 addBeamB(22);
304 }
305
306 // Case with gamma gamma incoming state.
307 else if (fluxType == "ggm") {
308 addBeamA(21);
309 addBeamA(22);
310 addBeamB(21);
311 addBeamB(22);
312 addPair(21, 22);
313 addPair(22, 21);
314 }
315
316 // Case with gamma gamma incoming state.
317 else if (fluxType == "gmgm") {
318 addBeamA(22);
319 addBeamB(22);
320 addPair(22, 22);
321 }
322
323 // Unrecognized fluxType is bad sign. Else done.
324 else {
325 infoPtr->errorMsg("Error in SigmaProcess::initFlux: "
326 "unrecognized inFlux type", fluxType);
327 return false;
328 }
329 return true;
330
331}
332
333//--------------------------------------------------------------------------
334
335// Convolute matrix-element expression(s) with parton flux and K factor.
336
337double SigmaProcess::sigmaPDF() {
338
339 // Evaluate and store the required parton densities.
340 for (int j = 0; j < sizeBeamA(); ++j)
341 inBeamA[j].pdf = beamAPtr->xfHard( inBeamA[j].id, x1Save, Q2FacSave);
342 for (int j = 0; j < sizeBeamB(); ++j)
343 inBeamB[j].pdf = beamBPtr->xfHard( inBeamB[j].id, x2Save, Q2FacSave);
344
345 // Loop over allowed incoming channels.
346 sigmaSumSave = 0.;
347 for (int i = 0; i < sizePair(); ++i) {
348
349 // Evaluate hard-scattering cross section. Include K factor.
350 inPair[i].pdfSigma = Kfactor
351 * sigmaHatWrap(inPair[i].idA, inPair[i].idB);
352
353 // Multiply by respective parton densities.
354 for (int j = 0; j < sizeBeamA(); ++j)
355 if (inPair[i].idA == inBeamA[j].id) {
356 inPair[i].pdfA = inBeamA[j].pdf;
357 inPair[i].pdfSigma *= inBeamA[j].pdf;
358 break;
359 }
360 for (int j = 0; j < sizeBeamB(); ++j)
361 if (inPair[i].idB == inBeamB[j].id) {
362 inPair[i].pdfB = inBeamB[j].pdf;
363 inPair[i].pdfSigma *= inBeamB[j].pdf;
364 break;
365 }
366
367 // Sum for all channels.
368 sigmaSumSave += inPair[i].pdfSigma;
369 }
370
371 // Done.
372 return sigmaSumSave;
373
374}
375
376//--------------------------------------------------------------------------
377
378// Select incoming parton channel and extract parton densities (resolved).
379
380void SigmaProcess::pickInState(int id1in, int id2in) {
381
382 // Multiparton interactions: partons already selected.
383 if (id1in != 0 && id2in != 0) {
384 id1 = id1in;
385 id2 = id2in;
386 }
387
388 // Pick channel. Extract channel flavours and pdf's.
389 double sigmaRand = sigmaSumSave * rndmPtr->flat();
390 for (int i = 0; i < sizePair(); ++i) {
391 sigmaRand -= inPair[i].pdfSigma;
392 if (sigmaRand <= 0.) {
393 id1 = inPair[i].idA;
394 id2 = inPair[i].idB;
395 pdf1Save = inPair[i].pdfA;
396 pdf2Save = inPair[i].pdfB;
397 break;
398 }
399 }
400
401}
402
403//--------------------------------------------------------------------------
404
405// Calculate incoming modified masses and four-vectors for matrix elements.
406
407bool SigmaProcess::setupForMEin() {
408
409 // Initially assume it will work out to set up modified kinematics.
410 bool allowME = true;
411
412 // Correct incoming c, b, mu and tau to be massive or not.
413 mME[0] = 0.;
414 int id1Tmp = abs(id1);
415 if (id1Tmp == 4) mME[0] = mcME;
416 if (id1Tmp == 5) mME[0] = mbME;
417 if (id1Tmp == 13) mME[0] = mmuME;
418 if (id1Tmp == 15) mME[0] = mtauME;
419 mME[1] = 0.;
420 int id2Tmp = abs(id2);
421 if (id2Tmp == 4) mME[1] = mcME;
422 if (id2Tmp == 5) mME[1] = mbME;
423 if (id2Tmp == 13) mME[1] = mmuME;
424 if (id2Tmp == 15) mME[1] = mtauME;
425
426 // If kinematically impossible return to massless case, but set error.
427 if (mME[0] + mME[1] >= mH) {
428 mME[0] = 0.;
429 mME[1] = 0.;
430 allowME = false;
431 }
432
433 // Do incoming two-body kinematics for massless or massive cases.
434 if (mME[0] == 0. && mME[1] == 0.) {
435 pME[0] = 0.5 * mH * Vec4( 0., 0., 1., 1.);
436 pME[1] = 0.5 * mH * Vec4( 0., 0., -1., 1.);
437 } else {
438 double e0 = 0.5 * (mH * mH + mME[0] * mME[0] - mME[1] * mME[1]) / mH;
439 double pz0 = sqrtpos(e0 * e0 - mME[0] * mME[0]);
440 pME[0] = Vec4( 0., 0., pz0, e0);
441 pME[1] = Vec4( 0., 0., -pz0, mH - e0);
442 }
443
444 // Done.
445 return allowME;
446
447}
448
449//--------------------------------------------------------------------------
450
451// Evaluate weight for W decay distribution in t -> W b -> f fbar b.
452
453double SigmaProcess::weightTopDecay( Event& process, int iResBeg,
454 int iResEnd) {
455
456 // If not pair W d/s/b and mother t then return unit weight.
457 if (iResEnd - iResBeg != 1) return 1.;
458 int iW1 = iResBeg;
459 int iB2 = iResBeg + 1;
460 int idW1 = process[iW1].idAbs();
461 int idB2 = process[iB2].idAbs();
462 if (idW1 != 24) {
463 swap(iW1, iB2);
464 swap(idW1, idB2);
465 }
466 if (idW1 != 24 || (idB2 != 1 && idB2 != 3 && idB2 != 5)) return 1.;
467 int iT = process[iW1].mother1();
468 if (iT <= 0 || process[iT].idAbs() != 6) return 1.;
469
470 // Find sign-matched order of W decay products.
471 int iF = process[iW1].daughter1();
472 int iFbar = process[iW1].daughter2();
473 if (iFbar - iF != 1) return 1.;
474 if (process[iT].id() * process[iF].id() < 0) swap(iF, iFbar);
475
476 // Weight and maximum weight.
477 double wt = (process[iT].p() * process[iFbar].p())
478 * (process[iF].p() * process[iB2].p());
479 double wtMax = ( pow4(process[iT].m()) - pow4(process[iW1].m()) ) / 8.;
480
481 // Done.
482 return wt / wtMax;
483
484}
485
486//--------------------------------------------------------------------------
487
488// Evaluate weight for Z0/W+- decay distributions in H -> Z0/W+ Z0/W- -> 4f
489// and H -> gamma Z0 -> gamma f fbar.
490
491double SigmaProcess::weightHiggsDecay( Event& process, int iResBeg,
492 int iResEnd) {
493
494 // If not pair Z0 Z0, W+ W- or gamma Z0 then return unit weight.
495 if (iResEnd - iResBeg != 1) return 1.;
496 int iZW1 = iResBeg;
497 int iZW2 = iResBeg + 1;
498 int idZW1 = process[iZW1].id();
499 int idZW2 = process[iZW2].id();
500 if (idZW1 < 0 || idZW2 == 22) {
501 swap(iZW1, iZW2);
502 swap(idZW1, idZW2);
503 }
504 if ( (idZW1 != 23 || idZW2 != 23) && (idZW1 != 24 || idZW2 != -24)
505 && (idZW1 != 22 || idZW2 != 23) ) return 1.;
506
507 // If mother is not Higgs then return unit weight.
508 int iH = process[iZW1].mother1();
509 if (iH <= 0) return 1.;
510 int idH = process[iH].id();
511 if (idH != 25 && idH != 35 && idH !=36) return 1.;
512
513 // H -> gamma Z0 -> gamma f fbar is 1 + cos^2(theta) in Z rest frame.
514 if (idZW1 == 22) {
515 int i5 = process[iZW2].daughter1();
516 int i6 = process[iZW2].daughter2();
517 double pgmZ = process[iZW1].p() * process[iZW2].p();
518 double pgm5 = process[iZW1].p() * process[i5].p();
519 double pgm6 = process[iZW1].p() * process[i6].p();
520 return (pow2(pgm5) + pow2(pgm6)) / pow2(pgmZ);
521 }
522
523 // Parameters depend on Higgs type: H0(H_1), H^0(H_2) or A^0(H_3).
524 int higgsParity = higgsH1parity;
525 double higgsEta = higgsH1eta;
526 if (idH == 35) {
527 higgsParity = higgsH2parity;
528 higgsEta = higgsH2eta;
529 } else if (idH == 36) {
530 higgsParity = higgsA3parity;
531 higgsEta = higgsA3eta;
532 }
533
534 // Option with isotropic decays.
535 if (higgsParity == 0) return 1.;
536
537 // Maximum and initial weight.
538 double wtMax = pow4(process[iH].m());
539 double wt = wtMax;
540
541 // Find sign-matched order of Z0/W+- decay products.
542 int i3 = process[iZW1].daughter1();
543 int i4 = process[iZW1].daughter2();
544 if (process[i3].id() < 0) swap( i3, i4);
545 int i5 = process[iZW2].daughter1();
546 int i6 = process[iZW2].daughter2();
547 if (process[i5].id() < 0) swap( i5, i6);
548
549 // Evaluate four-vector products and find masses..
550 double p35 = 2. * process[i3].p() * process[i5].p();
551 double p36 = 2. * process[i3].p() * process[i6].p();
552 double p45 = 2. * process[i4].p() * process[i5].p();
553 double p46 = 2. * process[i4].p() * process[i6].p();
554 double p34 = 2. * process[i3].p() * process[i4].p();
555 double p56 = 2. * process[i5].p() * process[i6].p();
556 double mZW1 = process[iZW1].m();
557 double mZW2 = process[iZW2].m();
558
559 // For mixed CP states need epsilon product and gauge boson masses.
560 double epsilonProd = 0.;
561 if (higgsParity == 3) {
562 double p[4][4];
563 for (int i = 0; i < 4; ++i) {
564 int ii = i3;
565 if (i == 1) ii = i4;
566 if (i == 2) ii = i5;
567 if (i == 3) ii = i6;
568 p[i][0] = process[ii].e();
569 p[i][1] = process[ii].px();
570 p[i][2] = process[ii].py();
571 p[i][3] = process[ii].pz();
572 }
573 epsilonProd
574 = p[0][0]*p[1][1]*p[2][2]*p[3][3] - p[0][0]*p[1][1]*p[2][3]*p[3][2]
575 - p[0][0]*p[1][2]*p[2][1]*p[3][3] + p[0][0]*p[1][2]*p[2][3]*p[3][1]
576 + p[0][0]*p[1][3]*p[2][1]*p[3][2] - p[0][0]*p[1][3]*p[2][2]*p[3][1]
577 - p[0][1]*p[1][0]*p[2][2]*p[3][3] + p[0][1]*p[1][0]*p[2][3]*p[3][2]
578 + p[0][1]*p[1][2]*p[2][0]*p[3][3] - p[0][1]*p[1][2]*p[2][3]*p[3][0]
579 - p[0][1]*p[1][3]*p[2][0]*p[3][2] + p[0][1]*p[1][3]*p[2][2]*p[3][0]
580 + p[0][2]*p[1][0]*p[2][1]*p[3][3] - p[0][2]*p[1][0]*p[2][3]*p[3][1]
581 - p[0][2]*p[1][1]*p[2][0]*p[3][3] + p[0][2]*p[1][1]*p[2][3]*p[3][0]
582 + p[0][2]*p[1][3]*p[2][0]*p[3][1] - p[0][2]*p[1][3]*p[2][1]*p[3][0]
583 - p[0][3]*p[1][0]*p[2][1]*p[3][2] + p[0][3]*p[1][0]*p[2][2]*p[3][1]
584 + p[0][3]*p[1][1]*p[2][0]*p[3][2] - p[0][3]*p[1][1]*p[2][2]*p[3][0]
585 - p[0][3]*p[1][2]*p[2][0]*p[3][1] + p[0][3]*p[1][2]*p[2][1]*p[3][0];
586 }
587
588 // Z0 Z0 decay: vector and axial couplings of two fermion pairs.
589 if (idZW1 == 23) {
590 double vf1 = couplingsPtr->vf(process[i3].idAbs());
591 double af1 = couplingsPtr->af(process[i3].idAbs());
592 double vf2 = couplingsPtr->vf(process[i5].idAbs());
593 double af2 = couplingsPtr->af(process[i5].idAbs());
594 double va12asym = 4. * vf1 * af1 * vf2 * af2
595 / ( (vf1*vf1 + af1*af1) * (vf2*vf2 + af2*af2) );
596 double etaMod = higgsEta / pow2( particleDataPtr->m0(23) );
597
598 // Normal CP-even decay.
599 if (higgsParity == 1) wt = 8. * (1. + va12asym) * p35 * p46
600 + 8. * (1. - va12asym) * p36 * p45;
601
602 // CP-odd decay (normal for A0(H_3)).
603 else if (higgsParity == 2) wt = ( pow2(p35 + p46)
604 + pow2(p36 + p45) - 2. * p34 * p56
605 - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
606 + va12asym * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) )
607 / (1. + va12asym);
608
609 // Mixed CP states.
610 else wt = 32. * ( 0.25 * ( (1. + va12asym) * p35 * p46
611 + (1. - va12asym) * p36 * p45 ) - 0.5 * etaMod * epsilonProd
612 * ( (1. + va12asym) * (p35 + p46) - (1. - va12asym) * (p36 + p45) )
613 + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
614 - 2. * pow2(p35 * p46 - p36 * p45)
615 + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
616 + va12asym * p34 * p56 * (p35 + p36 - p45 - p46)
617 * (p35 + p45 - p36 - p46) ) ) / ( 1. + 2. * etaMod * mZW1 * mZW2
618 + 2. * pow2(etaMod * mZW1 * mZW2) * (1. + va12asym) );
619
620 // W+ W- decay.
621 } else if (idZW1 == 24) {
622 double etaMod = higgsEta / pow2( particleDataPtr->m0(24) );
623
624 // Normal CP-even decay.
625 if (higgsParity == 1) wt = 16. * p35 * p46;
626
627 // CP-odd decay (normal for A0(H_3)).
628 else if (higgsParity == 2) wt = 0.5 * ( pow2(p35 + p46)
629 + pow2(p36 + p45) - 2. * p34 * p56
630 - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56)
631 + (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) );
632
633 // Mixed CP states.
634 else wt = 32. * ( 0.25 * 2. * p35 * p46
635 - 0.5 * etaMod * epsilonProd * 2. * (p35 + p46)
636 + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56)
637 - 2. * pow2(p35 * p46 - p36 * p45)
638 + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45))
639 + p34 * p56 * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ) )
640 / ( 1. * 2. * etaMod * mZW1 * mZW2 + 2. * pow2(etaMod * mZW1 * mZW2) );
641 }
642
643 // Done.
644 return wt / wtMax;
645
646}
647
648//==========================================================================
649
650// The Sigma1Process class.
651// Base class for resolved 2 -> 1 cross sections; derived from SigmaProcess.
652
653//--------------------------------------------------------------------------
654
655// Wrapper to sigmaHat, to (a) store current incoming flavours,
656// (b) convert from GeV^-2 to mb where required, and
657// (c) convert from |M|^2 to d(sigmaHat)/d(tHat) where required.
658
659double Sigma1Process::sigmaHatWrap(int id1in, int id2in) {
660
661 id1 = id1in;
662 id2 = id2in;
663 double sigmaTmp = sigmaHat();
664 if (convertM2()) {
665 sigmaTmp /= 2. * sH;
666 // Convert 2 * pi * delta(p^2 - m^2) to Breit-Wigner with same area.
667 int idTmp = resonanceA();
668 double mTmp = particleDataPtr->m0(idTmp);
669 double GamTmp = particleDataPtr->mWidth(idTmp);
670 sigmaTmp *= 2. * mTmp * GamTmp / ( pow2(sH - mTmp * mTmp)
671 + pow2(mTmp * GamTmp) );
672 }
673 if (convert2mb()) sigmaTmp *= CONVERT2MB;
674 return sigmaTmp;
675
676}
677
678//--------------------------------------------------------------------------
679
680// Input and complement kinematics for resolved 2 -> 1 process.
681
682void Sigma1Process::store1Kin( double x1in, double x2in, double sHin) {
683
684 // Default value only sensible for these processes.
685 swapTU = false;
686
687 // Incoming parton momentum fractions and sHat.
688 x1Save = x1in;
689 x2Save = x2in;
690 sH = sHin;
691 mH = sqrt(sH);
692 sH2 = sH * sH;
693
694 // Different options for renormalization scale, but normally sHat.
695 Q2RenSave = renormMultFac * sH;
696 if (renormScale1 == 2) Q2RenSave = renormFixScale;
697
698 // Different options for factorization scale, but normally sHat.
699 Q2FacSave = factorMultFac * sH;
700 if (factorScale1 == 2) Q2FacSave = factorFixScale;
701
702 // Evaluate alpha_strong and alpha_EM.
703 alpS = couplingsPtr->alphaS(Q2RenSave);
704 alpEM = couplingsPtr->alphaEM(Q2RenSave);
705
706}
707
708//--------------------------------------------------------------------------
709
710// Calculate modified masses and four-vectors for matrix elements.
711
712bool Sigma1Process::setupForME() {
713
714 // Common initial-state handling.
715 bool allowME = setupForMEin();
716
717 // Final state trivial here.
718 mME[2] = mH;
719 pME[2] = Vec4( 0., 0., 0., mH);
720
721 // Done.
722 return allowME;
723
724}
725
726//==========================================================================
727
728// The Sigma2Process class.
729// Base class for resolved 2 -> 2 cross sections; derived from SigmaProcess.
730
731//--------------------------------------------------------------------------
732
733// Input and complement kinematics for resolved 2 -> 2 process.
734
735void Sigma2Process::store2Kin( double x1in, double x2in, double sHin,
736 double tHin, double m3in, double m4in, double runBW3in, double runBW4in) {
737
738 // Default ordering of particles 3 and 4.
739 swapTU = false;
740
741 // Incoming parton momentum fractions.
742 x1Save = x1in;
743 x2Save = x2in;
744
745 // Incoming masses and their squares.
746 bool masslessKin = (id3Mass() == 0) && (id4Mass() == 0);
747 if (masslessKin) {
748 m3 = 0.;
749 m4 = 0.;
750 } else {
751 m3 = m3in;
752 m4 = m4in;
753 }
754 mSave[3] = m3;
755 mSave[4] = m4;
756 s3 = m3 * m3;
757 s4 = m4 * m4;
758
759 // Standard Mandelstam variables and their squares.
760 sH = sHin;
761 tH = tHin;
762 uH = (masslessKin) ? -(sH + tH) : s3 + s4 - (sH + tH);
763 mH = sqrt(sH);
764 sH2 = sH * sH;
765 tH2 = tH * tH;
766 uH2 = uH * uH;
767
768 // The nominal Breit-Wigner factors with running width.
769 runBW3 = runBW3in;
770 runBW4 = runBW4in;
771
772 // Calculate squared transverse momentum.
773 pT2 = (masslessKin) ? tH * uH / sH : (tH * uH - s3 * s4) / sH;
774
775 // Special case: pick scale as if 2 -> 1 process in disguise.
776 if (isSChannel()) {
777
778 // Different options for renormalization scale, but normally sHat.
779 Q2RenSave = renormMultFac * sH;
780 if (renormScale1 == 2) Q2RenSave = renormFixScale;
781
782 // Different options for factorization scale, but normally sHat.
783 Q2FacSave = factorMultFac * sH;
784 if (factorScale1 == 2) Q2FacSave = factorFixScale;
785
786 // Normal case with "true" 2 -> 2.
787 } else {
788
789 // Different options for renormalization scale.
790 if (masslessKin) Q2RenSave = (renormScale2 < 4) ? pT2 : sH;
791 else if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
792 else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
793 else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
794 else Q2RenSave = sH;
795 Q2RenSave *= renormMultFac;
796 if (renormScale2 == 5) Q2RenSave = renormFixScale;
797
798 // Different options for factorization scale.
799 if (masslessKin) Q2FacSave = (factorScale2 < 4) ? pT2 : sH;
800 else if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
801 else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
802 else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
803 else Q2FacSave = sH;
804 Q2FacSave *= factorMultFac;
805 if (factorScale2 == 5) Q2FacSave = factorFixScale;
806 }
807
808 // Evaluate alpha_strong and alpha_EM.
809 alpS = couplingsPtr->alphaS(Q2RenSave);
810 alpEM = couplingsPtr->alphaEM(Q2RenSave);
811
812}
813
814//--------------------------------------------------------------------------
815
816// As above, special kinematics for multiparton interactions.
817
818void Sigma2Process::store2KinMPI( double x1in, double x2in,
819 double sHin, double tHin, double uHin, double alpSin, double alpEMin,
820 bool needMasses, double m3in, double m4in) {
821
822 // Default ordering of particles 3 and 4.
823 swapTU = false;
824
825 // Incoming x values.
826 x1Save = x1in;
827 x2Save = x2in;
828
829 // Standard Mandelstam variables and their squares.
830 sH = sHin;
831 tH = tHin;
832 uH = uHin;
833 mH = sqrt(sH);
834 sH2 = sH * sH;
835 tH2 = tH * tH;
836 uH2 = uH * uH;
837
838 // Strong and electroweak couplings.
839 alpS = alpSin;
840 alpEM = alpEMin;
841
842 // Assume vanishing masses. (Will be modified in final kinematics.)
843 m3 = 0.;
844 s3 = 0.;
845 m4 = 0.;
846 s4 = 0.;
847 sHBeta = sH;
848
849 // Scattering angle.
850 cosTheta = (tH - uH) / sH;
851 sinTheta = 2. * sqrtpos( tH * uH ) / sH;
852
853 // In some cases must use masses and redefine meaning of tHat and uHat.
854 if (needMasses) {
855 m3 = m3in;
856 s3 = m3 * m3;
857 m4 = m4in;
858 s4 = m4 * m4;
859 sHMass = sH - s3 - s4;
860 sHBeta = sqrtpos(sHMass*sHMass - 4. * s3 * s4);
861 tH = -0.5 * (sHMass - sHBeta * cosTheta);
862 uH = -0.5 * (sHMass + sHBeta * cosTheta);
863 tH2 = tH * tH;
864 uH2 = uH * uH;
865 }
866
867 // pT2 with masses (at this stage) included.
868 pT2Mass = 0.25 * sHBeta * pow2(sinTheta);
869
870}
871
872//--------------------------------------------------------------------------
873
874// Perform kinematics for a multiparton interaction, including a rescattering.
875
876bool Sigma2Process::final2KinMPI( int i1Res, int i2Res, Vec4 p1Res, Vec4 p2Res,
877 double m1Res, double m2Res) {
878
879 // Have to set flavours and colours.
880 setIdColAcol();
881
882 // Check that masses of outgoing particles not too big.
883 m3 = particleDataPtr->m0(idSave[3]);
884 m4 = particleDataPtr->m0(idSave[4]);
885 mH = sqrt(sH);
886 if (m3 + m4 + MASSMARGIN > mH) return false;
887 s3 = m3 * m3;
888 s4 = m4 * m4;
889
890 // Do kinematics of the production; without or with masses.
891 double e1In = 0.5 * mH;
892 double e2In = e1In;
893 double pzIn = e1In;
894 if (i1Res > 0 || i2Res > 0) {
895 double s1 = m1Res * m1Res;
896 double s2 = m2Res * m2Res;
897 e1In = 0.5 * (sH + s1 - s2) / mH;
898 e2In = 0.5 * (sH + s2 - s1) / mH;
899 pzIn = sqrtpos( e1In*e1In - s1 );
900 }
901
902 // Do kinematics of the decay.
903 double e3 = 0.5 * (sH + s3 - s4) / mH;
904 double e4 = 0.5 * (sH + s4 - s3) / mH;
905 double pAbs = sqrtpos( e3*e3 - s3 );
906 phi = 2. * M_PI * rndmPtr->flat();
907 double pZ = pAbs * cosTheta;
908 pTFin = pAbs * sinTheta;
909 double pX = pTFin * sin(phi);
910 double pY = pTFin * cos(phi);
911 double scale = 0.5 * mH * sinTheta;
912
913 // Fill particle info.
914 int status1 = (i1Res == 0) ? -31 : -34;
915 int status2 = (i2Res == 0) ? -31 : -34;
916 parton[1] = Particle( idSave[1], status1, 0, 0, 3, 4,
917 colSave[1], acolSave[1], 0., 0., pzIn, e1In, m1Res, scale);
918 parton[2] = Particle( idSave[2], status2, 0, 0, 3, 4,
919 colSave[2], acolSave[2], 0., 0., -pzIn, e2In, m2Res, scale);
920 parton[3] = Particle( idSave[3], 33, 1, 2, 0, 0,
921 colSave[3], acolSave[3], pX, pY, pZ, e3, m3, scale);
922 parton[4] = Particle( idSave[4], 33, 1, 2, 0, 0,
923 colSave[4], acolSave[4], -pX, -pY, -pZ, e4, m4, scale);
924
925 // Boost particles from subprocess rest frame to event rest frame.
926 // Normal multiparton interaction: only longitudinal boost.
927 if (i1Res == 0 && i2Res == 0) {
928 double betaZ = (x1Save - x2Save) / (x1Save + x2Save);
929 for (int i = 1; i <= 4; ++i) parton[i].bst(0., 0., betaZ);
930 // Rescattering: generic rotation and boost required.
931 } else {
932 RotBstMatrix M;
933 M.fromCMframe( p1Res, p2Res);
934 for (int i = 1; i <= 4; ++i) parton[i].rotbst(M);
935 }
936
937 // Done.
938 return true;
939
940}
941
942//--------------------------------------------------------------------------
943
944// Calculate modified masses and four-vectors for matrix elements.
945
946bool Sigma2Process::setupForME() {
947
948 // Common initial-state handling.
949 bool allowME = setupForMEin();
950
951 // Correct outgoing c, b, mu and tau to be massive or not.
952 mME[2] = m3;
953 int id3Tmp = abs(id3Mass());
954 if (id3Tmp == 4) mME[2] = mcME;
955 if (id3Tmp == 5) mME[2] = mbME;
956 if (id3Tmp == 13) mME[2] = mmuME;
957 if (id3Tmp == 15) mME[2] = mtauME;
958 mME[3] = m4;
959 int id4Tmp = abs(id4Mass());
960 if (id4Tmp == 4) mME[3] = mcME;
961 if (id4Tmp == 5) mME[3] = mbME;
962 if (id4Tmp == 13) mME[3] = mmuME;
963 if (id4Tmp == 15) mME[3] = mtauME;
964
965 // If kinematically impossible turn to massless case, but set error.
966 if (mME[2] + mME[3] >= mH) {
967 mME[2] = 0.;
968 mME[3] = 0.;
969 allowME = false;
970 }
971
972 // Calculate scattering angle in subsystem rest frame.
973 double sH34 = sqrtpos( pow2(sH - s3 - s4) - 4. * s3 * s4);
974 double cThe = (tH - uH) / sH34;
975 double sThe = sqrtpos(1. - cThe * cThe);
976
977 // Setup massive kinematics with preserved scattering angle.
978 double s3ME = pow2(mME[2]);
979 double s4ME = pow2(mME[3]);
980 double sH34ME = sqrtpos( pow2(sH - s3ME - s4ME) - 4. * s3ME * s4ME);
981 double pAbsME = 0.5 * sH34ME / mH;
982
983 // Normally allowed with unequal (or vanishing) masses.
984 if (id3Tmp == 0 || id3Tmp != id4Tmp) {
985 pME[2] = Vec4( pAbsME * sThe, 0., pAbsME * cThe,
986 0.5 * (sH + s3ME - s4ME) / mH);
987 pME[3] = Vec4( -pAbsME * sThe, 0., -pAbsME * cThe,
988 0.5 * (sH + s4ME - s3ME) / mH);
989
990 // For equal (anti)particles (e.g. W+ W-) use averaged mass.
991 } else {
992 mME[2] = sqrtpos(0.5 * (s3ME + s4ME) - 0.25 * pow2(s3ME - s4ME) / sH);
993 mME[3] = mME[2];
994 pME[2] = Vec4( pAbsME * sThe, 0., pAbsME * cThe, 0.5 * mH);
995 pME[3] = Vec4( -pAbsME * sThe, 0., -pAbsME * cThe, 0.5 * mH);
996 }
997
998 // Done.
999 return allowME;
1000
1001}
1002
1003//==========================================================================
1004
1005// The Sigma3Process class.
1006// Base class for resolved 2 -> 3 cross sections; derived from SigmaProcess.
1007
1008//--------------------------------------------------------------------------
1009
1010// Input and complement kinematics for resolved 2 -> 3 process.
1011
1012void Sigma3Process::store3Kin( double x1in, double x2in, double sHin,
1013 Vec4 p3cmIn, Vec4 p4cmIn, Vec4 p5cmIn, double m3in, double m4in,
1014 double m5in, double runBW3in, double runBW4in, double runBW5in) {
1015
1016 // Default ordering of particles 3 and 4 - not relevant here.
1017 swapTU = false;
1018
1019 // Incoming parton momentum fractions.
1020 x1Save = x1in;
1021 x2Save = x2in;
1022
1023 // Incoming masses and their squares.
1024 if (id3Mass() == 0 && id4Mass() == 0 && id5Mass() == 0) {
1025 m3 = 0.;
1026 m4 = 0.;
1027 m5 = 0.;
1028 } else {
1029 m3 = m3in;
1030 m4 = m4in;
1031 m5 = m5in;
1032 }
1033 mSave[3] = m3;
1034 mSave[4] = m4;
1035 mSave[5] = m5;
1036 s3 = m3 * m3;
1037 s4 = m4 * m4;
1038 s5 = m5 * m5;
1039
1040 // Standard Mandelstam variables and four-momenta in rest frame.
1041 sH = sHin;
1042 mH = sqrt(sH);
1043 sH2 = sH * sH;
1044 p3cm = p3cmIn;
1045 p4cm = p4cmIn;
1046 p5cm = p5cmIn;
1047
1048 // The nominal Breit-Wigner factors with running width.
1049 runBW3 = runBW3in;
1050 runBW4 = runBW4in;
1051 runBW5 = runBW5in;
1052
1053 // Special case: pick scale as if 2 -> 1 process in disguise.
1054 if (isSChannel()) {
1055
1056 // Different options for renormalization scale, but normally sHat.
1057 Q2RenSave = renormMultFac * sH;
1058 if (renormScale1 == 2) Q2RenSave = renormFixScale;
1059
1060 // Different options for factorization scale, but normally sHat.
1061 Q2FacSave = factorMultFac * sH;
1062 if (factorScale1 == 2) Q2RenSave = factorFixScale;
1063
1064 // "Normal" 2 -> 3 processes, i.e. not vector boson fusion.
1065 } else if ( idTchan1() != 23 && idTchan1() != 24 && idTchan2() != 23
1066 && idTchan2() != 24 ) {
1067 double mT3S = s3 + p3cm.pT2();
1068 double mT4S = s4 + p4cm.pT2();
1069 double mT5S = s5 + p5cm.pT2();
1070
1071 // Different options for renormalization scale.
1072 if (renormScale3 == 1) Q2RenSave = min( mT3S, min(mT4S, mT5S) );
1073 else if (renormScale3 == 2) Q2RenSave = sqrt( mT3S * mT4S * mT5S
1074 / max( mT3S, max(mT4S, mT5S) ) );
1075 else if (renormScale3 == 3) Q2RenSave = pow( mT3S * mT4S * mT5S,
1076 1./3. );
1077 else if (renormScale3 == 4) Q2RenSave = (mT3S + mT4S + mT5S) / 3.;
1078 else Q2RenSave = sH;
1079 Q2RenSave *= renormMultFac;
1080 if (renormScale3 == 6) Q2RenSave = renormFixScale;
1081
1082 // Different options for factorization scale.
1083 if (factorScale3 == 1) Q2FacSave = min( mT3S, min(mT4S, mT5S) );
1084 else if (factorScale3 == 2) Q2FacSave = sqrt( mT3S * mT4S * mT5S
1085 / max( mT3S, max(mT4S, mT5S) ) );
1086 else if (factorScale3 == 3) Q2FacSave = pow( mT3S * mT4S * mT5S,
1087 1./3. );
1088 else if (factorScale3 == 4) Q2FacSave = (mT3S + mT4S + mT5S) / 3.;
1089 else Q2FacSave = sH;
1090 Q2FacSave *= factorMultFac;
1091 if (factorScale3 == 6) Q2FacSave = factorFixScale;
1092
1093 // Vector boson fusion 2 -> 3 processes; recoils in positions 4 and 5.
1094 } else {
1095 double sV4 = pow2( particleDataPtr->m0(idTchan1()) );
1096 double sV5 = pow2( particleDataPtr->m0(idTchan2()) );
1097 double mT3S = s3 + p3cm.pT2();
1098 double mTV4S = sV4 + p4cm.pT2();
1099 double mTV5S = sV5 + p5cm.pT2();
1100
1101 // Different options for renormalization scale.
1102 if (renormScale3VV == 1) Q2RenSave = max( sV4, sV5);
1103 else if (renormScale3VV == 2) Q2RenSave = sqrt( mTV4S * mTV5S );
1104 else if (renormScale3VV == 3) Q2RenSave = pow( mT3S * mTV4S * mTV5S,
1105 1./3. );
1106 else if (renormScale3VV == 4) Q2RenSave = (mT3S * mTV4S * mTV5S) / 3.;
1107 else Q2RenSave = sH;
1108 Q2RenSave *= renormMultFac;
1109 if (renormScale3VV == 6) Q2RenSave = renormFixScale;
1110
1111 // Different options for factorization scale.
1112 if (factorScale3VV == 1) Q2FacSave = max( sV4, sV5);
1113 else if (factorScale3VV == 2) Q2FacSave = sqrt( mTV4S * mTV5S );
1114 else if (factorScale3VV == 3) Q2FacSave = pow( mT3S * mTV4S * mTV5S,
1115 1./3. );
1116 else if (factorScale3VV == 4) Q2FacSave = (mT3S * mTV4S * mTV5S) / 3.;
1117 else Q2FacSave = sH;
1118 Q2FacSave *= factorMultFac;
1119 if (factorScale3VV == 6) Q2FacSave = factorFixScale;
1120 }
1121
1122 // Evaluate alpha_strong and alpha_EM.
1123 alpS = couplingsPtr->alphaS(Q2RenSave);
1124 alpEM = couplingsPtr->alphaEM(Q2RenSave);
1125
1126}
1127
1128//--------------------------------------------------------------------------
1129
1130// Calculate modified masses and four-vectors for matrix elements.
1131
1132bool Sigma3Process::setupForME() {
1133
1134 // Common initial-state handling.
1135 bool allowME = setupForMEin();
1136
1137 // Correct outgoing c, b, mu and tau to be massive or not.
1138 mME[2] = m3;
1139 int id3Tmp = abs(id3Mass());
1140 if (id3Tmp == 4) mME[2] = mcME;
1141 if (id3Tmp == 5) mME[2] = mbME;
1142 if (id3Tmp == 13) mME[2] = mmuME;
1143 if (id3Tmp == 15) mME[2] = mtauME;
1144 mME[3] = m4;
1145 int id4Tmp = abs(id4Mass());
1146 if (id4Tmp == 4) mME[3] = mcME;
1147 if (id4Tmp == 5) mME[3] = mbME;
1148 if (id4Tmp == 13) mME[3] = mmuME;
1149 if (id4Tmp == 15) mME[3] = mtauME;
1150 mME[4] = m5;
1151 int id5Tmp = abs(id5Mass());
1152 if (id5Tmp == 4) mME[4] = mcME;
1153 if (id5Tmp == 5) mME[4] = mbME;
1154 if (id5Tmp == 13) mME[4] = mmuME;
1155 if (id5Tmp == 15) mME[4] = mtauME;
1156
1157 // If kinematically impossible turn to massless case, but set error.
1158 if (mME[2] + mME[3] + mME[4] >= mH) {
1159 mME[2] = 0.;
1160 mME[3] = 0.;
1161 mME[4] = 0.;
1162 allowME = false;
1163 }
1164
1165 // Form new average masses if identical particles.
1166 if (id3Tmp != 0 && id4Tmp == id3Tmp && id5Tmp == id3Tmp) {
1167 double mAvg = (mME[2] + mME[3] + mME[4]) / 3.;
1168 mME[2] = mAvg;
1169 mME[3] = mAvg;
1170 mME[4] = mAvg;
1171 } else if (id3Tmp != 0 && id4Tmp == id3Tmp) {
1172 mME[2] = sqrtpos(0.5 * (pow2(mME[2]) + pow2(mME[3]))
1173 - 0.25 * pow2(pow2(mME[2]) - pow2(mME[3])) / sH);
1174 mME[3] = mME[2];
1175 } else if (id3Tmp != 0 && id5Tmp == id3Tmp) {
1176 mME[2] = sqrtpos(0.5 * (pow2(mME[2]) + pow2(mME[4]))
1177 - 0.25 * pow2(pow2(mME[2]) - pow2(mME[4])) / sH);
1178 mME[4] = mME[2];
1179 } else if (id4Tmp != 0 && id5Tmp == id4Tmp) {
1180 mME[3] = sqrtpos(0.5 * (pow2(mME[3]) + pow2(mME[4]))
1181 - 0.25 * pow2(pow2(mME[3]) - pow2(mME[4])) / sH);
1182 mME[4] = mME[2];
1183 }
1184
1185 // Iterate rescaled three-momenta until convergence.
1186 double m2ME3 = pow2(mME[2]);
1187 double m2ME4 = pow2(mME[3]);
1188 double m2ME5 = pow2(mME[4]);
1189 double p2ME3 = p3cm.pAbs2();
1190 double p2ME4 = p4cm.pAbs2();
1191 double p2ME5 = p5cm.pAbs2();
1192 double p2sum = p2ME3 + p2ME4 + p2ME5;
1193 double eME3 = sqrt(m2ME3 + p2ME3);
1194 double eME4 = sqrt(m2ME4 + p2ME4);
1195 double eME5 = sqrt(m2ME5 + p2ME5);
1196 double esum = eME3 + eME4 + eME5;
1197 double p2rat = p2ME3 / eME3 + p2ME4 / eME4 + p2ME5 / eME5;
1198 int iStep = 0;
1199 while ( abs(esum - mH) > COMPRELERR * mH && iStep < NCOMPSTEP ) {
1200 ++iStep;
1201 double compFac = 1. + 2. * (mH - esum) / p2rat;
1202 p2ME3 *= compFac;
1203 p2ME4 *= compFac;
1204 p2ME5 *= compFac;
1205 eME3 = sqrt(m2ME3 + p2ME3);
1206 eME4 = sqrt(m2ME4 + p2ME4);
1207 eME5 = sqrt(m2ME5 + p2ME5);
1208 esum = eME3 + eME4 + eME5;
1209 p2rat = p2ME3 / eME3 + p2ME4 / eME4 + p2ME5 / eME5;
1210 }
1211
1212 // If failed convergence set error flag.
1213 if (abs(esum - mH) > COMPRELERR * mH) allowME = false;
1214
1215 // Set up accepted kinematics.
1216 double totFac = sqrt( (p2ME3 + p2ME4 + p2ME5) / p2sum);
1217 pME[2] = totFac * p3cm;
1218 pME[2].e( eME3);
1219 pME[3] = totFac * p4cm;
1220 pME[3].e( eME4);
1221 pME[4] = totFac * p5cm;
1222 pME[4].e( eME5);
1223
1224 // Done.
1225 return allowME;
1226
1227}
1228
1229//==========================================================================
1230
1231// The SigmaLHAProcess class.
1232// Wrapper for Les Houches Accord external input; derived from SigmaProcess.
1233// Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
1234
1235//--------------------------------------------------------------------------
1236
1237// Evaluate weight for decay angles.
1238
1239double SigmaLHAProcess::weightDecay( Event& process, int iResBeg,
1240 int iResEnd) {
1241
1242 // Do nothing if decays present already at input.
1243 if (iResBeg < process.savedSizeValue()) return 1.;
1244
1245 // Identity of mother of decaying reseonance(s).
1246 int idMother = process[process[iResBeg].mother1()].idAbs();
1247
1248 // For Higgs decay hand over to standard routine.
1249 if (idMother == 25 || idMother == 35 || idMother == 36)
1250 return weightHiggsDecay( process, iResBeg, iResEnd);
1251
1252 // For top decay hand over to standard routine.
1253 if (idMother == 6)
1254 return weightTopDecay( process, iResBeg, iResEnd);
1255
1256 // Else done.
1257 return 1.;
1258
1259}
1260
1261//--------------------------------------------------------------------------
1262
1263// Set scale, alpha_strong and alpha_EM when not set.
1264
1265void SigmaLHAProcess::setScale() {
1266
1267 // If scale has not been set, then to set.
1268 double scaleLHA = lhaUpPtr->scale();
1269 if (scaleLHA < 0.) {
1270
1271 // Final-state partons and their invariant mass.
1272 vector<int> iFin;
1273 Vec4 pFinSum;
1274 for (int i = 3; i < lhaUpPtr->sizePart(); ++i)
1275 if (lhaUpPtr->mother1(i) == 1) {
1276 iFin.push_back(i);
1277 pFinSum += Vec4( lhaUpPtr->px(i), lhaUpPtr->py(i),
1278 lhaUpPtr->pz(i), lhaUpPtr->e(i) );
1279 }
1280 int nFin = iFin.size();
1281 sH = pFinSum * pFinSum;
1282 mH = sqrt(sH);
1283 sH2 = sH * sH;
1284
1285 // If 1 final-state particle then use Sigma1Process logic.
1286 if (nFin == 1) {
1287 Q2RenSave = renormMultFac * sH;
1288 if (renormScale1 == 2) Q2RenSave = renormFixScale;
1289 Q2FacSave = factorMultFac * sH;
1290 if (factorScale1 == 2) Q2FacSave = factorFixScale;
1291
1292 // If 2 final-state particles then use Sigma2Process logic.
1293 } else if (nFin == 2) {
1294 double s3 = pow2(lhaUpPtr->m(iFin[0]));
1295 double s4 = pow2(lhaUpPtr->m(iFin[1]));
1296 double pT2 = pow2(lhaUpPtr->px(iFin[0])) + pow2(lhaUpPtr->py(iFin[0]));
1297 if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
1298 else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
1299 else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
1300 else Q2RenSave = sH;
1301 Q2RenSave *= renormMultFac;
1302 if (renormScale2 == 5) Q2RenSave = renormFixScale;
1303 if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
1304 else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
1305 else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
1306 else Q2FacSave = sH;
1307 Q2FacSave *= factorMultFac;
1308 if (factorScale2 == 5) Q2FacSave = factorFixScale;
1309
1310 // If 3 or more final-state particles then use Sigma3Process logic.
1311 } else {
1312 double mTSlow = sH;
1313 double mTSmed = sH;
1314 double mTSprod = 1.;
1315 double mTSsum = 0.;
1316 for (int i = 0; i < nFin; ++i) {
1317 double mTSnow = pow2(lhaUpPtr->m(iFin[i]))
1318 + pow2(lhaUpPtr->px(iFin[i])) + pow2(lhaUpPtr->py(iFin[i]));
1319 if (mTSnow < mTSlow) {mTSmed = mTSlow; mTSlow = mTSnow;}
1320 else if (mTSnow < mTSmed) mTSmed = mTSnow;
1321 mTSprod *= mTSnow;
1322 mTSsum += mTSnow;
1323 }
1324 if (renormScale3 == 1) Q2RenSave = mTSlow;
1325 else if (renormScale3 == 2) Q2RenSave = sqrt(mTSlow * mTSmed);
1326 else if (renormScale3 == 3) Q2RenSave = pow(mTSprod, 1. / nFin);
1327 else if (renormScale3 == 4) Q2RenSave = mTSsum / nFin;
1328 else Q2RenSave = sH;
1329 Q2RenSave *= renormMultFac;
1330 if (renormScale3 == 6) Q2RenSave = renormFixScale;
1331 if (factorScale3 == 1) Q2FacSave = mTSlow;
1332 else if (factorScale3 == 2) Q2FacSave = sqrt(mTSlow * mTSmed);
1333 else if (factorScale3 == 3) Q2FacSave = pow(mTSprod, 1. / nFin);
1334 else if (factorScale3 == 4) Q2FacSave = mTSsum / nFin;
1335 else Q2FacSave = sH;
1336 Q2FacSave *= factorMultFac;
1337 if (factorScale3 == 6) Q2FacSave = factorFixScale;
1338 }
1339 }
1340
1341 // If alpha_strong and alpha_EM have not been set, then set them.
1342 if (lhaUpPtr->alphaQCD() < 0.001) {
1343 double Q2RenNow = (scaleLHA < 0.) ? Q2RenSave : pow2(scaleLHA);
1344 alpS = couplingsPtr->alphaS(Q2RenNow);
1345 }
1346 if (lhaUpPtr->alphaQED() < 0.001) {
1347 double Q2RenNow = (scaleLHA < 0.) ? Q2RenSave : pow2(scaleLHA);
1348 alpEM = couplingsPtr->alphaEM(Q2RenNow);
1349 }
1350
1351}
1352
1353//--------------------------------------------------------------------------
1354
1355// Obtain number of final-state partons from LHA object.
1356
1357int SigmaLHAProcess::nFinal() const {
1358
1359 // At initialization size unknown, so return 0.
1360 if (lhaUpPtr->sizePart() <= 0) return 0;
1361
1362 // Sum up all particles that has first mother = 1.
1363 int nFin = 0;
1364 for (int i = 3; i < lhaUpPtr->sizePart(); ++i)
1365 if (lhaUpPtr->mother1(i) == 1) ++nFin;
1366 return nFin;
1367
1368}
1369
1370//==========================================================================
1371
1372} // end namespace Pythia8