]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8170/src/SusyResonanceWidths.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Stupid bug fix in new superlight mode (from Zurich airport)
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / SusyResonanceWidths.cxx
CommitLineData
63ba5337 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