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