StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SusyResonanceWidths.cc
1 // SusyResonanceWidths.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand
3 // Authors: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
5 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 
7 // Function definitions (not found in the header) for the
8 // SUSY Resonance widths classes.
9 
10 #include "Pythia8/SusyResonanceWidths.h"
11 #include "Pythia8/SusyWidthFunctions.h"
12 #include "Pythia8/SusyCouplings.h"
13 #include "Pythia8/ParticleData.h"
14 #include "Pythia8/PythiaComplex.h"
15 
16 namespace Pythia8 {
17 
18 //==========================================================================
19 
20 // The SUSYResonanceWidths Class
21 // Derived class for SUSY resonances
22 
23 const bool SUSYResonanceWidths::DBSUSY = false;
24 
25 //--------------------------------------------------------------------------
26 
27 bool SUSYResonanceWidths::initBSM(){
28 
29  if (coupSUSYPtr->isSUSY) {
30  return true;
31  }
32 
33  return false;
34 }
35 
36 //--------------------------------------------------------------------------
37 
38 bool SUSYResonanceWidths::allowCalc(){
39 
40  // Check if decay calculations at all possible
41  if ( !coupSUSYPtr->isSUSY ) return false;
42  if ( (idRes == 45 || idRes == 46 || idRes == 1000045)
43  && !coupSUSYPtr->isNMSSM ) return false;
44 
45  if (settingsPtr->flag("SLHA:useDecayTable") ) {
46 
47  // Next check if decay table was read in via SLHA and takes precedence
48  for ( int iDec = 0; iDec < int((coupSUSYPtr->slhaPtr)->decays.size());
49  ++iDec)
50  if ( (coupSUSYPtr->slhaPtr)->decays[iDec].getId() == abs(idRes) ) {
51  if (DBSUSY) cout<<"Using external decay table for:"<<idRes<<endl;
52  return false;
53  }
54 
55  }
56 
57  // Else we should do the calculation; set available channels
58  bool done = getChannels(idRes);
59  stringstream idStream;
60  idStream << "ID = " << idRes ;
61  if (!done) infoPtr->errorMsg("Error in SusyResonanceWidths::allowcalc: "
62  "unable to reset decay table.", idStream.str(), true);
63  return done;
64 }
65 
66 //==========================================================================
67 
68 // The ResonanceSquark class
69 // Derived class for Squark resonances
70 
71 //--------------------------------------------------------------------------
72 
73 // Set up channels
74 
75 bool ResonanceSquark::getChannels(int idPDG){
76 
77  idPDG = abs(idPDG);
78 
79  int ksusy = 1000000;
80  if (idPDG < ksusy) return false;
81  if(idPDG % ksusy >= 7 || idPDG % ksusy < 1) return false;
82 
83  ParticleDataEntry* squarkEntryPtr
84  = particleDataPtr->particleDataEntryPtr(idPDG);
85 
86  // Delete any decay channels read
87  squarkEntryPtr->clearChannels();
88 
89  if(idPDG % 2 == 0) { // up-type squarks
90 
91  /*
92  // Gravitino
93  squarkEntryPtr->addChannel(1, 0.0, 103, 1000039, 2);
94  squarkEntryPtr->addChannel(1, 0.0, 0, 1000039, 4);
95  squarkEntryPtr->addChannel(1, 0.0, 0, 1000039, 6);
96  */
97 
98  // Gaugino - quark
99  squarkEntryPtr->addChannel(1, 0.0, 0, 1000024, 3);
100  squarkEntryPtr->addChannel(1, 0.0, 0, 1000024, 5);
101  squarkEntryPtr->addChannel(1, 0.0, 0, 1000037, 1);
102  squarkEntryPtr->addChannel(1, 0.0, 0, 1000037, 3);
103  squarkEntryPtr->addChannel(1, 0.0, 0, 1000037, 5);
104  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 2);
105  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 4);
106  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 6);
107  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 2);
108  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 4);
109  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 6);
110  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 2);
111  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 4);
112  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 6);
113  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 2);
114  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 4);
115  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 6);
116 
117  // Squark - Gauge boson
118  squarkEntryPtr->addChannel(1, 0.0, 0, 1000001, -24);
119  squarkEntryPtr->addChannel(1, 0.0, 0, 1000003, -24);
120  squarkEntryPtr->addChannel(1, 0.0, 0, 1000005, -24);
121  squarkEntryPtr->addChannel(1, 0.0, 0, 2000001, -24);
122  squarkEntryPtr->addChannel(1, 0.0, 0, 2000003, -24);
123  squarkEntryPtr->addChannel(1, 0.0, 0, 2000005, -24);
124  squarkEntryPtr->addChannel(1, 0.0, 0, 1000001, -37);
125  squarkEntryPtr->addChannel(1, 0.0, 0, 1000003, -37);
126  squarkEntryPtr->addChannel(1, 0.0, 0, 1000005, -37);
127  squarkEntryPtr->addChannel(1, 0.0, 0, 2000001, -37);
128  squarkEntryPtr->addChannel(1, 0.0, 0, 2000003, -37);
129  squarkEntryPtr->addChannel(1, 0.0, 0, 2000005, -37);
130 
131  // Gluino - quark
132  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 2);
133  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 4);
134  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 6);
135 
136  // lepton - quark via LQD
137  squarkEntryPtr->addChannel(1, 0.0, 0, -11, 1);
138  squarkEntryPtr->addChannel(1, 0.0, 0, -11, 3);
139  squarkEntryPtr->addChannel(1, 0.0, 0, -11, 5);
140  squarkEntryPtr->addChannel(1, 0.0, 0, -13, 1);
141  squarkEntryPtr->addChannel(1, 0.0, 0, -13, 3);
142  squarkEntryPtr->addChannel(1, 0.0, 0, -13, 5);
143  squarkEntryPtr->addChannel(1, 0.0, 0, -15, 1);
144  squarkEntryPtr->addChannel(1, 0.0, 0, -15, 3);
145  squarkEntryPtr->addChannel(1, 0.0, 0, -15, 5);
146 
147  // quark - quark via UDD
148  squarkEntryPtr->addChannel(1, 0.0, 0, -1 ,-3);
149  squarkEntryPtr->addChannel(1, 0.0, 0, -1 ,-5);
150  squarkEntryPtr->addChannel(1, 0.0, 0, -3 ,-5);
151 
152 
153  } else { //down-type squarks
154 
155  // Gaugino - quark
156  squarkEntryPtr->addChannel(1, 0.0, 0, -1000024, 2);
157  squarkEntryPtr->addChannel(1, 0.0, 0, -1000037, 2);
158  squarkEntryPtr->addChannel(1, 0.0, 0, -1000024, 4);
159  squarkEntryPtr->addChannel(1, 0.0, 0, -1000037, 4);
160  squarkEntryPtr->addChannel(1, 0.0, 0, -1000024, 6);
161  squarkEntryPtr->addChannel(1, 0.0, 0, -1000037, 6);
162  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 1);
163  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 3);
164  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 5);
165  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 1);
166  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 3);
167  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 5);
168  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 1);
169  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 3);
170  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 5);
171  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 1);
172  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 3);
173  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 5);
174 
175  // Squark - Gauge boson
176  squarkEntryPtr->addChannel(1, 0.0, 0, 1000002, -24);
177  squarkEntryPtr->addChannel(1, 0.0, 0, 1000004, -24);
178  squarkEntryPtr->addChannel(1, 0.0, 0, 1000006, -24);
179  squarkEntryPtr->addChannel(1, 0.0, 0, 2000002, -24);
180  squarkEntryPtr->addChannel(1, 0.0, 0, 2000004, -24);
181  squarkEntryPtr->addChannel(1, 0.0, 0, 2000006, -24);
182  squarkEntryPtr->addChannel(1, 0.0, 0, 1000002, -37);
183  squarkEntryPtr->addChannel(1, 0.0, 0, 1000004, -37);
184  squarkEntryPtr->addChannel(1, 0.0, 0, 1000006, -37);
185  squarkEntryPtr->addChannel(1, 0.0, 0, 2000002, -37);
186  squarkEntryPtr->addChannel(1, 0.0, 0, 2000004, -37);
187  squarkEntryPtr->addChannel(1, 0.0, 0, 2000006, -37);
188 
189  // Gluino - quark
190  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 1);
191  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 2);
192  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 5);
193 
194  // lepton - quark via LQD
195  squarkEntryPtr->addChannel(1, 0.0, 0, -12, 1);
196  squarkEntryPtr->addChannel(1, 0.0, 0, -12, 3);
197  squarkEntryPtr->addChannel(1, 0.0, 0, -12, 5);
198  squarkEntryPtr->addChannel(1, 0.0, 0, -14, 1);
199  squarkEntryPtr->addChannel(1, 0.0, 0, -14, 3);
200  squarkEntryPtr->addChannel(1, 0.0, 0, -14, 5);
201  squarkEntryPtr->addChannel(1, 0.0, 0, -16, 1);
202  squarkEntryPtr->addChannel(1, 0.0, 0, -16, 3);
203  squarkEntryPtr->addChannel(1, 0.0, 0, -16, 5);
204  squarkEntryPtr->addChannel(1, 0.0, 0, 12 ,1);
205  squarkEntryPtr->addChannel(1, 0.0, 0, 11 ,2);
206  squarkEntryPtr->addChannel(1, 0.0, 0, 12, 3);
207  squarkEntryPtr->addChannel(1, 0.0, 0, 11, 4);
208  squarkEntryPtr->addChannel(1, 0.0, 0, 12, 5);
209  squarkEntryPtr->addChannel(1, 0.0, 0, 11, 6);
210  squarkEntryPtr->addChannel(1, 0.0, 0, 14, 1);
211  squarkEntryPtr->addChannel(1, 0.0, 0, 13, 2);
212  squarkEntryPtr->addChannel(1, 0.0, 0, 14, 3);
213  squarkEntryPtr->addChannel(1, 0.0, 0, 13, 4);
214  squarkEntryPtr->addChannel(1, 0.0, 0, 14, 5);
215  squarkEntryPtr->addChannel(1, 0.0, 0, 13, 6);
216  squarkEntryPtr->addChannel(1, 0.0, 0, 16, 1);
217  squarkEntryPtr->addChannel(1, 0.0, 0, 15, 2);
218  squarkEntryPtr->addChannel(1, 0.0, 0, 16, 3);
219  squarkEntryPtr->addChannel(1, 0.0, 0, 15, 4);
220  squarkEntryPtr->addChannel(1, 0.0, 0, 16, 5);
221  squarkEntryPtr->addChannel(1, 0.0, 0, 15, 6);
222 
223 
224  // quark - quark via UDD
225  squarkEntryPtr->addChannel(1, 0.0, 0, -2, -1);
226  squarkEntryPtr->addChannel(1, 0.0, 0, -2, -3);
227  squarkEntryPtr->addChannel(1, 0.0, 0, -2, -5);
228  squarkEntryPtr->addChannel(1, 0.0, 0, -4, -1);
229  squarkEntryPtr->addChannel(1, 0.0, 0, -4, -3);
230  squarkEntryPtr->addChannel(1, 0.0, 0, -4, -5);
231  squarkEntryPtr->addChannel(1, 0.0, 0, -6, -1);
232  squarkEntryPtr->addChannel(1, 0.0, 0, -6, -3);
233  squarkEntryPtr->addChannel(1, 0.0, 0, -6, -5);
234  }
235 
236  return true;
237 
238 }
239 
240 //--------------------------------------------------------------------------
241 
242 // Initialize constants.
243 
244 void ResonanceSquark::initConstants() {
245 
246  // Locally stored properties and couplings.
247  s2W = coupSUSYPtr->sin2W;
248 
249 }
250 
251 //--------------------------------------------------------------------------
252 
253 // Calculate various common prefactors for the current mass.
254 
255 void ResonanceSquark::calcPreFac(bool) {
256 
257  // Common coupling factors.
258  alpS = coupSMPtr->alphaS(mHat * mHat );
259  alpEM = coupSMPtr->alphaEM(mHat * mHat);
260  preFac = 1.0 / (s2W * pow(mHat,3));
261  ps *= mHat * mHat;
262 
263 }
264 
265 //--------------------------------------------------------------------------
266 
267 // Calculate width for currently considered channel.
268 
269 void ResonanceSquark::calcWidth(bool) {
270 
271  // Squark type -- in u_i/d_i and generation
272  int ksusy = 1000000;
273  bool idown = (abs(idRes)%2 == 0 ? false : true);
274  int isq = (abs(idRes)/ksusy == 2) ?
275  (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
276  // int isqgen = (abs(idRes)%10 + 1)/2;
277 
278  // Check that mass is above threshold.
279  if (ps == 0.) return;
280  else{
281  // Two-body decays
282  kinFac = (mHat * mHat - mf1 * mf1 - mf2 * mf2);
283 
284  double fac = 0.0 , wid = 0.0;
285 
286  //RPV decays
287  //Case 1a: UDD-type
288  if (id1Abs < 7 && id2Abs < 7){
289 
290  // Quark generations
291  int iq1 = (id1Abs + 1)/2;
292  int iq2 = (id2Abs + 1)/2;
293 
294  // Check for RPV UDD couplings
295  if (!coupSUSYPtr->isUDD) {widNow = 0; return;}
296 
297  // ~q -> q_i + q_j
298 
299  fac = 2.0 * kinFac / (16.0 * M_PI * pow(mHat,3));
300  wid = 0.0;
301  if (idown) {
302  if ((id1Abs+id2Abs)%2 == 1){
303  if (id1Abs%2==1)
304  for (int isq2 = 1; isq2 < 4; isq2++)
305  wid += norm(coupSUSYPtr->rvUDD[iq2][iq1][isq2]
306  * coupSUSYPtr->Rdsq[isq][isq2+3]);
307  else
308  for (int isq2 = 1; isq2 < 4; isq2++)
309  wid += norm(coupSUSYPtr->rvUDD[iq1][iq2][isq2]
310  * coupSUSYPtr->Rdsq[isq][isq2+3]);
311  }
312  }
313  else {
314  if ((id1Abs+id2Abs)%2 != 0) widNow = 0.0;
315  else
316  for (int isq2 = 1; isq2 < 4; isq2++)
317  wid += norm(coupSUSYPtr->rvUDD[isq2][iq1][iq2]
318  * coupSUSYPtr->Rusq[isq][isq2+3]);
319  }
320  }
321 
322  //Case 1b: LQD-type
323  else if (id1Abs < 17 && id2Abs < 7){
324  if (!coupSUSYPtr->isLQD) {widNow = 0; return;}
325 
326  int ilep = (id1Abs - 9)/2;
327  int iq = (id2Abs + 1)/2;
328 
329  fac = kinFac / (16.0 * M_PI * pow(mHat,3));
330  wid = 0.0;
331  if (idown){
332  if (iq%2 == 0){
333  // q is up-type; ~q is right-handed down type
334  for (int isq2=1; isq2<3; isq2++)
335  wid += norm(coupSUSYPtr->Rdsq[isq][isq2+3]
336  * coupSUSYPtr->rvLQD[ilep][iq][isq2]);
337  }else{
338  //q is down type; ~q left-handed down-type
339  for (int isq2=1; isq2<3; isq2++)
340  wid += norm(coupSUSYPtr->Rdsq[isq][isq2]
341  * coupSUSYPtr->rvLQD[ilep][isq2][isq2]);
342  }
343  }
344  else{
345  if (iq%2 == 0) {widNow = 0.0; return;}
346  // q is down type; ~q is left-handed up-type
347  for (int isq2=1; isq2<3; isq2++)
348  wid += norm(coupSUSYPtr->Rusq[isq][isq2]
349  * coupSUSYPtr->rvLQD[ilep][isq2][iq]);
350  }
351  }
352 
353  //Case 2: quark + gaugino
354  else if (id1Abs > ksusy && id2Abs < 7) {
355 
356  int iq = (id2Abs + 1)/2;
357 
358  // ~q -> ~g + q
359  if (id1Abs == 1000021 && idRes%10 == id2Abs) {
360  // Removed factor of s2W in denominator: strong process -- no EW
361  fac = 2.0 * alpS / (3.0 * pow3(mHat));
362  if (idown)
363  wid = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])
364  + norm(coupSUSYPtr->RsddG[isq][iq]))
365  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddG[isq][iq]
366  * conj(coupSUSYPtr->RsddG[isq][iq]));
367  else
368  wid = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])
369  + norm(coupSUSYPtr->RsuuG[isq][iq]))
370  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]
371  * conj(coupSUSYPtr->RsuuG[isq][iq]));
372  }
373  else
374  for (int i=1; i<6 ; i++){
375  // ~q -> ~chi0 + q
376  if (coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
377  fac = alpEM * preFac / (2.0 * (1 - s2W));
378  if (idown)
379  wid = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][i])
380  + norm(coupSUSYPtr->RsddX[isq][iq][i]))
381  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddX[isq][iq][i]
382  * conj(coupSUSYPtr->RsddX[isq][iq][i]));
383  else
384  wid = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][i])
385  + norm(coupSUSYPtr->RsuuX[isq][iq][i]))
386  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][i]
387  * conj(coupSUSYPtr->RsuuX[isq][iq][i]));
388  }
389 
390  // ~q -> chi- + q
391  else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs
392  && idRes%2 != id2Abs%2){
393 
394  fac = alpEM * preFac / (4.0 * (1 - s2W));
395  if (idown)
396  wid = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][i])
397  + norm(coupSUSYPtr->RsduX[isq][iq][i]))
398  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsduX[isq][iq][i]
399  * conj(coupSUSYPtr->RsduX[isq][iq][i]));
400  else
401  wid = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][i])
402  + norm(coupSUSYPtr->RsudX[isq][iq][i]))
403  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsudX[isq][iq][i]
404  * conj(coupSUSYPtr->RsudX[isq][iq][i]));
405  }
406  }
407  }
408 
409  //Case 3: ~q_i -> ~q_j + Z/W
410  else if (id1Abs > ksusy && id1Abs%100 < 7
411  && (id2Abs == 23 || id2Abs == 24)){
412 
413  // factor of lambda^(3/2) = ps^(3/2) ;
414  fac = alpEM * preFac/(16.0 * pow2(particleDataPtr->m0(id2Abs))
415  * (1.0 - s2W)) * pow2(ps) ;
416 
417  int isq2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
418  isq2 = min(6, isq2);
419 
420  if (id2Abs == 23 && id1Abs%2 == idRes%2){
421  if (idown)
422  wid = norm(coupSUSYPtr->LsdsdZ[isq][isq2]
423  + coupSUSYPtr->RsdsdZ[isq][isq2]);
424  else
425  wid = norm(coupSUSYPtr->LsusuZ[isq][isq2]
426  + coupSUSYPtr->RsusuZ[isq][isq2]);
427  }
428  else if (id2Abs == 24 && id1Abs%2 != idRes%2){
429  if (idown)
430  wid = norm(coupSUSYPtr->LsusdW[isq2][isq]);
431  else
432  wid = norm(coupSUSYPtr->LsusdW[isq][isq2]);
433  }
434  }
435 
436  // TODO: Case ~q_i -> ~q_j + h/H
437  widNow = fac * wid * ps * pow2(mHat);
438  if (DBSUSY) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs
439  <<" Width: "<<widNow<<endl;
440  return;
441  }
442 
443 }
444 
445 //==========================================================================
446 
447 // The ResonanceGluino class
448 // Derived class for Gluino resonances
449 
450 //--------------------------------------------------------------------------
451 
452 // Set up Channels
453 
454 bool ResonanceGluino::getChannels(int idPDG){
455 
456  idPDG = abs(idPDG);
457  if (idPDG != 1000021) return false;
458 
459  ParticleDataEntry* gluinoEntryPtr
460  = particleDataPtr->particleDataEntryPtr(idPDG);
461 
462  // Delete any decay channels read
463  gluinoEntryPtr->clearChannels();
464 
465  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000001, -1);
466  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000001, 1);
467  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000001 ,-3);
468  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000001, 3);
469  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000001 ,-5);
470  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000001, 5);
471  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000001 ,-1);
472  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000001, 1);
473  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000001 ,-3);
474  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000001, 3);
475  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000001 ,-5);
476  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000001, 5);
477  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000002 ,-2);
478  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000002, 2);
479  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000002 ,-4);
480  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000002, 4);
481  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000002 ,-6);
482  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000002, 6);
483  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000002 ,-2);
484  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000002, 2);
485  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000002 ,-4);
486  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000002, 4);
487  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000002 ,-6);
488  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000002, 6);
489  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000003 ,-1);
490  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000003, 1);
491  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000003 ,-3);
492  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000003, 3);
493  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000003 ,-5);
494  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000003, 5);
495  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000003 ,-1);
496  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000003, 1);
497  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000003 ,-3);
498  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000003, 3);
499  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000003 ,-5);
500  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000003, 5);
501  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000004 ,-2);
502  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000004, 2);
503  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000004 ,-4);
504  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000004, 4);
505  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000004 ,-6);
506  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000004, 6);
507  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000004 ,-2);
508  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000004, 2);
509  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000004 ,-4);
510  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000004, 4);
511  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000004 ,-6);
512  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000004, 6);
513  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000005 ,-1);
514  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000005, 1);
515  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000005 ,-3);
516  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000005, 3);
517  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000005 ,-5);
518  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000005, 5);
519  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000005 ,-1);
520  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000005, 1);
521  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000005 ,-3);
522  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000005, 3);
523  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000005 ,-5);
524  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000005, 5);
525  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000006 ,-6);
526  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000006, 6);
527  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000006 ,-2);
528  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000006, 2);
529  gluinoEntryPtr->addChannel(1, 0.0, 0, 1000006 ,-4);
530  gluinoEntryPtr->addChannel(1, 0.0, 0, -1000006, 4);
531  gluinoEntryPtr->addChannel(1, 0.0, 0, 2000006 ,-6);
532  gluinoEntryPtr->addChannel(1, 0.0, 0, -2000006, 6);
533 
534  return true;
535 }
536 
537 //--------------------------------------------------------------------------
538 
539 // Initialize constants.
540 
541 void ResonanceGluino::initConstants() {
542  return;
543 }
544 
545 //--------------------------------------------------------------------------
546 
547 // Calculate various common prefactors for the current mass.
548 
549 void ResonanceGluino::calcPreFac(bool) {
550 
551  // Common coupling factors.
552  alpS = coupSMPtr->alphaS(mHat * mHat);
553  preFac = alpS /( 8.0 * pow(mHat,3));
554  return;
555 
556 }
557 
558 
559 //--------------------------------------------------------------------------
560 
561 // Calculate width for currently considered channel.
562 
563 void ResonanceGluino::calcWidth(bool) {
564 
565  widNow = 0.0;
566  if (ps == 0.) return;
567  kinFac = (mHat * mHat - mf1 * mf1 + mf2 * mf2);
568 
569  if (id1Abs > 1000000 && (id1Abs % 100) < 7 && id2Abs < 7) {
570 
571  int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
572  : (abs(id1Abs)%10+1)/2;
573  bool idown = id2Abs%2;
574  int iq = (id2Abs + 1)/2;
575 
576  // ~g -> ~q + q
577  if (idown){
578  widNow = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])
579  + norm(coupSUSYPtr->RsddG[isq][iq]))
580  + 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddG[isq][iq]
581  * conj(coupSUSYPtr->RsddG[isq][iq]));
582  }
583  else{
584  widNow = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])
585  + norm(coupSUSYPtr->RsuuG[isq][iq]))
586  + 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]
587  * conj(coupSUSYPtr->RsuuG[isq][iq]));
588  }
589  widNow = widNow * preFac * ps * pow2(mHat);
590  if (DBSUSY) {
591  cout<<"Gluino:: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
592  cout<<scientific<<widNow<<endl;
593  }
594  return;
595  }
596 }
597 
598 //==========================================================================
599 
600 // Class ResonanceNeut
601 // Derived class for Neutralino Resonances
602 //
603 //--------------------------------------------------------------------------
604 
605 // Set up Channels
606 
607 bool ResonanceNeut::getChannels(int idPDG){
608 
609  idPDG = abs(idPDG);
610 
611  int iNeut = coupSUSYPtr->typeNeut(idPDG);
612  if (iNeut < 1) return false;
613 
614  ParticleDataEntry* neutEntryPtr
615  = particleDataPtr->particleDataEntryPtr(idPDG);
616 
617  // Delete any decay channels read
618  neutEntryPtr->clearChannels();
619 
620  // RPV channels
621 
622  neutEntryPtr->addChannel(1, 0.0, 0, -12, -13, 11);
623  neutEntryPtr->addChannel(1, 0.0, 0, 12, 13, -11);
624  neutEntryPtr->addChannel(1, 0.0, 0, -12, -13, 13);
625  neutEntryPtr->addChannel(1, 0.0, 0, 12, 13, -13);
626  neutEntryPtr->addChannel(1, 0.0, 0, -12, -13, 15);
627  neutEntryPtr->addChannel(1, 0.0, 0, 12, 13, -15);
628  neutEntryPtr->addChannel(1, 0.0, 0, -12, -15, 11);
629  neutEntryPtr->addChannel(1, 0.0, 0, 12, 15, -11);
630  neutEntryPtr->addChannel(1, 0.0, 0, -12, -15, 13);
631  neutEntryPtr->addChannel(1, 0.0, 0, 12, 15, -13);
632  neutEntryPtr->addChannel(1, 0.0, 0, -12, -15, 15);
633  neutEntryPtr->addChannel(1, 0.0, 0, 12, 15, -15);
634  neutEntryPtr->addChannel(1, 0.0, 0, -14, -11, 11);
635  neutEntryPtr->addChannel(1, 0.0, 0, 14, 11, -11);
636  neutEntryPtr->addChannel(1, 0.0, 0, -14, -11, 13);
637  neutEntryPtr->addChannel(1, 0.0, 0, 14, 11, -13);
638  neutEntryPtr->addChannel(1, 0.0, 0, -14, -11, 15);
639  neutEntryPtr->addChannel(1, 0.0, 0, 14, 11, -15);
640  neutEntryPtr->addChannel(1, 0.0, 0, -14, -15, 11);
641  neutEntryPtr->addChannel(1, 0.0, 0, 14, 15, -11);
642  neutEntryPtr->addChannel(1, 0.0, 0, -14, -15, 13);
643  neutEntryPtr->addChannel(1, 0.0, 0, 14, 15, -13);
644  neutEntryPtr->addChannel(1, 0.0, 0, -14, -15, 15);
645  neutEntryPtr->addChannel(1, 0.0, 0, 14, 15, -15);
646  neutEntryPtr->addChannel(1, 0.0, 0, -16, -11, 11);
647  neutEntryPtr->addChannel(1, 0.0, 0, 16, 11, -11);
648  neutEntryPtr->addChannel(1, 0.0, 0, -16, -11, 13);
649  neutEntryPtr->addChannel(1, 0.0, 0, 16, 11, -13);
650  neutEntryPtr->addChannel(1, 0.0, 0, -16, -11, 15);
651  neutEntryPtr->addChannel(1, 0.0, 0, 16, 11, -15);
652  neutEntryPtr->addChannel(1, 0.0, 0, -16, -13, 11);
653  neutEntryPtr->addChannel(1, 0.0, 0, 16, 13, -11);
654  neutEntryPtr->addChannel(1, 0.0, 0, -16, -13, 13);
655  neutEntryPtr->addChannel(1, 0.0, 0, 16, 13, -13);
656  neutEntryPtr->addChannel(1, 0.0, 0, -16, -13, 15);
657  neutEntryPtr->addChannel(1, 0.0, 0, 16, 13, -15);
658  neutEntryPtr->addChannel(1, 0.0, 0, -12, -1, 1);
659  neutEntryPtr->addChannel(1, 0.0, 0, 12, 1, -1);
660  neutEntryPtr->addChannel(1, 0.0, 0, -11, -2, 1);
661  neutEntryPtr->addChannel(1, 0.0, 0, 11, 2, -1);
662  neutEntryPtr->addChannel(1, 0.0, 0, -12, -1, 3);
663  neutEntryPtr->addChannel(1, 0.0, 0, 12, 1, -3);
664  neutEntryPtr->addChannel(1, 0.0, 0, -11, -2, 3);
665  neutEntryPtr->addChannel(1, 0.0, 0, 11, 2, -3);
666  neutEntryPtr->addChannel(1, 0.0, 0, -12, -1, 5);
667  neutEntryPtr->addChannel(1, 0.0, 0, 12, 1, -5);
668  neutEntryPtr->addChannel(1, 0.0, 0, -11, -2, 5);
669  neutEntryPtr->addChannel(1, 0.0, 0, 11, 2, -5);
670  neutEntryPtr->addChannel(1, 0.0, 0, -12, -3, 1);
671  neutEntryPtr->addChannel(1, 0.0, 0, 12, 3, -1);
672  neutEntryPtr->addChannel(1, 0.0, 0, -11, -4, 1);
673  neutEntryPtr->addChannel(1, 0.0, 0, 11, 4, -1);
674  neutEntryPtr->addChannel(1, 0.0, 0, -12, -3, 3);
675  neutEntryPtr->addChannel(1, 0.0, 0, 12, 3, -3);
676  neutEntryPtr->addChannel(1, 0.0, 0, -11, -4, 3);
677  neutEntryPtr->addChannel(1, 0.0, 0, 11, 4, -3);
678  neutEntryPtr->addChannel(1, 0.0, 0, -12, -3, 5);
679  neutEntryPtr->addChannel(1, 0.0, 0, 12, 3, -5);
680  neutEntryPtr->addChannel(1, 0.0, 0, -11, -4, 5);
681  neutEntryPtr->addChannel(1, 0.0, 0, 11, 4, -5);
682  neutEntryPtr->addChannel(1, 0.0, 0, -12, -5, 1);
683  neutEntryPtr->addChannel(1, 0.0, 0, 12, 5, -1);
684  neutEntryPtr->addChannel(1, 0.0, 0, -11, -6, 1);
685  neutEntryPtr->addChannel(1, 0.0, 0, 11, 6, -1);
686  neutEntryPtr->addChannel(1, 0.0, 0, -12, -5, 3);
687  neutEntryPtr->addChannel(1, 0.0, 0, 12, 5, -3);
688  neutEntryPtr->addChannel(1, 0.0, 0, -11, -6, 3);
689  neutEntryPtr->addChannel(1, 0.0, 0, 11, 6, -3);
690  neutEntryPtr->addChannel(1, 0.0, 0, 12 ,-5 ,5);
691  neutEntryPtr->addChannel(1, 0.0, 0, 12, 5, -5);
692  neutEntryPtr->addChannel(1, 0.0, 0, -11, -6, 5);
693  neutEntryPtr->addChannel(1, 0.0, 0, 11, 6, -5);
694  neutEntryPtr->addChannel(1, 0.0, 0, -14, -1, 1);
695  neutEntryPtr->addChannel(1, 0.0, 0, 14, 1, -1);
696  neutEntryPtr->addChannel(1, 0.0, 0, -13, -2, 1);
697  neutEntryPtr->addChannel(1, 0.0, 0, 13, 2, -1);
698  neutEntryPtr->addChannel(1, 0.0, 0, -14, -1, 3);
699  neutEntryPtr->addChannel(1, 0.0, 0, 14, 1, -3);
700  neutEntryPtr->addChannel(1, 0.0, 0, -13, -2, 3);
701  neutEntryPtr->addChannel(1, 0.0, 0, 13, 2, -3);
702  neutEntryPtr->addChannel(1, 0.0, 0, -14, -1, 5);
703  neutEntryPtr->addChannel(1, 0.0, 0, 14, 1, -5);
704  neutEntryPtr->addChannel(1, 0.0, 0, -13, -2, 5);
705  neutEntryPtr->addChannel(1, 0.0, 0, 13, 2, -5);
706  neutEntryPtr->addChannel(1, 0.0, 0, -14, -3, 1);
707  neutEntryPtr->addChannel(1, 0.0, 0, 14, 3, -1);
708  neutEntryPtr->addChannel(1, 0.0, 0, -13, -4, 1);
709  neutEntryPtr->addChannel(1, 0.0, 0, 13, 4, -1);
710  neutEntryPtr->addChannel(1, 0.0, 0, -14, -3, 3);
711  neutEntryPtr->addChannel(1, 0.0, 0, 14, 3, -3);
712  neutEntryPtr->addChannel(1, 0.0, 0, -13, -4, 3);
713  neutEntryPtr->addChannel(1, 0.0, 0, 13, 4, -3);
714  neutEntryPtr->addChannel(1, 0.0, 0, -14, -3, 5);
715  neutEntryPtr->addChannel(1, 0.0, 0, 14, 3, -5);
716  neutEntryPtr->addChannel(1, 0.0, 0, -13, -4, 5);
717  neutEntryPtr->addChannel(1, 0.0, 0, 13, 4, -5);
718  neutEntryPtr->addChannel(1, 0.0, 0, -14, -5, 1);
719  neutEntryPtr->addChannel(1, 0.0, 0, 14, 5, -1);
720  neutEntryPtr->addChannel(1, 0.0, 0, -13, -6, 1);
721  neutEntryPtr->addChannel(1, 0.0, 0, 13, 6, -1);
722  neutEntryPtr->addChannel(1, 0.0, 0, -14, -5, 3);
723  neutEntryPtr->addChannel(1, 0.0, 0, 14, 5, -3);
724  neutEntryPtr->addChannel(1, 0.0, 0, -13, -6, 3);
725  neutEntryPtr->addChannel(1, 0.0, 0, 13, 6, -3);
726  neutEntryPtr->addChannel(1, 0.0, 0, -14, -5, 5);
727  neutEntryPtr->addChannel(1, 0.0, 0, 14, 5, -5);
728  neutEntryPtr->addChannel(1, 0.0, 0, -13, -6, 5);
729  neutEntryPtr->addChannel(1, 0.0, 0, 13, 6, -5);
730  neutEntryPtr->addChannel(1, 0.0, 0, -16, -1, 1);
731  neutEntryPtr->addChannel(1, 0.0, 0, 16, 1, -1);
732  neutEntryPtr->addChannel(1, 0.0, 0, -15, -2, 1);
733  neutEntryPtr->addChannel(1, 0.0, 0, 15, 2, -1);
734  neutEntryPtr->addChannel(1, 0.0, 0, -16, -1, 3);
735  neutEntryPtr->addChannel(1, 0.0, 0, 16, 1, -3);
736  neutEntryPtr->addChannel(1, 0.0, 0, -15, -2, 3);
737  neutEntryPtr->addChannel(1, 0.0, 0, 15, 2, -3);
738  neutEntryPtr->addChannel(1, 0.0, 0, -16, -1, 5);
739  neutEntryPtr->addChannel(1, 0.0, 0, 16, 1, -5);
740  neutEntryPtr->addChannel(1, 0.0, 0, -15, -2, 5);
741  neutEntryPtr->addChannel(1, 0.0, 0, 15, 2, -5);
742  neutEntryPtr->addChannel(1, 0.0, 0, -16, -3, 1);
743  neutEntryPtr->addChannel(1, 0.0, 0, 16, 3, -1);
744  neutEntryPtr->addChannel(1, 0.0, 0, -15, -4, 1);
745  neutEntryPtr->addChannel(1, 0.0, 0, 15, 4, -1);
746  neutEntryPtr->addChannel(1, 0.0, 0, -16, -3, 3);
747  neutEntryPtr->addChannel(1, 0.0, 0, 16, 3, -3);
748  neutEntryPtr->addChannel(1, 0.0, 0, -15, -4, 3);
749  neutEntryPtr->addChannel(1, 0.0, 0, 15, 4, -3);
750  neutEntryPtr->addChannel(1, 0.0, 0, -16, -3, 5);
751  neutEntryPtr->addChannel(1, 0.0, 0, 16, 3, -5);
752  neutEntryPtr->addChannel(1, 0.0, 0, -15, -4, 5);
753  neutEntryPtr->addChannel(1, 0.0, 0, 15, 4, -5);
754  neutEntryPtr->addChannel(1, 0.0, 0, -16, -5, 1);
755  neutEntryPtr->addChannel(1, 0.0, 0, 16, 5, -1);
756  neutEntryPtr->addChannel(1, 0.0, 0, -15, -6, 1);
757  neutEntryPtr->addChannel(1, 0.0, 0, 15, 6, -1);
758  neutEntryPtr->addChannel(1, 0.0, 0, -16, -5, 3);
759  neutEntryPtr->addChannel(1, 0.0, 0, 16, 5, -3);
760  neutEntryPtr->addChannel(1, 0.0, 0, -15, -6, 3);
761  neutEntryPtr->addChannel(1, 0.0, 0, 15, 6, -3);
762  neutEntryPtr->addChannel(1, 0.0, 0, -16, -5, 5);
763  neutEntryPtr->addChannel(1, 0.0, 0, 16, 5, -5);
764  neutEntryPtr->addChannel(1, 0.0, 0, -15, -6, 5);
765  neutEntryPtr->addChannel(1, 0.0, 0, 15, 6, -5);
766  neutEntryPtr->addChannel(1, 0.0, 0, -2 ,-1 ,-3);
767  neutEntryPtr->addChannel(1, 0.0, 0, 2, 1, 3);
768  neutEntryPtr->addChannel(1, 0.0, 0, -2, -1, -5);
769  neutEntryPtr->addChannel(1, 0.0, 0, 2, 1, 5);
770  neutEntryPtr->addChannel(1, 0.0, 0, -2, -3, -5);
771  neutEntryPtr->addChannel(1, 0.0, 0, 2, 3, 5);
772  neutEntryPtr->addChannel(1, 0.0, 0, -4, -1, -3);
773  neutEntryPtr->addChannel(1, 0.0, 0, 4, 1, 3);
774  neutEntryPtr->addChannel(1, 0.0, 0, -4, -1, -5);
775  neutEntryPtr->addChannel(1, 0.0, 0, 4, 1, 5);
776  neutEntryPtr->addChannel(1, 0.0, 0, -4, -3, -5);
777  neutEntryPtr->addChannel(1, 0.0, 0, 4, 3, 5);
778  neutEntryPtr->addChannel(1, 0.0, 0, -6, -1, -3);
779  neutEntryPtr->addChannel(1, 0.0, 0, 6, 1, 3);
780  neutEntryPtr->addChannel(1, 0.0, 0, -6, -1, -5);
781  neutEntryPtr->addChannel(1, 0.0, 0, 6, 1, 5);
782  neutEntryPtr->addChannel(1, 0.0, 0, -6, -3, -5);
783  neutEntryPtr->addChannel(1, 0.0, 0, 6, 3, 5);
784 
785  if (iNeut > 1) {
786 
787  neutEntryPtr->addChannel(1, 0.0, 0, 1000022, 22);
788  neutEntryPtr->addChannel(1, 0.0, 0, 1000022, 23);
789  neutEntryPtr->addChannel(1, 0.0, 0, 1000022, 25);
790  neutEntryPtr->addChannel(1, 0.0, 0, 1000022, 35);
791  neutEntryPtr->addChannel(1, 0.0, 0, 1000022, 36);
792 
793  if ( iNeut > 2) {
794  neutEntryPtr->addChannel(1, 0.0, 0, 1000023, 22);
795  neutEntryPtr->addChannel(1, 0.0, 0, 1000023, 23);
796  neutEntryPtr->addChannel(1, 0.0, 0, 1000023, 25);
797  neutEntryPtr->addChannel(1, 0.0, 0, 1000023, 35);
798  neutEntryPtr->addChannel(1, 0.0, 0, 1000023, 36);
799  }
800 
801  if (iNeut > 3) {
802  neutEntryPtr->addChannel(1, 0.0, 0, 1000025, 22);
803  neutEntryPtr->addChannel(1, 0.0, 0, 1000025, 23);
804  neutEntryPtr->addChannel(1, 0.0, 0, 1000025, 25);
805  neutEntryPtr->addChannel(1, 0.0, 0, 1000025, 35);
806  neutEntryPtr->addChannel(1, 0.0, 0, 1000025, 36);
807  }
808 
809  if (iNeut > 4) {
810  neutEntryPtr->addChannel(1, 0.0, 0, 1000035, 22);
811  neutEntryPtr->addChannel(1, 0.0, 0, 1000035, 23);
812  neutEntryPtr->addChannel(1, 0.0, 0, 1000035, 25);
813  neutEntryPtr->addChannel(1, 0.0, 0, 1000035, 35);
814  neutEntryPtr->addChannel(1, 0.0, 0, 1000035, 36);
815  }
816 
817  neutEntryPtr->addChannel(1, 0.0, 0, 1000024, -24);
818  neutEntryPtr->addChannel(1, 0.0, 0, -1000024, 24);
819  neutEntryPtr->addChannel(1, 0.0, 0, 1000037, -24);
820  neutEntryPtr->addChannel(1, 0.0, 0, -1000037, 24);
821  neutEntryPtr->addChannel(1, 0.0, 0, 1000024, -37);
822  neutEntryPtr->addChannel(1, 0.0, 0, -1000024, 37);
823  neutEntryPtr->addChannel(1, 0.0, 0, 1000037, -37);
824  neutEntryPtr->addChannel(1, 0.0, 0, -1000037, 37);
825  neutEntryPtr->addChannel(1, 0.0, 0, 1000011, -11);
826  neutEntryPtr->addChannel(1, 0.0, 0, -1000011, 11);
827  neutEntryPtr->addChannel(1, 0.0, 0, 2000011, -11);
828  neutEntryPtr->addChannel(1, 0.0, 0, -2000011, 11);
829  neutEntryPtr->addChannel(1, 0.0, 0, 1000012, -12);
830  neutEntryPtr->addChannel(1, 0.0, 0, -1000012, 12);
831  neutEntryPtr->addChannel(1, 0.0, 0, 1000013, -13);
832  neutEntryPtr->addChannel(1, 0.0, 0, -1000013, 13);
833  neutEntryPtr->addChannel(1, 0.0, 0, 2000013, -13);
834  neutEntryPtr->addChannel(1, 0.0, 0, -2000013, 13);
835  neutEntryPtr->addChannel(1, 0.0, 0, 1000014, -14);
836  neutEntryPtr->addChannel(1, 0.0, 0, -1000014, 14);
837  neutEntryPtr->addChannel(1, 0.0, 0, 1000015, -15);
838  neutEntryPtr->addChannel(1, 0.0, 0, -1000015, 15);
839  neutEntryPtr->addChannel(1, 0.0, 0, 2000015, -15);
840  neutEntryPtr->addChannel(1, 0.0, 0, -2000015, 15);
841  neutEntryPtr->addChannel(1, 0.0, 0, 1000016, -16);
842  neutEntryPtr->addChannel(1, 0.0, 0, -1000016, 16);
843  neutEntryPtr->addChannel(1, 0.0, 0, 1000001, -1);
844  neutEntryPtr->addChannel(1, 0.0, 0, -1000001, 1);
845  neutEntryPtr->addChannel(1, 0.0, 0, 1000001, -3);
846  neutEntryPtr->addChannel(1, 0.0, 0, -1000001, 3);
847  neutEntryPtr->addChannel(1, 0.0, 0, 1000001, -5);
848  neutEntryPtr->addChannel(1, 0.0, 0, -1000001, 5);
849  neutEntryPtr->addChannel(1, 0.0, 0, 2000001, -1);
850  neutEntryPtr->addChannel(1, 0.0, 0, -2000001, 1);
851  neutEntryPtr->addChannel(1, 0.0, 0, 2000001, -3);
852  neutEntryPtr->addChannel(1, 0.0, 0, -2000001, 3);
853  neutEntryPtr->addChannel(1, 0.0, 0, 2000001, -5);
854  neutEntryPtr->addChannel(1, 0.0, 0, -2000001, 5);
855  neutEntryPtr->addChannel(1, 0.0, 0, 1000002, -2);
856  neutEntryPtr->addChannel(1, 0.0, 0, -1000002, 2);
857  neutEntryPtr->addChannel(1, 0.0, 0, 1000002, -4);
858  neutEntryPtr->addChannel(1, 0.0, 0, -1000002, 4);
859  neutEntryPtr->addChannel(1, 0.0, 0, 1000002, -6);
860  neutEntryPtr->addChannel(1, 0.0, 0, -1000002, 6);
861  neutEntryPtr->addChannel(1, 0.0, 0, 2000002, -2);
862  neutEntryPtr->addChannel(1, 0.0, 0, -2000002, 2);
863  neutEntryPtr->addChannel(1, 0.0, 0, 2000002, -4);
864  neutEntryPtr->addChannel(1, 0.0, 0, -2000002, 4);
865  neutEntryPtr->addChannel(1, 0.0, 0, 2000002, -6);
866  neutEntryPtr->addChannel(1, 0.0, 0, -2000002, 6);
867  neutEntryPtr->addChannel(1, 0.0, 0, 1000003, -1);
868  neutEntryPtr->addChannel(1, 0.0, 0, -1000003, 1);
869  neutEntryPtr->addChannel(1, 0.0, 0, 1000003, -3);
870  neutEntryPtr->addChannel(1, 0.0, 0, -1000003, 3);
871  neutEntryPtr->addChannel(1, 0.0, 0, 1000003, -5);
872  neutEntryPtr->addChannel(1, 0.0, 0, -1000003, 5);
873  neutEntryPtr->addChannel(1, 0.0, 0, 2000003, -1);
874  neutEntryPtr->addChannel(1, 0.0, 0, -2000003, 1);
875  neutEntryPtr->addChannel(1, 0.0, 0, 2000003, -3);
876  neutEntryPtr->addChannel(1, 0.0, 0, -2000003, 3);
877  neutEntryPtr->addChannel(1, 0.0, 0, 2000003, -5);
878  neutEntryPtr->addChannel(1, 0.0, 0, -2000003, 5);
879  neutEntryPtr->addChannel(1, 0.0, 0, 1000004, -2);
880  neutEntryPtr->addChannel(1, 0.0, 0, -1000004, 2);
881  neutEntryPtr->addChannel(1, 0.0, 0, 1000004, -4);
882  neutEntryPtr->addChannel(1, 0.0, 0, -1000004, 4);
883  neutEntryPtr->addChannel(1, 0.0, 0, 1000004, -6);
884  neutEntryPtr->addChannel(1, 0.0, 0, -1000004, 6);
885  neutEntryPtr->addChannel(1, 0.0, 0, 2000004, -2);
886  neutEntryPtr->addChannel(1, 0.0, 0, -2000004, 2);
887  neutEntryPtr->addChannel(1, 0.0, 0, 2000004, -4);
888  neutEntryPtr->addChannel(1, 0.0, 0, -2000004, 4);
889  neutEntryPtr->addChannel(1, 0.0, 0, 2000004, -6);
890  neutEntryPtr->addChannel(1, 0.0, 0, -2000004, 6);
891  neutEntryPtr->addChannel(1, 0.0, 0, 1000005, -1);
892  neutEntryPtr->addChannel(1, 0.0, 0, -1000005, 1);
893  neutEntryPtr->addChannel(1, 0.0, 0, 1000005, -3);
894  neutEntryPtr->addChannel(1, 0.0, 0, -1000005, 3);
895  neutEntryPtr->addChannel(1, 0.0, 0, 1000005, -5);
896  neutEntryPtr->addChannel(1, 0.0, 0, -1000005, 5);
897  neutEntryPtr->addChannel(1, 0.0, 0, 2000005, -1);
898  neutEntryPtr->addChannel(1, 0.0, 0, -2000005, 1);
899  neutEntryPtr->addChannel(1, 0.0, 0, 2000005, -3);
900  neutEntryPtr->addChannel(1, 0.0, 0, -2000005, 3);
901  neutEntryPtr->addChannel(1, 0.0, 0, 2000005, -5);
902  neutEntryPtr->addChannel(1, 0.0, 0, -2000005, 5);
903  neutEntryPtr->addChannel(1, 0.0, 0, 1000006, -6);
904  neutEntryPtr->addChannel(1, 0.0, 0, -1000006, 6);
905  neutEntryPtr->addChannel(1, 0.0, 0, 1000006, -2);
906  neutEntryPtr->addChannel(1, 0.0, 0, -1000006, 2);
907  neutEntryPtr->addChannel(1, 0.0, 0, 1000006, -4);
908  neutEntryPtr->addChannel(1, 0.0, 0, -1000006, 4);
909  neutEntryPtr->addChannel(1, 0.0, 0, 2000006, -6);
910  neutEntryPtr->addChannel(1, 0.0, 0, -2000006, 6);
911 
912  // Modes involving right-handed sneutrinos are not included by default,
913  // but can be added by hand, by uncommenting the following lines.
914  // neutEntryPtr->addChannel(1, 0.0, 0, 2000012, -12);
915  // neutEntryPtr->addChannel(1, 0.0, 0, -2000012, 12);
916  // neutEntryPtr->addChannel(1, 0.0, 0, 2000014, -14);
917  // neutEntryPtr->addChannel(1, 0.0, 0, -2000014, 14);
918  // neutEntryPtr->addChannel(1, 0.0, 0, 2000016, -16);
919  // neutEntryPtr->addChannel(1, 0.0, 0, -2000016, 16);
920 
921 
922 
923  }
924  return true;
925 }
926 
927 //--------------------------------------------------------------------------
928 
929 void ResonanceNeut::initConstants() {
930 
931  s2W = coupSUSYPtr->sin2W;
932 
933  // Initialize functions for calculating 3-body widths
934  // psi.setPointers(particleDataPtr,coupSUSYPtr,infoPtr);
935  // phi.setPointers(particleDataPtr,coupSUSYPtr,infoPtr);
936  // upsil.setPointers(particleDataPtr,coupSUSYPtr,infoPtr);
937 
938 }
939 
940 //--------------------------------------------------------------------------
941 
942 // Calculate various common prefactors for the current mass.
943 
944 void ResonanceNeut::calcPreFac(bool) {
945 
946  // Common coupling factors.
947  alpEM = coupSMPtr->alphaEM(mHat * mHat);
948  preFac = alpEM / (8.0 * s2W * pow(mHat,3));
949  return;
950 
951 }
952 
953 //--------------------------------------------------------------------------
954 
955 // Calculate width for currently considered channel.
956 void ResonanceNeut::calcWidth(bool){
957 
958  widNow = 0.0;
959 
960  if (ps == 0.) return;
961  else if (mult ==2){
962  // Two-body decays
963 
964  kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
965  kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
966  + pow2(mHat) * pow2(mf2) + pow2(mf1) * pow2(mf2)
967  - 2.0 * pow2(mHat) * pow2(mf1);
968 
969  // Stable lightest neutralino
970  if (idRes == 1000022) return;
971 
972  double fac = 0.0;
973  int iNeut1 = coupSUSYPtr->typeNeut(idRes);
974  int iNeut2 = coupSUSYPtr->typeNeut(id1Abs);
975  int iChar1 = coupSUSYPtr->typeChar(id1Abs);
976 
977  if (iNeut2>0 && id2Abs == 23){
978  // ~chi0_i -> chi0_j + Z
979  fac = kinFac2 * (norm(coupSUSYPtr->OLpp[iNeut1][iNeut2])
980  + norm(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
981  fac -= 12.0 * mHat * mf1 * pow2(mf2)
982  * real(coupSUSYPtr->OLpp[iNeut1][iNeut2]
983  * conj(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
984  fac /= pow2(mf2) * (1.0 - s2W);
985  }
986  else if (iChar1>0 && id2Abs==24){
987  // ~chi0_i -> chi+_j + W- (or c.c.)
988 
989  fac = kinFac2 * (norm(coupSUSYPtr->OL[iNeut1][iChar1])
990  + norm(coupSUSYPtr->OR[iNeut1][iChar1]));
991  fac -= 12.0 * mHat * mf1 * pow2(mf2)
992  * real(coupSUSYPtr->OL[iNeut1][iChar1]
993  * conj(coupSUSYPtr->OR[iNeut1][iChar1]));
994  fac /= pow2(mf2);
995  }
996  else if (id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
997  // ~chi0_k -> ~q + q
998  bool idown = (id1Abs%2 == 1);
999  int iq = (id2Abs + 1 )/ 2;
1000  int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1001  : (abs(id1Abs)%10+1)/2;
1002 
1003  if (idown){
1004  fac = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][iNeut1])
1005  + norm(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
1006  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddX[isq][iq][iNeut1]
1007  * conj(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
1008  }
1009  else{
1010  fac = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][iNeut1])
1011  + norm(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
1012  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][iNeut1]
1013  * conj(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
1014  }
1015  // Extra multiplicative factor of 3 over sleptons
1016  fac *= 6.0/(1 - s2W);
1017  }
1018  else if (id1Abs > 2000010 && id1Abs%2 == 0 ) {
1019  // Check for right-handed neutralinos.
1020  widNow = 0;
1021  }
1022  else if (id1Abs > 1000000 && id1Abs%100 > 10 && id1Abs%100 < 17
1023  && id2Abs < 17){
1024 
1025  // ~chi0_k -> ~l + l
1026  bool idown = id2Abs%2;
1027  int il = (id2Abs - 9)/ 2;
1028  int isl = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1029  : (abs(id1Abs)%10+1)/2;
1030 
1031  if (idown){
1032  fac = kinFac * (norm(coupSUSYPtr->LsllX[isl][il][iNeut1])
1033  + norm(coupSUSYPtr->RsllX[isl][il][iNeut1]));
1034  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsllX[isl][il][iNeut1]
1035  * conj(coupSUSYPtr->RsllX[isl][il][iNeut1]));
1036  }
1037  else{
1038  fac = kinFac * (norm(coupSUSYPtr->LsvvX[isl][il][iNeut1]));
1039  }
1040  fac *= 2.0/(1 - s2W);
1041  }
1042  // TODO: Decays in higgs
1043  // Final width for 2-body decays
1044  widNow = fac * preFac * ps * pow2(mHat);
1045  if (DBSUSY) {
1046  cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
1047  cout<<scientific<<widNow<<endl;
1048  }
1049  }
1050  else {
1051  //RPV 3-body decays
1052 
1053  //Case: UDD-type (TO BE RE-DONE to fix bug)
1054  return;
1055 
1056  }
1057 
1058  // Normalisation. Extra factor of 2 to account for the fact that
1059  // d_i, d_j will always be ordered in ascending order. N_c! = 6
1060  widNow *= 12.0 /(pow3(2.0 * M_PI * mHat) * 32.0);
1061 
1062  return;
1063 }
1064 
1065 //==========================================================================
1066 
1067 // Class ResonanceChar
1068 // Derived class for Neutralino Resonances
1069 // Decays into higgses/sleptons not yet implemented
1070 
1071 //--------------------------------------------------------------------------
1072 
1073 // Set up Channels
1074 
1075 bool ResonanceChar::getChannels(int idPDG){
1076 
1077  idPDG = abs(idPDG);
1078  int iChar = coupSUSYPtr->typeChar(idPDG);
1079  if (iChar < 1) return false;
1080 
1081  ParticleDataEntry* charEntryPtr
1082  = particleDataPtr->particleDataEntryPtr(idPDG);
1083 
1084  // Delete any decay channels read
1085  charEntryPtr->clearChannels();
1086 
1087  charEntryPtr->addChannel(1, 0.0, 0, 1000022, 24);
1088  charEntryPtr->addChannel(1, 0.0, 0, 1000023, 24);
1089  charEntryPtr->addChannel(1, 0.0, 0, 1000025, 24);
1090  charEntryPtr->addChannel(1, 0.0, 0, 1000035, 24);
1091  charEntryPtr->addChannel(1, 0.0, 0, 1000022, 37);
1092  charEntryPtr->addChannel(1, 0.0, 0, 1000023, 37);
1093  charEntryPtr->addChannel(1, 0.0, 0, 1000025, 37);
1094  charEntryPtr->addChannel(1, 0.0, 0, 1000035, 37);
1095  charEntryPtr->addChannel(1, 0.0, 0, 1000012, -11);
1096  charEntryPtr->addChannel(1, 0.0, 0, -1000011, 12);
1097  charEntryPtr->addChannel(1, 0.0, 0, -2000011, 12);
1098  charEntryPtr->addChannel(1, 0.0, 0, 1000014, -13);
1099  charEntryPtr->addChannel(1, 0.0, 0, -1000013, 14);
1100  charEntryPtr->addChannel(1, 0.0, 0, -2000013, 14);
1101  charEntryPtr->addChannel(1, 0.0, 0, 1000016, -15);
1102  charEntryPtr->addChannel(1, 0.0, 0, -1000015, 16);
1103  charEntryPtr->addChannel(1, 0.0, 0, -2000015, 16);
1104  charEntryPtr->addChannel(1, 0.0, 0, 1000002, -1);
1105  charEntryPtr->addChannel(1, 0.0, 0, 1000002, -3);
1106  charEntryPtr->addChannel(1, 0.0, 0, 1000002, -5);
1107  charEntryPtr->addChannel(1, 0.0, 0, 2000002, -1);
1108  charEntryPtr->addChannel(1, 0.0, 0, 2000002, -3);
1109  charEntryPtr->addChannel(1, 0.0, 0, 2000002, -5);
1110  charEntryPtr->addChannel(1, 0.0, 0, -1000001, 2);
1111  charEntryPtr->addChannel(1, 0.0, 0, -1000001, 4);
1112  charEntryPtr->addChannel(1, 0.0, 0, -1000001, 6);
1113  charEntryPtr->addChannel(1, 0.0, 0, -2000001, 2);
1114  charEntryPtr->addChannel(1, 0.0, 0, -2000001, 4);
1115  charEntryPtr->addChannel(1, 0.0, 0, -2000001, 6);
1116  charEntryPtr->addChannel(1, 0.0, 0, 1000004, -1);
1117  charEntryPtr->addChannel(1, 0.0, 0, 1000004, -3);
1118  charEntryPtr->addChannel(1, 0.0, 0, 1000004, -5);
1119  charEntryPtr->addChannel(1, 0.0, 0, 2000004, -1);
1120  charEntryPtr->addChannel(1, 0.0, 0, 2000004, -3);
1121  charEntryPtr->addChannel(1, 0.0, 0, 2000004, -5);
1122  charEntryPtr->addChannel(1, 0.0, 0, -1000003, 2);
1123  charEntryPtr->addChannel(1, 0.0, 0, -1000003, 4);
1124  charEntryPtr->addChannel(1, 0.0, 0, -1000003, 6);
1125  charEntryPtr->addChannel(1, 0.0, 0, -2000003, 2);
1126  charEntryPtr->addChannel(1, 0.0, 0, -2000003, 4);
1127  charEntryPtr->addChannel(1, 0.0, 0, -2000003, 6);
1128  charEntryPtr->addChannel(1, 0.0, 0, 1000006, -1);
1129  charEntryPtr->addChannel(1, 0.0, 0, 1000006, -3);
1130  charEntryPtr->addChannel(1, 0.0, 0, 1000006, -5);
1131  charEntryPtr->addChannel(1, 0.0, 0, 2000006, -1);
1132  charEntryPtr->addChannel(1, 0.0, 0, 2000006, -3);
1133  charEntryPtr->addChannel(1, 0.0, 0, 2000006, -5);
1134  charEntryPtr->addChannel(1, 0.0, 0, -1000005, 2);
1135  charEntryPtr->addChannel(1, 0.0, 0, -1000005, 4);
1136  charEntryPtr->addChannel(1, 0.0, 0, -1000005, 6);
1137  charEntryPtr->addChannel(1, 0.0, 0, -2000005, 2);
1138  charEntryPtr->addChannel(1, 0.0, 0, -2000005, 4);
1139  charEntryPtr->addChannel(1, 0.0, 0, -2000005, 6);
1140 
1141  if (iChar > 1) {
1142  charEntryPtr->addChannel(1, 0.0, 0, 1000024, 23);
1143  charEntryPtr->addChannel(1, 0.0, 0, 1000024, 25);
1144  charEntryPtr->addChannel(1, 0.0, 0, 1000024, 35);
1145  charEntryPtr->addChannel(1, 0.0, 0, 1000024, 36);
1146  }
1147 
1148  // Modes involving right-handed sneutrinos are not included by default,
1149  // but can be added by hand, by uncommenting the following lines.
1150  // charEntryPtr->addChannel(1, 0.0, 0, 2000012, -11);
1151  // charEntryPtr->addChannel(1, 0.0, 0, 2000014, -13);
1152  // charEntryPtr->addChannel(1, 0.0, 0, 2000016, -15);
1153 
1154  return true;
1155 
1156 }
1157 
1158 //--------------------------------------------------------------------------
1159 
1160 void ResonanceChar::initConstants() {
1161 
1162  s2W = coupSUSYPtr->sin2W;
1163  return;
1164 
1165 }
1166 
1167 //--------------------------------------------------------------------------
1168 
1169 // Calculate various common prefactors for the current mass.
1170 
1171 void ResonanceChar::calcPreFac(bool) {
1172 
1173  alpEM = coupSMPtr->alphaEM(mHat * mHat);
1174  preFac = alpEM / (8.0 * s2W * pow(mHat,3));
1175  return;
1176 
1177 }
1178 
1179 //--------------------------------------------------------------------------
1180 
1181 // Calculate width for currently considered channel.
1182 
1183 void ResonanceChar::calcWidth(bool) {
1184 
1185  widNow = 0.0;
1186  if (ps == 0.) return;
1187 
1188  if (mult ==2){
1189  double fac = 0.0;
1190  kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
1191  kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
1192  + pow2(mHat) * pow2(mf2) + pow2(mf1)
1193  * pow2(mf2) - 2.0 * pow2(mHat) * pow2(mf1);
1194 
1195  int idChar1 = coupSUSYPtr->typeChar(idRes);
1196  int idChar2 = coupSUSYPtr->typeChar(id1Abs);
1197  int idNeut1 = coupSUSYPtr->typeNeut(id1Abs);
1198 
1199  if (idChar2>0 && id2Abs == 23){
1200  // ~chi_i -> chi_j + Z
1201  fac = kinFac2 * (norm(coupSUSYPtr->OLp[idChar1][idChar2])
1202  + norm(coupSUSYPtr->ORp[idChar1][idChar2]));
1203  fac -= 12.0 * mHat * mf1 * pow2(mf2)
1204  * real(coupSUSYPtr->OLp[idChar1][idChar2]
1205  * conj(coupSUSYPtr->ORp[idChar1][idChar2]));
1206  fac /= pow2(mf2) * (1.0 - s2W);
1207  }
1208  else if (idNeut1>0 && id2Abs==24){
1209  // ~chi_i -> chi0_j + W- (or c.c.)
1210 
1211  fac = kinFac2 * (norm(coupSUSYPtr->OL[idNeut1][idChar1])
1212  + norm(coupSUSYPtr->OR[idNeut1][idChar1]));
1213  fac -= 12.0 * mHat * mf1 * pow2(mf2)
1214  * real(coupSUSYPtr->OL[idNeut1][idChar1]
1215  * conj(coupSUSYPtr->OR[idNeut1][idChar1]));
1216  fac /= pow2(mf2);
1217  }
1218  else if (id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
1219  // ~chi0_k -> ~q + q
1220  bool idown = (id1Abs%2 == 1);
1221  int iq = (id2Abs + 1 )/ 2;
1222  int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1223  : (abs(id1Abs)%10+1)/2;
1224 
1225  if (idown){
1226  fac = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][idChar1])
1227  + norm(coupSUSYPtr->RsduX[isq][iq][idChar1]));
1228  fac += 4.0 * mHat * mf2
1229  * real(coupSUSYPtr->LsduX[isq][iq][idChar1]
1230  * conj(coupSUSYPtr->RsduX[isq][iq][idChar1]));
1231  }
1232  else{
1233  fac = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][idChar1])
1234  + norm(coupSUSYPtr->RsudX[isq][iq][idChar1]));
1235  fac += 4.0 * mHat * mf2
1236  * real(coupSUSYPtr->LsudX[isq][iq][idChar1]
1237  * conj(coupSUSYPtr->RsudX[isq][iq][idChar1]));
1238  }
1239  fac *= 6.0/(1 - s2W);
1240  }
1241  else if (id1Abs > 2000010 && id1Abs%2 == 0 ) {
1242  // Check for right-handed neutralinos.
1243  widNow = 0;
1244  }
1245  else if (id1Abs > 1000000 && id1Abs%100 > 10 && id1Abs%100 < 17
1246  && id2Abs < 17){
1247  // ~chi+_k -> ~l + l
1248  bool idown = id2Abs%2;
1249  int il = (id2Abs - 9)/ 2;
1250  int isl = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1251  : (abs(id1Abs)%10+1)/2;
1252 
1253  if (idown){
1254  fac = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][idChar1])
1255  + norm(coupSUSYPtr->RslvX[isl][il][idChar1]));
1256  fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][idChar1]
1257  * conj(coupSUSYPtr->RslvX[isl][il][idChar1]));
1258  }
1259  else{
1260  fac = kinFac * (norm(coupSUSYPtr->LsvlX[isl][il][idChar1]));
1261  }
1262  fac *= 2.0/(1 - s2W);
1263  }
1264 
1265  // TODO: Decays in higgs
1266  widNow = fac * preFac * ps * pow2(mHat);
1267  if (DBSUSY) {
1268  cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
1269  cout<<scientific<<widNow<<endl;
1270  }
1271  }else{
1272  //TODO: Implement Chargino 3-body decays
1273  }
1274  return;
1275 }
1276 
1277 //==========================================================================
1278 // The ResonanceSlepton class
1279 // Derived class for Slepton (and sneutrino) resonances
1280 
1281 //--------------------------------------------------------------------------
1282 
1283 // Set up Channels
1284 
1285 bool ResonanceSlepton::getChannels(int idPDG){
1286 
1287  idPDG = abs(idPDG);
1288 
1289  int ksusy = 1000000;
1290  if (idPDG < ksusy) return false;
1291  if(idPDG % ksusy < 7 || idPDG % ksusy > 17) return false;
1292 
1293  ParticleDataEntry* slepEntryPtr
1294  = particleDataPtr->particleDataEntryPtr(idPDG);
1295 
1296  // Delete any decay channels read
1297  slepEntryPtr->clearChannels();
1298 
1299  if(idPDG % 2 == 1) {
1300 
1301  slepEntryPtr->addChannel(1, 0.0, 0, -1000024, 16);
1302  slepEntryPtr->addChannel(1, 0.0, 0, -1000037, 16);
1303  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, 15);
1304  slepEntryPtr->addChannel(1, 0.0, 0, 1000023, 15);
1305  slepEntryPtr->addChannel(1, 0.0, 0, 1000025, 15);
1306  slepEntryPtr->addChannel(1, 0.0, 0, 1000035, 15);
1307  slepEntryPtr->addChannel(1, 0.0, 0, 1000016, -24);
1308  slepEntryPtr->addChannel(1, 0.0, 0, 2000016, -24);
1309  slepEntryPtr->addChannel(1, 0.0, 0, 1000016, -37);
1310  slepEntryPtr->addChannel(1, 0.0, 0, 2000016, -37);
1311  slepEntryPtr->addChannel(1, 0.0, 0, 12, 13);
1312  slepEntryPtr->addChannel(1, 0.0, 0, 12, 15);
1313  slepEntryPtr->addChannel(1, 0.0, 0, 14, 11);
1314  slepEntryPtr->addChannel(1, 0.0, 0, 14, 15);
1315  slepEntryPtr->addChannel(1, 0.0, 0, 16, 11);
1316  slepEntryPtr->addChannel(1, 0.0, 0, 16, 13);
1317  slepEntryPtr->addChannel(1, 0.0, 0, -12, 11);
1318  slepEntryPtr->addChannel(1, 0.0, 0, -12, 13);
1319  slepEntryPtr->addChannel(1, 0.0, 0, -12, 15);
1320  slepEntryPtr->addChannel(1, 0.0, 0, -14, 11);
1321  slepEntryPtr->addChannel(1, 0.0, 0, -14, 13);
1322  slepEntryPtr->addChannel(1, 0.0, 0, -14, 15);
1323  slepEntryPtr->addChannel(1, 0.0, 0, -2, 1);
1324  slepEntryPtr->addChannel(1, 0.0, 0, -2, 3);
1325  slepEntryPtr->addChannel(1, 0.0, 0, -2, 5);
1326  slepEntryPtr->addChannel(1, 0.0, 0, -4, 1);
1327  slepEntryPtr->addChannel(1, 0.0, 0, -4, 3);
1328  slepEntryPtr->addChannel(1, 0.0, 0, -4, 5);
1329  slepEntryPtr->addChannel(1, 0.0, 0, -6, 1);
1330  slepEntryPtr->addChannel(1, 0.0, 0, -6, 3);
1331  slepEntryPtr->addChannel(1, 0.0, 0, -6, 5);
1332  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, -211, 16);
1333  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, -213, 16);
1334  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, -9000211, 16);
1335  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, 16, 12, 11);
1336  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, 16, 14, 13);
1337 
1338  } else {
1339  slepEntryPtr->addChannel(1, 0.0, 0, 1000024, 15);
1340  slepEntryPtr->addChannel(1, 0.0, 0, 1000037, 15);
1341  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, 16);
1342  slepEntryPtr->addChannel(1, 0.0, 0, 1000023, 16);
1343  slepEntryPtr->addChannel(1, 0.0, 0, 1000025, 16);
1344  slepEntryPtr->addChannel(1, 0.0, 0, 1000035, 16);
1345  slepEntryPtr->addChannel(1, 0.0, 0, 1000015, 24);
1346  slepEntryPtr->addChannel(1, 0.0, 0, 2000015, 24);
1347  slepEntryPtr->addChannel(1, 0.0, 0, 1000015, 37);
1348  slepEntryPtr->addChannel(1, 0.0, 0, 2000015, 37);
1349  slepEntryPtr->addChannel(1, 0.0, 0, -11, 11);
1350  slepEntryPtr->addChannel(1, 0.0, 0, -11, 13);
1351  slepEntryPtr->addChannel(1, 0.0, 0, -11, 15);
1352  slepEntryPtr->addChannel(1, 0.0, 0, -13, 11);
1353  slepEntryPtr->addChannel(1, 0.0, 0, -13, 13);
1354  slepEntryPtr->addChannel(1, 0.0, 0, -13, 15);
1355  slepEntryPtr->addChannel(1, 0.0, 0, -1, 1);
1356  slepEntryPtr->addChannel(1, 0.0, 0, -1, 3);
1357  slepEntryPtr->addChannel(1, 0.0, 0, -1, 5);
1358  slepEntryPtr->addChannel(1, 0.0, 0, -3, 1);
1359  slepEntryPtr->addChannel(1, 0.0, 0, -3, 3);
1360  slepEntryPtr->addChannel(1, 0.0, 0, -3, 5);
1361  slepEntryPtr->addChannel(1, 0.0, 0, -5, 1);
1362  slepEntryPtr->addChannel(1, 0.0, 0, -5, 3);
1363  slepEntryPtr->addChannel(1, 0.0, 0, -5, 5);
1364  }
1365 
1366  return true;
1367 }
1368 
1369 //--------------------------------------------------------------------------
1370 
1371 // Initialize constants.
1372 
1373 void ResonanceSlepton::initConstants() {
1374 
1375  // Locally stored properties and couplings.
1376  s2W = coupSUSYPtr->sin2W;
1377 
1378  // Initialize functions for calculating 3-body widths
1379  stauWidths.setPointers(infoPtr);
1380 
1381 }
1382 
1383 //--------------------------------------------------------------------------
1384 
1385 // Calculate various common prefactors for the current mass.
1386 
1387 void ResonanceSlepton::calcPreFac(bool) {
1388 
1389  // Common coupling factors.
1390  alpEM = coupSMPtr->alphaEM(mHat * mHat);
1391  preFac = 1.0 / (s2W * pow(mHat,3));
1392 
1393 }
1394 
1395 //--------------------------------------------------------------------------
1396 
1397 // Calculate width for currently considered channel.
1398 
1399 void ResonanceSlepton::calcWidth(bool) {
1400 
1401  // Slepton type -- in u_i/d_i and generation
1402  int ksusy = 1000000;
1403  int isl = (abs(idRes)/ksusy == 2) ? (abs(idRes)%10+1)/2 + 3
1404  : (abs(idRes)%10+1)/2;
1405  int il = (id2Abs-9)/2;
1406  bool islep = abs(idRes)%2;
1407 
1408  // Check that mass is above threshold.
1409  if (ps == 0.) return;
1410  widNow = 0.0;
1411 
1412  double fac = 0.0 , wid = 0.0;
1413 
1414  if (mult == 2) {
1415  // Two-body decays
1416 
1417  kinFac = (mHat * mHat - mf1 * mf1 - mf2 * mf2);
1418  fac = kinFac / (16.0 * M_PI * pow(mHat,3));
1419 
1420  // Case 1: RPV decays
1421  if (id1Abs < 17 && id2Abs < 17) {
1422 
1423  wid = 0;
1424  int il2 = (id1Abs - 9)/2;
1425 
1426  //Case 1a: RPV LLE
1427  if(id1Abs > 10 && id2Abs > 10) {
1428  if (!coupSUSYPtr->isLLE) { widNow = 0.0; return;}
1429 
1430  if (!islep){ // sneutrino
1431  for (int isl2=1; isl2<3; isl2++)
1432  wid += norm(coupSUSYPtr->Rsv[isl][isl2]
1433  * coupSUSYPtr->rvLLE[il][isl2][il2]);
1434  } else {
1435  for (int isl2=1; isl2<3; isl2++)
1436  wid += norm(coupSUSYPtr->Rsl[isl][isl2+3]
1437  * coupSUSYPtr->rvLLE[isl2][il][il2]);
1438  }
1439  }
1440  //Case 1b: RPV LQD
1441  if(id1Abs < 10 && id2Abs < 10) {
1442  if (!coupSUSYPtr->isLQD) { widNow = 0.0; return;}
1443  if (!islep){ // sneutrino
1444  for (int isl2=1; isl2<3; isl2++)
1445  wid += norm(coupSUSYPtr->Rsv[isl][isl2]
1446  * coupSUSYPtr->rvLQD[isl2][id1Abs][id2Abs]);
1447  wid *= 3.0; // colour factor
1448 
1449  } else {
1450  for (int isl2=1; isl2<3; isl2++)
1451  wid += norm(coupSUSYPtr->Rsl[isl][isl2+3]
1452  * coupSUSYPtr->rvLLE[isl2][id1Abs][id2Abs]);
1453  wid *= 3.0; // colour factor
1454  }
1455  }
1456  }
1457 
1458  //Case 2: lepton + gaugino
1459 
1460  if (id1Abs > ksusy && id2Abs > 10 && id2Abs < 17) {
1461  for (int i=1; i<6 ; i++){
1462  // ~ell/~nu -> ~chi0 + ell/nu
1463  if (coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
1464  fac = alpEM * preFac / (2.0 * (1 - s2W));
1465  if (islep)
1466  wid = kinFac * (norm(coupSUSYPtr->LsllX[isl][il][i])
1467  + norm(coupSUSYPtr->RsllX[isl][il][i]))
1468  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsllX[isl][il][i]
1469  * conj(coupSUSYPtr->RsllX[isl][il][i]));
1470  else
1471  wid = kinFac * (norm(coupSUSYPtr->LsvvX[isl][il][i])
1472  + norm(coupSUSYPtr->RsvvX[isl][il][i]))
1473  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsvvX[isl][il][i]
1474  * conj(coupSUSYPtr->RsvvX[isl][il][i]));
1475  }
1476 
1477  // ~ell/~nu -> ~chi- + nu/ell
1478  else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs
1479  && idRes%2 != id2Abs%2){
1480 
1481  fac = alpEM * preFac / (4.0 * (1 - s2W));
1482 
1483  if (islep)
1484  wid = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][i])
1485  + norm(coupSUSYPtr->RslvX[isl][il][i]))
1486  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][i]
1487  * conj(coupSUSYPtr->RslvX[isl][il][i]));
1488  else
1489  wid = kinFac * (norm(coupSUSYPtr->LsvvX[isl][il][i])
1490  + norm(coupSUSYPtr->RsvvX[isl][il][i]))
1491  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsvvX[isl][il][i]
1492  * conj(coupSUSYPtr->RsvvX[isl][il][i]));
1493  }
1494  }
1495  }
1496 
1497  //Case 3: ~l_i -> ~l_j + Z/W
1498  else if (id1Abs > ksusy+10 && id1Abs%100 < 17
1499  && (id2Abs == 23 || id2Abs == 24)){
1500 
1501  // factor of lambda^(3/2) = ps^3;
1502  fac = alpEM * preFac/(16.0 * pow2(mf2) * (1.0 - s2W)) * pow2(ps) ;
1503 
1504  int isl2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
1505  isl2 = min(6, isl2);
1506 
1507  if (id2Abs == 23 && id1Abs%2 == idRes%2){
1508  if (islep)
1509  wid = norm(coupSUSYPtr->LslslZ[isl][isl2]
1510  + coupSUSYPtr->RslslZ[isl][isl2]);
1511  else
1512  wid = norm(coupSUSYPtr->LsvsvZ[isl][isl2]
1513  + coupSUSYPtr->RsvsvZ[isl][isl2]);
1514  }
1515  else if (id2Abs == 24 && id1Abs%2 != idRes%2){
1516  if (islep)
1517  wid = norm(coupSUSYPtr->LslsvW[isl2][isl]);
1518  else
1519  wid = norm(coupSUSYPtr->LslsvW[isl][isl2]);
1520  }
1521  }
1522 
1523  // TODO: Case ~l_i -> ~l_j + h/H
1524 
1525  widNow = fac * wid * ps * pow2(mHat);
1526 
1527  if (DBSUSY) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs
1528  <<" Width: "<<widNow<<endl;
1529  }
1530  else { // Case 4: special many-body stau decays
1531  // Case 4a: ~tau -> pi0 nu_tau ~chi0
1532  // Case 4b: ~tau -> rho/a1 nu_tau ~chi0
1533  // Case 4c: ~tau -> l nu_l nu_tau ~chi0
1534 
1535 
1536  // Check we are in the compressed regime
1537  if (mRes - particleDataPtr->m0(1000022) > particleDataPtr->m0(15)) return;
1538 
1539 
1540  // Check that there is a stau component
1541  double staufac = norm(coupSUSYPtr->Rsl[isl][3])
1542  + norm(coupSUSYPtr->Rsl[isl][6]);
1543 
1544  if (staufac < 1.0e-6) return;
1545  if(id2Abs > 17) {
1546  int idIn = id2Abs == 1000022 ? id1Abs : id2Abs;
1547  widNow = stauWidths.getWidth(idRes, idIn);
1548  }
1549  else
1550  widNow = stauWidths.getWidth(idRes, id3Abs);
1551 
1552  widNow *= staufac;
1553 
1554  }
1555 
1556  return;
1557 
1558 }
1559 
1560 //==========================================================================
1561 
1562 } //end namespace Pythia8
Generated by   doxygen 1.8.5