source: HiSusy/trunk/Pythia8/pythia8170/src/SigmaProcess.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES
File size: 45.4 KB
Line 
1// SigmaProcess.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2012 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
Note: See TracBrowser for help on using the repository browser.