StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SusyResonanceWidths.cc
1 // SusyResonanceWidths.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 Torbjorn Sjostrand
3 // Authors: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
5 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 
7 // Function definitions (not found in the header) for the
8 // SUSY Resonance widths classes.
9 
10 #include "Pythia8/SusyResonanceWidths.h"
11 #include "Pythia8/SusyCouplings.h"
12 #include "Pythia8/ParticleData.h"
13 #include "Pythia8/PythiaComplex.h"
14 
15 namespace Pythia8 {
16 
17 //==========================================================================
18 
19 // WidthFunctions Class
20 
21 // Contains functions to be integrated for calculating the 3-body
22 // decay widths
23 
24 //--------------------------------------------------------------------------
25 
26 void WidthFunction::init( ParticleData* particleDataPtrIn,
27  CoupSUSY* coupSUSYPtrIn) {
28 
29  particleDataPtr = particleDataPtrIn;
30  coupSUSYPtr = coupSUSYPtrIn;
31 
32 }
33 
34 //--------------------------------------------------------------------------
35 
36 void WidthFunction::setInternal2(int idResIn, int id1In, int id2In,
37  int id3In, int idIntIn) {
38 
39  // Res -> 1,2,3
40  idRes = idResIn;
41  id1 = id1In;
42  id2 = id2In;
43  id3 = id3In;
44  idInt = idIntIn;
45 
46  mRes = particleDataPtr->m0(idRes);
47  m1 = particleDataPtr->m0(id1);
48  m2 = particleDataPtr->m0(id2);
49  m3 = particleDataPtr->m0(id3);
50 
51  // Internal propagator
52  mInt = particleDataPtr->m0(idInt);
53  gammaInt = particleDataPtr->mWidth(idInt);
54 
55  return;
56 }
57 
58 //--------------------------------------------------------------------------
59 
60 double WidthFunction::function(double) {
61 
62  cout<<"Warning using dummy width function"<<endl;
63  return 0.;
64 }
65 
66 //--------------------------------------------------------------------------
67 
68 double WidthFunction::function(double,double) {
69 
70  cout<<"Warning using dummy width function"<<endl;
71  return 0.;
72 }
73 
74 //==========================================================================
75 
76 // Psi, Upsilon and Phi classes.
77 
78 //--------------------------------------------------------------------------
79 
80 void Psi::setInternal (int idResIn, int id1In, int id2In, int id3In,
81  int idIntIn, int) {
82 
83  setInternal2(idResIn, id1In, id2In, id3In, idIntIn);
84 
85  mInt = particleDataPtr->m0(idInt);
86  gammaInt = particleDataPtr->mWidth(idInt);
87  iX = coupSUSYPtr->typeNeut(idRes);
88  iQ = (id3+1)/2;
89  iSq = (idInt>1000000)? 3 + (idInt%10+1)/2 : (idInt%10+1)/2;
90  isSqDown = (idInt % 2 == 1)? true : false;
91  m1 = particleDataPtr->m0(id1);
92  m2 = particleDataPtr->m0(id2);
93  m3 = particleDataPtr->m0(id3);
94  return;
95 }
96 
97 //--------------------------------------------------------------------------
98 
99 void Upsilon::setInternal (int idResIn, int id1In, int id2In, int id3In,
100  int idIntIn, int idInt2In) {
101 
102  setInternal2(idResIn, id1In, id2In, id3In, idIntIn);
103 
104  idInt2 = idInt2In;
105  mInt = particleDataPtr->m0(idInt);
106  gammaInt = particleDataPtr->mWidth(idInt);
107  mInt2 = particleDataPtr->m0(idInt2);
108  gammaInt2 = particleDataPtr->mWidth(idInt2);
109 
110  iX = coupSUSYPtr->typeNeut(idRes);
111  iQ = (id3+1)/2;
112  iSq = (idInt>1000000)? 3 + (idInt%10+1)/2 : (idInt%10+1)/2;
113  iSq2 = (idInt2>1000000)? 3 + (idInt2%10+1)/2 : (idInt2%10+1)/2;
114  isSqDown = (idIntIn % 2 == 1)? true : false;
115  m1 = particleDataPtr->m0(id1);
116  m2 = particleDataPtr->m0(id2);
117  m3 = particleDataPtr->m0(id3);
118  return;
119 }
120 
121 //--------------------------------------------------------------------------
122 
123 void Phi::setInternal (int idResIn, int id1In, int id2In, int id3In,
124  int idIntIn, int idInt2In) {
125 
126  setInternal2(idResIn, id1In, id2In, id3In, idIntIn);
127 
128  idInt2 = idInt2In;
129  mInt = particleDataPtr->m0(idInt);
130  gammaInt = particleDataPtr->mWidth(idInt);
131  mInt2 = particleDataPtr->m0(idInt2);
132  gammaInt2 = particleDataPtr->mWidth(idInt2);
133 
134  iX = coupSUSYPtr->typeNeut(idRes);
135  iQ = (id3+1)/2;
136  iSq = (idInt>1000000)? 3 + (idInt%10+1)/2 : (idInt%10+1)/2;
137  iSq2 = (idInt2>1000000)? 3 + (idInt2%10+1)/2 : (idInt2%10+1)/2;
138  isSqDown = (idIntIn % 2 == 1)? true : false;
139  m1 = particleDataPtr->m0(id1);
140  m2 = particleDataPtr->m0(id2);
141  m3 = particleDataPtr->m0(id3);
142  return;
143 }
144 
145 //--------------------------------------------------------------------------
146 
147 double Psi::function(double m12sq) {
148 
149  double R, factor1, factor2, value;
150 
151  // Check that the propagators are offshell
152  if (m12sq > pow2(mInt) || abs(m12sq - pow2(mInt)) < gammaInt) return 0;
153 
154  R = 1.0/(pow2(m12sq-pow2(mInt)) + pow2(gammaInt*mInt));
155  if (isSqDown){
156  factor1 = (norm(coupSUSYPtr->LsddX[iSq][iQ][iX])
157  + norm(coupSUSYPtr->RsddX[iSq][iQ][iX]))*
158  (mRes*mRes + m3*m3 - m12sq);
159  factor2 = 4.0 * real(coupSUSYPtr->LsddX[iSq][iQ][iX]
160  * conj(coupSUSYPtr->RsddX[iSq][iQ][iX])) * m3 * mRes;
161 
162  } else {
163  factor1 = (norm(coupSUSYPtr->LsuuX[iSq][iQ][iX])
164  + norm(coupSUSYPtr->RsuuX[iSq][iQ][iX]))*
165  (mRes*mRes + m3*m3 - m12sq);
166  factor2 = 4.0 * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
167  * conj(coupSUSYPtr->RsuuX[iSq][iQ][iX])) * m3 * mRes;
168  }
169 
170  value = R * (m12sq - m1*m1 - m2*m2) * (factor1+factor2);
171 
172  return value;
173 }
174 
175 
176 //--------------------------------------------------------------------------
177 
178 double Upsilon::function(double m12sq) {
179 
180  double R1,R2, S, factor1, factor2, value;
181 
182  // Check that the propagators are offshell
183  if (m12sq > pow2(mInt) || abs(m12sq - pow2(mInt)) < gammaInt) return 0;
184  if (m12sq > pow2(mInt2) || abs(m12sq - pow2(mInt2)) < gammaInt2) return 0;
185 
186  R1 = 1.0/(pow2(m12sq-pow2(mInt)) + pow2(gammaInt*mInt));
187  R2 = 1.0/(pow2(m12sq-pow2(mInt2)) + pow2(gammaInt2*mInt2));
188  S = R1 * R2 * ((m12sq - pow2(mInt)) * (m12sq - pow2(mInt2))
189  + gammaInt * mInt * gammaInt2 * mInt2);
190 
191  if (isSqDown){
192  factor1 = real(coupSUSYPtr->LsddX[iSq][iQ][iX]
193  * conj(coupSUSYPtr->LsddX[iSq2][iQ][iX]))
194  + real(coupSUSYPtr->RsddX[iSq][iQ][iX]
195  * conj(coupSUSYPtr->RsddX[iSq2][iQ][iX]));
196  factor2 = real(coupSUSYPtr->LsddX[iSq][iQ][iX]
197  * conj(coupSUSYPtr->RsddX[iSq2][iQ][iX]))
198  + real(coupSUSYPtr->RsddX[iSq][iQ][iX]
199  * conj(coupSUSYPtr->LsddX[iSq2][iQ][iX]));
200  }else{
201  factor1 = real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
202  * conj(coupSUSYPtr->LsuuX[iSq2][iQ][iX]))
203  + real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
204  * conj(coupSUSYPtr->RsuuX[iSq2][iQ][iX]));
205  factor2 = real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
206  * conj(coupSUSYPtr->RsuuX[iSq2][iQ][iX]))
207  + real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
208  * conj(coupSUSYPtr->LsuuX[iSq2][iQ][iX]));
209  }
210 
211  value = S * (m12sq - pow2(m1) - pow2(m2)) *
212  ( factor1 * (pow2(mRes) + pow2(m3) - m12sq) + 2.0 * factor2 * m3 * mRes);
213 
214  // cout<<"I1: "<<idInt<<" I2:"<<idInt2<<" factor1: "<<factor1
215  // <<" factor2:"<<factor2<<" value:"<<value<<endl;
216 
217  return value;
218 
219 }
220 
221 //--------------------------------------------------------------------------
222 
223 double Phi::function(double m12sqIn) {
224 
225  m12sq = m12sqIn;
226  // Check that the propagators are offshell
227  if (m12sq > pow2(mInt) || abs(m12sq - pow2(mInt)) < gammaInt) return 0;
228 
229  double m23max, m23min, E2, E3;
230 
231  E2 = (m12sq - pow2(m1) - pow2(m2))/(2.0 * sqrt(m12sq));
232  E3 = (pow2(mRes) - m12sq - pow2(m3))/(2.0 * sqrt(m12sq));
233  m23max = pow2(E2+E3) - (sqrt(E2*E2 - m2*m2) - sqrt(E3*E3 - m3*m3)) ;
234  m23min = pow2(E2+E3) - (sqrt(E2*E2 - m2*m2) + sqrt(E3*E3 - m3*m3)) ;
235 
236  if (E2 < m2 || E3 < m3){
237  cout <<"Error in Phi:function"<<endl;
238  return 0.;
239  }
240 
241  double integral2 = integrateGauss(m23min,m23max,1.0e-4);
242  return integral2;
243 }
244 
245 //--------------------------------------------------------------------------
246 
247 double Phi::function2(double m23sq) {
248 
249  // Check that the propagators are offshell
250  if (m23sq > pow2(mInt2) || abs(m23sq - pow2(mInt2)) < gammaInt2) return 0;
251 
252  double R1, R2, S, factor1, factor2, factor3, factor4, value, fac;
253  int iQ2;
254 
255  R1 = 1.0/(pow2(m12sq-pow2(mInt)) + pow2(gammaInt*mInt));
256  R2 = 1.0/(pow2(m12sq-pow2(mInt2)) + pow2(gammaInt2*mInt2));
257  S = R1 * R2 * ((m12sq - pow2(mInt)) * (m12sq - pow2(mInt2))
258  + gammaInt * mInt * gammaInt2 * mInt2);
259 
260  fac = 1.0;
261 
262  if (isSqDown){
263  // Only factor is when both d_i and d_j are near.
264  iQ2 = (id1%2 == 1)? (id1+1)/2 : (id2+1)/2;
265 
266  if (mRes > mInt2 + particleDataPtr->m0(iQ2)) fac = 0.;
267 
268  factor1 = m1 * m3 * real(coupSUSYPtr->LsddX[iSq][iQ][iX]
269  * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
270  * (m12sq + m23sq - pow2(m1) - pow2(m3));
271  factor2 = m1 * mRes * real(coupSUSYPtr->RsddX[iSq][iQ][iX]
272  * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
273  * (m23sq - pow2(m2) - pow2(m3));
274  factor3 = m3 * mRes * real(coupSUSYPtr->LsddX[iSq][iQ][iX]
275  * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
276  * (m12sq - pow2(m1) - pow2(m2));
277  factor4 = m3 * mRes * real(coupSUSYPtr->RsddX[iSq][iQ][iX]
278  * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
279  * (m12sq * m23sq - pow2(m1 * m3) - pow2(m2 * mRes));
280 
281  }else{
282  // Factor A: u and d_1
283  iQ2 = (id1+1)/2;
284 
285  if (mRes > mInt2 + particleDataPtr->m0(iQ2)) fac = 0.;
286 
287  factor1 = m1 * m3 * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
288  * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
289  * (m12sq + m23sq - pow2(m1) - pow2(m3));
290  factor2 = m1 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
291  * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
292  * (m23sq - pow2(m2) - pow2(m3));
293  factor3 = m3 * mRes * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
294  * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
295  * (m12sq - pow2(m1) - pow2(m2));
296  factor4 = m3 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
297  * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
298  * (m12sq * m23sq - pow2(m1 * m3) - pow2(m2 * mRes));
299 
300 
301  // Factor B: u and d_2; change 1 <=> 2
302  iQ2 = (id2+1)/2;
303 
304  if (mRes > mInt2 + particleDataPtr->m0(iQ2)) fac = 0.;
305 
306  factor1 += m2 * m3 * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
307  * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
308  * (m12sq + m23sq - pow2(m2) - pow2(m3));
309  factor2 += m2 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
310  * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
311  * (m23sq - pow2(m1) - pow2(m3));
312  factor3 += m3 * mRes * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
313  * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
314  * (m12sq - pow2(m2) - pow2(m1));
315  factor4 += m3 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
316  * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
317  * (m12sq * m23sq - pow2(m2 * m3) - pow2(m1 * mRes));
318  }
319 
320  value = S * (factor1 + factor2 + factor3 + factor4);
321 
322  // cout<<"I1: "<<idInt<<" I2:"<<idInt2<<" factor1: "<<factor1
323  // <<" factor2:"<<factor2<<" value:"<<value<<endl;
324 
325  return (fac * value);
326 }
327 
328 //--------------------------------------------------------------------------
329 
330 double Phi::integrateGauss(double xlo, double xhi, double tol) {
331 
332  //8-point unweighted
333  static double x8[4]={0.96028985649753623,
334  0.79666647741362674,
335  0.52553240991632899,
336  0.18343464249564980};
337  static double w8[4]={0.10122853629037626,
338  0.22238103445337447,
339  0.31370664587788729,
340  0.36268378337836198};
341  //16-point unweighted
342  static double x16[8]={0.98940093499164993,
343  0.94457502307323258,
344  0.86563120238783174,
345  0.75540440835500303,
346  0.61787624440264375,
347  0.45801677765722739,
348  0.28160355077925891,
349  0.09501250983763744};
350  static double w16[8]={0.027152459411754095,
351  0.062253523938647893,
352  0.095158511682492785,
353  0.12462897125553387,
354  0.14959598881657673,
355  0.16915651939500254,
356  0.18260341504492359,
357  0.18945061045506850};
358  //boundary checks
359  if (xlo == xhi) {
360  cerr<<"xlo = xhi"<<endl;
361  return 0.0;
362  }
363  if ( xlo > xhi ) {
364  cerr<<" (integrateGauss:) -> xhi < xlo"<<endl;
365  return 0.0;
366  }
367  //initialize
368  double sum = 0.0;
369  double c = 0.001/abs(xhi-xlo);
370  double zlo = xlo;
371  double zhi = xhi;
372 
373  bool nextbin = true;
374 
375  while ( nextbin ) {
376 
377  double zmi = 0.5*(zhi+zlo); // midpoint
378  double zmr = 0.5*(zhi-zlo); // midpoint, relative to zlo
379 
380  //calculate 8-point and 16-point quadratures
381  double s8=0.0;
382  for (int i=0;i<4;i++) {
383  double dz = zmr * x8[i];
384  s8 += w8[i]*(function2(zmi+dz) + function2(zmi-dz));
385  }
386  s8 *= zmr;
387  double s16=0.0;
388  for (int i=0;i<8;i++) {
389  double dz = zmr * x16[i];
390  s16 += w16[i]*(function2(zmi+dz) + function2(zmi-dz));
391  }
392  s16 *= zmr;
393  if (abs(s16-s8) < tol*(1+abs(s16))) {
394  //precision in this bin OK, add to cumulative and go to next
395  nextbin=true;
396  sum += s16;
397  //next bin: LO = end of current, HI = end of integration region.
398  zlo=zhi;
399  zhi=xhi;
400  if ( zlo == zhi ) nextbin = false;
401  } else {
402  //precision in this bin not OK, subdivide.
403  if (1.0 + c*abs(zmr) == 1.0) {
404  cerr << " (integrateGauss:) too high accuracy required"<<endl;
405  sum = 0.0 ;
406  break;
407  }
408  zhi=zmi;
409  nextbin=true;
410  }
411  }
412  return sum;
413 }
414 
415 //==========================================================================
416 
417 // The SUSYResonanceWidths Class
418 // Derived class for SUSY resonances
419 
420 const bool SUSYResonanceWidths::DBSUSY = false;
421 
422 //--------------------------------------------------------------------------
423 
424 bool SUSYResonanceWidths::initBSM(){
425 
426  if (couplingsPtr->isSUSY) {
427  coupSUSYPtr = (CoupSUSY *) couplingsPtr;
428  return true;
429  }
430 
431  return false;
432 }
433 
434 
435 //--------------------------------------------------------------------------
436 
437 bool SUSYResonanceWidths::allowCalc(){
438 
439  // Check if decay calculations at all possible
440  if ( !couplingsPtr->isSUSY ) return false;
441 
442  // Next check if SLHA decay tables are ignored (= always do calculation)
443  if ( !settingsPtr->flag("SLHA:useDecayTable") ) return true;
444 
445  // Next check if decay table was read in via SLHA and takes precedence
446  for ( int iDec = 1; iDec < int((coupSUSYPtr->slhaPtr)->decays.size());
447  ++iDec)
448  if ( (coupSUSYPtr->slhaPtr)->decays[iDec].getId() == abs(idRes) ) {
449  if (DBSUSY) cout<<"Using external decay table for:"<<idRes<<endl;
450  return false;
451  }
452 
453  // Else we should do the calculation
454  return true;
455 }
456 
457 //==========================================================================
458 
459 // The ResonanceSquark class
460 // Derived class for Squark resonances
461 
462 //--------------------------------------------------------------------------
463 
464 // Initialize constants.
465 
466 void ResonanceSquark::initConstants() {
467 
468  // Locally stored properties and couplings.
469  s2W = coupSUSYPtr->sin2W;
470 
471 }
472 
473 //--------------------------------------------------------------------------
474 
475 // Calculate various common prefactors for the current mass.
476 
477 void ResonanceSquark::calcPreFac(bool) {
478 
479  // Common coupling factors.
480  alpS = coupSUSYPtr->alphaS(mHat * mHat );
481  alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
482  preFac = 1.0 / (s2W * pow(mHat,3));
483  ps *= mHat * mHat;
484 
485 }
486 
487 //--------------------------------------------------------------------------
488 
489 // Calculate width for currently considered channel.
490 
491 void ResonanceSquark::calcWidth(bool) {
492 
493  // Squark type -- in u_i/d_i and generation
494  int ksusy = 1000000;
495  bool idown = (abs(idRes)%2 == 0 ? false : true);
496  int isq = (abs(idRes)/ksusy == 2) ?
497  (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
498  // int isqgen = (abs(idRes)%10 + 1)/2;
499 
500  // Check that mass is above threshold.
501  if (ps == 0.) return;
502  else{
503  // Two-body decays
504  kinFac = (mHat * mHat - mf1 * mf1 - mf2 * mf2);
505 
506  double fac = 0.0 , wid = 0.0;
507 
508  //RPV decays
509  //Case 1a: UDD-type
510  if (id1Abs < 7 && id2Abs < 7){
511 
512  // Quark generations
513  int iq1 = (id1Abs + 1)/2;
514  int iq2 = (id2Abs + 1)/2;
515 
516  // Check for RPV UDD couplings
517  if (!coupSUSYPtr->isUDD) {widNow = 0; return;}
518 
519  // ~q -> q_i + q_j
520 
521  fac = 2.0 * kinFac / (16.0 * M_PI * pow(mHat,3));
522  wid = 0.0;
523  if (idown) {
524  if ((id1Abs+id2Abs)%2 == 1){
525  if (id1Abs%2==1)
526  for (int isq2 = 1; isq2 < 4; isq2++)
527  wid += norm(coupSUSYPtr->rvUDD[iq2][iq1][isq2]
528  * coupSUSYPtr->Rdsq[isq][isq2+3]);
529  else
530  for (int isq2 = 1; isq2 < 4; isq2++)
531  wid += norm(coupSUSYPtr->rvUDD[iq1][iq2][isq2]
532  * coupSUSYPtr->Rdsq[isq][isq2+3]);
533  }
534  }
535  else {
536  if ((id1Abs+id2Abs)%2 != 0) widNow = 0.0;
537  else
538  for (int isq2 = 1; isq2 < 4; isq2++)
539  wid += norm(coupSUSYPtr->rvUDD[isq2][iq1][iq2]
540  * coupSUSYPtr->Rusq[isq][isq2+3]);
541  }
542  }
543 
544  //Case 1b: LQD-type
545  else if (id1Abs < 17 && id2Abs < 7){
546  if (!coupSUSYPtr->isLQD) {widNow = 0; return;}
547 
548  int ilep = (id1Abs - 9)/2;
549  int iq = (id2Abs + 1)/2;
550 
551  fac = kinFac / (16.0 * M_PI * pow(mHat,3));
552  wid = 0.0;
553  if (idown){
554  if (iq%2 == 0){
555  // q is up-type; ~q is right-handed down type
556  for (int isq2=1; isq2<3; isq2++)
557  wid += norm(coupSUSYPtr->Rdsq[isq][isq2+3]
558  * coupSUSYPtr->rvLQD[ilep][iq][isq2]);
559  }else{
560  //q is down type; ~q left-handed down-type
561  for (int isq2=1; isq2<3; isq2++)
562  wid += norm(coupSUSYPtr->Rdsq[isq][isq2]
563  * coupSUSYPtr->rvLQD[ilep][isq2][isq2]);
564  }
565  }
566  else{
567  if (iq%2 == 0) {widNow = 0.0; return;}
568  // q is down type; ~q is left-handed up-type
569  for (int isq2=1; isq2<3; isq2++)
570  wid += norm(coupSUSYPtr->Rusq[isq][isq2]
571  * coupSUSYPtr->rvLQD[ilep][isq2][iq]);
572  }
573  }
574 
575  //Case 2: quark + gaugino
576  else if (id1Abs > ksusy && id2Abs < 7) {
577 
578  int iq = (id2Abs + 1)/2;
579 
580  // ~q -> ~g + q
581  if (id1Abs == 1000021 && idRes%10 == id2Abs) {
582  // Removed factor of s2W in denominator: strong process -- no EW
583  fac = 2.0 * alpS / (3.0 * pow3(mHat));
584  if (idown)
585  wid = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])
586  + norm(coupSUSYPtr->RsddG[isq][iq]))
587  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddG[isq][iq]
588  * conj(coupSUSYPtr->RsddG[isq][iq]));
589  else
590  wid = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])
591  + norm(coupSUSYPtr->RsuuG[isq][iq]))
592  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]
593  * conj(coupSUSYPtr->RsuuG[isq][iq]));
594  }
595  else
596  for (int i=1; i<6 ; i++){
597  // ~q -> ~chi0 + q
598  if (coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
599  fac = alpEM * preFac / (2.0 * (1 - s2W));
600  if (idown)
601  wid = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][i])
602  + norm(coupSUSYPtr->RsddX[isq][iq][i]))
603  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddX[isq][iq][i]
604  * conj(coupSUSYPtr->RsddX[isq][iq][i]));
605  else
606  wid = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][i])
607  + norm(coupSUSYPtr->RsuuX[isq][iq][i]))
608  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][i]
609  * conj(coupSUSYPtr->RsuuX[isq][iq][i]));
610  }
611 
612  // ~q -> chi- + q
613  else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs
614  && idRes%2 != id2Abs%2){
615 
616  fac = alpEM * preFac / (4.0 * (1 - s2W));
617  if (idown)
618  wid = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][i])
619  + norm(coupSUSYPtr->RsduX[isq][iq][i]))
620  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsduX[isq][iq][i]
621  * conj(coupSUSYPtr->RsduX[isq][iq][i]));
622  else
623  wid = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][i])
624  + norm(coupSUSYPtr->RsudX[isq][iq][i]))
625  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsudX[isq][iq][i]
626  * conj(coupSUSYPtr->RsudX[isq][iq][i]));
627  }
628  }
629  }
630 
631  //Case 3: ~q_i -> ~q_j + Z/W
632  else if (id1Abs > ksusy && id1Abs%100 < 7
633  && (id2Abs == 23 || id2Abs == 24)){
634 
635  // factor of lambda^(3/2) = ps^(3/2) ;
636  fac = alpEM * preFac/(16.0 * pow2(particleDataPtr->m0(id2Abs))
637  * (1.0 - s2W)) * pow2(ps) ;
638 
639  int isq2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
640 
641  if (id2Abs == 23 && id1Abs%2 == idRes%2){
642  if (idown)
643  wid = norm(coupSUSYPtr->LsdsdZ[isq][isq2]
644  + coupSUSYPtr->RsdsdZ[isq][isq2]);
645  else
646  wid = norm(coupSUSYPtr->LsusuZ[isq][isq2]
647  + coupSUSYPtr->RsusuZ[isq][isq2]);
648  }
649  else if (id2Abs == 24 && id1Abs%2 != idRes%2){
650  if (idown)
651  wid = norm(coupSUSYPtr->LsusdW[isq2][isq]);
652  else
653  wid = norm(coupSUSYPtr->LsusdW[isq][isq2]);
654  }
655  }
656 
657  // TODO: Case ~q_i -> ~q_j + h/H
658  widNow = fac * wid * ps * pow2(mHat);
659  if (DBSUSY) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs
660  <<" Width: "<<widNow<<endl;
661  return;
662  }
663 
664 }
665 
666 //==========================================================================
667 
668 // The ResonanceGluino class
669 // Derived class for Gluino resonances
670 
671 //--------------------------------------------------------------------------
672 
673 // Initialize constants.
674 
675 void ResonanceGluino::initConstants() {
676  return;
677 }
678 
679 //--------------------------------------------------------------------------
680 
681 // Calculate various common prefactors for the current mass.
682 
683 void ResonanceGluino::calcPreFac(bool) {
684 
685  // Common coupling factors.
686  alpS = coupSUSYPtr->alphaS(mHat * mHat);
687  preFac = alpS /( 8.0 * pow(mHat,3));
688  return;
689 
690 }
691 
692 //--------------------------------------------------------------------------
693 
694 // Calculate width for currently considered channel.
695 
696 void ResonanceGluino::calcWidth(bool) {
697 
698  widNow = 0.0;
699  if (ps == 0.) return;
700  kinFac = (mHat * mHat - mf1 * mf1 + mf2 * mf2);
701 
702  if (id1Abs > 1000000 && (id1Abs % 100) < 7 && id2Abs < 7) {
703 
704  int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
705  : (abs(id1Abs)%10+1)/2;
706  bool idown = id2Abs%2;
707  int iq = (id2Abs + 1)/2;
708 
709  // ~g -> ~q + q
710  if (idown){
711  widNow = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])
712  + norm(coupSUSYPtr->RsddG[isq][iq]))
713  + 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddG[isq][iq]
714  * conj(coupSUSYPtr->RsddG[isq][iq]));
715  }
716  else{
717  widNow = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])
718  + norm(coupSUSYPtr->RsuuG[isq][iq]))
719  + 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]
720  * conj(coupSUSYPtr->RsuuG[isq][iq]));
721  }
722  widNow = widNow * preFac * ps * pow2(mHat);
723  if (DBSUSY) {
724  cout<<"Gluino:: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
725  cout<<scientific<<widNow<<endl;
726  }
727  return;
728  }
729 }
730 
731 //==========================================================================
732 
733 // Class ResonanceNeut
734 // Derived class for Neutralino Resonances
735 //
736 //--------------------------------------------------------------------------
737 
738 void ResonanceNeut::initConstants() {
739 
740  s2W = coupSUSYPtr->sin2W;
741 
742  // Initialize functions for calculating 3-body widths
743  psi.init(particleDataPtr,coupSUSYPtr);
744  phi.init(particleDataPtr,coupSUSYPtr);
745  upsil.init(particleDataPtr,coupSUSYPtr);
746 
747 }
748 
749 //--------------------------------------------------------------------------
750 
751 // Calculate various common prefactors for the current mass.
752 
753 void ResonanceNeut::calcPreFac(bool) {
754 
755  // Common coupling factors.
756  alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
757  preFac = alpEM / (8.0 * s2W * pow(mHat,3));
758  return;
759 
760 }
761 
762 //--------------------------------------------------------------------------
763 
764 // Calculate width for currently considered channel.
765 void ResonanceNeut::calcWidth(bool){
766 
767  widNow = 0.0;
768 
769  if (ps == 0.) return;
770  else if (mult ==2){
771  // Two-body decays
772 
773  kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
774  kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
775  + pow2(mHat) * pow2(mf2) + pow2(mf1) * pow2(mf2)
776  - 2.0 * pow2(mHat) * pow2(mf1);
777 
778  // Stable lightest neutralino
779  if (idRes == 1000022) return;
780 
781  double fac = 0.0;
782  int iNeut1 = coupSUSYPtr->typeNeut(idRes);
783  int iNeut2 = coupSUSYPtr->typeNeut(id1Abs);
784  int iChar1 = coupSUSYPtr->typeChar(id1Abs);
785 
786  if (iNeut2>0 && id2Abs == 23){
787  // ~chi0_i -> chi0_j + Z
788  fac = kinFac2 * (norm(coupSUSYPtr->OLpp[iNeut1][iNeut2])
789  + norm(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
790  fac -= 12.0 * mHat * mf1 * pow2(mf2)
791  * real(coupSUSYPtr->OLpp[iNeut1][iNeut2]
792  * conj(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
793  fac /= pow2(mf2) * (1.0 - s2W);
794  }
795  else if (iChar1>0 && id2Abs==24){
796  // ~chi0_i -> chi+_j + W- (or c.c.)
797 
798  fac = kinFac2 * (norm(coupSUSYPtr->OL[iNeut1][iChar1])
799  + norm(coupSUSYPtr->OR[iNeut1][iChar1]));
800  fac -= 12.0 * mHat * mf1 * pow2(mf2)
801  * real(coupSUSYPtr->OL[iNeut1][iChar1]
802  * conj(coupSUSYPtr->OR[iNeut1][iChar1]));
803  fac /= pow2(mf2);
804  }
805  else if (id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
806  // ~chi0_k -> ~q + q
807  bool idown = (id1Abs%2 == 1);
808  int iq = (id2Abs + 1 )/ 2;
809  int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
810  : (abs(id1Abs)%10+1)/2;
811 
812  if (idown){
813  fac = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][iNeut1])
814  + norm(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
815  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddX[isq][iq][iNeut1]
816  * conj(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
817  }
818  else{
819  fac = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][iNeut1])
820  + norm(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
821  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][iNeut1]
822  * conj(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
823  }
824  // Extra multiplicative factor of 3 over sleptons
825  fac *= 6.0/(1 - s2W);
826  }
827  else if (id1Abs > 2000010 && id1Abs%2 == 0 ) {
828  // Check for right-handed neutralinos.
829  widNow = 0;
830  }
831  else if (id1Abs > 1000000 && id1Abs%100 > 10 && id1Abs%100 < 17
832  && id2Abs < 17){
833  // ~chi0_k -> ~l + l
834  bool idown = id2Abs%2;
835  int il = (id2Abs - 9)/ 2;
836  int isl = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
837  : (abs(id1Abs)%10+1)/2;
838 
839  if (idown){
840  fac = kinFac * (norm(coupSUSYPtr->LsllX[isl][il][iNeut1])
841  + norm(coupSUSYPtr->RsllX[isl][il][iNeut1]));
842  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsllX[isl][il][iNeut1]
843  * conj(coupSUSYPtr->RsllX[isl][il][iNeut1]));
844  }
845  else{
846  fac = kinFac * (norm(coupSUSYPtr->LsvvX[isl][il][iNeut1]));
847  }
848  fac *= 2.0/(1 - s2W);
849  }
850  // TODO: Decays in higgs
851  // Final width for 2-body decays
852  widNow = fac * preFac * ps * pow2(mHat);
853  if (DBSUSY) {
854  cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
855  cout<<scientific<<widNow<<endl;
856  }
857  }
858  else {
859  //RPV 3-body decays
860  //Case: UDD-type
861 
862  if (!coupSUSYPtr->isUDD) return;
863 
864  if (id1Abs < 7 && id2Abs < 7 && id3Abs < 7){
865 
866  // Check that quarks compatible with UDD structure
867  if ((id1Abs+id2Abs+id3Abs)%2 == 1) return;
868  double rvfac,m12min,m12max,integral;
869  int idInt;
870 
871  // Loop over mass eigenstate in actual propagator
872  for (int idIntRes = 1; idIntRes <= 6; idIntRes++){
873  // Right handed field in the UDD coupling
874  for (int iSq = 1; iSq <= 3; iSq++){
875  double m1, m2, m3, mixfac1(0.), mixfac2(0.), mixfac3(0.);
876  int itemp1,itemp2,itemp3;
877  // up/down-type internal lines
878  for (int itype = 1; itype<=3; itype++){
879  //itype = 1: up
880  //itype = 2: down1
881  //itype = 3: down2
882  if (itype ==1 ) idInt = coupSUSYPtr->idSup(idIntRes);
883  else idInt = coupSUSYPtr->idSdown(idIntRes);
884  if (id1Abs%2 == 0){
885  if (itype == 1){
886  itemp3 = id1Abs;
887  itemp1 = id2Abs;
888  itemp2 = id3Abs;
889  rvfac = pow2(
890  coupSUSYPtr->rvUDD[iSq][(id2Abs+1)/2][(id3Abs+1)/2]);
891  }else if (itype ==2){
892  itemp3 = id2Abs;
893  itemp1 = id1Abs;
894  itemp2 = id3Abs;
895  rvfac = pow2(
896  coupSUSYPtr->rvUDD[(id1Abs+1)/2][iSq][(id3Abs+1)/2]);
897  } else{
898  itemp3 = id3Abs;
899  itemp1 = id1Abs;
900  itemp2 = id2Abs;
901  rvfac = pow2(
902  coupSUSYPtr->rvUDD[(id1Abs+1)/2][(id2Abs+1)/2][iSq]);
903  }
904  }else if (id2Abs%2 == 0){
905  if (itype==1){
906  itemp3 = id2Abs;
907  itemp1 = id1Abs;
908  itemp2 = id3Abs;
909  rvfac = pow2(
910  coupSUSYPtr->rvUDD[iSq][(id1Abs+1)/2][(id3Abs+1)/2]);
911  }else if (itype ==2){
912  itemp3 = id1Abs;
913  itemp1 = id2Abs;
914  itemp2 = id3Abs;
915  rvfac = pow2(
916  coupSUSYPtr->rvUDD[(id2Abs+1)/2][iSq][(id3Abs+1)/2]);
917  } else{
918  itemp3 = id3Abs;
919  itemp1 = id2Abs;
920  itemp2 = id1Abs;
921  rvfac = pow2(
922  coupSUSYPtr->rvUDD[(id2Abs+1)/2][(id2Abs+1)/2][iSq]);
923  }
924  }else{
925  if (itype==1){
926  itemp3 = id3Abs;
927  itemp1 = id1Abs;
928  itemp2 = id2Abs;
929  rvfac = pow2(
930  coupSUSYPtr->rvUDD[iSq][(id1Abs+1)/2][(id2Abs+1)/2]);
931  }else if (itype ==2){
932  itemp3 = id2Abs;
933  itemp1 = id1Abs;
934  itemp2 = id3Abs;
935  rvfac = pow2(
936  coupSUSYPtr->rvUDD[(id3Abs+1)/2][iSq][(id2Abs+1)/2]);
937  } else{
938  itemp3 = id3Abs;
939  itemp1 = id1Abs;
940  itemp2 = id2Abs;
941  rvfac = pow2(
942  coupSUSYPtr->rvUDD[(id3Abs+1)/2][(id1Abs+1)/2][iSq]);
943  }
944 
945  }
946 
947  m1 = particleDataPtr->m0(itemp1);
948  m2 = particleDataPtr->m0(itemp2);
949  m3 = particleDataPtr->m0(itemp3);
950 
951  m12min = pow2(m1+m2);
952  m12max = pow2(mHat-m3);
953 
954  // Ignore mode when 2-body decay is possible
955  if (mRes > particleDataPtr->m0(idInt)
956  + particleDataPtr->m0(itemp3)) continue;
957 
958  // Single diagram squared terms
959  psi.setInternal(idRes, itemp1, itemp2, itemp3, idInt, 0);
960  // Mixing with R-states
961  if (itype == 1)
962  mixfac1 = norm(coupSUSYPtr->Rusq[idIntRes][iSq+3]);
963  else
964  mixfac1 = norm(coupSUSYPtr->Rdsq[idIntRes][iSq+3]);
965 
966  if (abs(rvfac * mixfac1) > 1.0e-8) {
967  integral = integrateGauss(&psi,m12min,m12max,1.0e-4);
968  widNow += rvfac * mixfac1 * integral;
969  // if (DBSUSY || idRes == 1000023)
970  // cout << scientific << setw(10) <<"Psi: intRes: "<<idInt
971  // <<" integral:"<<integral<<" mixfac:"<<mixfac1
972  // <<" widNow:"<<widNow<<endl;
973  }
974 
975  // Mixing of diagrams with different internal squarks
976  // of same isospin
977  for (int idIntRes2 = 1; idIntRes2 <= 6; idIntRes2++){
978  if (idIntRes2 == idIntRes) continue;
979  int idInt2;
980  if (itype == 1 ){
981  idInt2 = coupSUSYPtr->idSup(idIntRes2);
982  mixfac2 = 2.0 * real(coupSUSYPtr->Rusq[idIntRes][iSq+3]
983  * conj(coupSUSYPtr->Rusq[idIntRes2][iSq+3]));
984  } else {
985  idInt2 = coupSUSYPtr->idSdown(idIntRes2);
986  mixfac2 = 2.0 * real(coupSUSYPtr->Rdsq[idIntRes][iSq+3]
987  * conj(coupSUSYPtr->Rdsq[idIntRes2][iSq+3]));
988  }
989 
990  // Ignore mode when 2-body decay is possible
991  if (mRes > particleDataPtr->m0(idInt2)
992  + particleDataPtr->m0(itemp3)) continue;
993 
994  upsil.setInternal(idRes,itemp1, itemp2,itemp3,idInt,idInt2);
995  if (abs(rvfac * mixfac2) > 0.0) {
996  integral = integrateGauss(&upsil,m12min,m12max,1.0e-4);
997  widNow += rvfac * mixfac2 * integral;
998  // if (DBSUSY || idRes == 1000023)
999  // cout << scientific << setw(10) <<"Upsilon: intRes: "
1000  // <<idInt<<" intRes2:"<<idInt2<<" integral:"<<integral
1001  // <<" mixfac:"<<mixfac2<<" widNow:"<<widNow<<endl;
1002  }
1003  }
1004 
1005  // Interference between two diagrams with quarks
1006  // of different isospin
1007 
1008  for (int idIntRes2 = 1; idIntRes2 <= 6; idIntRes2++){
1009  if (itype != 1 && idIntRes2 == idIntRes) continue;
1010  int idInt2;
1011 
1012  for (int iSq2 = 1; iSq2 <= 3; iSq2++){
1013  if (itype == 1 ){
1014  idInt2 = coupSUSYPtr->idSdown(idIntRes2);
1015  mixfac3 = 2.0 * real(coupSUSYPtr->Rusq[idIntRes][iSq+3]
1016  * conj(coupSUSYPtr->Rdsq[idIntRes2][iSq2+3]));
1017  } else {
1018  idInt2 = coupSUSYPtr->idSdown(idIntRes2);
1019  mixfac3 = 2.0 * real(coupSUSYPtr->Rdsq[idIntRes][iSq+3]
1020  * conj(coupSUSYPtr->Rdsq[idIntRes2][iSq2+3]));
1021  }
1022 
1023  if (abs(rvfac * mixfac3) > 0.0) {
1024  phi.setInternal(idRes,itemp1, itemp2,itemp3,idInt,idInt2);
1025  //integral = 0;
1026  //if (idIntRes == 2 && iSq2 ==4)
1027  integral = integrateGauss(&phi,m12min,m12max,1.0e-4);
1028  widNow -= rvfac * mixfac2 * integral;
1029  }
1030  }
1031  }
1032  }
1033  }
1034  }
1035  }
1036  // Normalisation. Extra factor of 2 to account for the fact that
1037  // d_i, d_j will always be ordered in ascending order. N_c! = 6
1038  widNow *= 12.0 /(pow3(2.0 * M_PI * mHat) * 32.0);
1039  }
1040 
1041  return;
1042 }
1043 
1044 //==========================================================================
1045 
1046 // Class ResonanceChar
1047 // Derived class for Neutralino Resonances
1048 // Decays into higgses/sleptons not yet implemented
1049 
1050 //--------------------------------------------------------------------------
1051 
1052 void ResonanceChar::initConstants() {
1053 
1054  s2W = coupSUSYPtr->sin2W;
1055  return;
1056 
1057 }
1058 
1059 //--------------------------------------------------------------------------
1060 
1061 // Calculate various common prefactors for the current mass.
1062 
1063 void ResonanceChar::calcPreFac(bool) {
1064 
1065  alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
1066  preFac = alpEM / (8.0 * s2W * pow(mHat,3));
1067  return;
1068 
1069 }
1070 
1071 //--------------------------------------------------------------------------
1072 
1073 // Calculate width for currently considered channel.
1074 
1075 void ResonanceChar::calcWidth(bool) {
1076 
1077  widNow = 0.0;
1078  if (ps == 0.) return;
1079 
1080  if (mult ==2){
1081  double fac = 0.0;
1082  kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
1083  kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
1084  + pow2(mHat) * pow2(mf2) + pow2(mf1)
1085  * pow2(mf2) - 2.0 * pow2(mHat) * pow2(mf1);
1086 
1087  int idChar1 = coupSUSYPtr->typeChar(idRes);
1088  int idChar2 = coupSUSYPtr->typeChar(id1Abs);
1089  int idNeut1 = coupSUSYPtr->typeNeut(id1Abs);
1090 
1091  if (idChar2>0 && id2Abs == 23){
1092  // ~chi_i -> chi_j + Z
1093  fac = kinFac2 * (norm(coupSUSYPtr->OLp[idChar1][idChar2])
1094  + norm(coupSUSYPtr->ORp[idChar1][idChar2]));
1095  fac -= 12.0 * mHat * mf1 * pow2(mf2)
1096  * real(coupSUSYPtr->OLp[idChar1][idChar2]
1097  * conj(coupSUSYPtr->ORp[idChar1][idChar2]));
1098  fac /= pow2(mf2) * (1.0 - s2W);
1099  }
1100  else if (idNeut1>0 && id2Abs==24){
1101  // ~chi_i -> chi0_j + W- (or c.c.)
1102 
1103  fac = kinFac2 * (norm(coupSUSYPtr->OL[idNeut1][idChar1])
1104  + norm(coupSUSYPtr->OR[idNeut1][idChar1]));
1105  fac -= 12.0 * mHat * mf1 * pow2(mf2)
1106  * real(coupSUSYPtr->OL[idNeut1][idChar1]
1107  * conj(coupSUSYPtr->OR[idNeut1][idChar1]));
1108  fac /= pow2(mf2);
1109  }
1110  else if (id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
1111  // ~chi0_k -> ~q + q
1112  bool idown = (id1Abs%2 == 1);
1113  int iq = (id2Abs + 1 )/ 2;
1114  int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1115  : (abs(id1Abs)%10+1)/2;
1116 
1117  if (idown){
1118  fac = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][idChar1])
1119  + norm(coupSUSYPtr->RsduX[isq][iq][idChar1]));
1120  fac += 4.0 * mHat * mf2
1121  * real(coupSUSYPtr->LsduX[isq][iq][idChar1]
1122  * conj(coupSUSYPtr->RsduX[isq][iq][idChar1]));
1123  }
1124  else{
1125  fac = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][idChar1])
1126  + norm(coupSUSYPtr->RsudX[isq][iq][idChar1]));
1127  fac += 4.0 * mHat * mf2
1128  * real(coupSUSYPtr->LsudX[isq][iq][idChar1]
1129  * conj(coupSUSYPtr->RsudX[isq][iq][idChar1]));
1130  }
1131  fac *= 6.0/(1 - s2W);
1132  }
1133  else if (id1Abs > 2000010 && id1Abs%2 == 0 ) {
1134  // Check for right-handed neutralinos.
1135  widNow = 0;
1136  }
1137  else if (id1Abs > 1000000 && id1Abs%100 > 10 && id1Abs%100 < 17
1138  && id2Abs < 17){
1139  // ~chi+_k -> ~l + l
1140  bool idown = id2Abs%2;
1141  int il = (id2Abs - 9)/ 2;
1142  int isl = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1143  : (abs(id1Abs)%10+1)/2;
1144 
1145  if (idown){
1146  fac = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][idChar1])
1147  + norm(coupSUSYPtr->RslvX[isl][il][idChar1]));
1148  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][idChar1]
1149  * conj(coupSUSYPtr->RslvX[isl][il][idChar1]));
1150  }
1151  else{
1152  fac = kinFac * (norm(coupSUSYPtr->LsvlX[isl][il][idChar1]));
1153  }
1154  fac *= 2.0/(1 - s2W);
1155  }
1156 
1157  // TODO: Decays in higgs
1158  widNow = fac * preFac * ps * pow2(mHat);
1159  if (DBSUSY) {
1160  cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
1161  cout<<scientific<<widNow<<endl;
1162  }
1163  }else{
1164  //TODO: Implement Chargino 3-body decays
1165  }
1166  return;
1167 }
1168 
1169 
1170 //==========================================================================
1171 // The ResonanceSlepton class
1172 // Derived class for Slepton (and sneutrino) resonances
1173 
1174 //--------------------------------------------------------------------------
1175 
1176 // Initialize constants.
1177 
1178 void ResonanceSlepton::initConstants() {
1179 
1180  // Locally stored properties and couplings.
1181  s2W = coupSUSYPtr->sin2W;
1182 
1183 }
1184 
1185 //--------------------------------------------------------------------------
1186 
1187 // Calculate various common prefactors for the current mass.
1188 
1189 void ResonanceSlepton::calcPreFac(bool) {
1190 
1191  // Common coupling factors.
1192  alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
1193  preFac = 1.0 / (s2W * pow(mHat,3));
1194 
1195 }
1196 
1197 //--------------------------------------------------------------------------
1198 
1199 // Calculate width for currently considered channel.
1200 
1201 void ResonanceSlepton::calcWidth(bool) {
1202 
1203  // Slepton type -- in u_i/d_i and generation
1204  int ksusy = 1000000;
1205  int isl = (abs(idRes)/ksusy == 2) ? (abs(idRes)%10+1)/2 + 3
1206  : (abs(idRes)%10+1)/2;
1207  int il = (id2Abs-9)/2;
1208  bool islep = idRes%2;
1209 
1210  // Check that mass is above threshold.
1211  if (ps == 0.) return;
1212  else{
1213  // Two-body decays
1214  kinFac = (mHat * mHat - mf1 * mf1 - mf2 * mf2);
1215 
1216  double fac = 0.0 , wid = 0.0;
1217 
1218  //Case 1: RPV: To be implemented
1219  //Case 2: slepton + gaugino
1220 
1221  if (id1Abs > ksusy && id2Abs > 10 && id2Abs < 17) {
1222  for (int i=1; i<6 ; i++){
1223  // ~ell/~nu -> ~chi0 + ell/nu
1224  if (coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
1225  fac = alpEM * preFac / (2.0 * (1 - s2W));
1226  if (islep)
1227  wid = kinFac * (norm(coupSUSYPtr->LsllX[isl][il][i])
1228  + norm(coupSUSYPtr->RsllX[isl][il][i]))
1229  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsllX[isl][il][i]
1230  * conj(coupSUSYPtr->RsllX[isl][il][i]));
1231  else
1232  wid = kinFac * (norm(coupSUSYPtr->LsvvX[isl][il][i])
1233  + norm(coupSUSYPtr->RsvvX[isl][il][i]))
1234  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsvvX[isl][il][i]
1235  * conj(coupSUSYPtr->RsvvX[isl][il][i]));
1236  }
1237 
1238  // ~ell/~nu -> ~chi- + nu/ell
1239  else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs
1240  && idRes%2 != id2Abs%2){
1241 
1242  fac = alpEM * preFac / (4.0 * (1 - s2W));
1243  if (islep)
1244  wid = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][i])
1245  + norm(coupSUSYPtr->RslvX[isl][il][i]))
1246  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][i]
1247  * conj(coupSUSYPtr->RslvX[isl][il][i]));
1248  else
1249  wid = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][i])
1250  + norm(coupSUSYPtr->RslvX[isl][il][i]))
1251  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][i]
1252  * conj(coupSUSYPtr->RslvX[isl][il][i]));
1253  }
1254  }
1255  }
1256 
1257  //Case 3: ~l_i -> ~l_j + Z/W
1258  else if (id1Abs > ksusy+10 && id1Abs%100 < 17
1259  && (id2Abs == 23 || id2Abs == 24)){
1260 
1261  // factor of lambda^(3/2) = ps^3;
1262  fac = alpEM * preFac/(16.0 * pow2(mf2) * (1.0 - s2W)) * pow2(ps) ;
1263 
1264  int isl2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
1265 
1266  if (id2Abs == 23 && id1Abs%2 == idRes%2){
1267  if (islep)
1268  wid = norm(coupSUSYPtr->LslslZ[isl][isl2]
1269  + coupSUSYPtr->RslslZ[isl][isl2]);
1270  else
1271  wid = norm(coupSUSYPtr->LsvsvZ[isl][isl2]
1272  + coupSUSYPtr->RsvsvZ[isl][isl2]);
1273  }
1274  else if (id2Abs == 24 && id1Abs%2 != idRes%2){
1275  if (islep)
1276  wid = norm(coupSUSYPtr->LslsvW[isl2][isl]);
1277  else
1278  wid = norm(coupSUSYPtr->LslsvW[isl][isl2]);
1279  }
1280  }
1281 
1282  // TODO: Case ~l_i -> ~l_j + h/H
1283 
1284 
1285  widNow = fac * wid * ps * pow2(mHat);
1286  if (DBSUSY) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs
1287  <<" Width: "<<widNow<<endl;
1288  return;
1289  }
1290 
1291 }
1292 
1293 //==========================================================================
1294 
1295 // Gaussian Integrator for 3-body decay widths
1296 
1297 double SUSYResonanceWidths::integrateGauss(WidthFunction* widthFn,
1298  double xlo, double xhi, double tol) {
1299 
1300  //8-point unweighted
1301  static double x8[4]={0.96028985649753623,
1302  0.79666647741362674,
1303  0.52553240991632899,
1304  0.18343464249564980};
1305  static double w8[4]={0.10122853629037626,
1306  0.22238103445337447,
1307  0.31370664587788729,
1308  0.36268378337836198};
1309  //16-point unweighted
1310  static double x16[8]={0.98940093499164993,
1311  0.94457502307323258,
1312  0.86563120238783174,
1313  0.75540440835500303,
1314  0.61787624440264375,
1315  0.45801677765722739,
1316  0.28160355077925891,
1317  0.09501250983763744};
1318  static double w16[8]={0.027152459411754095,
1319  0.062253523938647893,
1320  0.095158511682492785,
1321  0.12462897125553387,
1322  0.14959598881657673,
1323  0.16915651939500254,
1324  0.18260341504492359,
1325  0.18945061045506850};
1326  //boundary checks
1327  if (xlo == xhi) {
1328  cerr<<"xlo = xhi"<<endl;
1329  return 0.0;
1330  }
1331  if (xlo > xhi) {
1332  cerr<<" (integrateGauss:) -> xhi < xlo"<<endl;
1333  return 0.0;
1334  }
1335  //initialize
1336  double sum = 0.0;
1337  double c = 0.001/abs(xhi-xlo);
1338  double zlo = xlo;
1339  double zhi = xhi;
1340 
1341  bool nextbin = true;
1342 
1343  while ( nextbin ) {
1344 
1345  double zmi = 0.5*(zhi+zlo); // midpoint
1346  double zmr = 0.5*(zhi-zlo); // midpoint, relative to zlo
1347 
1348  //calculate 8-point and 16-point quadratures
1349  double s8=0.0;
1350  for (int i=0;i<4;i++) {
1351  double dz = zmr * x8[i];
1352  s8 += w8[i]*(widthFn->function(zmi+dz) + widthFn->function(zmi-dz));
1353  }
1354  s8 *= zmr;
1355  double s16=0.0;
1356  for (int i=0;i<8;i++) {
1357  double dz = zmr * x16[i];
1358  s16 += w16[i]*(widthFn->function(zmi+dz) + widthFn->function(zmi-dz));
1359  }
1360  s16 *= zmr;
1361  if (abs(s16-s8) < tol*(1+abs(s16))) {
1362  //precision in this bin OK, add to cumulative and go to next
1363  nextbin=true;
1364  sum += s16;
1365  //next bin: LO = end of current, HI = end of integration region.
1366  zlo=zhi;
1367  zhi=xhi;
1368  if ( zlo == zhi ) nextbin = false;
1369  } else {
1370  //precision in this bin not OK, subdivide.
1371  if (1.0 + c*abs(zmr) == 1.0) {
1372  cerr << " (integrateGauss:) too high accuracy required"<<endl;
1373  sum = 0.0 ;
1374  break;
1375  }
1376  zhi=zmi;
1377  nextbin=true;
1378  }
1379  }
1380 
1381  return sum;
1382 }
1383 
1384 //======================================================================
1385 
1386 } //end namespace Pythia8
1387 
1388 
Generated by   doxygen 1.8.5