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) 2018 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 (couplingsPtr->isSUSY) {
30  coupSUSYPtr = (CoupSUSY *) couplingsPtr;
31  return true;
32  }
33 
34  return false;
35 }
36 
37 //--------------------------------------------------------------------------
38 
39 bool SUSYResonanceWidths::allowCalc(){
40 
41  // Check if decay calculations at all possible
42  if ( !couplingsPtr->isSUSY ) return false;
43  if ( (idRes == 45 || idRes == 46 || idRes == 1000045)
44  && !coupSUSYPtr->isNMSSM ) return false;
45 
46  if (settingsPtr->flag("SLHA:useDecayTable") ) {
47 
48  // Next check if decay table was read in via SLHA and takes precedence
49  for ( int iDec = 0; iDec < int((coupSUSYPtr->slhaPtr)->decays.size());
50  ++iDec)
51  if ( (coupSUSYPtr->slhaPtr)->decays[iDec].getId() == abs(idRes) ) {
52  if (DBSUSY) cout<<"Using external decay table for:"<<idRes<<endl;
53  return false;
54  }
55 
56  }
57 
58  // Else we should do the calculation; set available channels
59  bool done = getChannels(idRes);
60  stringstream idStream;
61  idStream << "ID = " << idRes ;
62  if (!done) infoPtr->errorMsg("Error in SusyResonanceWidths::allowcalc: "
63  "unable to reset decay table.", idStream.str(), true);
64  return done;
65 }
66 
67 //==========================================================================
68 
69 // The ResonanceSquark class
70 // Derived class for Squark resonances
71 
72 //--------------------------------------------------------------------------
73 
74 // Set up channels
75 
76 bool ResonanceSquark::getChannels(int idPDG){
77 
78  idPDG = abs(idPDG);
79 
80  int ksusy = 1000000;
81  if (idPDG < ksusy) return false;
82  if(idPDG % ksusy >= 7 || idPDG % ksusy < 1) return false;
83 
84  ParticleDataEntry* squarkEntryPtr
85  = particleDataPtr->particleDataEntryPtr(idPDG);
86 
87  // Delete any decay channels read
88  squarkEntryPtr->clearChannels();
89 
90  if(idPDG % 2 == 0) { // up-type squarks
91 
92  /*
93  // Gravitino
94  squarkEntryPtr->addChannel(1, 0.0, 103, 1000039, 2);
95  squarkEntryPtr->addChannel(1, 0.0, 0, 1000039, 4);
96  squarkEntryPtr->addChannel(1, 0.0, 0, 1000039, 6);
97  */
98 
99  // Gaugino - quark
100  squarkEntryPtr->addChannel(1, 0.0, 0, 1000024, 3);
101  squarkEntryPtr->addChannel(1, 0.0, 0, 1000024, 5);
102  squarkEntryPtr->addChannel(1, 0.0, 0, 1000037, 1);
103  squarkEntryPtr->addChannel(1, 0.0, 0, 1000037, 3);
104  squarkEntryPtr->addChannel(1, 0.0, 0, 1000037, 5);
105  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 2);
106  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 4);
107  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 6);
108  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 2);
109  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 4);
110  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 6);
111  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 2);
112  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 4);
113  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 6);
114  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 2);
115  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 4);
116  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 6);
117 
118  // Squark - Gauge boson
119  squarkEntryPtr->addChannel(1, 0.0, 0, 1000001, -24);
120  squarkEntryPtr->addChannel(1, 0.0, 0, 1000003, -24);
121  squarkEntryPtr->addChannel(1, 0.0, 0, 1000005, -24);
122  squarkEntryPtr->addChannel(1, 0.0, 0, 2000001, -24);
123  squarkEntryPtr->addChannel(1, 0.0, 0, 2000003, -24);
124  squarkEntryPtr->addChannel(1, 0.0, 0, 2000005, -24);
125  squarkEntryPtr->addChannel(1, 0.0, 0, 1000001, -37);
126  squarkEntryPtr->addChannel(1, 0.0, 0, 1000003, -37);
127  squarkEntryPtr->addChannel(1, 0.0, 0, 1000005, -37);
128  squarkEntryPtr->addChannel(1, 0.0, 0, 2000001, -37);
129  squarkEntryPtr->addChannel(1, 0.0, 0, 2000003, -37);
130  squarkEntryPtr->addChannel(1, 0.0, 0, 2000005, -37);
131 
132  // Gluino - quark
133  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 2);
134  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 4);
135  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 6);
136 
137  // lepton - quark via LQD
138  squarkEntryPtr->addChannel(1, 0.0, 0, -11, 1);
139  squarkEntryPtr->addChannel(1, 0.0, 0, -11, 3);
140  squarkEntryPtr->addChannel(1, 0.0, 0, -11, 5);
141  squarkEntryPtr->addChannel(1, 0.0, 0, -13, 1);
142  squarkEntryPtr->addChannel(1, 0.0, 0, -13, 3);
143  squarkEntryPtr->addChannel(1, 0.0, 0, -13, 5);
144  squarkEntryPtr->addChannel(1, 0.0, 0, -15, 1);
145  squarkEntryPtr->addChannel(1, 0.0, 0, -15, 3);
146  squarkEntryPtr->addChannel(1, 0.0, 0, -15, 5);
147 
148  // quark - quark via UDD
149  squarkEntryPtr->addChannel(1, 0.0, 0, -1 ,-3);
150  squarkEntryPtr->addChannel(1, 0.0, 0, -1 ,-5);
151  squarkEntryPtr->addChannel(1, 0.0, 0, -3 ,-5);
152 
153 
154  } else { //down-type squarks
155 
156  // Gaugino - quark
157  squarkEntryPtr->addChannel(1, 0.0, 0, -1000024, 2);
158  squarkEntryPtr->addChannel(1, 0.0, 0, -1000037, 2);
159  squarkEntryPtr->addChannel(1, 0.0, 0, -1000024, 4);
160  squarkEntryPtr->addChannel(1, 0.0, 0, -1000037, 4);
161  squarkEntryPtr->addChannel(1, 0.0, 0, -1000024, 6);
162  squarkEntryPtr->addChannel(1, 0.0, 0, -1000037, 6);
163  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 1);
164  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 3);
165  squarkEntryPtr->addChannel(1, 0.0, 0, 1000022, 5);
166  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 1);
167  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 3);
168  squarkEntryPtr->addChannel(1, 0.0, 0, 1000023, 5);
169  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 1);
170  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 3);
171  squarkEntryPtr->addChannel(1, 0.0, 0, 1000025, 5);
172  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 1);
173  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 3);
174  squarkEntryPtr->addChannel(1, 0.0, 0, 1000035, 5);
175 
176  // Squark - Gauge boson
177  squarkEntryPtr->addChannel(1, 0.0, 0, 1000002, -24);
178  squarkEntryPtr->addChannel(1, 0.0, 0, 1000004, -24);
179  squarkEntryPtr->addChannel(1, 0.0, 0, 1000006, -24);
180  squarkEntryPtr->addChannel(1, 0.0, 0, 2000002, -24);
181  squarkEntryPtr->addChannel(1, 0.0, 0, 2000004, -24);
182  squarkEntryPtr->addChannel(1, 0.0, 0, 2000006, -24);
183  squarkEntryPtr->addChannel(1, 0.0, 0, 1000002, -37);
184  squarkEntryPtr->addChannel(1, 0.0, 0, 1000004, -37);
185  squarkEntryPtr->addChannel(1, 0.0, 0, 1000006, -37);
186  squarkEntryPtr->addChannel(1, 0.0, 0, 2000002, -37);
187  squarkEntryPtr->addChannel(1, 0.0, 0, 2000004, -37);
188  squarkEntryPtr->addChannel(1, 0.0, 0, 2000006, -37);
189 
190  // Gluino - quark
191  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 1);
192  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 2);
193  squarkEntryPtr->addChannel(1, 0.0, 0, 1000021, 5);
194 
195  // lepton - quark via LQD
196  squarkEntryPtr->addChannel(1, 0.0, 0, -12, 1);
197  squarkEntryPtr->addChannel(1, 0.0, 0, -12, 3);
198  squarkEntryPtr->addChannel(1, 0.0, 0, -12, 5);
199  squarkEntryPtr->addChannel(1, 0.0, 0, -14, 1);
200  squarkEntryPtr->addChannel(1, 0.0, 0, -14, 3);
201  squarkEntryPtr->addChannel(1, 0.0, 0, -14, 5);
202  squarkEntryPtr->addChannel(1, 0.0, 0, -16, 1);
203  squarkEntryPtr->addChannel(1, 0.0, 0, -16, 3);
204  squarkEntryPtr->addChannel(1, 0.0, 0, -16, 5);
205  squarkEntryPtr->addChannel(1, 0.0, 0, 12 ,1);
206  squarkEntryPtr->addChannel(1, 0.0, 0, 11 ,2);
207  squarkEntryPtr->addChannel(1, 0.0, 0, 12, 3);
208  squarkEntryPtr->addChannel(1, 0.0, 0, 11, 4);
209  squarkEntryPtr->addChannel(1, 0.0, 0, 12, 5);
210  squarkEntryPtr->addChannel(1, 0.0, 0, 11, 6);
211  squarkEntryPtr->addChannel(1, 0.0, 0, 14, 1);
212  squarkEntryPtr->addChannel(1, 0.0, 0, 13, 2);
213  squarkEntryPtr->addChannel(1, 0.0, 0, 14, 3);
214  squarkEntryPtr->addChannel(1, 0.0, 0, 13, 4);
215  squarkEntryPtr->addChannel(1, 0.0, 0, 14, 5);
216  squarkEntryPtr->addChannel(1, 0.0, 0, 13, 6);
217  squarkEntryPtr->addChannel(1, 0.0, 0, 16, 1);
218  squarkEntryPtr->addChannel(1, 0.0, 0, 15, 2);
219  squarkEntryPtr->addChannel(1, 0.0, 0, 16, 3);
220  squarkEntryPtr->addChannel(1, 0.0, 0, 15, 4);
221  squarkEntryPtr->addChannel(1, 0.0, 0, 16, 5);
222  squarkEntryPtr->addChannel(1, 0.0, 0, 15, 6);
223 
224 
225  // quark - quark via UDD
226  squarkEntryPtr->addChannel(1, 0.0, 0, -2, -1);
227  squarkEntryPtr->addChannel(1, 0.0, 0, -2, -3);
228  squarkEntryPtr->addChannel(1, 0.0, 0, -2, -5);
229  squarkEntryPtr->addChannel(1, 0.0, 0, -4, -1);
230  squarkEntryPtr->addChannel(1, 0.0, 0, -4, -3);
231  squarkEntryPtr->addChannel(1, 0.0, 0, -4, -5);
232  squarkEntryPtr->addChannel(1, 0.0, 0, -6, -1);
233  squarkEntryPtr->addChannel(1, 0.0, 0, -6, -3);
234  squarkEntryPtr->addChannel(1, 0.0, 0, -6, -5);
235  }
236 
237  return true;
238 
239 }
240 
241 //--------------------------------------------------------------------------
242 
243 // Initialize constants.
244 
245 void ResonanceSquark::initConstants() {
246 
247  // Locally stored properties and couplings.
248  s2W = coupSUSYPtr->sin2W;
249 
250 }
251 
252 //--------------------------------------------------------------------------
253 
254 // Calculate various common prefactors for the current mass.
255 
256 void ResonanceSquark::calcPreFac(bool) {
257 
258  // Common coupling factors.
259  alpS = coupSUSYPtr->alphaS(mHat * mHat );
260  alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
261  preFac = 1.0 / (s2W * pow(mHat,3));
262  ps *= mHat * mHat;
263 
264 }
265 
266 //--------------------------------------------------------------------------
267 
268 // Calculate width for currently considered channel.
269 
270 void ResonanceSquark::calcWidth(bool) {
271 
272  // Squark type -- in u_i/d_i and generation
273  int ksusy = 1000000;
274  bool idown = (abs(idRes)%2 == 0 ? false : true);
275  int isq = (abs(idRes)/ksusy == 2) ?
276  (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
277  // int isqgen = (abs(idRes)%10 + 1)/2;
278 
279  // Check that mass is above threshold.
280  if (ps == 0.) return;
281  else{
282  // Two-body decays
283  kinFac = (mHat * mHat - mf1 * mf1 - mf2 * mf2);
284 
285  double fac = 0.0 , wid = 0.0;
286 
287  //RPV decays
288  //Case 1a: UDD-type
289  if (id1Abs < 7 && id2Abs < 7){
290 
291  // Quark generations
292  int iq1 = (id1Abs + 1)/2;
293  int iq2 = (id2Abs + 1)/2;
294 
295  // Check for RPV UDD couplings
296  if (!coupSUSYPtr->isUDD) {widNow = 0; return;}
297 
298  // ~q -> q_i + q_j
299 
300  fac = 2.0 * kinFac / (16.0 * M_PI * pow(mHat,3));
301  wid = 0.0;
302  if (idown) {
303  if ((id1Abs+id2Abs)%2 == 1){
304  if (id1Abs%2==1)
305  for (int isq2 = 1; isq2 < 4; isq2++)
306  wid += norm(coupSUSYPtr->rvUDD[iq2][iq1][isq2]
307  * coupSUSYPtr->Rdsq[isq][isq2+3]);
308  else
309  for (int isq2 = 1; isq2 < 4; isq2++)
310  wid += norm(coupSUSYPtr->rvUDD[iq1][iq2][isq2]
311  * coupSUSYPtr->Rdsq[isq][isq2+3]);
312  }
313  }
314  else {
315  if ((id1Abs+id2Abs)%2 != 0) widNow = 0.0;
316  else
317  for (int isq2 = 1; isq2 < 4; isq2++)
318  wid += norm(coupSUSYPtr->rvUDD[isq2][iq1][iq2]
319  * coupSUSYPtr->Rusq[isq][isq2+3]);
320  }
321  }
322 
323  //Case 1b: LQD-type
324  else if (id1Abs < 17 && id2Abs < 7){
325  if (!coupSUSYPtr->isLQD) {widNow = 0; return;}
326 
327  int ilep = (id1Abs - 9)/2;
328  int iq = (id2Abs + 1)/2;
329 
330  fac = kinFac / (16.0 * M_PI * pow(mHat,3));
331  wid = 0.0;
332  if (idown){
333  if (iq%2 == 0){
334  // q is up-type; ~q is right-handed down type
335  for (int isq2=1; isq2<3; isq2++)
336  wid += norm(coupSUSYPtr->Rdsq[isq][isq2+3]
337  * coupSUSYPtr->rvLQD[ilep][iq][isq2]);
338  }else{
339  //q is down type; ~q left-handed down-type
340  for (int isq2=1; isq2<3; isq2++)
341  wid += norm(coupSUSYPtr->Rdsq[isq][isq2]
342  * coupSUSYPtr->rvLQD[ilep][isq2][isq2]);
343  }
344  }
345  else{
346  if (iq%2 == 0) {widNow = 0.0; return;}
347  // q is down type; ~q is left-handed up-type
348  for (int isq2=1; isq2<3; isq2++)
349  wid += norm(coupSUSYPtr->Rusq[isq][isq2]
350  * coupSUSYPtr->rvLQD[ilep][isq2][iq]);
351  }
352  }
353 
354  //Case 2: quark + gaugino
355  else if (id1Abs > ksusy && id2Abs < 7) {
356 
357  int iq = (id2Abs + 1)/2;
358 
359  // ~q -> ~g + q
360  if (id1Abs == 1000021 && idRes%10 == id2Abs) {
361  // Removed factor of s2W in denominator: strong process -- no EW
362  fac = 2.0 * alpS / (3.0 * pow3(mHat));
363  if (idown)
364  wid = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])
365  + norm(coupSUSYPtr->RsddG[isq][iq]))
366  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddG[isq][iq]
367  * conj(coupSUSYPtr->RsddG[isq][iq]));
368  else
369  wid = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])
370  + norm(coupSUSYPtr->RsuuG[isq][iq]))
371  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]
372  * conj(coupSUSYPtr->RsuuG[isq][iq]));
373  }
374  else
375  for (int i=1; i<6 ; i++){
376  // ~q -> ~chi0 + q
377  if (coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
378  fac = alpEM * preFac / (2.0 * (1 - s2W));
379  if (idown)
380  wid = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][i])
381  + norm(coupSUSYPtr->RsddX[isq][iq][i]))
382  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddX[isq][iq][i]
383  * conj(coupSUSYPtr->RsddX[isq][iq][i]));
384  else
385  wid = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][i])
386  + norm(coupSUSYPtr->RsuuX[isq][iq][i]))
387  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][i]
388  * conj(coupSUSYPtr->RsuuX[isq][iq][i]));
389  }
390 
391  // ~q -> chi- + q
392  else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs
393  && idRes%2 != id2Abs%2){
394 
395  fac = alpEM * preFac / (4.0 * (1 - s2W));
396  if (idown)
397  wid = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][i])
398  + norm(coupSUSYPtr->RsduX[isq][iq][i]))
399  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsduX[isq][iq][i]
400  * conj(coupSUSYPtr->RsduX[isq][iq][i]));
401  else
402  wid = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][i])
403  + norm(coupSUSYPtr->RsudX[isq][iq][i]))
404  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsudX[isq][iq][i]
405  * conj(coupSUSYPtr->RsudX[isq][iq][i]));
406  }
407  }
408  }
409 
410  //Case 3: ~q_i -> ~q_j + Z/W
411  else if (id1Abs > ksusy && id1Abs%100 < 7
412  && (id2Abs == 23 || id2Abs == 24)){
413 
414  // factor of lambda^(3/2) = ps^(3/2) ;
415  fac = alpEM * preFac/(16.0 * pow2(particleDataPtr->m0(id2Abs))
416  * (1.0 - s2W)) * pow2(ps) ;
417 
418  int isq2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
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 = coupSUSYPtr->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 = coupSUSYPtr->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 = coupSUSYPtr->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, 111, 16);
1333  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, 113, 16);
1334  slepEntryPtr->addChannel(1, 0.0, 0, 1000022, 900111, 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(particleDataPtr,coupSUSYPtr,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 = coupSUSYPtr->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  if (islep)
1483  wid = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][i])
1484  + norm(coupSUSYPtr->RslvX[isl][il][i]))
1485  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][i]
1486  * conj(coupSUSYPtr->RslvX[isl][il][i]));
1487  else
1488  wid = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][i])
1489  + norm(coupSUSYPtr->RslvX[isl][il][i]))
1490  - 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][i]
1491  * conj(coupSUSYPtr->RslvX[isl][il][i]));
1492  }
1493  }
1494  }
1495 
1496  //Case 3: ~l_i -> ~l_j + Z/W
1497  else if (id1Abs > ksusy+10 && id1Abs%100 < 17
1498  && (id2Abs == 23 || id2Abs == 24)){
1499 
1500  // factor of lambda^(3/2) = ps^3;
1501  fac = alpEM * preFac/(16.0 * pow2(mf2) * (1.0 - s2W)) * pow2(ps) ;
1502 
1503  int isl2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
1504 
1505  if (id2Abs == 23 && id1Abs%2 == idRes%2){
1506  if (islep)
1507  wid = norm(coupSUSYPtr->LslslZ[isl][isl2]
1508  + coupSUSYPtr->RslslZ[isl][isl2]);
1509  else
1510  wid = norm(coupSUSYPtr->LsvsvZ[isl][isl2]
1511  + coupSUSYPtr->RsvsvZ[isl][isl2]);
1512  }
1513  else if (id2Abs == 24 && id1Abs%2 != idRes%2){
1514  if (islep)
1515  wid = norm(coupSUSYPtr->LslsvW[isl2][isl]);
1516  else
1517  wid = norm(coupSUSYPtr->LslsvW[isl][isl2]);
1518  }
1519  }
1520 
1521  // TODO: Case ~l_i -> ~l_j + h/H
1522 
1523  widNow = fac * wid * ps * pow2(mHat);
1524 
1525  if (DBSUSY) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs
1526  <<" Width: "<<widNow<<endl;
1527  }
1528  else { // Case 4: special many-body stau decays
1529  // Case 4a: ~tau -> pi0 nu_tau ~chi0
1530  // Case 4b: ~tau -> rho/a1 nu_tau ~chi0
1531  // Case 4c: ~tau -> l nu_l nu_tau ~chi0
1532 
1533  // Check that there is a stau component
1534  double staufac = norm(coupSUSYPtr->Rsl[isl][3])
1535  + norm(coupSUSYPtr->Rsl[isl][6]);
1536  if (staufac < 1.0e-6 || abs(mRes - particleDataPtr->m0(15)) > 0.0 ) return;
1537  if(id2Abs > 17)
1538  widNow = stauWidths.getWidth(idRes, id2Abs);
1539  else
1540  widNow = stauWidths.getWidth(idRes, id3Abs);
1541 
1542  widNow *= staufac;
1543 
1544  }
1545 
1546  return;
1547 
1548 }
1549 
1550 //==========================================================================
1551 
1552 } //end namespace Pythia8
Generated by   doxygen 1.8.5