- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / src / SigmaNewGaugeBosons.cxx
1 // SigmaNewGaugeBosons.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 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 // leptoquark simulation classes. 
8
9 #include "SigmaNewGaugeBosons.h"
10
11 namespace Pythia8 {
12
13
14 //==========================================================================
15
16 // Sigma1ffbarZprimeWprime class.
17 // Collects common methods for f fbar -> Z'/W' -> WW/WZ -> 4 fermions.
18 // Copied from SigmaEW for gauge-boson-pair production.
19
20 //--------------------------------------------------------------------------
21
22 // Calculate and store internal products.
23
24 void Sigma1ffbarZprimeWprime::setupProd( Event& process, int i1, int i2, 
25   int i3, int i4, int i5, int i6) {
26
27   // Store incoming and outgoing momenta,
28   pRot[1] = process[i1].p();
29   pRot[2] = process[i2].p();
30   pRot[3] = process[i3].p();
31   pRot[4] = process[i4].p();
32   pRot[5] = process[i5].p();
33   pRot[6] = process[i6].p();
34
35   // Do random rotation to avoid accidental zeroes in HA expressions. 
36   bool smallPT = false;
37   do {
38     smallPT = false;
39     double thetaNow = acos(2. * rndmPtr->flat() - 1.);
40     double phiNow   = 2. * M_PI * rndmPtr->flat();
41     for (int i = 1; i <= 6; ++i) { 
42       pRot[i].rot( thetaNow, phiNow);
43       if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
44     }
45   } while (smallPT); 
46
47   // Calculate internal products.
48   for (int i = 1; i < 6; ++i) {
49     for (int j = i + 1; j <= 6; ++j) { 
50       hA[i][j] = 
51           sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz()) 
52         / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() ) 
53         - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz()) 
54         / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() ); 
55       hC[i][j] = conj( hA[i][j] );
56       if (i <= 2) {
57         hA[i][j] *= complex( 0., 1.);
58         hC[i][j] *= complex( 0., 1.);
59       }
60       hA[j][i] = - hA[i][j]; 
61       hC[j][i] = - hC[i][j]; 
62     }
63   }
64
65 }
66
67 //--------------------------------------------------------------------------
68
69 // Evaluate the F function of Gunion and Kunszt.
70
71 complex Sigma1ffbarZprimeWprime::fGK(int j1, int j2, int j3, int j4, 
72   int j5, int j6) {
73  
74   return 4. * hA[j1][j3] * hC[j2][j6] 
75          * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] ); 
76
77 }
78
79 //--------------------------------------------------------------------------
80
81 // Evaluate the Xi function of Gunion and Kunszt.
82
83 double Sigma1ffbarZprimeWprime::xiGK( double tHnow, double uHnow, 
84   double s3now, double s4now) {
85
86   return - 4. * s3now * s4now + tHnow * (3. * tHnow + 4. * uHnow) 
87     + tHnow * tHnow * ( tHnow * uHnow / (s3now * s4now)
88     - 2. * (1. / s3now + 1./s4now) * (tHnow + uHnow)  
89     + 2. * (s3now / s4now + s4now / s3now) );
90
91 }
92
93 //--------------------------------------------------------------------------
94
95 // Evaluate the Xj function of Gunion and Kunszt.
96
97 double Sigma1ffbarZprimeWprime::xjGK( double tHnow, double uHnow,
98   double s3now, double s4now) {
99
100   return 8. * pow2(s3now + s4now) - 8. * (s3now + s4now) * (tHnow + uHnow)
101     - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow 
102     / (s3now * s4now) - 2. * (1. / s3now + 1. / s4now) * (tHnow + uHnow) 
103     + 2. * (s3now / s4now + s4now / s3now) );
104
105 }
106
107 //==========================================================================
108
109 // Sigma1ffbar2gmZZprime class.
110 // Cross section for f fbar -> gamma*/Z0/Z'0 (f is quark or lepton). 
111
112 //--------------------------------------------------------------------------
113
114 // Initialize process. 
115   
116 void Sigma1ffbar2gmZZprime::initProc() {
117
118   // Allow to pick only parts of full gamma*/Z0/Z'0 expression.
119   gmZmode     = settingsPtr->mode("Zprime:gmZmode");
120
121   // Store Z'0 mass and width for propagator. 
122   mRes        = particleDataPtr->m0(32);
123   GammaRes    = particleDataPtr->mWidth(32);
124   m2Res       = mRes*mRes;
125   GamMRat     = GammaRes / mRes;
126   sin2tW      = coupSMPtr->sin2thetaW();
127   cos2tW      = 1. - sin2tW;
128   thetaWRat   = 1. / (16. * sin2tW * cos2tW);
129
130   // Properties of Z0 resonance also needed.
131   mZ          = particleDataPtr->m0(23);
132   GammaZ      = particleDataPtr->mWidth(23);
133   m2Z         = mZ*mZ;
134   GamMRatZ    = GammaZ / mZ;
135
136   // Ensure that arrays initially are empty.
137   for (int i = 0; i < 20; ++i) afZp[i] = 0.;
138   for (int i = 0; i < 20; ++i) vfZp[i] = 0.;
139   
140   // Store first-generation axial and vector couplings.
141   afZp[1]     = settingsPtr->parm("Zprime:ad");
142   afZp[2]     = settingsPtr->parm("Zprime:au");
143   afZp[11]    = settingsPtr->parm("Zprime:ae");
144   afZp[12]    = settingsPtr->parm("Zprime:anue");
145   vfZp[1]     = settingsPtr->parm("Zprime:vd");
146   vfZp[2]     = settingsPtr->parm("Zprime:vu");
147   vfZp[11]    = settingsPtr->parm("Zprime:ve");
148   vfZp[12]    = settingsPtr->parm("Zprime:vnue");
149
150   // Second and third generation could be carbon copy of this...
151   if (settingsPtr->flag("Zprime:universality")) {
152     for (int i = 3; i <= 6; ++i) {
153       afZp[i]    = afZp[i-2];
154       vfZp[i]    = vfZp[i-2];
155       afZp[i+10] = afZp[i+8];
156       vfZp[i+10] = vfZp[i+8];
157     }
158
159   // ... or could have different couplings.   
160   } else {
161     afZp[3]   = settingsPtr->parm("Zprime:as");
162     afZp[4]   = settingsPtr->parm("Zprime:ac");
163     afZp[5]   = settingsPtr->parm("Zprime:ab");
164     afZp[6]   = settingsPtr->parm("Zprime:at");
165     afZp[13]  = settingsPtr->parm("Zprime:amu");
166     afZp[14]  = settingsPtr->parm("Zprime:anumu");
167     afZp[15]  = settingsPtr->parm("Zprime:atau");
168     afZp[16]  = settingsPtr->parm("Zprime:anutau");
169     vfZp[3]   = settingsPtr->parm("Zprime:vs");
170     vfZp[4]   = settingsPtr->parm("Zprime:vc");
171     vfZp[5]   = settingsPtr->parm("Zprime:vb");
172     vfZp[6]   = settingsPtr->parm("Zprime:vt");
173     vfZp[13]  = settingsPtr->parm("Zprime:vmu");
174     vfZp[14]  = settingsPtr->parm("Zprime:vnumu");
175     vfZp[15]  = settingsPtr->parm("Zprime:vtau");
176     vfZp[16]  = settingsPtr->parm("Zprime:vnutau");
177   }
178
179   // Coupling for Z' -> W+ W- and decay angular admixture.
180   coupZpWW    = settingsPtr->parm("Zprime:coup2WW");
181   anglesZpWW  = settingsPtr->parm("Zprime:anglesWW");
182
183   // Set pointer to particle properties and decay table.
184   particlePtr = particleDataPtr->particleDataEntryPtr(32);
185
186
187
188 //--------------------------------------------------------------------------
189
190 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
191
192 void Sigma1ffbar2gmZZprime::sigmaKin() { 
193
194   // Common coupling factors.
195   double colQ = 3. * (1. + alpS / M_PI);
196
197   // Reset quantities to sum. Declare variables in loop.
198   gamSum   = 0.;
199   gamZSum  = 0.;
200   ZSum     = 0.;
201   gamZpSum = 0.;
202   ZZpSum   = 0.;
203   ZpSum    = 0.;
204   int    idAbs, onMode;
205   double mf, mr, ps, kinFacA, kinFacV, ef, af, vf, apf, vpf,
206          ef2, efvf, vaf2, efvpf, vafvapf, vapf2, colf;
207
208   // Loop over all open Z'0 decay channels. 
209   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
210     onMode = particlePtr->channel(i).onMode();
211     if (onMode != 1 && onMode != 2) continue;
212     idAbs = abs( particlePtr->channel(i).product(0) );
213
214     // Contributions from three fermion generations.
215     if ( (idAbs > 0 && idAbs < 7) || ( idAbs > 10 && idAbs < 17) ) {
216       mf = particleDataPtr->m0(idAbs);
217
218       // Check that above threshold. 
219       if (mH > 2. * mf + MASSMARGIN) {
220         mr        = pow2(mf / mH);
221         ps        = sqrtpos(1. - 4. * mr); 
222
223         // Couplings of gamma^*/Z^0/Z'^0  to final flavour 
224         ef        = coupSMPtr->ef(idAbs);
225         af        = coupSMPtr->af(idAbs);
226         vf        = coupSMPtr->vf(idAbs);
227         apf       = afZp[idAbs];
228         vpf       = vfZp[idAbs];
229
230         // Combine couplings with kinematical factors.
231         kinFacA   = pow3(ps);
232         kinFacV   = ps * (1. + 2. * mr);
233         ef2       = ef * ef * kinFacV;
234         efvf      = ef * vf * kinFacV;
235         vaf2      = vf * vf * kinFacV + af * af * kinFacA; 
236         efvpf     = ef * vpf * kinFacV;
237         vafvapf   = vf * vpf * kinFacV + af * apf * kinFacA;
238         vapf2     = vpf * vpf * kinFacV + apf * apf * kinFacA; 
239
240         // Colour factor. Additionally secondary width for top.
241         colf      = (idAbs < 9) ? colQ : 1.;
242         if (idAbs == 6) colf *= particleDataPtr->resOpenFrac(6, -6);
243
244         // Store sum of combinations.
245         gamSum   += colf * ef2;
246         gamZSum  += colf * efvf;
247         ZSum     += colf * vaf2;
248         gamZpSum += colf * efvpf;
249         ZZpSum   += colf * vafvapf;
250         ZpSum    += colf * vapf2;
251       }
252
253     // Optional contribution from W+ W-.
254     } else if (idAbs == 24) {
255       mf = particleDataPtr->m0(idAbs);
256       if (mH > 2. * mf + MASSMARGIN) {
257         mr        = pow2(mf / mH);
258         ps        = sqrtpos(1. - 4. * mr); 
259         ZpSum    += pow2(coupZpWW * cos2tW) * pow3(ps)
260                   * (1. + 20. * mr + 12. * mr*mr)
261                   * particleDataPtr->resOpenFrac(24, -24);
262       }
263     }
264   }
265
266   // Calculate prefactors for gamma/Z0/Z'0 cross section terms.
267   double propZ  = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
268   double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
269   gamNorm   = 4. * M_PI * pow2(alpEM) / (3. * sH); 
270   gamZNorm  = gamNorm * 2. * thetaWRat * (sH - m2Z) * propZ;
271   ZNorm     = gamNorm * pow2(thetaWRat) * sH * propZ;
272   gamZpNorm = gamNorm * 2. * thetaWRat * (sH - m2Res) * propZp;
273   ZZpNorm   = gamNorm * 2. * pow2(thetaWRat) * ((sH - m2Z) * (sH - m2Res)
274               + sH * GamMRatZ * sH * GamMRat) * propZ * propZp;
275   ZpNorm    = gamNorm * pow2(thetaWRat) * sH * propZp;
276
277   // Optionally only keep some of gamma*, Z0 and Z' terms.
278   if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.; 
279     ZZpNorm = 0.; ZpNorm = 0.;}
280   if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.; 
281     ZZpNorm = 0.; ZpNorm = 0.;}
282   if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
283     gamZpNorm = 0.; ZZpNorm = 0.;}
284   if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
285   if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
286   if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}
287
288 }
289
290 //--------------------------------------------------------------------------
291
292 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
293
294 double Sigma1ffbar2gmZZprime::sigmaHat() { 
295
296   // Couplings to an incoming flavour.
297   int idAbs      = abs(id1); 
298   double ei      = coupSMPtr->ef(idAbs);
299   double ai      = coupSMPtr->af(idAbs);
300   double vi      = coupSMPtr->vf(idAbs);
301   double api     = afZp[idAbs];
302   double vpi     = vfZp[idAbs];
303   double ei2     = ei * ei;
304   double eivi    = ei * vi;
305   double vai2    = vi * vi + ai * ai; 
306   double eivpi   = ei * vpi;
307   double vaivapi = vi * vpi + ai * api;; 
308   double vapi2   = vpi * vpi + api * api;
309
310   // Combine gamma, interference and Z0 parts.
311   double sigma = ei2 * gamNorm * gamSum + eivi * gamZNorm * gamZSum
312                + vai2 * ZNorm * ZSum + eivpi * gamZpNorm * gamZpSum
313                + vaivapi * ZZpNorm * ZZpSum + vapi2 * ZpNorm * ZpSum;
314
315   // Colour factor. Answer.
316   if (idAbs < 9) sigma /= 3.;
317   return sigma;    
318
319 }
320
321 //--------------------------------------------------------------------------
322
323 // Select identity, colour and anticolour.
324
325 void Sigma1ffbar2gmZZprime::setIdColAcol() {
326
327   // Flavours trivial.
328   setId( id1, id2, 32);
329
330   // Colour flow topologies. Swap when antiquarks.
331   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
332   else              setColAcol( 0, 0, 0, 0, 0, 0);
333   if (id1 < 0) swapColAcol();
334
335 }
336
337 //--------------------------------------------------------------------------
338
339 // Evaluate weight for gamma*/Z0/Z'0 decay angle.
340   
341 double Sigma1ffbar2gmZZprime::weightDecay( Event& process, int iResBeg, 
342   int iResEnd) {
343
344   // Default values, in- and out-flavours in process.
345   double wt    = 1.;
346   double wtMax = 1.;
347   int idInAbs  = process[3].idAbs();
348   int idOutAbs = process[6].idAbs();
349
350   // Angular weight for outgoing fermion pair.
351   if (iResBeg == 5 && iResEnd == 5 &&
352     (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
353
354     // Couplings for in- and out-flavours.
355     double ei  = coupSMPtr->ef(idInAbs);
356     double vi  = coupSMPtr->vf(idInAbs);
357     double ai  = coupSMPtr->af(idInAbs);
358     double vpi = vfZp[idInAbs];
359     double api = afZp[idInAbs];
360     double ef  = coupSMPtr->ef(idOutAbs);
361     double vf  = coupSMPtr->vf(idOutAbs);
362     double af  = coupSMPtr->af(idOutAbs);
363     double vpf = vfZp[idOutAbs];
364     double apf = afZp[idOutAbs];
365
366     // Phase space factors. (One power of beta left out in formulae.)
367     double mr1 = pow2(process[6].m()) / sH;
368     double mr2 = pow2(process[7].m()) / sH;
369     double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
370     double mrAvg = 0.5 * (mr1 + mr2) - 0.25 * pow2(mr1 - mr2);
371
372     // Coefficients of angular expression.
373     double coefTran = ei*ei * gamNorm * ef*ef + ei * vi * gamZNorm * ef * vf
374       + (vi*vi + ai*ai) * ZNorm * (vf*vf + ps*ps * af*af)
375       + ei * vpi * gamZpNorm * ef * vpf
376       + (vi * vpi + ai * api) * ZZpNorm * (vf * vpf + ps*ps * af * apf) 
377       + (vpi*vpi + api*api) * ZpNorm * (vpf*vpf + ps*ps * apf*apf);
378     double coefLong = 4. * mrAvg * ( ei*ei * gamNorm * ef*ef 
379       + ei * vi * gamZNorm * ef * vf + (vi*vi + ai*ai) * ZNorm * vf*vf
380       + ei * vpi * gamZpNorm * ef * vpf 
381       + (vi * vpi + ai * api) * ZZpNorm * vf * vpf
382       + (vpi*vpi + api*api) * ZpNorm * vpf*vpf );
383     double coefAsym = ps * ( ei * ai * gamZNorm * ef * af 
384       + 4. * vi * ai * ZNorm * vf * af + ei * api * gamZpNorm * ef * apf 
385       + (vi * api + vpi * ai) * ZZpNorm * (vf * apf + vpf * af)
386       + 4. * vpi * api * ZpNorm * vpf * apf );
387
388     // Flip asymmetry for in-fermion + out-antifermion.
389     if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
390
391     // Reconstruct decay angle and weight for it.
392     double cosThe = (process[3].p() - process[4].p()) 
393       * (process[7].p() - process[6].p()) / (sH * ps);
394     wt    = coefTran * (1. + pow2(cosThe)) 
395        + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
396     wtMax = 2. * (coefTran + abs(coefAsym));
397   }
398
399   // Angular weight for Z' -> W+ W-.
400   else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
401     double mr1 = pow2(process[6].m()) / sH;
402     double mr2 = pow2(process[7].m()) / sH;
403     double ps  = sqrtpos(pow2(1. - mr1 -mr2) - 4. * mr1 * mr2);
404     double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2 
405       + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
406     double cFlat = -cCos2 + 0.5 * (mr1 + mr2) 
407       * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
408
409     // Reconstruct decay angle and weight for it.
410     double cosThe = (process[3].p() - process[4].p()) 
411       * (process[7].p() - process[6].p()) / (sH * ps);
412     wt    = cFlat + cCos2 * cosThe*cosThe; 
413     wtMax = cFlat + max(0., cCos2);
414   }
415
416   // Angular weight for f + fbar -> Z' -> W+ + W- -> 4 fermions.
417   else if (iResBeg == 6 && iResEnd == 7 && idOutAbs == 24) {
418  
419     // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
420     // with f' fbar' from W- and f" fbar" from W+.
421     int i1 = (process[3].id() < 0) ? 3 : 4;
422     int i2 = 7 - i1; 
423     int i3 = (process[8].id() > 0) ? 8 : 9;
424     int i4 = 17 - i3;
425     int i5 = (process[10].id() > 0) ? 10 : 11;
426     int i6 = 21 - i5;
427     if (process[6].id() > 0) {swap(i3, i5); swap(i4, i6);}
428   
429     // Decay distribution like in f fbar -> Z^* -> W+ W-. 
430     if (rndmPtr->flat() > anglesZpWW) {
431
432       // Set up four-products and internal products.
433       setupProd( process, i1, i2, i3, i4, i5, i6); 
434
435       // tHat and uHat of fbar f -> W- W+, and their squared masses.
436       int iNeg     = (process[6].id() < 0) ? 6 : 7;
437       int iPos     = 13 - iNeg;
438       double tHres = (process[i1].p() - process[iNeg].p()).m2Calc();
439       double uHres = (process[i1].p() - process[iPos].p()).m2Calc(); 
440       double s3now = process[iNeg].m2();
441       double s4now = process[iPos].m2();
442
443       // Kinematics combinations (norm(x) = |x|^2).
444       double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
445       double fGK253 = norm(fGK( 2, 1, 5, 6, 3, 4) - fGK( 2, 1, 3, 4, 5, 6) );
446       double xiT    = xiGK( tHres, uHres, s3now, s4now);
447       double xiU    = xiGK( uHres, tHres, s3now, s4now);
448       double xjTU   = xjGK( tHres, uHres, s3now, s4now);
449
450       //  Couplings of incoming (anti)fermion. Combine with kinematics.
451       int idAbs     = process[i1].idAbs();
452       double li     = 0.5 * (vfZp[idAbs] + afZp[idAbs]); 
453       double ri     = 0.5 * (vfZp[idAbs] - afZp[idAbs]);  
454       wt            = li*li * fGK135 + ri*ri * fGK253;
455       wtMax         = 4. * s3now * s4now * (li*li + ri*ri) 
456                     * (xiT + xiU - xjTU);
457    
458     // Decay distribution like in f fbar -> h^0 -> W+ W-. 
459     } else {
460       double p35  = 2. * process[i3].p() * process[i5].p(); 
461       double p46  = 2. * process[i4].p() * process[i6].p(); 
462       wt          = 16. * p35 * p46;
463       wtMax       = sH2;
464     }
465   }
466
467   // Angular weight in top decay by standard routine.
468   else if (process[process[iResBeg].mother1()].idAbs() == 6)
469     return weightTopDecay( process, iResBeg, iResEnd);
470
471
472   // Done.
473   return (wt / wtMax);
474
475 }
476
477 //==========================================================================
478
479 // Sigma1ffbar2Wprime class.
480 // Cross section for f fbar' -> W'+- (f is quark or lepton). 
481
482 //--------------------------------------------------------------------------
483
484 // Initialize process. 
485   
486 void Sigma1ffbar2Wprime::initProc() {
487
488   // Store W+- mass and width for propagator. 
489   mRes     = particleDataPtr->m0(34);
490   GammaRes = particleDataPtr->mWidth(34);
491   m2Res    = mRes*mRes;
492   GamMRat  = GammaRes / mRes;
493   thetaWRat = 1. / (12. * coupSMPtr->sin2thetaW());
494
495   // Axial and vector couplings of fermions.
496   aqWp      = settingsPtr->parm("Wprime:aq");
497   vqWp      = settingsPtr->parm("Wprime:vq");
498   alWp      = settingsPtr->parm("Wprime:al");
499   vlWp      = settingsPtr->parm("Wprime:vl");
500
501   // Coupling for W' -> W Z and decay angular admixture.
502   coupWpWZ    = settingsPtr->parm("Wprime:coup2WZ");
503   anglesWpWZ  = settingsPtr->parm("Wprime:anglesWZ");
504
505   // Set pointer to particle properties and decay table.
506   particlePtr = particleDataPtr->particleDataEntryPtr(34);
507
508
509
510 //--------------------------------------------------------------------------
511
512 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
513
514 void Sigma1ffbar2Wprime::sigmaKin() { 
515
516   // Set up Breit-Wigner. Cross section for W+ and W- separately.
517   double sigBW  = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
518   double preFac = alpEM * thetaWRat * mH;
519   sigma0Pos     = preFac * sigBW * particlePtr->resWidthOpen(34, mH);    
520   sigma0Neg     = preFac * sigBW * particlePtr->resWidthOpen(-34, mH);    
521
522 }
523
524 //--------------------------------------------------------------------------
525
526 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
527
528 double Sigma1ffbar2Wprime::sigmaHat() {
529
530   // Secondary width for W+ or W-. CKM and colour factors.
531   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
532   double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
533   if (abs(id1) < 7) sigma *= coupSMPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
534
535   // Couplings.
536   if (abs(id1) < 7) sigma *= 0.5 * (aqWp * aqWp + vqWp * vqWp);
537   else              sigma *= 0.5 * (alWp * alWp + vlWp * vlWp);        
538
539   // Answer.
540   return sigma;    
541
542 }
543
544 //--------------------------------------------------------------------------
545
546 // Select identity, colour and anticolour.
547
548 void Sigma1ffbar2Wprime::setIdColAcol() {
549
550   // Sign of outgoing W.
551   int sign          = 1 - 2 * (abs(id1)%2);
552   if (id1 < 0) sign = -sign;
553   setId( id1, id2, 34 * sign);
554
555   // Colour flow topologies. Swap when antiquarks.
556   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
557   else              setColAcol( 0, 0, 0, 0, 0, 0);
558   if (id1 < 0) swapColAcol();
559
560 }
561
562 //--------------------------------------------------------------------------
563
564 // Evaluate weight for W decay angle.
565   
566 double Sigma1ffbar2Wprime::weightDecay( Event& process, int iResBeg, 
567   int iResEnd) {
568
569   // Default values, in- and out-flavours in process.
570   double wt    = 1.;
571   double wtMax = 1.;
572   int idInAbs  = process[3].idAbs();
573   int idOutAbs = process[6].idAbs();
574
575   // Angular weight for outgoing fermion pair.
576   if (iResBeg == 5 && iResEnd == 5 &&
577     (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
578
579     // Couplings for in- and out-flavours.
580     double ai  = (idInAbs < 9) ? aqWp : alWp;
581     double vi  = (idInAbs < 9) ? vqWp : vlWp;
582     double af  = (idOutAbs < 9) ? aqWp : alWp;
583     double vf  = (idOutAbs < 9) ? vqWp : vlWp;
584
585     // Asymmetry expression.
586     double coefAsym = 8. * vi * ai * vf * af 
587       / ((vi*vi + ai*ai) * (vf*vf + af*af));
588
589     // Flip asymmetry for in-fermion + out-antifermion.
590     if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
591
592     // Phase space factors.
593     double mr1 = pow2(process[6].m()) / sH;
594     double mr2 = pow2(process[7].m()) / sH;
595     double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
596
597     // Reconstruct decay angle and weight for it.
598     double cosThe = (process[3].p() - process[4].p()) 
599       * (process[7].p() - process[6].p()) / (sH * ps);
600     wt    = 1. + coefAsym * cosThe + cosThe * cosThe;
601     wtMax = 2. + abs(coefAsym);
602   }
603
604   // Angular weight for W' -> W Z.
605   else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
606     double mr1 = pow2(process[6].m()) / sH;
607     double mr2 = pow2(process[7].m()) / sH;
608     double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
609     double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2 
610       + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
611     double cFlat = -cCos2 + 0.5 * (mr1 + mr2) 
612       * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
613
614     // Reconstruct decay angle and weight for it.
615     double cosThe = (process[3].p() - process[4].p()) 
616       * (process[7].p() - process[6].p()) / (sH * ps);
617     wt    = cFlat + cCos2 * cosThe*cosThe; 
618     wtMax = cFlat + max(0., cCos2);
619   }
620
621   // Angular weight for f + fbar -> W' -> W + Z -> 4 fermions.
622   else if (iResBeg == 6 && iResEnd == 7 
623     && (idOutAbs == 24 || idOutAbs == 23)) {
624  
625     // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
626     // with f' fbar' from W and f" fbar" from Z.
627     int i1 = (process[3].id() < 0) ? 3 : 4;
628     int i2 = 7 - i1; 
629     int i3 = (process[8].id() > 0) ? 8 : 9;
630     int i4 = 17 - i3;
631     int i5 = (process[10].id() > 0) ? 10 : 11;
632     int i6 = 21 - i5;
633     if (process[6].id() == 23) {swap(i3, i5); swap(i4, i6);}
634   
635     // Decay distribution like in f fbar -> Z^* -> W+ W-. 
636     if (rndmPtr->flat() > anglesWpWZ) {
637
638       // Set up four-products and internal products.
639       setupProd( process, i1, i2, i3, i4, i5, i6); 
640
641       // tHat and uHat of fbar f -> W Z, and their squared masses.
642       int iW       = (process[6].id() == 23) ? 7 : 6;
643       int iZ       = 13 - iW;
644       double tHres = (process[i1].p() - process[iW].p()).m2Calc();
645       double uHres = (process[i1].p() - process[iZ].p()).m2Calc(); 
646       double s3now = process[iW].m2();
647       double s4now = process[iZ].m2();
648
649       // Kinematics combinations (norm(x) = |x|^2).
650       double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
651       double fGK136 = norm(fGK( 1, 2, 3, 4, 6, 5) - fGK( 1, 2, 6, 5, 3, 4) );
652       double xiT    = xiGK( tHres, uHres, s3now, s4now);
653       double xiU    = xiGK( uHres, tHres, s3now, s4now);
654       double xjTU   = xjGK( tHres, uHres, s3now, s4now);
655
656       //  Couplings of outgoing fermion from Z. Combine with kinematics.
657       int idAbs     = process[i5].idAbs();
658       double lfZ    = coupSMPtr->lf(idAbs); 
659       double rfZ    = coupSMPtr->rf(idAbs); 
660       wt            = lfZ*lfZ * fGK135 + rfZ*rfZ * fGK136;
661       wtMax         = 4. * s3now * s4now * (lfZ*lfZ + rfZ*rfZ) 
662                     * (xiT + xiU - xjTU);
663    
664     // Decay distribution like in f fbar -> H^+- -> W+- Z0. 
665     } else {
666       double p35  = 2. * process[i3].p() * process[i5].p(); 
667       double p46  = 2. * process[i4].p() * process[i6].p(); 
668       wt          = 16. * p35 * p46;
669       wtMax       = sH2;
670     }
671   }
672
673   // Angular weight in top decay by standard routine.
674   else if (process[process[iResBeg].mother1()].idAbs() == 6)
675     return weightTopDecay( process, iResBeg, iResEnd);
676  
677   // Done.
678   return (wt / wtMax);
679
680 }
681
682
683 //==========================================================================
684
685 // Sigma1ffbar2Rhorizontal class.
686 // Cross section for f fbar' -> R^0 (f is a quark or lepton). 
687
688 //--------------------------------------------------------------------------
689
690 // Initialize process. 
691   
692 void Sigma1ffbar2Rhorizontal::initProc() {
693
694   // Store R^0 mass and width for propagator. 
695   mRes     = particleDataPtr->m0(41);
696   GammaRes = particleDataPtr->mWidth(41);
697   m2Res    = mRes*mRes;
698   GamMRat  = GammaRes / mRes;
699   thetaWRat = 1. / (12. * coupSMPtr->sin2thetaW());
700
701   // Set pointer to particle properties and decay table.
702   particlePtr = particleDataPtr->particleDataEntryPtr(41);
703
704
705
706 //--------------------------------------------------------------------------
707
708 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
709
710 void Sigma1ffbar2Rhorizontal::sigmaKin() { 
711
712   // Set up Breit-Wigner. Cross section for W+ and W- separately.
713   double sigBW  = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
714   double preFac = alpEM * thetaWRat * mH;
715   sigma0Pos     = preFac * sigBW * particlePtr->resWidthOpen(41, mH);    
716   sigma0Neg     = preFac * sigBW * particlePtr->resWidthOpen(-41, mH);    
717
718 }
719
720 //--------------------------------------------------------------------------
721
722 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
723
724 double Sigma1ffbar2Rhorizontal::sigmaHat() {
725
726   // Check for allowed flavour combinations, one generation apart.
727   if (id1 * id2 > 0 || abs(id1 + id2) != 2) return 0.;  
728
729   // Find whether R0 or R0bar. Colour factors.
730   double sigma = (id1 + id2 > 0) ? sigma0Pos : sigma0Neg;
731   if (abs(id1) < 7) sigma /= 3.;
732
733   // Answer.
734   return sigma;    
735
736 }
737
738 //--------------------------------------------------------------------------
739
740 // Select identity, colour and anticolour.
741
742 void Sigma1ffbar2Rhorizontal::setIdColAcol() {
743
744   // Outgoing R0 or R0bar.
745   id3 = (id1 +id2 > 0) ? 41 : -41;
746   setId( id1, id2, id3);
747
748   // Colour flow topologies. Swap when antiquarks.
749   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
750   else              setColAcol( 0, 0, 0, 0, 0, 0);
751   if (id1 < 0) swapColAcol();
752
753 }
754
755 //==========================================================================
756
757 } // end namespace Pythia8