using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SigmaLeftRightSym.cxx
1 // SigmaLeftRightSym.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 // left-right-symmetry simulation classes. 
8
9 #include "SigmaLeftRightSym.h"
10
11 namespace Pythia8 {
12
13 //**************************************************************************
14
15 // Sigma1ffbar2ZRight class.
16 // Cross section for f fbar -> Z_R^0 (righthanded gauge boson). 
17
18 //*********
19
20 // Initialize process. 
21   
22 void Sigma1ffbar2ZRight::initProc() {
23
24   // Store Z_R mass and width for propagator. 
25   idZR     = 9900023;
26   mRes     = ParticleDataTable::m0(idZR);
27   GammaRes = ParticleDataTable::mWidth(idZR);
28   m2Res    = mRes*mRes;
29   GamMRat  = GammaRes / mRes;
30   sin2tW   = CoupEW::sin2thetaW();
31
32   // Set pointer to particle properties and decay table.
33   ZRPtr    = ParticleDataTable::particleDataPtr(idZR);
34
35
36
37 //*********
38
39 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
40
41 void Sigma1ffbar2ZRight::sigmaKin() { 
42
43   // Set up Breit-Wigner. Width out only includes open channels. 
44   double sigBW    = 12. * M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );    
45   double widthOut = ZRPtr->resWidthOpen(idZR, mH);
46
47   // Prefactor for incoming widths. Combine. Done. 
48   double preFac   = alpEM * mH / ( 48. * sin2tW * (1. - sin2tW) 
49                   * (1. - 2. * sin2tW) ); 
50   sigma0          = preFac * sigBW * widthOut;    
51
52 }
53
54 //*********
55
56 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
57
58 double Sigma1ffbar2ZRight::sigmaHat() { 
59
60   // Vector and axial couplings of incoming fermion pair.
61   int    idAbs = abs(id1); 
62   double af = 0.;
63   double vf = 0.;
64   if (idAbs < 9 && idAbs%2 == 1) {
65     af      = -1. + 2. * sin2tW; 
66     vf      = -1. + 4. * sin2tW / 3.;
67   } else if (idAbs < 9) {   
68     af      = 1. - 2. * sin2tW; 
69     vf      = 1. - 8. * sin2tW / 3.;
70   } else if (idAbs < 19 && idAbs%2 == 1) {   
71     af      = -1. + 2. * sin2tW; 
72     vf      = -1. + 4. * sin2tW;
73   }
74
75   // Colour factor. Answer.
76   double sigma = (vf*vf + af*af) * sigma0;
77   if (idAbs < 9) sigma /= 3.;
78   return sigma;    
79
80 }
81
82 //*********
83
84 // Select identity, colour and anticolour.
85
86 void Sigma1ffbar2ZRight::setIdColAcol() {
87
88   // Flavours trivial.
89   setId( id1, id2, idZR);
90
91   // Colour flow topologies. Swap when antiquarks.
92   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
93   else              setColAcol( 0, 0, 0, 0, 0, 0);
94   if (id1 < 0) swapColAcol();
95
96 }
97
98 //*********
99
100 // Evaluate weight for Z_R decay angle.
101   
102 double Sigma1ffbar2ZRight::weightDecay( Event& process, int iResBeg, 
103   int iResEnd) {
104
105   // Identity of mother of decaying reseonance(s).
106   int idMother = process[process[iResBeg].mother1()].idAbs();
107
108   // For top decay hand over to standard routine.
109   if (idMother == 6) 
110     return weightTopDecay( process, iResBeg, iResEnd);
111
112   // Z_R should sit in entry 5.
113   if (iResBeg != 5 || iResEnd != 5) return 1.;
114
115   // Couplings for in- and out-flavours.
116   double ai, vi, af, vf;
117   int idInAbs   = process[3].idAbs();
118   if (idInAbs < 9 && idInAbs%2 == 1) {
119     ai          = -1. + 2. * sin2tW; 
120     vi          = -1. + 4. * sin2tW / 3.;
121   } else if (idInAbs < 9) {   
122     ai          = 1. - 2. * sin2tW; 
123     vi          = 1. - 8. * sin2tW / 3.;
124   } else {   
125     ai          = -1. + 2. * sin2tW; 
126     vi          = -1. + 4. * sin2tW;
127   }
128   int idOutAbs = process[6].idAbs();
129   if (idOutAbs < 9 && idOutAbs%2 == 1) {
130     af          = -1. + 2. * sin2tW; 
131     vf          = -1. + 4. * sin2tW / 3.;
132   } else if (idOutAbs < 9) {   
133     af          = 1. - 2. * sin2tW; 
134     vf          = 1. - 8. * sin2tW / 3.;
135   } else {   
136     af          = -1. + 2. * sin2tW; 
137     vf          = -1. + 4. * sin2tW;
138   }
139
140   // Phase space factors. Reconstruct decay angle.
141   double mr1    = pow2(process[6].m()) / sH;
142   double mr2    = pow2(process[7].m()) / sH;
143   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
144   double cosThe = (process[3].p() - process[4].p()) 
145     * (process[7].p() - process[6].p()) / (sH * betaf);
146
147   // Angular weight and its maximum.
148   double wt1    = (vi*vi + ai*ai) * (vf*vf + af*af * betaf*betaf);
149   double wt2    = (1. - betaf*betaf) * (vi*vi + ai*ai) * vf*vf;
150   double wt3    = betaf * 4. * vi * ai * vf * af;  
151   if (process[3].id() * process[6].id() < 0) wt3 = -wt3;
152   double wt     = wt1 * (1. + cosThe*cosThe) + wt2 * (1. - cosThe*cosThe)
153                 + 2. * wt3 * cosThe;
154   double wtMax  = 2. * (wt1 + abs(wt3));
155
156   // Done.
157   return wt / wtMax;
158
159 }
160
161 //**************************************************************************
162
163 // Sigma1ffbar2WRight class.
164 // Cross section for f fbar' -> W_R^+- (righthanded gauge boson). 
165
166 //*********
167
168 // Initialize process. 
169   
170 void Sigma1ffbar2WRight::initProc() {
171
172   // Store W_R^+- mass and width for propagator. 
173   idWR     = 9900024;
174   mRes     = ParticleDataTable::m0(idWR);
175   GammaRes = ParticleDataTable::mWidth(idWR);
176   m2Res    = mRes*mRes;
177   GamMRat  = GammaRes / mRes;
178   thetaWRat = 1. / (12. * CoupEW::sin2thetaW());
179
180   // Set pointer to particle properties and decay table.
181   particlePtr = ParticleDataTable::particleDataPtr(idWR);
182
183
184
185 //*********
186
187 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
188
189 void Sigma1ffbar2WRight::sigmaKin() { 
190
191   // Common coupling factors.
192   double colQ   = 3. * (1. + alpS / M_PI);
193
194   // Reset quantities to sum. Declare variables inside loop.
195   double widOutPos = 0.; 
196   double widOutNeg = 0.; 
197   int    id1Now, id2Now, id1Abs, id2Abs, id1Neg, id2Neg, onMode;
198   double widNow, widSecPos, widSecNeg, mf1, mf2, mr1, mr2, kinFac;
199
200   // Loop over all W_R^+- decay channels. 
201   for (int i = 0; i < particlePtr->decay.size(); ++i) {
202     widNow = 0.;
203     id1Now      = particlePtr->decay[i].product(0);
204     id2Now      = particlePtr->decay[i].product(1);
205     id1Abs      = abs(id1Now);
206     id2Abs      = abs(id2Now);
207
208     // Check that above threshold. Phase space.
209     mf1 = ParticleDataTable::m0(id1Abs);
210     mf2 = ParticleDataTable::m0(id2Abs);
211     if (mH > mf1 + mf2 + MASSMARGIN) {
212       mr1    = pow2(mf1 / mH);
213       mr2    = pow2(mf2 / mH);
214       kinFac = (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
215              * sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 ); 
216
217       // Combine kinematics with colour factor and CKM couplings.
218       widNow = kinFac;
219       if (id1Abs < 9) widNow *= colQ * VCKM::V2id(id1Abs, id2Abs);
220  
221       // Secondary width from top and righthanded neutrino decay.
222       id1Neg    = (id1Abs < 19) ? -id1Now : id1Abs; 
223       id2Neg    = (id2Abs < 19) ? -id2Now : id2Abs; 
224       widSecPos = ParticleDataTable::resOpenFrac(id1Now, id2Now); 
225       widSecNeg = ParticleDataTable::resOpenFrac(id1Neg, id2Neg); 
226
227       // Add weight for channels on for all, W_R^+ and W_R^-, respectively.
228       onMode = particlePtr->decay[i].onMode();
229       if (onMode == 1 || onMode == 2) widOutPos += widNow * widSecPos;
230       if (onMode == 1 || onMode == 3) widOutNeg += widNow * widSecNeg;
231
232     // End loop over fermions.
233     }
234   }
235
236   // Set up Breit-Wigner. Cross section for W_R^+ and W_R^- separately.
237   double sigBW = 12. * M_PI * pow2(alpEM * thetaWRat) * sH
238                / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
239   sigma0Pos = sigBW * widOutPos;    
240   sigma0Neg = sigBW * widOutNeg;    
241
242 }
243
244 //*********
245
246 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
247
248 double Sigma1ffbar2WRight::sigmaHat() {
249
250   // Secondary width for W_R^+ or W_R^-. CKM and colour factors.
251   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
252   double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
253   if (abs(id1) < 9) sigma *= VCKM::V2id(abs(id1), abs(id2)) / 3.;
254
255   // Answer.
256   return sigma;    
257
258 }
259
260 //*********
261
262 // Select identity, colour and anticolour.
263
264 void Sigma1ffbar2WRight::setIdColAcol() {
265
266   // Sign of outgoing W_R.
267   int sign          = (abs(id1)%2 == 0) ? 1 : -1;
268   if (id1 < 0) sign = -sign;
269   setId( id1, id2, idWR * sign);
270
271   // Colour flow topologies. Swap when antiquarks.
272   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
273   else              setColAcol( 0, 0, 0, 0, 0, 0);
274   if (id1 < 0) swapColAcol();
275
276 }
277
278 //*********
279
280 // Evaluate weight for W_R decay angle.
281   
282 double Sigma1ffbar2WRight::weightDecay( Event& process, int iResBeg, 
283   int iResEnd) {
284
285   // Identity of mother of decaying reseonance(s).
286   int idMother = process[process[iResBeg].mother1()].idAbs();
287
288   // For top decay hand over to standard routine.
289   if (idMother == 6) 
290     return weightTopDecay( process, iResBeg, iResEnd);
291
292   // W_R should sit in entry 5.
293   if (iResBeg != 5 || iResEnd != 5) return 1.;
294
295   // Phase space factors.
296   double mr1    = pow2(process[6].m()) / sH;
297   double mr2    = pow2(process[7].m()) / sH;
298   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
299    
300   // Sign of asymmetry.
301   double eps    = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
302
303   // Reconstruct decay angle and weight for it.
304   double cosThe = (process[3].p() - process[4].p()) 
305     * (process[7].p() - process[6].p()) / (sH * betaf);
306   double wtMax  = 4.;
307   double wt     = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2); 
308  
309   // Done.
310   return (wt / wtMax);
311
312 }
313
314 //**************************************************************************
315
316 // Sigma1ll2Hchgchg class.
317 // Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
318
319 //*********
320
321 // Initialize process. 
322   
323 void Sigma1ll2Hchgchg::initProc() {
324
325   // Set process properties: H_L^++-- or H_R^++--.
326   if (leftRight == 1) {
327     idHLR    = 9900041;
328     codeSave = 3121;
329     nameSave = "l l -> H_L^++--";
330   } else {
331     idHLR    = 9900042;
332     codeSave = 3141;
333     nameSave = "l l -> H_R^++--";
334   }
335
336   // Read in Yukawa matrix for couplings to a lepton pair.
337   yukawa[1][1]  = Settings::parm("LeftRightSymmmetry:coupHee");
338   yukawa[2][1]  = Settings::parm("LeftRightSymmmetry:coupHmue");
339   yukawa[2][2]  = Settings::parm("LeftRightSymmmetry:coupHmumu");
340   yukawa[3][1]  = Settings::parm("LeftRightSymmmetry:coupHtaue");
341   yukawa[3][2]  = Settings::parm("LeftRightSymmmetry:coupHtaumu");
342   yukawa[3][3]  = Settings::parm("LeftRightSymmmetry:coupHtautau");
343
344   // Store H_L/R mass and width for propagator. 
345   mRes     = ParticleDataTable::m0(idHLR);
346   GammaRes = ParticleDataTable::mWidth(idHLR);
347   m2Res    = mRes*mRes;
348   GamMRat  = GammaRes / mRes;
349
350   // Set pointer to particle properties and decay table.
351   particlePtr = ParticleDataTable::particleDataPtr(idHLR);
352
353
354
355 //*********
356
357 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
358
359 double Sigma1ll2Hchgchg::sigmaHat() {
360
361   // Initial state must consist of two identical-sign leptons.
362   if (id1 * id2 < 0) return 0.;
363   int id1Abs = abs(id1);
364   int id2Abs = abs(id2);  
365   if (id1Abs != 11 && id1Abs != 13 && id1Abs != 15) return 0.;
366   if (id2Abs != 11 && id2Abs != 13 && id2Abs != 15) return 0.;
367
368   // Set up Breit-Wigner, inwidth and outwidth.
369   double sigBW  = 8. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
370   double widIn  = pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) 
371                 * mH / (8. * M_PI);
372   int idSgn     = (id1 < 0) ? idHLR : -idHLR;
373   double widOut = particlePtr->resWidthOpen( idSgn, mH);
374
375   // Answer.
376   return widIn * sigBW * widOut;    
377
378 }
379
380 //*********
381
382 // Select identity, colour and anticolour.
383
384 void Sigma1ll2Hchgchg::setIdColAcol() {
385
386   // Sign of outgoing H_L/R.
387   int idSgn     = (id1 < 0) ? idHLR : -idHLR;
388   setId( id1, id2, idSgn);
389
390   // No colours whatsoever.
391   setColAcol( 0, 0, 0, 0, 0, 0);
392
393 }
394
395 //*********
396
397 // Evaluate weight for H_L/R sequential decay angles.
398   
399 double Sigma1ll2Hchgchg::weightDecay( Event& process, int iResBeg, 
400   int iResEnd) {
401
402   // Identity of mother of decaying reseonance(s).
403   int idMother = process[process[iResBeg].mother1()].idAbs();
404
405   // For top decay hand over to standard routine.
406   if (idMother == 6) 
407     return weightTopDecay( process, iResBeg, iResEnd);
408  
409   // Else isotropic decay.
410   return 1.;
411
412 }
413
414 //**************************************************************************
415
416 // Sigma2lgm2Hchgchgl class.
417 // Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
418
419 //*********
420
421 // Initialize process. 
422   
423 void Sigma2lgm2Hchgchgl::initProc() {
424
425   // Set process properties: H_L^++-- or H_R^++-- and e/mu/tau.
426   idHLR        = (leftRight == 1) ? 9900041 : 9900042;
427   codeSave     = (leftRight == 1) ? 3122 : 3142;
428   if (idLep == 13) codeSave += 2;
429   if (idLep == 15) codeSave += 4;
430   if      (codeSave == 3122) nameSave = "l^+- gamma -> H_L^++-- e^-+";
431   else if (codeSave == 3123) nameSave = "l^+- gamma -> H_L^++-- mu^-+";
432   else if (codeSave == 3124) nameSave = "l^+- gamma -> H_L^++-- tau^-+";
433   else if (codeSave == 3142) nameSave = "l^+- gamma -> H_R^++-- e^-+";
434   else if (codeSave == 3143) nameSave = "l^+- gamma -> H_R^++-- mu^-+";
435   else                       nameSave = "l^+- gamma -> H_R^++-- tau^-+";
436
437   // Read in relevantYukawa matrix for couplings to a lepton pair.
438   if (idLep == 11) {
439     yukawa[1]  = Settings::parm("LeftRightSymmmetry:coupHee");
440     yukawa[2]  = Settings::parm("LeftRightSymmmetry:coupHmue");
441     yukawa[3]  = Settings::parm("LeftRightSymmmetry:coupHtaue");
442   } else if (idLep == 13) { 
443     yukawa[1]  = Settings::parm("LeftRightSymmmetry:coupHmue");
444     yukawa[2]  = Settings::parm("LeftRightSymmmetry:coupHmumu");
445     yukawa[3]  = Settings::parm("LeftRightSymmmetry:coupHtaumu");
446   } else {
447     yukawa[1]  = Settings::parm("LeftRightSymmmetry:coupHtaue");
448     yukawa[2]  = Settings::parm("LeftRightSymmmetry:coupHtaumu");
449     yukawa[3]  = Settings::parm("LeftRightSymmmetry:coupHtautau");
450   }
451
452   // Secondary open width fractions.
453   openFracPos  = ParticleDataTable::resOpenFrac( idHLR);
454   openFracNeg  = ParticleDataTable::resOpenFrac(-idHLR);
455
456
457
458 //*********
459
460 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
461
462 double Sigma2lgm2Hchgchgl::sigmaHat() {
463
464   // Initial state must consist of a lepton and a photon.
465   int idIn     = (id2 == 22) ? id1 : id2;  
466   int idInAbs  = abs(idIn);
467   if (idInAbs != 11 && idInAbs != 13 && idInAbs != 15) return 0.;
468
469   // Incoming squared lepton mass.
470   double s1    = pow2( ParticleDataTable::m0(idInAbs) );
471   
472   // Kinematical expressions.
473   double smm1  = 8. * (sH + tH - s3) * (sH + tH - 2. * s3 - s1 - s4) 
474                / pow2(uH - s3);
475   double smm2  = 2. * ( (2. * s3 - 3. * s1) * s4 + (s1 - 2. * s4) * tH
476                - (tH - s4) * sH ) / pow2(tH - s4);
477   double smm3  = 2. * ( (2. * s3 - 3. * s4 + tH) * s1 
478                - (2. * s1 - s4 + tH) * sH ) / pow2(sH - s1);
479   double smm12 = 4. * ( (2. * s1 - s4 - 2. * s3 + tH) * sH 
480                + (tH - 3. * s3 - 3. * s4) * tH + (2. * s3 - 2. * s1 
481                + 3. * s4) * s3 ) / ( (uH - s3) * (tH - s4) );
482   double smm13 = -4. * ( (tH + s1 - 2. * s4) * tH - (s3 + 3. * s1 - 2. * s4)
483                * s3 + (s3 + 3. * s1 + tH) * sH - pow2(tH - s3 + sH) )
484                / ( (uH - s3) * (sH - s1) );
485   double smm23 = -4. * ( (s1 - s4 + s3) * tH - s3*s3 + s3 * (s1 + s4)
486                - 3. * s1 * s4 - (s1 - s4 - s3 + tH) * sH)
487                / ( (sH - s1) * (tH - s4) );
488   double sigma = alpEM * pow2(sH / (sH - s1) ) * (smm1 + smm2 + smm3
489                + smm12 + smm13 + smm23) / (4. * sH2);
490
491   // Lepton Yukawa and secondary widths.
492   sigma       *= pow2(yukawa[(idInAbs-9)/2]);
493   sigma       *= (idIn < 0) ? openFracPos : openFracNeg; 
494
495   // Answer.
496   return sigma;    
497
498 }
499
500 //*********
501
502 // Select identity, colour and anticolour.
503
504 void Sigma2lgm2Hchgchgl::setIdColAcol() {
505
506   // Sign of outgoing H_L/R.
507   int idIn     = (id2 == 22) ? id1 : id2;
508   int idSgn    = (idIn < 0) ? idHLR : -idHLR;
509   int idOut    = (idIn < 0) ? idLep : -idLep;
510   setId( id1, id2, idSgn, idOut);
511
512   // tHat is defined between incoming lepton and outgoing Higgs.
513   if (id1 == 22) swapTU = true;
514
515   // No colours whatsoever.
516   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
517
518 }
519
520 //*********
521
522 // Evaluate weight for H_L/R sequential decay angles.
523   
524 double Sigma2lgm2Hchgchgl::weightDecay( Event& process, int iResBeg, 
525   int iResEnd) {
526
527   // Identity of mother of decaying reseonance(s).
528   int idMother = process[process[iResBeg].mother1()].idAbs();
529
530   // For top decay hand over to standard routine.
531   if (idMother == 6) 
532     return weightTopDecay( process, iResBeg, iResEnd);
533  
534   // Else isotropic decay.
535   return 1.;
536
537 }
538
539 //**************************************************************************
540
541 // Sigma3ff2HchgchgfftWW class.
542 // Cross section for f_1 f_2 -> H0 f_3 f_4 (W+ W- fusion of SM Higgs).
543
544 //*********
545
546 // Initialize process. 
547   
548 void Sigma3ff2HchgchgfftWW::initProc() {
549
550   // Set process properties: H_L^++-- or H_R^++--.
551   if (leftRight == 1) {
552     idHLR      = 9900041;
553     codeSave   = 3125;
554     nameSave   = "f_1 f_2 -> H_L^++-- f_3 f_4 (W+- W+- fusion)";
555   } else {
556     idHLR      = 9900042;
557     codeSave   = 3145;
558     nameSave   = "f_1 f_2 -> H_R^++-- f_3 f_4 (W+- W+- fusion)";
559   }
560
561   // Common fixed mass and coupling factor.
562   double mW    = ParticleDataTable::m0(24); 
563   double mWR   = ParticleDataTable::m0(9900024);
564   mWS          = (leftRight == 1) ? pow2(mW) : pow2(mWR);
565   double gL    = Settings::parm("LeftRightSymmmetry:gL");
566   double gR    = Settings::parm("LeftRightSymmmetry:gR");
567   double vL    = Settings::parm("LeftRightSymmmetry:vL"); 
568   prefac       = (leftRight == 1) ? pow2(pow4(gL) * vL) 
569                                   : 2. * pow2(pow3(gR) * mWR); 
570   // Secondary open width fractions.
571   openFracPos  = ParticleDataTable::resOpenFrac( idHLR);
572   openFracNeg  = ParticleDataTable::resOpenFrac(-idHLR);
573
574
575
576 //*********
577
578 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
579
580 void Sigma3ff2HchgchgfftWW::sigmaKin() { 
581
582   // Required four-vector products.
583   double pp12  = 0.5 * sH;
584   double pp14  = 0.5 * mH * p4cm.pMinus();
585   double pp15  = 0.5 * mH * p5cm.pMinus();
586   double pp24  = 0.5 * mH * p5cm.pPlus();
587   double pp25  = 0.5 * mH * p5cm.pPlus();
588   double pp45  = p4cm * p5cm;
589
590   // Cross section: kinematics part. Combine with couplings.
591   double propT = 1. / ( (2. * pp14 + mWS) * (2. * pp25 + mWS) );  
592   double propU = 1. / ( (2. * pp24 + mWS) * (2. * pp15 + mWS) );
593   sigma0TU     = prefac * pp12 * pp45 * pow2(propT + propU); 
594   sigma0T      = prefac * pp12 * pp45 * 2. * pow2(propT); 
595
596 }
597
598 //*********
599
600 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
601
602 double Sigma3ff2HchgchgfftWW::sigmaHat() { 
603
604   // Do not allow creation of righthanded neutrinos for H_R.
605   int id1Abs   = abs(id1);
606   int id2Abs   = abs(id2);
607   if ( leftRight == 2 && (id1Abs > 10 || id2Abs > 10) ) return 0.; 
608
609   // Many flavour combinations not possible because of charge.
610   int chg1     = (( id1Abs%2 == 0 && id1 > 0) 
611                || (id1Abs%2 == 1 && id1 < 0) ) ? 1 : -1;
612   int chg2     = (( id2Abs%2 == 0 && id2 > 0) 
613                || (id2Abs%2 == 1 && id2 < 0) ) ? 1 : -1;
614   if (abs(chg1 + chg2) != 2) return 0.;
615
616   // Basic cross section. CKM factors for final states.
617   double sigma = (id2 == id1 && id1Abs > 10) ? sigma0TU : sigma0T;
618   sigma       *= VCKM::V2sum(id1Abs) * VCKM::V2sum(id2Abs);
619
620   // Secondary width for H0.
621   sigma       *= (chg1 + chg2 == 2) ? openFracPos : openFracNeg;
622
623   // Spin-state extra factor 2 per incoming neutrino.
624   if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.; 
625   if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.; 
626
627   // Answer.
628   return sigma;
629  
630 }
631
632 //*********
633
634 // Select identity, colour and anticolour.
635
636 void Sigma3ff2HchgchgfftWW::setIdColAcol() {
637
638   // Pick out-flavours by relative CKM weights.
639   int id1Abs   = abs(id1);
640   int id2Abs   = abs(id2);
641   id4          = VCKM::V2pick(id1);
642   id5          = VCKM::V2pick(id2);
643
644   // Find charge of Higgs .
645   id3 = (( id1Abs%2 == 0 && id1 > 0) || (id1Abs%2 == 1 && id1 < 0) ) 
646       ? idHLR : -idHLR;
647   setId( id1, id2, id3, id4, id5);
648
649   // Colour flow topologies. Swap when antiquarks.
650   if (id1Abs < 9 && id2Abs < 9 && id1*id2 > 0) 
651                        setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
652   else if (id1Abs < 9 && id2Abs < 9)
653                        setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2); 
654   else if (id1Abs < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0); 
655   else if (id2Abs < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0); 
656   else                 setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
657   if ( (id1Abs < 9 && id1 < 0) || (id1Abs > 10 && id2 < 0) ) 
658     swapColAcol();
659
660 }
661
662 //*********
663
664 // Evaluate weight for decay angles.
665
666 double Sigma3ff2HchgchgfftWW::weightDecay( Event& process, int iResBeg,
667   int iResEnd) {
668
669   // Identity of mother of decaying reseonance(s).
670   int idMother = process[process[iResBeg].mother1()].idAbs();
671
672   // For top decay hand over to standard routine.
673   if (idMother == 6) 
674     return weightTopDecay( process, iResBeg, iResEnd);
675
676   // Else done.
677   return 1.; 
678
679 }
680
681 //**************************************************************************
682
683 // Sigma2ffbar2HchgchgHchgchg class.
684 // Cross section for f fbar -> H_(L/R)^++ H_(L/R)^-- (doubly charged Higgs).
685
686 //*********
687
688 // Initialize process. 
689   
690 void Sigma2ffbar2HchgchgHchgchg::initProc() {
691
692   // Set process properties: H_L^++ H_L^-- or H_R^++ H_R^--.
693   if (leftRight == 1) {
694     idHLR      = 9900041;
695     codeSave   = 3126;
696     nameSave   = "f fbar -> H_L^++ H_L^--";
697   } else {
698     idHLR      = 9900042;
699     codeSave   = 3146;
700     nameSave   = "f fbar -> H_R^++ H_R^--";
701   }
702
703   // Read in Yukawa matrix for couplings to a lepton pair.
704   yukawa[1][1] = Settings::parm("LeftRightSymmmetry:coupHee");
705   yukawa[2][1] = Settings::parm("LeftRightSymmmetry:coupHmue");
706   yukawa[2][2] = Settings::parm("LeftRightSymmmetry:coupHmumu");
707   yukawa[3][1] = Settings::parm("LeftRightSymmmetry:coupHtaue");
708   yukawa[3][2] = Settings::parm("LeftRightSymmmetry:coupHtaumu");
709   yukawa[3][3] = Settings::parm("LeftRightSymmmetry:coupHtautau");
710
711   // Electroweak parameters.
712   mRes         = ParticleDataTable::m0(23);
713   GammaRes     = ParticleDataTable::mWidth(23);
714   m2Res        = mRes*mRes;
715   GamMRat      = GammaRes / mRes;
716   sin2tW       = CoupEW::sin2thetaW();
717   preFac       = (1. - 2. * sin2tW) / ( 8. * sin2tW * (1. - sin2tW) );
718
719   // Open fraction from secondary widths.
720   openFrac     = ParticleDataTable::resOpenFrac( idHLR, -idHLR);
721
722
723
724 //*********
725
726 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
727
728 double Sigma2ffbar2HchgchgHchgchg::sigmaHat() {
729
730   // Electroweak couplings to gamma^*/Z^0.
731   int    idAbs   = abs(id1);
732   double ei      = CoupEW::ef(idAbs);
733   double vi      = CoupEW::vf(idAbs); 
734   double ai      = CoupEW::af(idAbs); 
735
736   // Part via gamma^*/Z^0 propagator.No Z^0 coupling to H_R.
737   double resProp = 1. / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
738   double sigma   = 8. * pow2(alpEM) * ei*ei / sH2;
739   if (leftRight == 1) sigma += 8. * pow2(alpEM) 
740     * (2. * ei * vi * preFac * (sH - m2Res) * resProp / sH
741     + (vi * vi + ai * ai) * pow2(preFac) * resProp);      
742
743   // Part via t-channel lepton + interference; sum over possibilities. 
744   if (idAbs == 11 || idAbs == 13 || idAbs == 15) {
745     double yuk2Sum;
746     if (idAbs == 11) yuk2Sum 
747       = pow2(yukawa[1][1]) + pow2(yukawa[2][1]) + pow2(yukawa[3][1]);
748     else if (idAbs == 13) yuk2Sum 
749       = pow2(yukawa[2][1]) + pow2(yukawa[2][2]) + pow2(yukawa[3][2]);
750     else yuk2Sum 
751       = pow2(yukawa[3][1]) + pow2(yukawa[3][2]) + pow2(yukawa[3][3]);
752     yuk2Sum /= 4. * M_PI;
753     sigma   += 8. * alpEM * ei * yuk2Sum / (sH * tH) 
754       + 4. * pow2(yuk2Sum) / tH2;
755     if (leftRight == 1) sigma += 8. * alpEM * (vi + ai) * yuk2Sum 
756       * preFac * (sH - m2Res) * resProp / tH;   
757   }  
758
759   // Common kinematical factor. Colour factor.
760   sigma *= M_PI * (tH * uH - s3 * s4) / sH2;
761   if (idAbs < 9) sigma /= 3.;  
762
763   // Answer.
764   return sigma;    
765
766 }
767
768 //*********
769
770 // Select identity, colour and anticolour.
771
772 void Sigma2ffbar2HchgchgHchgchg::setIdColAcol() {
773
774   // Outgoing flavours trivial.
775   setId( id1, id2, idHLR, -idHLR);
776
777   // tHat is defined between incoming fermion and outgoing H--.
778   if (id1 > 0) swapTU = true;
779
780   // No colours at all or one flow topology. Swap if first is antiquark.
781   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
782   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
783   if (id1 < 0) swapColAcol();
784
785 }
786
787 //*********
788
789 // Evaluate weight for H_L/R sequential decay angles.
790   
791 double Sigma2ffbar2HchgchgHchgchg::weightDecay( Event& process, 
792   int iResBeg, int iResEnd) {
793
794   // Identity of mother of decaying reseonance(s).
795   int idMother = process[process[iResBeg].mother1()].idAbs();
796
797   // For top decay hand over to standard routine.
798   if (idMother == 6) 
799     return weightTopDecay( process, iResBeg, iResEnd);
800  
801   // Else isotropic decay.
802   return 1.;
803
804 }
805
806 //**************************************************************************
807
808 } // end namespace Pythia8