using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SigmaProcess.cxx
1 // SigmaProcess.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2008 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
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 //*********
30
31 // Perform simple initialization and store pointers.
32
33 void SigmaProcess::init(Info* infoPtrIn, BeamParticle* beamAPtrIn, 
34   BeamParticle* beamBPtrIn, AlphaStrong* alphaSPtrIn, 
35   AlphaEM* alphaEMPtrIn,SigmaTotal* sigmaTotPtrIn,
36   SusyLesHouches* slhaPtrIn) {
37
38   // Store pointers.
39   infoPtr     = infoPtrIn;
40   beamAPtr    = beamAPtrIn;
41   beamBPtr    = beamBPtrIn;
42   alphaSPtr   = alphaSPtrIn;
43   alphaEMPtr  = alphaEMPtrIn;
44   sigmaTotPtr = sigmaTotPtrIn;
45   slhaPtr     = slhaPtrIn;
46
47   // Read out some properties of beams to allow shorthand.
48   idA         = beamAPtr->id();
49   idB         = beamBPtr->id();
50   mA          = beamAPtr->m();
51   mB          = beamBPtr->m();
52   isLeptonA   = beamAPtr->isLepton();
53   isLeptonB   = beamBPtr->isLepton();
54   hasLeptonBeams = isLeptonA || isLeptonB;
55
56   // K factor, multiplying resolved processes. (But not here for MI.)
57   Kfactor         = Settings::parm("SigmaProcess:Kfactor");
58
59   // Maximum incoming quark flavour.
60   nQuarkIn        = Settings::mode("PDFinProcess:nQuarkIn");
61
62   // Renormalization scale choice.
63   renormScale1    = Settings::mode("SigmaProcess:renormScale1"); 
64   renormScale2    = Settings::mode("SigmaProcess:renormScale2"); 
65   renormScale3    = Settings::mode("SigmaProcess:renormScale3"); 
66   renormScale3VV  = Settings::mode("SigmaProcess:renormScale3VV"); 
67   renormMultFac   = Settings::parm("SigmaProcess:renormMultFac"); 
68   renormFixScale  = Settings::parm("SigmaProcess:renormFixScale"); 
69
70   // Factorization scale choice.
71   factorScale1    = Settings::mode("SigmaProcess:factorScale1"); 
72   factorScale2    = Settings::mode("SigmaProcess:factorScale2"); 
73   factorScale3    = Settings::mode("SigmaProcess:factorScale3"); 
74   factorScale3VV  = Settings::mode("SigmaProcess:factorScale3VV"); 
75   factorMultFac   = Settings::parm("SigmaProcess:factorMultFac"); 
76   factorFixScale  = Settings::parm("SigmaProcess:factorFixScale"); 
77
78   // CP violation parameters for the BSM Higgs sector.
79   higgsH1parity   = Settings::mode("HiggsH1:parity");
80   higgsH1eta      = Settings::parm("HiggsH1:etaParity");
81   higgsH2parity   = Settings::mode("HiggsH2:parity");
82   higgsH2eta      = Settings::parm("HiggsH2:etaParity");
83   higgsA3parity   = Settings::mode("HiggsA3:parity");
84   higgsA3eta      = Settings::parm("HiggsA3:etaParity");
85
86   // If BSM not switched on then H1 should have SM properties.
87   if (!Settings::flag("Higgs:useBSM")){
88     higgsH1parity = 1;
89     higgsH1eta    = 0.;
90   }
91
92 }
93
94 //*********
95
96 // Set up allowed flux of incoming partons.
97 // addBeam: set up PDF's that need to be evaluated for the two beams.
98 // addPair: set up pairs of incoming partons from the two beams.
99
100 bool SigmaProcess::initFlux() {
101
102   // Read in process-specific channel information.
103   string fluxType = inFlux();
104
105   // Case with g g incoming state.
106   if (fluxType == "gg") {
107     addBeamA(21);
108     addBeamB(21);
109     addPair(21, 21);
110   }
111
112   // Case with q g incoming state.
113   else if (fluxType == "qg") {
114     for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
115       int idNow = (i == 0) ? 21 : i;
116       addBeamA(idNow);
117       addBeamB(idNow);
118     }
119     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
120     if (idNow != 0) {
121       addPair(idNow, 21);
122       addPair(21, idNow);
123     }
124   }
125
126   // Case with q q', q qbar' or qbar qbar' incoming state.
127   else if (fluxType == "qq") {
128     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
129     if (idNow != 0) {
130       addBeamA(idNow);
131       addBeamB(idNow);
132     }
133     for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now) 
134     if (id1Now != 0) 
135     for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now) 
136     if (id2Now != 0) 
137       addPair(id1Now, id2Now);
138   }
139
140   // Case with q qbar incoming state.
141   else if (fluxType == "qqbarSame") {
142     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
143     if (idNow != 0) {
144       addBeamA(idNow);
145       addBeamB(idNow);
146     }
147     for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
148     if (idNow != 0) 
149       addPair(idNow, -idNow);
150   }
151
152   // Case with f f', f fbar', fbar fbar' incoming state.
153   else if (fluxType == "ff") {
154     // If beams are leptons then they are also the colliding partons.
155     if ( isLeptonA && isLeptonB ) {
156       addBeamA(idA);
157       addBeamB(idB);
158       addPair(idA, idB);
159     // First beam is lepton and second is hadron.
160     } else if ( isLeptonA ) {
161       addBeamA(idA);
162       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
163       if (idNow != 0) {
164         addBeamB(idNow);
165         addPair(idA, idNow);
166       }
167     // First beam is hadron and second is lepton.
168     } else if ( isLeptonB ) {
169       addBeamB(idB);
170       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
171       if (idNow != 0) {
172         addBeamA(idNow);
173         addPair(idNow, idB);
174       }
175     // Hadron beams gives quarks.
176     } else {
177       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
178       if (idNow != 0) {
179         addBeamA(idNow);
180         addBeamB(idNow);
181       }
182       for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now) 
183       if (id1Now != 0) 
184       for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now) 
185       if (id2Now != 0) 
186         addPair(id1Now, id2Now);
187     }
188   }
189
190   // Case with f fbar incoming state.
191   else if (fluxType == "ffbarSame") {
192     // If beams are antiparticle pair and leptons then also colliding partons.
193     if ( idA + idB == 0 && isLeptonA ) {
194       addBeamA(idA);
195       addBeamB(idB);
196       addPair(idA, idB);
197     // Else assume both to be hadrons, for better or worse.
198     } else {
199       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
200       if (idNow != 0) {
201         addBeamA(idNow);
202         addBeamB(idNow);
203       }
204       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
205       if (idNow != 0) 
206         addPair(idNow, -idNow);
207     }
208   }
209
210   // Case with f fbar' charged(+-1) incoming state.
211   else if (fluxType == "ffbarChg") {
212     // If beams are leptons then also colliding partons.
213     if ( isLeptonA && isLeptonB && abs( ParticleDataTable::chargeType(idA) 
214              + ParticleDataTable::chargeType(idB) ) == 3 ) {
215       addBeamA(idA);
216       addBeamB(idB);
217       addPair(idA, idB);
218     // Hadron beams gives quarks.
219     } else {
220       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
221       if (idNow != 0) {
222         addBeamA(idNow);
223         addBeamB(idNow);
224       }
225       for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now) 
226       if (id1Now != 0) 
227       for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now) 
228       if (id2Now != 0 && id1Now * id2Now < 0 
229         && (abs(id1Now) + abs(id2Now))%2 == 1) addPair(id1Now, id2Now);
230     }
231   }
232
233   // Case with f fbar' generic incoming state.
234   else if (fluxType == "ffbar") {
235     // If beams are leptons then also colliding partons.
236     if (isLeptonA && isLeptonB && idA * idB < 0) {
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         addPair(id1Now, id2Now);
252     }
253   }
254
255   // Case with f gamma incoming state.
256   else if (fluxType == "fgm") {
257     // Fermion from incoming side A.
258     if ( isLeptonA ) {
259       addBeamA(idA);
260       addPair(idA, 22);
261     } else {  
262       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
263       if (idNow != 0) {
264         addBeamA(idNow);
265         addPair(idNow, 22);
266       }
267     }
268     // Fermion from incoming side B.
269     if ( isLeptonB ) {
270       addBeamB( idB);
271       addPair(22, idB);
272     } else {  
273       for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) 
274       if (idNow != 0) {
275         addBeamB(idNow);
276         addPair(22, idNow);
277       }
278     }
279     // Photons in the beams.
280     addBeamA(22);
281     addBeamB(22);
282   }
283
284   // Case with gamma gamma incoming state.
285   else if (fluxType == "ggm") {
286     addBeamA(21);
287     addBeamA(22);
288     addBeamB(21);
289     addBeamB(22);
290     addPair(21, 22);
291     addPair(22, 21);
292   }
293
294   // Case with gamma gamma incoming state.
295   else if (fluxType == "gmgm") {
296     addBeamA(22);
297     addBeamB(22);
298     addPair(22, 22);
299   }
300
301   // Unrecognized fluxType is bad sign. Else done.
302   else {
303     infoPtr->errorMsg("Error in SigmaProcess::initFlux: "
304     "unrecognized inFlux type", fluxType);
305     return false;
306   }
307   return true;
308
309 }
310
311 //*********
312
313 // Convolute matrix-element expression(s) with parton flux and K factor.
314
315 double SigmaProcess::sigmaPDF() {
316
317   // Evaluate and store the required parton densities.
318   for (int j = 0; j < sizeBeamA(); ++j) 
319     inBeamA[j].pdf = beamAPtr->xfHard( inBeamA[j].id, x1Save, Q2FacSave); 
320   for (int j = 0; j < sizeBeamB(); ++j) 
321     inBeamB[j].pdf = beamBPtr->xfHard( inBeamB[j].id, x2Save, Q2FacSave); 
322
323   // Loop over allowed incoming channels.
324   sigmaSumSave = 0.;
325   for (int i = 0; i < sizePair(); ++i) {
326     
327     // Evaluate hard-scattering cross section. Include K factor.
328     inPair[i].pdfSigma = Kfactor 
329                        * sigmaHatWrap(inPair[i].idA, inPair[i].idB);
330     
331     // Multiply by respective parton densities.
332     for (int j = 0; j < sizeBeamA(); ++j) 
333     if (inPair[i].idA == inBeamA[j].id) {
334       inPair[i].pdfA      = inBeamA[j].pdf;
335       inPair[i].pdfSigma *= inBeamA[j].pdf;
336       break;
337     }
338     for (int j = 0; j < sizeBeamB(); ++j) 
339     if (inPair[i].idB == inBeamB[j].id) {
340       inPair[i].pdfB      = inBeamB[j].pdf;
341       inPair[i].pdfSigma *= inBeamB[j].pdf;
342       break;
343     }
344
345     // Sum for all channels.
346     sigmaSumSave += inPair[i].pdfSigma;
347   }
348  
349   // Done.
350   return sigmaSumSave;
351
352 }
353
354 //*********
355
356 // Select incoming parton channel and extract parton densities (resolved).
357
358 void SigmaProcess::pickInState(int id1in, int id2in) {
359
360   // Multiple interactions: partons already selected.
361   if (id1in != 0 && id2in != 0) {
362     id1 = id1in;
363     id2 = id2in;
364   }
365
366   // Pick channel. Extract channel flavours and pdf's.
367   double sigmaRand =  sigmaSumSave * Rndm::flat();
368   for (int i = 0; i < sizePair(); ++i) {
369     sigmaRand -= inPair[i].pdfSigma;
370     if (sigmaRand <= 0.) {
371       id1      = inPair[i].idA;
372       id2      = inPair[i].idB;
373       pdf1Save = inPair[i].pdfA; 
374       pdf2Save = inPair[i].pdfB; 
375       break;
376     }
377   }
378
379 }
380
381 //*********
382
383 // Evaluate weight for W decay distribution in t -> W b -> f fbar b.
384
385 double SigmaProcess::weightTopDecay( Event& process, int iResBeg,
386   int iResEnd) {
387
388   // If not pair W d/s/b and mother t then return unit weight.
389   if (iResEnd - iResBeg != 1) return 1.;
390   int iW1  = iResBeg;
391   int iB2  = iResBeg + 1;
392   int idW1 = process[iW1].idAbs();
393   int idB2 = process[iB2].idAbs();
394   if (idW1 != 24) {
395     swap(iW1, iB2); 
396     swap(idW1, idB2);
397   } 
398   if (idW1 != 24 || (idB2 != 1 && idB2 != 3 && idB2 != 5)) return 1.;
399   int iT   = process[iW1].mother1(); 
400   if (iT <= 0 || process[iT].idAbs() != 6) return 1.;
401
402   // Find sign-matched order of W decay products. 
403   int iF    = process[iW1].daughter1(); 
404   int iFbar = process[iW1].daughter2();
405   if (iFbar - iF != 1) return 1.; 
406   if (process[iT].id() * process[iF].id() < 0) swap(iF, iFbar);
407
408   // Weight and maximum weight.
409   double wt    = (process[iT].p() * process[iFbar].p()) 
410                * (process[iF].p() * process[iB2].p());
411   double wtMax = ( pow4(process[iT].m()) - pow4(process[iW1].m()) ) / 8.;  
412
413   // Done.
414   return wt / wtMax;
415
416 }
417
418
419 //*********
420
421 // Evaluate weight for Z0/W+- decay distributions in H -> Z0/W+ Z0/W- -> 4f.
422
423 double SigmaProcess::weightHiggsDecay( Event& process, int iResBeg, 
424   int iResEnd) {
425
426   // If not pair Z0 Z0 or W+ W- then return unit weight.
427   if (iResEnd - iResBeg != 1) return 1.;
428   int iZW1  = iResBeg;
429   int iZW2  = iResBeg + 1;
430   int idZW1 = process[iZW1].id();
431   int idZW2 = process[iZW2].id();
432   if (idZW1 < 0) {
433     swap(iZW1, iZW2); 
434     swap(idZW1, idZW2);
435   } 
436   if ( (idZW1 != 23 || idZW2 != 23) && (idZW1 != 24 || idZW2 != -24) )
437     return 1.;
438
439   // If mother is not Higgs then return unit weight.
440   int iH  = process[iZW1].mother1(); 
441   if (iH <= 0) return 1.;
442   int idH = process[iH].id();
443   if (idH != 25 && idH != 35 && idH !=36) return 1.;
444
445   // Parameters depend on Higgs type: H0(H_1), H^0(H_2) or A^0(H_3).
446   int    higgsParity = higgsH1parity; 
447   double higgsEta    = higgsH1eta;
448   if (idH == 35) {
449     higgsParity      = higgsH2parity;
450     higgsEta         = higgsH2eta;
451   } else if (idH == 36) {
452     higgsParity      = higgsA3parity;
453     higgsEta         = higgsA3eta;
454   }
455
456   // Option with isotropic decays.
457   if (higgsParity == 0) return 1.;
458
459   // Maximum and initial weight. 
460   double wtMax = pow4(process[iH].m());
461   double wt    = wtMax; 
462
463   // Find sign-matched order of Z0/W+- decay products. 
464   int i3 = process[iZW1].daughter1();
465   int i4 = process[iZW1].daughter2();
466   if (process[i3].id() < 0) swap( i3, i4); 
467   int i5 = process[iZW2].daughter1();
468   int i6 = process[iZW2].daughter2();
469   if (process[i5].id() < 0) swap( i5, i6); 
470
471   // Evaluate four-vector products and find masses..
472   double p35  = 2. * process[i3].p() * process[i5].p(); 
473   double p36  = 2. * process[i3].p() * process[i6].p(); 
474   double p45  = 2. * process[i4].p() * process[i5].p(); 
475   double p46  = 2. * process[i4].p() * process[i6].p(); 
476   double p34  = 2. * process[i3].p() * process[i4].p(); 
477   double p56  = 2. * process[i5].p() * process[i6].p(); 
478   double mZW1 = process[iZW1].m();
479   double mZW2 = process[iZW2].m();
480
481   // For mixed CP states need epsilon product and gauge boson masses.
482   double epsilonProd = 0.;
483   if (higgsParity == 3) {
484     double p[4][4];
485     for (int i = 0; i < 4; ++i) {
486       int         ii = i3;
487       if (i == 1) ii = i4;
488       if (i == 2) ii = i5;
489       if (i == 3) ii = i6;
490       p[i][0] = process[ii].e();
491       p[i][1] = process[ii].px();
492       p[i][2] = process[ii].py();
493       p[i][3] = process[ii].pz();
494     }     
495     epsilonProd 
496       = p[0][0]*p[1][1]*p[2][2]*p[3][3] - p[0][0]*p[1][1]*p[2][3]*p[3][2] 
497       - p[0][0]*p[1][2]*p[2][1]*p[3][3] + p[0][0]*p[1][2]*p[2][3]*p[3][1]
498       + p[0][0]*p[1][3]*p[2][1]*p[3][2] - p[0][0]*p[1][3]*p[2][2]*p[3][1]
499       - p[0][1]*p[1][0]*p[2][2]*p[3][3] + p[0][1]*p[1][0]*p[2][3]*p[3][2]
500       + p[0][1]*p[1][2]*p[2][0]*p[3][3] - p[0][1]*p[1][2]*p[2][3]*p[3][0]
501       - p[0][1]*p[1][3]*p[2][0]*p[3][2] + p[0][1]*p[1][3]*p[2][2]*p[3][0]
502       + p[0][2]*p[1][0]*p[2][1]*p[3][3] - p[0][2]*p[1][0]*p[2][3]*p[3][1]
503       - p[0][2]*p[1][1]*p[2][0]*p[3][3] + p[0][2]*p[1][1]*p[2][3]*p[3][0] 
504       + p[0][2]*p[1][3]*p[2][0]*p[3][1] - p[0][2]*p[1][3]*p[2][1]*p[3][0]
505       - p[0][3]*p[1][0]*p[2][1]*p[3][2] + p[0][3]*p[1][0]*p[2][2]*p[3][1] 
506       + p[0][3]*p[1][1]*p[2][0]*p[3][2] - p[0][3]*p[1][1]*p[2][2]*p[3][0] 
507       - p[0][3]*p[1][2]*p[2][0]*p[3][1] + p[0][3]*p[1][2]*p[2][1]*p[3][0];
508   }
509
510   // Z0 Z0 decay: vector and axial couplings of two fermion pairs.
511   if (idZW1 == 23) {
512     double vf1 = CoupEW::vf(process[i3].idAbs());
513     double af1 = CoupEW::af(process[i3].idAbs());
514     double vf2 = CoupEW::vf(process[i5].idAbs());
515     double af2 = CoupEW::af(process[i5].idAbs());
516     double va12asym = 4. * vf1 * af1 * vf2 * af2 
517       / ( (vf1*vf1 + af1*af1) * (vf2*vf2 + af2*af2) );
518     double etaMod = higgsEta / pow2( ParticleDataTable::m0(23) );
519     
520     // Normal CP-even decay.
521     if (higgsParity == 1) wt = 8. * (1. + va12asym) * p35 * p46 
522       + 8. * (1. - va12asym) * p36 * p45;
523
524     // CP-odd decay (normal for A0(H_3)).
525     else if (higgsParity == 2) wt = ( pow2(p35 + p46) 
526       + pow2(p36 + p45) - 2. * p34 * p56
527       - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56) 
528       + va12asym * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) )
529       / (1. +  va12asym);
530
531     // Mixed CP states. 
532     else wt = 32. * ( 0.25 * ( (1. + va12asym) * p35 * p46 
533       + (1. - va12asym) * p36 * p45 ) - 0.5 * etaMod * epsilonProd
534       * ( (1. + va12asym) * (p35 + p46) - (1. - va12asym) * (p36 + p45) )
535       + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56) 
536       - 2. * pow2(p35 * p46 - p36 * p45) 
537       + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45)) 
538       + va12asym * p34 * p56 * (p35 + p36 - p45 - p46) 
539       * (p35 + p45 - p36 - p46) ) ) / ( 1. + 2. * etaMod * mZW1 * mZW2 
540       + 2. * pow2(etaMod * mZW1 * mZW2) * (1. + va12asym) );
541
542   // W+ W- decay.
543   } else if (idZW1 == 24) {
544     double etaMod = higgsEta / pow2( ParticleDataTable::m0(24) );
545     
546     // Normal CP-even decay.
547     if (higgsParity == 1) wt = 16. * p35 * p46; 
548
549     // CP-odd decay (normal for A0(H_3)).
550     else if (higgsParity == 2) wt = 0.5 * ( pow2(p35 + p46) 
551       + pow2(p36 + p45) - 2. * p34 * p56  
552       - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56) 
553       + (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) );
554
555     // Mixed CP states. 
556     else wt = 32. * ( 0.25 * 2. * p35 * p46 
557       - 0.5 * etaMod * epsilonProd * 2. * (p35 + p46)
558       + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56) 
559       - 2. * pow2(p35 * p46 - p36 * p45) 
560       + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45)) 
561       + p34 * p56 * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ) ) 
562       / ( 1. * 2. * etaMod * mZW1 * mZW2 + 2. * pow2(etaMod * mZW1 * mZW2) );
563   }
564
565   // Done.
566   return wt / wtMax;
567
568 }
569
570 //**************************************************************************
571
572 // The Sigma1Process class.
573 // Base class for resolved 2 -> 1 cross sections; derived from SigmaProcess.
574
575 //*********
576
577 // Input and complement kinematics for resolved 2 -> 1 process. 
578
579 void Sigma1Process::store1Kin( double x1in, double x2in, double sHin) {
580
581   // Default value only sensible for these processes.
582   swapTU = false;
583
584   // Incoming parton momentum fractions and sHat.
585   x1Save = x1in;
586   x2Save = x2in;
587   sH     = sHin;
588   mH     = sqrt(sH);
589   sH2    = sH * sH;
590
591   // Different options for renormalization scale, but normally sHat.
592   Q2RenSave = renormMultFac * sH;
593   if (renormScale1 == 2) Q2RenSave = renormFixScale; 
594
595   // Different options for factorization scale, but normally sHat.
596   Q2FacSave = factorMultFac * sH;
597   if (factorScale1 == 2) Q2RenSave = factorFixScale; 
598
599   // Evaluate alpha_strong and alpha_EM.
600   alpS   = alphaSPtr->alphaS(Q2RenSave);  
601   alpEM  = alphaEMPtr->alphaEM(Q2RenSave);  
602
603 }
604
605 //**************************************************************************
606
607 // The Sigma2Process class.
608 // Base class for resolved 2 -> 2 cross sections; derived from SigmaProcess.
609
610 //*********
611
612 // Input and complement kinematics for resolved 2 -> 2 process. 
613
614 void Sigma2Process::store2Kin( double x1in, double x2in, double sHin, 
615   double tHin, double m3in, double m4in, double runBW3in, double runBW4in) {
616
617   // Default ordering of particles 3 and 4.
618   swapTU   = false;
619
620   // Incoming parton momentum fractions.
621   x1Save   = x1in;
622   x2Save   = x2in;
623
624   // Incoming masses and their squares.
625   bool masslessKin = (id3Mass() == 0) && (id4Mass() == 0);
626   if (masslessKin) {
627     m3     = 0.;
628     m4     = 0.;
629   } else {
630     m3     = m3in;
631     m4     = m4in;
632   }
633   mSave[3] = m3;
634   mSave[4] = m4;
635   s3       = m3 * m3;
636   s4       = m4 * m4;
637
638   // Standard Mandelstam variables and their squares.
639   sH       = sHin;
640   tH       = tHin;
641   uH       = (masslessKin) ? -(sH + tH) : s3 + s4 - (sH + tH); 
642   mH       = sqrt(sH);
643   sH2      = sH * sH;
644   tH2      = tH * tH;
645   uH2      = uH * uH;
646
647   // The nominal Breit-Wigner factors with running width.
648   runBW3   = runBW3in;
649   runBW4   = runBW4in; 
650
651   // Calculate squared transverse momentum.
652   pT2 = (masslessKin) ?  tH * uH / sH : (tH * uH - s3 * s4) / sH;
653
654   // Special case: pick scale as if 2 -> 1 process in disguise.
655   if (isSChannel()) {
656
657     // Different options for renormalization scale, but normally sHat.
658     Q2RenSave = renormMultFac * sH;
659     if (renormScale1 == 2) Q2RenSave = renormFixScale; 
660
661     // Different options for factorization scale, but normally sHat.
662     Q2FacSave = factorMultFac * sH;
663     if (factorScale1 == 2) Q2RenSave = factorFixScale; 
664
665   // Normal case with "true" 2 -> 2.  
666   } else { 
667
668     // Different options for renormalization scale.
669     if (masslessKin)            Q2RenSave = (renormScale2 < 4) ? pT2 : sH;
670     else if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4);
671     else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4));
672     else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4);
673     else                        Q2RenSave = sH;
674     Q2RenSave *= renormMultFac;
675     if      (renormScale2 == 5) Q2RenSave = renormFixScale; 
676
677     // Different options for factorization scale.
678     if (masslessKin)            Q2FacSave = (factorScale2 < 4) ? pT2 : sH;
679     else if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4);
680     else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4));
681     else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4);
682     else                        Q2FacSave = sH;
683     Q2FacSave *= factorMultFac;
684     if      (factorScale2 == 5) Q2FacSave = factorFixScale; 
685   }
686
687   // Evaluate alpha_strong and alpha_EM.
688   alpS = alphaSPtr->alphaS(Q2RenSave);  
689   alpEM = alphaEMPtr->alphaEM(Q2RenSave);  
690
691 }
692
693 //*********
694
695 // As above, special kinematics for multiple interactions. 
696
697 void Sigma2Process::store2KinMI( double x1in, double x2in,
698   double sHin, double tHin, double uHin, double alpSin, double alpEMin,
699   bool needMasses, double m3in, double m4in) {
700
701   // Default ordering of particles 3 and 4.
702   swapTU    = false;
703  
704   // Incoming x values.
705   x1Save    = x1in;
706   x2Save    = x2in;
707
708   // Standard Mandelstam variables and their squares.
709   sH        = sHin;
710   tH        = tHin;
711   uH        = uHin; 
712   mH        = sqrt(sH);
713   sH2       = sH * sH;
714   tH2       = tH * tH;
715   uH2       = uH * uH;
716
717   // Strong and electroweak couplings.
718   alpS      = alpSin;
719   alpEM     = alpEMin;
720
721   // Assume vanishing masses. (Will be modified in final kinematics.) 
722   m3        = 0.;
723   s3        = 0.;
724   m4        = 0.;
725   s4        = 0.;
726   sHBeta    = sH; 
727
728   // Scattering angle.
729   cosTheta  = (tH - uH) / sH;
730   sinTheta  = 2. * sqrtpos( tH * uH ) / sH;
731
732   // In some cases must use masses and redefine meaning of tHat and uHat.
733   if (needMasses) { 
734     m3      = m3in;
735     s3      = m3 * m3;
736     m4      = m4in;
737     s4      = m4 * m4;
738     sHMass  = sH - s3 - s4;
739     sHBeta  = sqrtpos(sHMass*sHMass - 4. * s3 * s4);   
740     tH      = -0.5 * (sHMass - sHBeta * cosTheta); 
741     uH      = -0.5 * (sHMass + sHBeta * cosTheta); 
742     tH2     = tH * tH;
743     uH2     = uH * uH;
744   }
745
746   // pT2 with masses (at this stage) included.
747   pT2Mass   = 0.25 * sHBeta * pow2(sinTheta);
748
749 }
750
751 //*********
752
753 // Perform kinematics for a Multiple Interaction.
754
755 bool Sigma2Process::final2KinMI() {
756
757   // Have to set flavours and colours.
758   setIdColAcol();
759
760   // Check that masses of outgoing particles not too big.
761   m3           = ParticleDataTable::m0(idSave[3]);
762   m4           = ParticleDataTable::m0(idSave[4]);
763   mH           = sqrt(sH);
764   if (m3 + m4 + MASSMARGIN > mH) return false;
765   s3           = m3 * m3;
766   s4           = m4 * m4;
767
768   // Do kinematics of the decay.
769   double eIn   = 0.5 * mH;
770   double e3    = 0.5 * (sH + s3 - s4) / mH;
771   double e4    = 0.5 * (sH + s4 - s3) / mH;
772   double pAbs  = sqrtpos( e3*e3 - s3 );
773   phi          = 2. * M_PI * Rndm::flat();
774   double pZ    = pAbs * cosTheta;
775   double pX    = pAbs * sinTheta * sin(phi);
776   double pY    = pAbs * sinTheta * cos(phi);
777   double scale = eIn * sinTheta;
778
779   // Fill particle info.
780   parton[1] = Particle( idSave[1], -31, 0, 0, 3, 4, colSave[1], acolSave[1],
781     0., 0., eIn, eIn, 0., scale);
782   parton[2] = Particle( idSave[2], -31, 0, 0, 3, 4, colSave[2], acolSave[2],
783     0., 0., -eIn, eIn, 0., scale);
784   parton[3] = Particle( idSave[3],  33, 1, 2, 0, 0, colSave[3], acolSave[3],
785     pX, pY, pZ, e3, m3, scale);
786   parton[4] = Particle( idSave[4],  33, 1, 2, 0, 0, colSave[4], acolSave[4],
787     -pX, -pY, -pZ, e4, m4, scale);
788
789   // Boost particles from subprocess rest frame to event rest frame.
790   double betaZ = (x1Save - x2Save) / (x1Save + x2Save);
791   for (int i = 1; i <= 4; ++i) parton[i].bst(0., 0., betaZ);
792
793   // Done.
794   return true;
795
796 }  
797
798 //**************************************************************************
799
800 // The Sigma3Process class.
801 // Base class for resolved 2 -> 3 cross sections; derived from SigmaProcess.
802
803 //*********
804
805 // Input and complement kinematics for resolved 2 -> 3 process. 
806
807 void Sigma3Process::store3Kin( double x1in, double x2in, double sHin, 
808   Vec4 p3cmIn, Vec4 p4cmIn, Vec4 p5cmIn, double m3in, double m4in, 
809   double m5in, double runBW3in, double runBW4in, double runBW5in) {
810
811   // Default ordering of particles 3 and 4 - not relevant here.
812   swapTU   = false;
813
814   // Incoming parton momentum fractions.
815   x1Save   = x1in;
816   x2Save   = x2in;
817
818   // Incoming masses and their squares.
819   if (id3Mass() == 0 && id4Mass() == 0 && id5Mass() == 0) {
820     m3     = 0.;
821     m4     = 0.;
822     m5     = 0.;
823   } else {
824     m3     = m3in;
825     m4     = m4in;
826     m5     = m5in;
827   }
828   mSave[3] = m3;
829   mSave[4] = m4;
830   mSave[5] = m5;
831   s3       = m3 * m3;
832   s4       = m4 * m4;
833   s5       = m5 * m5;
834
835   // Standard Mandelstam variables and four-momenta in rest frame.
836   sH       = sHin;
837   mH       = sqrt(sH);
838   p3cm     = p3cmIn;
839   p4cm     = p4cmIn;
840   p5cm     = p5cmIn;
841
842   // The nominal Breit-Wigner factors with running width.
843   runBW3   = runBW3in;
844   runBW4   = runBW4in; 
845   runBW5   = runBW5in; 
846
847   // Special case: pick scale as if 2 -> 1 process in disguise.
848   if (isSChannel()) {
849
850     // Different options for renormalization scale, but normally sHat.
851     Q2RenSave = renormMultFac * sH;
852     if (renormScale1 == 2) Q2RenSave = renormFixScale; 
853
854     // Different options for factorization scale, but normally sHat.
855     Q2FacSave = factorMultFac * sH;
856     if (factorScale1 == 2) Q2RenSave = factorFixScale; 
857
858   // "Normal" 2 -> 3 processes, i.e. not vector boson fusion.
859   } else if ( idTchan1() != 23 && idTchan1() != 24 && idTchan2() != 23 
860     && idTchan1() != 24 ) {
861     double mT3S = s3 + p3cm.pT2();
862     double mT4S = s4 + p4cm.pT2();
863     double mT5S = s5 + p5cm.pT2();
864     
865     // Different options for renormalization scale.
866     if      (renormScale3 == 1) Q2RenSave = min( mT3S, min(mT4S, mT5S) ); 
867     else if (renormScale3 == 2) Q2RenSave = sqrt( mT3S * mT4S * mT5S
868       / max( mT3S, max(mT4S, mT5S) ) );
869     else if (renormScale3 == 3) Q2RenSave = pow( mT3S * mT4S * mT5S, 
870                                             1./3. );
871     else if (renormScale3 == 4) Q2RenSave = (mT3S * mT4S * mT5S) / 3.;
872     else                        Q2RenSave = sH;
873     Q2RenSave *= renormMultFac;
874     if      (renormScale3 == 6) Q2RenSave = renormFixScale; 
875     
876     // Different options for factorization scale.
877     if      (factorScale3 == 1) Q2FacSave = min( mT3S, min(mT4S, mT5S) ); 
878     else if (factorScale3 == 2) Q2FacSave = sqrt( mT3S * mT4S * mT5S
879       / max( mT3S, max(mT4S, mT5S) ) );
880     else if (factorScale3 == 3) Q2FacSave = pow( mT3S * mT4S * mT5S, 
881                                             1./3. );
882     else if (factorScale3 == 4) Q2FacSave = (mT3S * mT4S * mT5S) / 3.;
883     else                        Q2FacSave = sH;
884     Q2RenSave *= factorMultFac;
885     if      (factorScale3 == 6) Q2FacSave = factorFixScale; 
886
887   // Vector boson fusion 2 -> 3 processes; recoils in positions 4 and 5.
888   } else {
889     double sV4   = pow2( ParticleDataTable::m0(idTchan1()) ); 
890     double sV5   = pow2( ParticleDataTable::m0(idTchan2()) ); 
891     double mT3S  = s3  + p3cm.pT2();
892     double mTV4S = sV4 + p4cm.pT2();
893     double mTV5S = sV5 + p5cm.pT2();
894
895     // Different options for renormalization scale.
896     if      (renormScale3VV == 1) Q2RenSave = max( sV4, sV5); 
897     else if (renormScale3VV == 2) Q2RenSave = sqrt( mTV4S * mTV5S );
898     else if (renormScale3VV == 3) Q2RenSave = pow( mT3S * mTV4S * mTV5S, 
899                                               1./3. );
900     else if (renormScale3VV == 4) Q2RenSave = (mT3S * mTV4S * mTV5S) / 3.;
901     else                          Q2RenSave = sH;
902     Q2RenSave *= renormMultFac;
903     if      (renormScale3VV == 6) Q2RenSave = renormFixScale; 
904     
905     // Different options for factorization scale.
906     if      (factorScale3VV == 1) Q2FacSave = max( sV4, sV5); 
907     else if (factorScale3VV == 2) Q2FacSave = sqrt( mTV4S * mTV5S );
908     else if (factorScale3VV == 3) Q2FacSave = pow( mT3S * mTV4S * mTV5S, 
909                                               1./3. );
910     else if (factorScale3VV == 4) Q2FacSave = (mT3S * mTV4S * mTV5S) / 3.;
911     else                          Q2FacSave = sH;
912     Q2RenSave *= factorMultFac;
913     if      (factorScale3VV == 6) Q2FacSave = factorFixScale; 
914   }
915
916   // Evaluate alpha_strong and alpha_EM.
917   alpS = alphaSPtr->alphaS(Q2RenSave);  
918   alpEM = alphaEMPtr->alphaEM(Q2RenSave);  
919
920 }
921
922 //**************************************************************************
923
924 // The SigmaLHAProcess class.
925 // Wrapper for Les Houches Accord external input; derived from SigmaProcess.
926 // Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks.
927
928 //*********
929
930 // Obtain number of final-state partons from LHA object.
931
932 int SigmaLHAProcess::nFinal() const {
933
934   // At initialization size unknown, so return 0.
935   if (lhaUpPtr->sizePart() <= 0) return 0;
936
937   // Sum up all particles that has first mother = 1.
938   int nFin = 0; 
939   for (int i = 3; i < lhaUpPtr->sizePart(); ++i) 
940     if (lhaUpPtr->mother1(i) == 1) ++nFin;
941   return nFin;
942
943 }
944
945 //**************************************************************************
946
947 } // end namespace Pythia8