Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / SigmaLeftRightSym.cxx
1 // SigmaLeftRightSym.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 // 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     = particleDataPtr->m0(idZR);
27   GammaRes = particleDataPtr->mWidth(idZR);
28   m2Res    = mRes*mRes;
29   GamMRat  = GammaRes / mRes;
30   sin2tW   = couplingsPtr->sin2thetaW();
31
32   // Set pointer to particle properties and decay table.
33   ZRPtr    = particleDataPtr->particleDataEntryPtr(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     = particleDataPtr->m0(idWR);
175   GammaRes = particleDataPtr->mWidth(idWR);
176   m2Res    = mRes*mRes;
177   GamMRat  = GammaRes / mRes;
178   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
179
180   // Set pointer to particle properties and decay table.
181   particlePtr = particleDataPtr->particleDataEntryPtr(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->sizeChannels(); ++i) {
202     id1Now      = particlePtr->channel(i).product(0);
203     id2Now      = particlePtr->channel(i).product(1);
204     id1Abs      = abs(id1Now);
205     id2Abs      = abs(id2Now);
206
207     // Check that above threshold. Phase space.
208     mf1 = particleDataPtr->m0(id1Abs);
209     mf2 = particleDataPtr->m0(id2Abs);
210     if (mH > mf1 + mf2 + MASSMARGIN) {
211       mr1    = pow2(mf1 / mH);
212       mr2    = pow2(mf2 / mH);
213       kinFac = (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
214              * sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 ); 
215
216       // Combine kinematics with colour factor and CKM couplings.
217       widNow = kinFac;
218       if (id1Abs < 9) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
219  
220       // Secondary width from top and righthanded neutrino decay.
221       id1Neg    = (id1Abs < 19) ? -id1Now : id1Abs; 
222       id2Neg    = (id2Abs < 19) ? -id2Now : id2Abs; 
223       widSecPos = particleDataPtr->resOpenFrac(id1Now, id2Now); 
224       widSecNeg = particleDataPtr->resOpenFrac(id1Neg, id2Neg); 
225
226       // Add weight for channels on for all, W_R^+ and W_R^-, respectively.
227       onMode = particlePtr->channel(i).onMode();
228       if (onMode == 1 || onMode == 2) widOutPos += widNow * widSecPos;
229       if (onMode == 1 || onMode == 3) widOutNeg += widNow * widSecNeg;
230
231     // End loop over fermions.
232     }
233   }
234
235   // Set up Breit-Wigner. Cross section for W_R^+ and W_R^- separately.
236   double sigBW = 12. * M_PI * pow2(alpEM * thetaWRat) * sH
237                / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
238   sigma0Pos = sigBW * widOutPos;    
239   sigma0Neg = sigBW * widOutNeg;    
240
241 }
242
243 //--------------------------------------------------------------------------
244
245 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
246
247 double Sigma1ffbar2WRight::sigmaHat() {
248
249   // Secondary width for W_R^+ or W_R^-. CKM and colour factors.
250   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
251   double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
252   if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
253
254   // Answer.
255   return sigma;    
256
257 }
258
259 //--------------------------------------------------------------------------
260
261 // Select identity, colour and anticolour.
262
263 void Sigma1ffbar2WRight::setIdColAcol() {
264
265   // Sign of outgoing W_R.
266   int sign          = (abs(id1)%2 == 0) ? 1 : -1;
267   if (id1 < 0) sign = -sign;
268   setId( id1, id2, idWR * sign);
269
270   // Colour flow topologies. Swap when antiquarks.
271   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
272   else              setColAcol( 0, 0, 0, 0, 0, 0);
273   if (id1 < 0) swapColAcol();
274
275 }
276
277 //--------------------------------------------------------------------------
278
279 // Evaluate weight for W_R decay angle.
280   
281 double Sigma1ffbar2WRight::weightDecay( Event& process, int iResBeg, 
282   int iResEnd) {
283
284   // Identity of mother of decaying reseonance(s).
285   int idMother = process[process[iResBeg].mother1()].idAbs();
286
287   // For top decay hand over to standard routine.
288   if (idMother == 6) 
289     return weightTopDecay( process, iResBeg, iResEnd);
290
291   // W_R should sit in entry 5.
292   if (iResBeg != 5 || iResEnd != 5) return 1.;
293
294   // Phase space factors.
295   double mr1    = pow2(process[6].m()) / sH;
296   double mr2    = pow2(process[7].m()) / sH;
297   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
298    
299   // Sign of asymmetry.
300   double eps    = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
301
302   // Reconstruct decay angle and weight for it.
303   double cosThe = (process[3].p() - process[4].p()) 
304     * (process[7].p() - process[6].p()) / (sH * betaf);
305   double wtMax  = 4.;
306   double wt     = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2); 
307  
308   // Done.
309   return (wt / wtMax);
310
311 }
312
313 //==========================================================================
314
315 // Sigma1ll2Hchgchg class.
316 // Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
317
318 //--------------------------------------------------------------------------
319
320 // Initialize process. 
321   
322 void Sigma1ll2Hchgchg::initProc() {
323
324   // Set process properties: H_L^++-- or H_R^++--.
325   if (leftRight == 1) {
326     idHLR    = 9900041;
327     codeSave = 3121;
328     nameSave = "l l -> H_L^++--";
329   } else {
330     idHLR    = 9900042;
331     codeSave = 3141;
332     nameSave = "l l -> H_R^++--";
333   }
334
335   // Read in Yukawa matrix for couplings to a lepton pair.
336   yukawa[1][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHee");
337   yukawa[2][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
338   yukawa[2][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
339   yukawa[3][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
340   yukawa[3][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
341   yukawa[3][3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
342
343   // Store H_L/R mass and width for propagator. 
344   mRes     = particleDataPtr->m0(idHLR);
345   GammaRes = particleDataPtr->mWidth(idHLR);
346   m2Res    = mRes*mRes;
347   GamMRat  = GammaRes / mRes;
348
349   // Set pointer to particle properties and decay table.
350   particlePtr = particleDataPtr->particleDataEntryPtr(idHLR);
351
352
353
354 //--------------------------------------------------------------------------
355
356 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
357
358 double Sigma1ll2Hchgchg::sigmaHat() {
359
360   // Initial state must consist of two identical-sign leptons.
361   if (id1 * id2 < 0) return 0.;
362   int id1Abs = abs(id1);
363   int id2Abs = abs(id2);  
364   if (id1Abs != 11 && id1Abs != 13 && id1Abs != 15) return 0.;
365   if (id2Abs != 11 && id2Abs != 13 && id2Abs != 15) return 0.;
366
367   // Set up Breit-Wigner, inwidth and outwidth.
368   double sigBW  = 8. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
369   double widIn  = pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) 
370                 * mH / (8. * M_PI);
371   int idSgn     = (id1 < 0) ? idHLR : -idHLR;
372   double widOut = particlePtr->resWidthOpen( idSgn, mH);
373
374   // Answer.
375   return widIn * sigBW * widOut;    
376
377 }
378
379 //--------------------------------------------------------------------------
380
381 // Select identity, colour and anticolour.
382
383 void Sigma1ll2Hchgchg::setIdColAcol() {
384
385   // Sign of outgoing H_L/R.
386   int idSgn     = (id1 < 0) ? idHLR : -idHLR;
387   setId( id1, id2, idSgn);
388
389   // No colours whatsoever.
390   setColAcol( 0, 0, 0, 0, 0, 0);
391
392 }
393
394 //--------------------------------------------------------------------------
395
396 // Evaluate weight for H_L/R sequential decay angles.
397   
398 double Sigma1ll2Hchgchg::weightDecay( Event& process, int iResBeg, 
399   int iResEnd) {
400
401   // Identity of mother of decaying reseonance(s).
402   int idMother = process[process[iResBeg].mother1()].idAbs();
403
404   // For top decay hand over to standard routine.
405   if (idMother == 6) 
406     return weightTopDecay( process, iResBeg, iResEnd);
407  
408   // Else isotropic decay.
409   return 1.;
410
411 }
412
413 //==========================================================================
414
415 // Sigma2lgm2Hchgchgl class.
416 // Cross section for l gamma -> H_L^++-- l or H_R^++-- l 
417 // (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]  = settingsPtr->parm("LeftRightSymmmetry:coupHee");
440     yukawa[2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
441     yukawa[3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
442   } else if (idLep == 13) { 
443     yukawa[1]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
444     yukawa[2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
445     yukawa[3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
446   } else {
447     yukawa[1]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
448     yukawa[2]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
449     yukawa[3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
450   }
451
452   // Secondary open width fractions.
453   openFracPos  = particleDataPtr->resOpenFrac( idHLR);
454   openFracNeg  = particleDataPtr->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( particleDataPtr->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 -> H_(L/R)^++-- f_3 f_4 (W+- W+- fusion).
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    = particleDataPtr->m0(24); 
563   double mWR   = particleDataPtr->m0(9900024);
564   mWS          = (leftRight == 1) ? pow2(mW) : pow2(mWR);
565   double gL    = settingsPtr->parm("LeftRightSymmmetry:gL");
566   double gR    = settingsPtr->parm("LeftRightSymmmetry:gR");
567   double vL    = settingsPtr->parm("LeftRightSymmmetry:vL"); 
568   prefac       = (leftRight == 1) ? pow2(pow4(gL) * vL) 
569                                   : 2. * pow2(pow3(gR) * mWR); 
570   // Secondary open width fractions.
571   openFracPos  = particleDataPtr->resOpenFrac( idHLR);
572   openFracNeg  = particleDataPtr->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.pNeg();
585   double pp15  = 0.5 * mH * p5cm.pNeg();
586   double pp24  = 0.5 * mH * p4cm.pPos();
587   double pp25  = 0.5 * mH * p5cm.pPos();
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       *= couplingsPtr->V2CKMsum(id1Abs) 
619                * couplingsPtr->V2CKMsum(id2Abs);
620
621   // Secondary width for H0.
622   sigma       *= (chg1 + chg2 == 2) ? openFracPos : openFracNeg;
623
624   // Spin-state extra factor 2 per incoming neutrino.
625   if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.; 
626   if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.; 
627
628   // Answer.
629   return sigma;
630  
631 }
632
633 //--------------------------------------------------------------------------
634
635 // Select identity, colour and anticolour.
636
637 void Sigma3ff2HchgchgfftWW::setIdColAcol() {
638
639   // Pick out-flavours by relative CKM weights.
640   int id1Abs   = abs(id1);
641   int id2Abs   = abs(id2);
642   id4          = couplingsPtr->V2CKMpick(id1);
643   id5          = couplingsPtr->V2CKMpick(id2);
644
645   // Find charge of Higgs .
646   id3 = (( id1Abs%2 == 0 && id1 > 0) || (id1Abs%2 == 1 && id1 < 0) ) 
647       ? idHLR : -idHLR;
648   setId( id1, id2, id3, id4, id5);
649
650   // Colour flow topologies. Swap when antiquarks.
651   if (id1Abs < 9 && id2Abs < 9 && id1*id2 > 0) 
652                        setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
653   else if (id1Abs < 9 && id2Abs < 9)
654                        setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2); 
655   else if (id1Abs < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0); 
656   else if (id2Abs < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0); 
657   else                 setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
658   if ( (id1Abs < 9 && id1 < 0) || (id1Abs > 10 && id2 < 0) ) 
659     swapColAcol();
660
661 }
662
663 //--------------------------------------------------------------------------
664
665 // Evaluate weight for decay angles.
666
667 double Sigma3ff2HchgchgfftWW::weightDecay( Event& process, int iResBeg,
668   int iResEnd) {
669
670   // Identity of mother of decaying reseonance(s).
671   int idMother = process[process[iResBeg].mother1()].idAbs();
672
673   // For top decay hand over to standard routine.
674   if (idMother == 6) 
675     return weightTopDecay( process, iResBeg, iResEnd);
676
677   // Else done.
678   return 1.; 
679
680 }
681
682 //==========================================================================
683
684 // Sigma2ffbar2HchgchgHchgchg class.
685 // Cross section for f fbar -> H_(L/R)^++ H_(L/R)^-- (doubly charged Higgs).
686
687 //--------------------------------------------------------------------------
688
689 // Initialize process. 
690   
691 void Sigma2ffbar2HchgchgHchgchg::initProc() {
692
693   // Set process properties: H_L^++ H_L^-- or H_R^++ H_R^--.
694   if (leftRight == 1) {
695     idHLR      = 9900041;
696     codeSave   = 3126;
697     nameSave   = "f fbar -> H_L^++ H_L^--";
698   } else {
699     idHLR      = 9900042;
700     codeSave   = 3146;
701     nameSave   = "f fbar -> H_R^++ H_R^--";
702   }
703
704   // Read in Yukawa matrix for couplings to a lepton pair.
705   yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
706   yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
707   yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
708   yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
709   yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
710   yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
711
712   // Electroweak parameters.
713   mRes         = particleDataPtr->m0(23);
714   GammaRes     = particleDataPtr->mWidth(23);
715   m2Res        = mRes*mRes;
716   GamMRat      = GammaRes / mRes;
717   sin2tW       = couplingsPtr->sin2thetaW();
718   preFac       = (1. - 2. * sin2tW) / ( 8. * sin2tW * (1. - sin2tW) );
719
720   // Open fraction from secondary widths.
721   openFrac     = particleDataPtr->resOpenFrac( idHLR, -idHLR);
722
723
724
725 //--------------------------------------------------------------------------
726
727 // Evaluate sigmaHat(sHat), including incoming flavour dependence. 
728
729 double Sigma2ffbar2HchgchgHchgchg::sigmaHat() {
730
731   // Electroweak couplings to gamma^*/Z^0.
732   int    idAbs   = abs(id1);
733   double ei      = couplingsPtr->ef(idAbs);
734   double vi      = couplingsPtr->vf(idAbs); 
735   double ai      = couplingsPtr->af(idAbs); 
736
737   // Part via gamma^*/Z^0 propagator. No Z^0 coupling to H_R.
738   double resProp = 1. / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
739   double sigma   = 8. * pow2(alpEM) * ei*ei / sH2;
740   if (leftRight == 1) sigma += 8. * pow2(alpEM) 
741     * (2. * ei * vi * preFac * (sH - m2Res) * resProp / sH
742     + (vi * vi + ai * ai) * pow2(preFac) * resProp);      
743
744   // Part via t-channel lepton + interference; sum over possibilities. 
745   if (idAbs == 11 || idAbs == 13 || idAbs == 15) {
746     double yuk2Sum;
747     if (idAbs == 11) yuk2Sum 
748       = pow2(yukawa[1][1]) + pow2(yukawa[2][1]) + pow2(yukawa[3][1]);
749     else if (idAbs == 13) yuk2Sum 
750       = pow2(yukawa[2][1]) + pow2(yukawa[2][2]) + pow2(yukawa[3][2]);
751     else yuk2Sum 
752       = pow2(yukawa[3][1]) + pow2(yukawa[3][2]) + pow2(yukawa[3][3]);
753     yuk2Sum /= 4. * M_PI;
754     sigma   += 8. * alpEM * ei * yuk2Sum / (sH * tH) 
755       + 4. * pow2(yuk2Sum) / tH2;
756     if (leftRight == 1) sigma += 8. * alpEM * (vi + ai) * yuk2Sum 
757       * preFac * (sH - m2Res) * resProp / tH;   
758   }  
759
760   // Common kinematical factor. Colour factor.
761   sigma *= M_PI * (tH * uH - s3 * s4) / sH2;
762   if (idAbs < 9) sigma /= 3.;  
763
764   // Answer.
765   return sigma;    
766
767 }
768
769 //--------------------------------------------------------------------------
770
771 // Select identity, colour and anticolour.
772
773 void Sigma2ffbar2HchgchgHchgchg::setIdColAcol() {
774
775   // Outgoing flavours trivial.
776   setId( id1, id2, idHLR, -idHLR);
777
778   // tHat is defined between incoming fermion and outgoing H--.
779   if (id1 > 0) swapTU = true;
780
781   // No colours at all or one flow topology. Swap if first is antiquark.
782   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
783   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
784   if (id1 < 0) swapColAcol();
785
786 }
787
788 //--------------------------------------------------------------------------
789
790 // Evaluate weight for H_L/R sequential decay angles.
791   
792 double Sigma2ffbar2HchgchgHchgchg::weightDecay( Event& process, 
793   int iResBeg, int iResEnd) {
794
795   // Identity of mother of decaying reseonance(s).
796   int idMother = process[process[iResBeg].mother1()].idAbs();
797
798   // For top decay hand over to standard routine.
799   if (idMother == 6) 
800     return weightTopDecay( process, iResBeg, iResEnd);
801  
802   // Else isotropic decay.
803   return 1.;
804
805 }
806
807 //==========================================================================
808
809 } // end namespace Pythia8