source: HiSusy/trunk/Pythia8/pythia8170/src/SusyResonanceWidths.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES
File size: 47.5 KB
Line 
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
14namespace 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
25void WidthFunction::init( ParticleData* particleDataPtrIn, 
26  CoupSUSY* coupSUSYPtrIn) {
27
28  particleDataPtr = particleDataPtrIn;
29  coupSUSYPtr = coupSUSYPtrIn;
30
31}
32
33//--------------------------------------------------------------------------
34
35void 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
59double WidthFunction::function(double) {
60
61  cout<<"Warning using dummy width function"<<endl;
62  return 0.;
63}
64
65//--------------------------------------------------------------------------
66
67double 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
79void 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
98void 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
122void 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
146double 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
177double 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
222double 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
246double 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
329double 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
419const bool SUSYResonanceWidths::DEBUG = false;
420
421//--------------------------------------------------------------------------
422
423bool 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
636void 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
648void 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
659void 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
843void 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
854void 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
864void 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
908void 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.
922void  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.
932void  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
1215void 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.
1225void  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.
1234void  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
1333void 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
1344void 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
1355void 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
1451double 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
Note: See TracBrowser for help on using the repository browser.