]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8145/src/SusyResonanceWidths.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Use Output directive instead of the old OutputFile and OUtputArchive. Save fileinfo...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8145 / src / SusyResonanceWidths.cxx
CommitLineData
9419eeef 1// ResonanceWidths.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2010
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6// Function definitions (not found in the header) for
7// the ResonanceWidths class and classes derived from it.
8
9#include "SusyResonanceWidths.h"
10#include "SusyCouplings.h"
11#include "PythiaComplex.h"
12
13namespace Pythia8 {
14
15//==========================================================================
16
17// The SUSYResonanceWidths Class
18// Derived class for SUSY resonances
19
20const bool SUSYResonanceWidths::DEBUG = false;
21
22//--------------------------------------------------------------------------
23
24bool SUSYResonanceWidths::init(Info* infoPtrIn, Settings* settingsPtrIn,
25 ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
26
27 // Save pointers.
28 infoPtr = infoPtrIn;
29 settingsPtr = settingsPtrIn;
30 particleDataPtr = particleDataPtrIn;
31 coupSUSYPtr = (couplingsPtrIn->isSUSY ? (CoupSUSY *) couplingsPtrIn: 0 );
32
33 // No initialization necessary for SM-only
34 if(!couplingsPtrIn->isSUSY) return true;
35
36 // Minimal decaying-resonance width. Minimal phase space for meMode = 103.
37 minWidth = settingsPtr->parm("ResonanceWidths:minWidth");
38 minThreshold = settingsPtr->parm("ResonanceWidths:minThreshold");
39
40 // Pointer to particle species.
41 particlePtr = particleDataPtr->particleDataEntryPtr(idRes);
42 if (particlePtr == 0) {
43 infoPtr->errorMsg("Error in ResonanceWidths::init:"
44 " unknown resonance identity code");
45 return false;
46 }
47
48 // Generic particles should not have meMode < 100, but allow
49 // some exceptions: not used Higgses and not used Technicolor.
50 if (idRes == 35 || idRes == 36 || idRes == 37
51 || idRes/1000000 == 3) isGeneric = false;
52
53 // Resonance properties: antiparticle, mass, width
54 hasAntiRes = particlePtr->hasAnti();
55 mRes = particlePtr->m0();
56 GammaRes = particlePtr->mWidth();
57 m2Res = mRes*mRes;
58
59 // For very narrow resonances assign fictitious small width.
60 if (GammaRes < minWidth) GammaRes = 0.1 * minWidth;
61 GamMRat = GammaRes / mRes;
62
63 // Secondary decay chains by default all on.
64 openPos = 1.;
65 openNeg = 1.;
66
67 // Allow option where on-shell width is forced to current value.
68 doForceWidth = particlePtr->doForceWidth();
69 forceFactor = 1.;
70
71 // Check that resonance OK.
72 if (particlePtr == 0) infoPtr->errorMsg("Error in ResonanceWidths::init:"
73 " unknown resonance identity code");
74
75 // Calculate various common prefactors for the current mass.
76 mHat = mRes;
77
78 // Initialize constants used for a resonance.
79
80 initConstants();
81 calcPreFac(true);
82
83 // Reset quantities to sum. Declare variables inside loop.
84 double widTot = 0.;
85 double widPos = 0.;
86 double widNeg = 0.;
87 int idNow, idAnti;
88 double openSecPos, openSecNeg;
89
90 // Loop over all decay channels. Basic properties of channel.
91 for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
92 iChannel = i;
93 onMode = particlePtr->channel(i).onMode();
94 meMode = particlePtr->channel(i).meMode();
95 mult = particlePtr->channel(i).multiplicity();
96 widNow = 0.;
97
98 // Warn if not relevant meMode.
99 if ( meMode < 0 || meMode > 103 || (isGeneric && meMode < 100) ) {
100 infoPtr->errorMsg("Error in ResonanceWidths::init:"
101 " resonance meMode not acceptable");
102 }
103
104 // Check if decay table was read in via SLHA
105 bool hasDecayTable = false;
106 for(unsigned int iDec = 1; iDec < (coupSUSYPtr->slhaPtr)->decays.size(); iDec++)
107 hasDecayTable = ((coupSUSYPtr->slhaPtr)->decays[iDec].getId() == abs(idRes));
108
109 // Calculation of SUSY particle widths
110 if (meMode == 103 && GammaRes > 0. &&
111 (!settingsPtr->flag("SLHA:useDecayTable") || !hasDecayTable)) {
112 // Read out information on channel: primarily use first two.
113 id1 = particlePtr->channel(i).product(0);
114 id2 = particlePtr->channel(i).product(1);
115 id1Abs = abs(id1);
116 id2Abs = abs(id2);
117
118 // Order first two in descending order of absolute values.
119 if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
120
121 // Allow for third product to be treated in derived classes.
122 if (mult > 2) {
123 id3 = particlePtr->channel(i).product(2);
124 id3Abs = abs(id3);
125
126 // Also order third into descending order of absolute values.
127 if (id3Abs > id2Abs) {swap( id2, id3); swap( id2Abs, id3Abs);}
128 if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
129 }
130
131 // Read out masses. Calculate two-body phase space.
132 mf1 = particleDataPtr->m0(id1Abs);
133 mf2 = particleDataPtr->m0(id2Abs);
134 mr1 = pow2(mf1 / mHat);
135 mr2 = pow2(mf2 / mHat);
136 ps = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
137 : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
138 if (mult > 2) {
139 mf3 = particleDataPtr->m0(id3Abs);
140 mr3 = pow2(mf3 / mHat);
141 }
142
143 // Let derived class calculate width for channel provided.
144 calcWidth(true);
145 }
146
147 // Channels with meMode >= 100 are calculated based on stored values.
148 else widNow = GammaRes * particlePtr->channel(i).bRatio();
149
150 // Find secondary open fractions of partial width.
151 openSecPos = 1.;
152 openSecNeg = 1.;
153 if (widNow > 0.) for (int j = 0; j < mult; ++j) {
154 idNow = particlePtr->channel(i).product(j);
155 idAnti = (particleDataPtr->hasAnti(idNow)) ? -idNow : idNow;
156 // Secondary widths not yet initialized for heavier states,
157 // so have to assume unit open fraction there.
158 if (idNow == 23 || abs(idNow) == 24
159 || particleDataPtr->m0(abs(idNow)) < mRes) {
160 openSecPos *= particleDataPtr->resOpenFrac(idNow);
161 openSecNeg *= particleDataPtr->resOpenFrac(idAnti);
162 }
163 }
164
165 // Store partial widths and secondary open fractions.
166 particlePtr->channel(i).onShellWidth(widNow);
167 particlePtr->channel(i).openSec( idRes, openSecPos);
168 particlePtr->channel(i).openSec(-idRes, openSecNeg);
169
170 // Update sum over all channnels and over open channels only.
171 widTot += widNow;
172 if (onMode == 1 || onMode == 2) widPos += widNow * openSecPos;
173 if (onMode == 1 || onMode == 3) widNeg += widNow * openSecNeg;
174 }
175
176 // If no decay channels are open then set particle stable and done.
177 if (widTot < minWidth) {
178 particlePtr->setMayDecay(false, false);
179 particlePtr->setMWidth(0., false);
180 for (int i = 0; i < particlePtr->sizeChannels(); ++i)
181 particlePtr->channel(i).bRatio( 0., false);
182 return true;
183 }
184
185 // Normalize branching ratios to unity.
186 double bRatio;
187 for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
188 bRatio = particlePtr->channel(i).onShellWidth() / widTot;
189 particlePtr->channel(i).bRatio( bRatio, false);
190 }
191
192 // Optionally force total width by rescaling of all partial ones.
193 if (doForceWidth) {
194 forceFactor = GammaRes / widTot;
195 for (int i = 0; i < particlePtr->sizeChannels(); ++i)
196 particlePtr->channel(i).onShellWidthFactor( forceFactor);
197 }
198
199 // Else update newly calculated partial width.
200 else {
201 particlePtr->setMWidth(widTot, false);
202 GammaRes = widTot;
203 }
204
205 // Updated width-to-mass ratio. Secondary widths for open.
206 GamMRat = GammaRes / mRes;
207 openPos = widPos / widTot;
208 openNeg = widNeg / widTot;
209
210 // Done.
211 return true;
212
213}
214
215
216
217//==========================================================================
218
219// The ResonanceSquark class
220// Derived class for Squark resonances
221
222//--------------------------------------------------------------------------
223
224// Initialize constants.
225
226void ResonanceSquark::initConstants() {
227
228 // Locally stored properties and couplings.
229 alpS = coupSUSYPtr->alphaS(mHat * mHat );
230 alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
231 s2W = coupSUSYPtr->sin2W;
232}
233
234//--------------------------------------------------------------------------
235
236// Calculate various common prefactors for the current mass.
237
238void ResonanceSquark::calcPreFac(bool) {
239
240 // Common coupling factors.
241 preFac = 1.0 / (s2W * pow(mHat,3));
242
243}
244
245//--------------------------------------------------------------------------
246
247// Calculate width for currently considered channel.
248
249void ResonanceSquark::calcWidth(bool) {
250
251 // Squark type -- in u_i/d_i and generation
252 int ksusy = 1000000;
253 bool idown = (abs(idRes)%2 == 0 ? false : true);
254 int isq = (abs(idRes)/ksusy == 2) ?
255 (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
256 int isqgen = (abs(idRes)%10 + 1)/2;
257
258 // Check that above threshold.
259 if (ps == 0.) return;
260 kinFac = (mHat * mHat - mf1 * mf1 -mf2 * mf2);
261 lambda = lam(mHat*mHat, mf1*mf1, mf2*mf2);
262
263 double fac = 0.0 , wid = 0.0;
264
265 //Case 1: RPV decay
266 if(id1Abs < 7 && id2Abs < 7){
267
268 //Temporary till colour assignments sorted out
269 //widNow = 0.0;
270 //return;
271
272 // Quark generations
273 int iq1 = (id1Abs + 1)/2;
274 int iq2 = (id2Abs + 1)/2;
275
276 // Check for RPV UDD couplings
277 if(!coupSUSYPtr->isUDD) {widNow = 0; return;}
278
279 // ~q -> q_i + q_j
280
281 fac = mHat;
282 if(idown) {
283 if ((id1Abs+id2Abs)%2 == 1){
284 if(id1Abs%2==1)
285 for(int isq2 = 4; isq2 < 7; isq2++)
286 wid = coupSUSYPtr->rvUDD[iq2][iq1][isqgen] *
287 norm(coupSUSYPtr->Rdsq[isq][isq2]);
288 else
289 for(int isq2 = 4; isq2 < 7; isq2++)
290 wid = coupSUSYPtr->rvUDD[iq1][iq2][isqgen] *
291 norm(coupSUSYPtr->Rdsq[isq][isq2]);
292 }
293 }
294 else {
295 if ((id1Abs+id2Abs)%2 != 0) widNow = 0.0;
296 else
297 for(int isq2 = 4; isq2 < 7; isq2++)
298 wid = coupSUSYPtr->rvUDD[isq][iq1][iq2] *
299 norm(coupSUSYPtr->Rusq[isq][isq2]);
300 }
301 }
302
303 //Case 2: quark + gaugino (higgsino)
304 else if (id1Abs > ksusy && id2Abs < 7) {
305
306 int iq = (id2Abs + 1)/2;
307
308 // ~q -> ~g + q
309 if(id1Abs == 1000021 && idRes%10 == id2Abs) {
310 fac = 2.0 / 3.0 * alpS * preFac * sqrt(lambda);
311 if(idown)
312 wid = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])+norm(coupSUSYPtr->RsddG[isq][iq]))
313 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsddG[isq][iq]*coupSUSYPtr->RsddG[isq][iq]);
314 else
315 wid = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])+norm(coupSUSYPtr->RsuuG[isq][iq]))
316 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]*coupSUSYPtr->RsuuG[isq][iq]);
317
318 }
319 else
320 for(int i=1; i<6 ; i++){
321 // ~q -> ~chi0 + q
322 if(coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
323 fac = alpEM * preFac * sqrt(lambda)/ (4.0 * (1 - s2W));
324 if(idown)
325 wid = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][i]) + norm(coupSUSYPtr->RsddX[isq][iq][i]))
326 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsddX[isq][iq][i]*coupSUSYPtr->RsddX[isq][iq][i]);
327 else
328 wid = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][i]) + norm(coupSUSYPtr->RsuuX[isq][iq][i]))
329 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][i]*coupSUSYPtr->RsuuX[isq][iq][i]);
330 }
331
332 // ~q -> chi- + q
333 else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs && idRes%2 != id2Abs%2){
334
335 fac = alpEM * preFac * sqrt(lambda)/ (4.0 * (1 - s2W));
336 if(idown)
337 wid = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][i]) + norm(coupSUSYPtr->RsduX[isq][iq][i]))
338 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsduX[isq][iq][i]*coupSUSYPtr->RsduX[isq][iq][i]);
339 else
340 wid = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][i]) + norm(coupSUSYPtr->RsudX[isq][iq][i]))
341 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsudX[isq][iq][i]*coupSUSYPtr->RsudX[isq][iq][i]);
342 }
343 }
344 }
345
346 //Case 3: ~q_i -> ~q_j + Z/W
347 else if (id1Abs > ksusy && id1Abs%100 < 7 && (id2Abs == 23 || id2Abs == 24)){
348
349 fac = alpEM * preFac/(16.0 * pow2(particleDataPtr->m0(id2Abs)) * (1.0 - s2W))
350 * pow(lambda, 1.5);
351
352 int isq2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
353
354 if(id2Abs == 23 && id1Abs%2 == idRes%2){
355 if(idown)
356 wid = norm(coupSUSYPtr->LsdsdZ[isq][isq2] + coupSUSYPtr->RsdsdZ[isq][isq2]);
357 else
358 wid = norm(coupSUSYPtr->LsusuZ[isq][isq2] + coupSUSYPtr->RsusuZ[isq][isq2]);
359 }
360 else if (id2Abs == 24 && id1Abs%2 != idRes%2){
361 if(idown)
362 wid = norm(coupSUSYPtr->LsusdW[isq2][isq]);
363 else
364 wid = norm(coupSUSYPtr->LsusdW[isq][isq2]);
365 }
366 }
367
368 widNow = fac * wid * ps;
369 if(DEBUG) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: "<<widNow<<endl;
370 return;
371
372}
373
374//==========================================================================
375
376// The ResonanceGluino class
377// Derived class for Gluino resonances
378
379//--------------------------------------------------------------------------
380
381// Initialize constants.
382
383void ResonanceGluino::initConstants() {
384
385 // Locally stored properties and couplings.
386 alpS = coupSUSYPtr->alphaS(mHat * mHat );
387 return;
388}
389
390//--------------------------------------------------------------------------
391
392// Calculate various common prefactors for the current mass.
393
394void ResonanceGluino::calcPreFac(bool) {
395 // Common coupling factors.
396 preFac = alpS /( 8.0 * pow(mHat,3));
397 return;
398}
399
400//--------------------------------------------------------------------------
401
402// Calculate width for currently considered channel.
403
404void ResonanceGluino::calcWidth(bool) {
405
406
407 widNow = 0.0;
408 if(ps == 0.) return;
409 kinFac = (mHat * mHat - mf1 * mf1 + mf2 * mf2);
410 lambda = lam(mHat*mHat, mf1*mf1 , mf2*mf2);
411
412 if(id1Abs > 1000000 && (id1Abs % 100) < 7 && id2Abs < 7) {
413
414 int isq = (abs(id1Abs)/1000000 == 2) ?
415 (abs(id1Abs)%10+1)/2 + 3: (abs(id1Abs)%10+1)/2;
416 bool idown = (id2Abs%2 == 1);
417 int iq = (id2Abs + 1)/2;
418
419 // ~g -> ~q + q
420 if(idown){
421 widNow = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])+norm(coupSUSYPtr->RsddG[isq][iq]))
422 + 4.0 * mHat * mf2 * real( coupSUSYPtr->LsddG[isq][iq]
423 * conj(coupSUSYPtr->RsddG[isq][iq]));
424 }
425 else{
426 widNow = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])+norm(coupSUSYPtr->RsuuG[isq][iq]))
427 + 4.0 * mHat * mf2 * real( coupSUSYPtr->LsuuG[isq][iq]
428 * conj(coupSUSYPtr->RsuuG[isq][iq]));
429
430 }
431 widNow = widNow * preFac * ps * sqrt(lambda);
432 if(DEBUG) {
433 cout<<"Gluino:: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
434 cout<<scientific<<widNow<<endl;
435 }
436 return;
437 }
438}
439
440//==========================================================================
441
442// Class ResonanceNeut
443// Derived class for Neutralino Resonances
444
445//--------------------------------------------------------------------------
446
447
448void ResonanceNeut::initConstants(){
449
450 alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
451 s2W = coupSUSYPtr->sin2W;
452}
453
454//--------------------------------------------------------------------------
455
456// Calculate various common prefactors for the current mass.
457void ResonanceNeut::calcPreFac(bool){
458
459 // Common coupling factors.
460 preFac = alpEM / (8.0 * s2W * pow(mHat,3));
461 return;
462}
463
464//--------------------------------------------------------------------------
465
466// Calculate width for currently considered channel.
467void ResonanceNeut::calcWidth(bool){
468
469 widNow = 0.0;
470
471 if(ps ==0.) return;
472 kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
473 kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
474 + pow2(mHat) * pow2(mf2) + pow2(mf1) * pow2(mf2) - pow2(mHat) * pow2(mf1);
475 lambda = lam(mHat*mHat, mf1*mf1, mf2*mf2);
476
477 // Stable lightest neutralino
478 if(idRes == 1000022) return;
479
480 double fac = 0.0;
481 int iNeut1 = typeNeut(idRes);
482 int iNeut2 = typeNeut(id1Abs);
483 int iChar1 = typeChar(id1Abs);
484
485 if(iNeut2>0 && id2Abs == 23){
486 // ~chi0_i -> chi0_j + Z
487 fac = kinFac2 * (norm(coupSUSYPtr->OLpp[iNeut1][iNeut2]) + norm(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
488 fac -= 12.0 * mHat * mf1 * pow2(mf2) *
489 real(coupSUSYPtr->OLpp[iNeut1][iNeut2] * conj(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
490 fac /= pow2(mf2) * (1.0 - s2W);
491 }
492 else if(iChar1>0 && id2Abs==24){
493 // ~chi0_i -> chi+_j + W- (or c.c.)
494
495 fac = kinFac2 * (norm(coupSUSYPtr->OL[iNeut1][iChar1]) + norm(coupSUSYPtr->OR[iNeut1][iChar1]));
496 fac -= 12.0 * mHat * mf1 * pow2(mf2) *
497 real(coupSUSYPtr->OL[iNeut1][iChar1] * conj(coupSUSYPtr->OR[iNeut1][iChar1]));
498 fac /= pow2(mf2);
499 }
500 else if(id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
501 // ~chi0_k -> ~q + q
502 bool idown = (id1Abs%2 == 1);
503 int iq = (id2Abs + 1 )/ 2;
504 int isq = (abs(idRes)/1000000 == 2) ?
505 (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
506
507 if(idown){
508 fac = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][iNeut1]) + norm(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
509 fac += 4.0 * mHat * mf2 *
510 real(coupSUSYPtr->LsddX[isq][iq][iNeut1] * conj(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
511 }
512 else{
513 fac = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][iNeut1]) + norm(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
514 fac += 4.0 * mHat * mf2 * sqrt(lambda) *
515 real(coupSUSYPtr->LsuuX[isq][iq][iNeut1] * conj(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
516 }
517 fac *= 2.0/(1 - s2W);
518 }
519
520 widNow = fac * preFac * ps * sqrt(lambda);
521 if(DEBUG) {
522 cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
523 cout<<scientific<<widNow<<endl;
524 }
525 return;
526}
527
528//==========================================================================
529
530// Class ResonanceChar
531// Derived class for Neutralino Resonances
532
533//--------------------------------------------------------------------------
534
535
536void ResonanceChar::initConstants(){
537
538 alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
539 s2W = coupSUSYPtr->sin2W;
540 return;
541}
542
543//--------------------------------------------------------------------------
544
545// Calculate various common prefactors for the current mass.
546void ResonanceChar::calcPreFac(bool){
547
548 preFac = alpEM / (8.0 * s2W * pow(mHat,3));
549 return;
550}
551
552//--------------------------------------------------------------------------
553
554// Calculate width for currently considered channel.
555void ResonanceChar::calcWidth(bool){
556
557 widNow = 0.0;
558 if(ps == 0.) return;
559
560 double fac = 0.0;
561 kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
562 kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
563 + pow2(mHat) * pow2(mf2) + pow2(mf1) * pow2(mf2) - pow2(mHat) * pow2(mf1);
564 lambda = lam(mHat*mHat , mf1*mf1 , mf2*mf2);
565
566 int idChar1 = typeChar(idRes);
567 int idChar2 = typeChar(id1Abs);
568 int idNeut1 = typeNeut(id1Abs);
569
570 if(idChar2>0 && id2Abs == 23){
571 // ~chi_i -> chi_j + Z
572 fac = kinFac2 * (norm(coupSUSYPtr->OLp[idChar1][idChar2])
573 + norm(coupSUSYPtr->ORp[idChar1][idChar2]));
574 fac -= 12.0 * mHat * mf1 * pow2(mf2) *
575 real(coupSUSYPtr->OLp[idChar1][idChar2]
576 * conj(coupSUSYPtr->ORp[idChar1][idChar2]));
577 fac /= pow2(mf2) * (1.0 - s2W);
578 }
579 else if(idNeut1>0 && id2Abs==24){
580 // ~chi_i -> chi0_j + W- (or c.c.)
581
582 fac = kinFac2 * (norm(coupSUSYPtr->OL[idNeut1][idChar1]) + norm(coupSUSYPtr->OR[idNeut1][idChar1]));
583 fac -= 12.0 * mHat * mf1 * pow2(mf2) *
584 real(coupSUSYPtr->OL[idNeut1][idChar1] * conj(coupSUSYPtr->OR[idNeut1][idChar1]));
585 fac /= pow2(mf2);
586 }
587 else if(id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
588 // ~chi0_k -> ~q + q
589 bool idown = (id1Abs%2 == 1);
590 int iq = (id2Abs + 1 )/ 2;
591 int isq = (abs(idRes)/1000000 == 2) ?
592 (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
593
594 if(idown){
595 fac = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][idChar1]) + norm(coupSUSYPtr->RsduX[isq][iq][idChar1]));
596 fac += 4.0 * mHat * mf2 *
597 real(coupSUSYPtr->LsduX[isq][iq][idChar1] * conj(coupSUSYPtr->RsduX[isq][iq][idChar1]));
598 }
599 else{
600 fac = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][idChar1]) + norm(coupSUSYPtr->RsudX[isq][iq][idChar1]));
601 fac += 4.0 * mHat * mf2 *
602 real(coupSUSYPtr->LsudX[isq][iq][idChar1] * conj(coupSUSYPtr->RsudX[isq][iq][idChar1]));
603 }
604 fac *= 2.0/(1 - s2W);
605 }
606
607 widNow = fac * preFac * ps * sqrt(lambda) ;
608 if(DEBUG) {
609 cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
610 cout<<scientific<<widNow<<endl;
611 }
612 return;
613}
614
615
616//==========================================================================
617
618//Return neutralino code; zero if not a neutralino
619
620int SUSYResonanceWidths::typeNeut(int idPDG) {
621 int type = 0;
622 int idAbs = abs(idPDG);
623 if(idAbs == 1000022) type = 1;
624 else if(idAbs == 1000023) type = 2;
625 else if(idAbs == 1000025) type = 3;
626 else if(idAbs == 1000035) type = 4;
627 else if(coupSUSYPtr->isNMSSM && idAbs == 1000045) type = 5;
628 return type;
629
630}
631
632
633//--------------------------------------------------------------------------
634
635//Check whether particle is a Chargino
636
637int SUSYResonanceWidths::typeChar(int idPDG) {
638 int type = 0;
639 if(abs(idPDG) == 1000024) type = 1;
640 else if (abs(idPDG) == 1000037)type = 2;
641 return type;
642}
643
644//--------------------------------------------------------------------------
645
646// Function for Kallen function
647
648double SUSYResonanceWidths::lam(double x, double y, double z){
649
650 double val = x*x + y*y + z*z - 2.0* (x*y + y*z + z*x);
651 return val;
652
653}
654
655
656//--------------------------------------------------------------------------
657
658
659
660} //end namespace Pythia8
661
662