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