gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
MSSM.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
27 
31 
33 
34 namespace Gambit
35 {
36  namespace DarkBit
37  {
38 
40  //
41  // Backend point initialization
42  //
44 
54  //
56  // General catalog for annihilation/decay process definition
57  //
59 
61  double DSgamma3bdy(double(*IBfunc)(int&,double&,double&),
62  int IBch, double Eg, double
63  E1, double M_DM, double m_1, double m_2)
64  {
65  double E2 = 2*M_DM - Eg - E1;
66  double p12 = E1*E1-m_1*m_1;
67  double p22 = E2*E2-m_2*m_2;
68  double p22min = Eg*Eg+p12-2*Eg*sqrt(p12);
69  double p22max = Eg*Eg+p12+2*Eg*sqrt(p12);
70  // Check if the process is kinematically allowed
71  if((E1 < m_1) || (E2 < m_2) || (p22<p22min) || (p22>p22max))
72  {
73  return 0;
74  }
75  double x = Eg/M_DM; // x = E_gamma/mx
76  double y = (m_2*m_2 + 4*M_DM * (M_DM - E2) ) / (4*M_DM*M_DM);
77  // Here, y = (p+k)^2/(4 mx^2), where p denotes the W+ momentum and k the
78  // photon momentum. (note that the expressions above and below only
79  // apply to the v->0 limit)
80 
81  double result = IBfunc(IBch,x,y);
82 
83  #ifdef DARKBIT_DEBUG
84  std::cout << " x, y = " << x << ", " << y << std::endl;
85  std::cout << " E, E1, E2 = " << Eg << ", " << E1 << ", "
86  << E2 << std::endl;
87  std::cout << " mDM, m1, m2 = " << M_DM << ", " << m_1 << ", "
88  << m_2 << std::endl;
89  std::cout << " IBfunc = " << result << std::endl;
90  #endif
91 
92  // M_DM^-2 is from the Jacobi determinant
93  return std::max(0., result) / (M_DM*M_DM);
94  }
95 
96 
101  {
102  using namespace Pipes::TH_ProcessCatalog_DS5_MSSM;
103  using std::vector;
104  using std::string;
105 
106  std::string DMid = *Dep::DarkMatter_ID;
107  if ( DMid != "~chi0_1" )
108  {
110  "TH_ProcessCatalog_DS5_MSSM requires DMid to be ~chi0_1.");
111  }
112 
113  // Instantiate new empty ProcessCatalog
114  TH_ProcessCatalog catalog;
115 
116 
118  // Import particle masses
120 
121  // Import based on Spectrum objects
122  const Spectrum& matched_spectra = *Dep::MSSM_spectrum;
123  const SubSpectrum& spec = matched_spectra.get_HE();
124  const SubSpectrum& SM = matched_spectra.get_LE();
125  const SMInputs& SMI = matched_spectra.get_SMInputs();
126 
127  // Get SM masses
128  auto getSMmass = [&](str Name, int spinX2)
129  {
130  catalog.particleProperties.insert(
131  std::pair<std::string, TH_ParticleProperty>(
132  Name , TH_ParticleProperty(SM.get(Par::Pole_Mass,Name), spinX2)));
133  };
134 
135  getSMmass("e-_1", 1);
136  getSMmass("e+_1", 1);
137  getSMmass("e-_2", 1);
138  getSMmass("e+_2", 1);
139  getSMmass("e-_3", 1);
140  getSMmass("e+_3", 1);
141  getSMmass("Z0", 2);
142  getSMmass("W+", 2);
143  getSMmass("W-", 2);
144  getSMmass("g", 2);
145  getSMmass("gamma", 2);
146  getSMmass("d_3", 1);
147  getSMmass("dbar_3", 1);
148  getSMmass("u_3", 1);
149  getSMmass("ubar_3", 1);
150 
151  // Pole masses not available for the light quarks.
152  auto addParticle = [&](str Name, double Mass, int spinX2)
153  {
154  catalog.particleProperties.insert(
155  std::pair<std::string, TH_ParticleProperty>(
156  Name , TH_ParticleProperty(Mass, spinX2)));
157  };
158 
159  addParticle("d_1" , SMI.mD, 1); // md(2 GeV)^MS-bar, not pole mass
160  addParticle("dbar_1", SMI.mD, 1); // md(2 GeV)^MS-bar, not pole mass
161  addParticle("u_1" , SMI.mU, 1); // mu(2 GeV)^MS-bar, not pole mass
162  addParticle("ubar_1", SMI.mU, 1); // mu(2 GeV)^MS-bar, not pole mass
163  addParticle("d_2" , SMI.mS, 1); // ms(2 GeV)^MS-bar, not pole mass
164  addParticle("dbar_2", SMI.mS, 1); // ms(2 GeV)^MS-bar, not pole mass
165  addParticle("u_2" , SMI.mCmC,1); // mc(mc)^MS-bar, not pole mass
166  addParticle("ubar_2", SMI.mCmC,1); // mc(mc)^MS-bar, not pole mass
167  // Masses for neutrino flavour eigenstates. Set to zero.
168  addParticle("nu_e", 0.0, 1);
169  addParticle("nubar_e", 0.0, 1);
170  addParticle("nu_mu", 0.0, 1);
171  addParticle("nubar_mu", 0.0, 1);
172  addParticle("nu_tau", 0.0, 1);
173  addParticle("nubar_tau",0.0, 1);
174 
175  addParticle("pi0", meson_masses.pi0, 0);
176  addParticle("pi+", meson_masses.pi_plus, 0);
177  addParticle("pi-", meson_masses.pi_minus, 0);
178  addParticle("eta", meson_masses.eta, 0);
179  addParticle("rho0", meson_masses.rho0, 1);
180  addParticle("rho+", meson_masses.rho_plus, 1);
181  addParticle("rho-", meson_masses.rho_minus, 1);
182  addParticle("omega", meson_masses.omega, 1);
183 
184 
185  // Get MSSM masses
186  auto getMSSMmass = [&](str Name, int spinX2)
187  {
188  catalog.particleProperties.insert(
189  std::pair<std::string, TH_ParticleProperty> (
190  Name , TH_ParticleProperty(std::abs(spec.get(Par::Pole_Mass,Name)), spinX2)));
191  };
192 
193  getMSSMmass("H+" , 0);
194  getMSSMmass("H-" , 0);
195  getMSSMmass("h0_1" , 0);
196  getMSSMmass("h0_2" , 0);
197  getMSSMmass("A0" , 0);
198  getMSSMmass("~chi0_1" , 1);
199  getMSSMmass("~d_1" , 0);
200  getMSSMmass("~dbar_1" , 0);
201  getMSSMmass("~u_1" , 0);
202  getMSSMmass("~ubar_1" , 0);
203  getMSSMmass("~d_2" , 0);
204  getMSSMmass("~dbar_2" , 0);
205  getMSSMmass("~u_2" , 0);
206  getMSSMmass("~ubar_2" , 0);
207  getMSSMmass("~d_3" , 0);
208  getMSSMmass("~dbar_3" , 0);
209  getMSSMmass("~u_3" , 0);
210  getMSSMmass("~ubar_3" , 0);
211  getMSSMmass("~d_4" , 0);
212  getMSSMmass("~dbar_4" , 0);
213  getMSSMmass("~u_4" , 0);
214  getMSSMmass("~ubar_4" , 0);
215  getMSSMmass("~d_5" , 0);
216  getMSSMmass("~dbar_5" , 0);
217  getMSSMmass("~u_5" , 0);
218  getMSSMmass("~ubar_5" , 0);
219  getMSSMmass("~d_6" , 0);
220  getMSSMmass("~dbar_6" , 0);
221  getMSSMmass("~u_6" , 0);
222  getMSSMmass("~ubar_6" , 0);
223  //getMSSMmass("~e_1" , 0);
224  //getMSSMmass("~ebar_1" , 0);
225  //getMSSMmass("~e-_1" , 0);
226  getMSSMmass("~e+_1" , 0);
227  getMSSMmass("~e-_1" , 0);
228  getMSSMmass("~e+_2" , 0);
229  getMSSMmass("~e-_2" , 0);
230  getMSSMmass("~e+_3" , 0);
231  getMSSMmass("~e-_3" , 0);
232  getMSSMmass("~e+_4" , 0);
233  getMSSMmass("~e-_4" , 0);
234  getMSSMmass("~e+_5" , 0);
235  getMSSMmass("~e-_5" , 0);
236  getMSSMmass("~e+_6" , 0);
237  getMSSMmass("~e-_6" , 0);
238  getMSSMmass("~nu_1" , 0);
239  getMSSMmass("~nubar_1", 0);
240  getMSSMmass("~nu_2" , 0);
241  getMSSMmass("~nubar_2", 0);
242  getMSSMmass("~nu_3" , 0);
243  getMSSMmass("~nubar_3", 0);
244 
245 
247  // Import two-body annihilation process
249 
250  // Set of possible final state particles. Used to determine which decays to import.
251  std::set<string> annFinalStates;
252 
253  // Declare DM annihilation process
254  TH_Process process(DMid, DMid);
255 
256  // Explicitly state that Neutralino DM is self-conjugate
257  process.isSelfConj = true;
258 
259  double M_DM = catalog.getParticleProperty(DMid).mass;
260 
261  // lambda for setting up 2-body annihilations (chi chi -> X X) from results in DS
262  auto setup_ds_process = [&](int index, str p1, str p2, double prefac)
263  {
264  /* Check if process is kinematically allowed */
265  double m_1 = catalog.getParticleProperty(p1).mass;
266  double m_2 = catalog.getParticleProperty(p2).mass;
267  if(m_1 + m_2 < 2*M_DM)
268  {
269  /* Set cross-section */
270  double sigma = BEreq::dssigmav(index);
271  /* Create associated kinematical functions (just dependent on vrel)
272  * here: s-wave, vrel independent 1-dim constant function */
273  daFunk::Funk kinematicFunction = daFunk::cnst(sigma*prefac, "v");
274  /* Create channel identifier string */
275  std::vector<std::string> finalStates;
276  finalStates.push_back(p1);
277  finalStates.push_back(p2);
278  /* Create channel and push it into channel list of process */
279  process.channelList.push_back(TH_Channel(finalStates, kinematicFunction));
280  annFinalStates.insert(p1);
281  annFinalStates.insert(p1);
282  }
283  };
284 
285  setup_ds_process(1 , "h0_2", "h0_2", 1 );
286  setup_ds_process(2 , "h0_2", "h0_1", 1 );
287  setup_ds_process(3 , "h0_1", "h0_1", 1 );
288  setup_ds_process(4 , "A0", "A0", 1 );
289  setup_ds_process(5 , "h0_2", "A0", 1 );
290  setup_ds_process(6 , "h0_1", "A0", 1 );
291  setup_ds_process(7 , "H+", "H-", 1 );
292  setup_ds_process(8 , "h0_2", "Z0", 1 );
293  setup_ds_process(9 , "h0_1", "Z0", 1 );
294  setup_ds_process(10, "A0", "Z0", 1 );
295  // Prefactor 0.5 since W+H- and W-H+ are summed in DS
296  setup_ds_process(11, "W+", "H-", 0.5 );
297  setup_ds_process(11, "W-", "H+", 0.5 );
298  setup_ds_process(12, "Z0", "Z0", 1 );
299  setup_ds_process(13, "W+", "W-", 1 );
300  setup_ds_process(14, "nu_e", "nubar_e",1 );
301  setup_ds_process(15, "e+_1", "e-_1", 1 );
302  setup_ds_process(16, "nu_mu", "nubar_mu",1 );
303  setup_ds_process(17, "e+_2", "e-_2", 1 );
304  setup_ds_process(18, "nu_tau", "nubar_tau",1 );
305  setup_ds_process(19, "e+_3", "e-_3", 1 );
306  setup_ds_process(20, "u_1", "ubar_1", 1 );
307  setup_ds_process(21, "d_1", "dbar_1", 1 );
308  setup_ds_process(22, "u_2", "ubar_2", 1 );
309  setup_ds_process(23, "d_2", "dbar_2", 1 );
310  setup_ds_process(24, "u_3", "ubar_3", 1 );
311  setup_ds_process(25, "d_3", "dbar_3", 1 );
312  setup_ds_process(26, "g", "g", 1 );
313  setup_ds_process(28, "gamma", "gamma", 1 );
314  // setup_ds_process(29, "Z0", gamma, 1 );
315 
316 
318  // Import three-body annihilation process
320 
322 
323  // lambda for setting up 3-body decays with gammas
324  auto setup_ds_process_gamma3body = [&](int IBch, str p1, str p2, double (*IBfunc)(int&, double&, double&), int index, double prefac)
325  {
326  /* Check if process is kinematically allowed */
327  double m_1 = catalog.getParticleProperty(str_flav_to_mass(p1)).mass;
328  double m_2 = catalog.getParticleProperty(str_flav_to_mass(p2)).mass;
329  if(m_1 + m_2 < 2*M_DM)
330  {
331  double sv = prefac*BEreq::dssigmav(index);
332  daFunk::Funk kinematicFunction =
334  IBfunc, IBch, daFunk::var("E"), daFunk::var("E1"),
335  M_DM, m_1, m_2);
336  /* Create channel identifier string */
337  std::vector<std::string> finalStates;
338  finalStates.push_back("gamma");
339  finalStates.push_back(str_flav_to_mass(p1));
340  finalStates.push_back(str_flav_to_mass(p2));
341  /* Create channel and push it into channel list of process */
342  process.channelList.push_back(TH_Channel(finalStates, kinematicFunction));
343  annFinalStates.insert(str_flav_to_mass(p1));
344  annFinalStates.insert(str_flav_to_mass(p2));
345  };
346  };
347 
349  if ( not runOptions->getValueOrDef<bool>(false, "ignore_three_body") )
350  {
351  // Set DarkSUSY DM mass parameter used in 3-body decays
352  BEreq::IBintvars->ibcom_mx = catalog.getParticleProperty(DMid).mass;
353 
354  setup_ds_process_gamma3body(1, "W+", "W-", BEreq::dsIBwwdxdy.pointer(), 13, 1 );
355  // Prefactor 0.5 since W+H- and W-H+ are summed in DS
356  setup_ds_process_gamma3body(2, "W+", "H-", BEreq::dsIBwhdxdy.pointer(), 11, 0.5);
357  // Prefactor 0.5 since W+H- and W-H+ are summed in DS
358  setup_ds_process_gamma3body(2, "W-", "H+", BEreq::dsIBwhdxdy.pointer(), 11, 0.5);
359  setup_ds_process_gamma3body(3, "H+", "H-", BEreq::dsIBhhdxdy.pointer(), 0, 1 );
360  setup_ds_process_gamma3body(4, "e+", "e-", BEreq::dsIBffdxdy.pointer(), 15, 1 );
361  setup_ds_process_gamma3body(5, "mu+", "mu-", BEreq::dsIBffdxdy.pointer(), 17, 1 );
362  setup_ds_process_gamma3body(6, "tau+", "tau-", BEreq::dsIBffdxdy.pointer(), 19, 1 );
363  setup_ds_process_gamma3body(7, "u", "ubar", BEreq::dsIBffdxdy.pointer(), 20, 1 );
364  setup_ds_process_gamma3body(8, "d", "dbar", BEreq::dsIBffdxdy.pointer(), 21, 1 );
365  setup_ds_process_gamma3body(9, "c", "cbar", BEreq::dsIBffdxdy.pointer(), 22, 1 );
366  setup_ds_process_gamma3body(10, "s", "sbar", BEreq::dsIBffdxdy.pointer(), 23, 1 );
367  setup_ds_process_gamma3body(11, "t", "tbar", BEreq::dsIBffdxdy.pointer(), 24, 1 );
368  setup_ds_process_gamma3body(12, "b", "bbar", BEreq::dsIBffdxdy.pointer(), 25, 1 );
369  };
370 
371 
373  // Import Decay information
375 
376  // Import based on decay table from DecayBit
377  const DecayTable* tbl = &(*Dep::decay_rates);
378 
379  // Set of imported decays - avoids double imports
380  std::set<string> importedDecays;
381 
382  // Option minBranching <double>: Minimum branching fraction of included
383  // processes (default 0.)
384  double minBranching = runOptions->getValueOrDef<double>(0.0,
385  "ProcessCatalog_MinBranching");
386 
387  // Exclude also ttbar final states
388  auto excludeDecays = daFunk::vec<std::string>("Z0", "W+", "W-", "u_3", "ubar_3", "e+_2", "e-_2", "e+_3", "e-_3");
389 
390  // Import relevant decays
392  if(annFinalStates.count("H+") == 1)
393  ImportDecays("H+", catalog, importedDecays, tbl, minBranching, excludeDecays);
394  if(annFinalStates.count("H-") == 1)
395  ImportDecays("H-", catalog, importedDecays, tbl, minBranching, excludeDecays);
396  if(annFinalStates.count("h0_1") == 1)
397  ImportDecays("h0_1", catalog, importedDecays, tbl, minBranching, excludeDecays);
398  if(annFinalStates.count("h0_2") == 1)
399  ImportDecays("h0_2", catalog, importedDecays, tbl, minBranching, excludeDecays);
400  if(annFinalStates.count("A0") == 1)
401  ImportDecays("A0", catalog, importedDecays, tbl, minBranching, excludeDecays);
402 
403  // Add process to process list
404  catalog.processList.push_back(process);
405 
406  // Validate
407  catalog.validate();
408 
409  // Return the finished process catalog
410  result = catalog;
411  }
412 
413 
418  {
419  using namespace Pipes::TH_ProcessCatalog_DS_MSSM;
420  using std::vector;
421  using std::string;
422 
423  std::string DMid = *Dep::DarkMatter_ID;
424  if ( DMid != "~chi0_1" )
425  {
427  "TH_ProcessCatalog_DS_MSSM requires DMid to be ~chi0_1.");
428  }
429 
430  // Instantiate new empty ProcessCatalog
431  TH_ProcessCatalog catalog;
432 
433 
435  // Import particle masses
437 
438  // Import based on Spectrum objects
439  const Spectrum& matched_spectra = *Dep::MSSM_spectrum;
440  const SubSpectrum& spec = matched_spectra.get_HE();
441  const SubSpectrum& SM = matched_spectra.get_LE();
442  const SMInputs& SMI = matched_spectra.get_SMInputs();
443 
444  // Get SM masses
445  auto getSMmass = [&](str Name, int spinX2)
446  {
447  catalog.particleProperties.insert(
448  std::pair<std::string, TH_ParticleProperty>(
449  Name , TH_ParticleProperty(SM.get(Par::Pole_Mass,Name), spinX2)));
450  };
451 
452  getSMmass("e-_1", 1);
453  getSMmass("e+_1", 1);
454  getSMmass("e-_2", 1);
455  getSMmass("e+_2", 1);
456  getSMmass("e-_3", 1);
457  getSMmass("e+_3", 1);
458  getSMmass("Z0", 2);
459  getSMmass("W+", 2);
460  getSMmass("W-", 2);
461  getSMmass("g", 2);
462  getSMmass("gamma", 2);
463  getSMmass("d_3", 1);
464  getSMmass("dbar_3", 1);
465  getSMmass("u_3", 1);
466  getSMmass("ubar_3", 1);
467 
468  // Pole masses not available for the light quarks.
469  auto addParticle = [&](str Name, double Mass, int spinX2)
470  {
471  catalog.particleProperties.insert(
472  std::pair<std::string, TH_ParticleProperty>(
473  Name , TH_ParticleProperty(Mass, spinX2)));
474  };
475 
476  addParticle("d_1" , SMI.mD, 1); // md(2 GeV)^MS-bar, not pole mass
477  addParticle("dbar_1", SMI.mD, 1); // md(2 GeV)^MS-bar, not pole mass
478  addParticle("u_1" , SMI.mU, 1); // mu(2 GeV)^MS-bar, not pole mass
479  addParticle("ubar_1", SMI.mU, 1); // mu(2 GeV)^MS-bar, not pole mass
480  addParticle("d_2" , SMI.mS, 1); // ms(2 GeV)^MS-bar, not pole mass
481  addParticle("dbar_2", SMI.mS, 1); // ms(2 GeV)^MS-bar, not pole mass
482  addParticle("u_2" , SMI.mCmC,1); // mc(mc)^MS-bar, not pole mass
483  addParticle("ubar_2", SMI.mCmC,1); // mc(mc)^MS-bar, not pole mass
484  // Masses for neutrino flavour eigenstates. Set to zero.
485  addParticle("nu_e", 0.0, 1);
486  addParticle("nubar_e", 0.0, 1);
487  addParticle("nu_mu", 0.0, 1);
488  addParticle("nubar_mu", 0.0, 1);
489  addParticle("nu_tau", 0.0, 1);
490  addParticle("nubar_tau",0.0, 1);
491 
492  addParticle("pi0", meson_masses.pi0, 0);
493  addParticle("pi+", meson_masses.pi_plus, 0);
494  addParticle("pi-", meson_masses.pi_minus, 0);
495  addParticle("eta", meson_masses.eta, 0);
496  addParticle("rho0", meson_masses.rho0, 1);
497  addParticle("rho+", meson_masses.rho_plus, 1);
498  addParticle("rho-", meson_masses.rho_minus, 1);
499  addParticle("omega", meson_masses.omega, 1);
500 
501  // Get MSSM masses
502  auto getMSSMmass = [&](str Name, int spinX2)
503  {
504  catalog.particleProperties.insert(
505  std::pair<std::string, TH_ParticleProperty>(
506  Name , TH_ParticleProperty(std::abs(spec.get(Par::Pole_Mass,Name)), spinX2)));
507  };
508 
509  getMSSMmass("H+" , 0);
510  getMSSMmass("H-" , 0);
511  getMSSMmass("h0_1" , 0);
512  getMSSMmass("h0_2" , 0);
513  getMSSMmass("A0" , 0);
514  getMSSMmass("~chi0_1" , 1);
515  getMSSMmass("~d_1" , 0);
516  getMSSMmass("~dbar_1" , 0);
517  getMSSMmass("~u_1" , 0);
518  getMSSMmass("~ubar_1" , 0);
519  getMSSMmass("~d_2" , 0);
520  getMSSMmass("~dbar_2" , 0);
521  getMSSMmass("~u_2" , 0);
522  getMSSMmass("~ubar_2" , 0);
523  getMSSMmass("~d_3" , 0);
524  getMSSMmass("~dbar_3" , 0);
525  getMSSMmass("~u_3" , 0);
526  getMSSMmass("~ubar_3" , 0);
527  getMSSMmass("~d_4" , 0);
528  getMSSMmass("~dbar_4" , 0);
529  getMSSMmass("~u_4" , 0);
530  getMSSMmass("~ubar_4" , 0);
531  getMSSMmass("~d_5" , 0);
532  getMSSMmass("~dbar_5" , 0);
533  getMSSMmass("~u_5" , 0);
534  getMSSMmass("~ubar_5" , 0);
535  getMSSMmass("~d_6" , 0);
536  getMSSMmass("~dbar_6" , 0);
537  getMSSMmass("~u_6" , 0);
538  getMSSMmass("~ubar_6" , 0);
539  //getMSSMmass("~e_1" , 0);
540  //getMSSMmass("~ebar_1" , 0);
541  //getMSSMmass("~e-_1" , 0);
542  getMSSMmass("~e+_1" , 0);
543  getMSSMmass("~e-_1" , 0);
544  getMSSMmass("~e+_2" , 0);
545  getMSSMmass("~e-_2" , 0);
546  getMSSMmass("~e+_3" , 0);
547  getMSSMmass("~e-_3" , 0);
548  getMSSMmass("~e+_4" , 0);
549  getMSSMmass("~e-_4" , 0);
550  getMSSMmass("~e+_5" , 0);
551  getMSSMmass("~e-_5" , 0);
552  getMSSMmass("~e+_6" , 0);
553  getMSSMmass("~e-_6" , 0);
554  getMSSMmass("~nu_1" , 0);
555  getMSSMmass("~nubar_1", 0);
556  getMSSMmass("~nu_2" , 0);
557  getMSSMmass("~nubar_2", 0);
558  getMSSMmass("~nu_3" , 0);
559  getMSSMmass("~nubar_3", 0);
560 
562  // Import two-body annihilation process
564 
565  // Set of possible final state particles. Used to determine which decays to import.
566  std::set<string> annFinalStates;
567 
568  // Declare DM annihilation process
569  TH_Process process(DMid, DMid);
570 
571  // Explicitly state that Neutralino DM is self-conjugate
572  process.isSelfConj = true;
573 
574  double M_DM = catalog.getParticleProperty(DMid).mass;
575 
576  // lambda for setting up 2-body annihilations (chi chi -> X X) from results in DS
577  auto setup_DS6_process = [&](int pdg1, int pdg2, str p1, str p2, double prefac)
578  {
579  /* Check if process is kinematically allowed */
580  double m_1 = catalog.getParticleProperty(p1).mass;
581  double m_2 = catalog.getParticleProperty(p2).mass;
582  if(m_1 + m_2 < 2*M_DM)
583  {
584  /* Set cross-section */
585  double sigma = BEreq::dssigmav0(pdg1,pdg2);
586  /* Create associated kinematical functions (just dependent on vrel)
587  * here: s-wave, vrel independent 1-dim constant function */
588  daFunk::Funk kinematicFunction = daFunk::cnst(sigma*prefac, "v");
589  /* Create channel identifier string */
590  std::vector<std::string> finalStates;
591  finalStates.push_back(p1);
592  finalStates.push_back(p2);
593  /* Create channel and push it into channel list of process */
594  process.channelList.push_back(TH_Channel(finalStates, kinematicFunction));
595  annFinalStates.insert(p1);
596  annFinalStates.insert(p2);
597  }
598  };
599 
600  setup_DS6_process( 35, 35, "h0_2", "h0_2", 1 );
601  setup_DS6_process( 35, 25, "h0_2", "h0_1", 1 );
602  setup_DS6_process( 25, 25, "h0_1", "h0_1", 1 );
603  setup_DS6_process( 36, 36, "A0", "A0", 1 );
604  setup_DS6_process( 35, 36, "h0_2", "A0", 1 );
605  setup_DS6_process( 25, 36, "h0_1", "A0", 1 );
606  setup_DS6_process( 37, -37, "H+", "H-", 1 );
607  setup_DS6_process( 35, 23, "h0_2", "Z0", 1 );
608  setup_DS6_process( 25, 23, "h0_1", "Z0", 1 );
609  setup_DS6_process( 36, 23, "A0", "Z0", 1 );
610  // Prefactor 0.5 since W+H- and W-H+ are summed in DS
611  setup_DS6_process( 24, -37, "W+", "H-", 0.5 );
612  setup_DS6_process(-24, 37, "W-", "H+", 0.5 );
613  setup_DS6_process( 23, 23, "Z0", "Z0", 1 );
614  setup_DS6_process( 24, -24, "W+", "W-", 1 );
615  setup_DS6_process( 12, -12, "nu_e", "nubar_e",1 );
616  setup_DS6_process( 11, -11, "e+_1", "e-_1", 1 );
617  setup_DS6_process( 14, -14, "nu_mu", "nubar_mu",1 );
618  setup_DS6_process( 13, -13, "e+_2", "e-_2", 1 );
619  setup_DS6_process( 16, -16, "nu_tau", "nubar_tau",1 );
620  setup_DS6_process( 15, -15, "e+_3", "e-_3", 1 );
621  setup_DS6_process( 2, - 2, "u_1", "ubar_1", 1 );
622  setup_DS6_process( 1, - 1, "d_1", "dbar_1", 1 );
623  setup_DS6_process( 4, - 4, "u_2", "ubar_2", 1 );
624  setup_DS6_process( 3, - 3, "d_2", "dbar_2", 1 );
625  setup_DS6_process( 6, - 6, "u_3", "ubar_3", 1 );
626  setup_DS6_process( 5, - 5, "d_3", "dbar_3", 1 );
627  setup_DS6_process( 21, 21, "g", "g", 1 );
628  setup_DS6_process( 22, 22, "gamma", "gamma", 1 );
629  setup_DS6_process( 23, 22, "Z0", "gamma", 1 );
630 
631 
633  // Import three-body annihilation process
635 
637 
638  //PS: commented out for now, as this can't be a backend function in its current form.
639  //BEreq::registerMassesForIB(catalog.particleProperties);
640 
641  // Macro for setting up 3-body decays with gammas
642 
643  auto setup_DS6_process_gamma3body = [&](int IBch, str p1, str p2, double (*IBfunc)(int&, double&, double&), int sv_pdg1, int sv_pdg2, double prefac)
644  {
645  /* Check if process is kinematically allowed */
646  double m_1 = catalog.getParticleProperty(str_flav_to_mass(p1)).mass;
647  double m_2 = catalog.getParticleProperty(str_flav_to_mass(p2)).mass;
648  if(m_1 + m_2 < 2*M_DM)
649  {
650  double sv;
651  if(sv_pdg1==0 && sv_pdg2==0)
652  {
653  sv = prefac*BEreq::dssigmav0tot();
654  }
655  else
656  {
657  sv = prefac*BEreq::dssigmav0(sv_pdg1,sv_pdg2);
658  }
659  daFunk::Funk kinematicFunction = daFunk::cnst(sv,"v")*daFunk::func(DSgamma3bdy,
660  IBfunc, IBch, daFunk::var("E"), daFunk::var("E1"), M_DM, m_1, m_2);
661  /* Create channel identifier string */
662  std::vector<std::string> finalStates;
663  finalStates.push_back("gamma");
664  finalStates.push_back(str_flav_to_mass(p1));
665  finalStates.push_back(str_flav_to_mass(p2));
666  /* Create channel and push it into channel list of process */
667  process.channelList.push_back(TH_Channel(finalStates, kinematicFunction));
668  annFinalStates.insert(str_flav_to_mass(p1));
669  annFinalStates.insert(str_flav_to_mass(p2));
670  }
671  };
672 
674  if ( not runOptions->getValueOrDef<bool>(false, "ignore_three_body") )
675  {
676  // Set DarkSUSY DM mass parameter used in 3-body decays
677  BEreq::IBintvars->ibcom_mx = catalog.getParticleProperty(DMid).mass;
678 
679  setup_DS6_process_gamma3body( 1, "W+", "W-", BEreq::dsIBwwdxdy.pointer(),24, -24, 1 );
680  // Prefactor 0.5 since W+H- and W-H+ are summed in DS
681  setup_DS6_process_gamma3body( 2, "W+", "H-", BEreq::dsIBwhdxdy.pointer(),24, -37, 0.5);
682  // Prefactor 0.5 since W+H- and W-H+ are summed in DS
683  setup_DS6_process_gamma3body( 2, "W-", "H+", BEreq::dsIBwhdxdy.pointer(),37, -24, 0.5);
684  setup_DS6_process_gamma3body( 3, "H+", "H-", BEreq::dsIBhhdxdy.pointer(), 0, 0, 1 );
685 
686  setup_DS6_process_gamma3body( 4, "e+", "e-", BEreq::dsIBffdxdy.pointer(),11, -11, 1 );
687  setup_DS6_process_gamma3body( 5, "mu+", "mu-", BEreq::dsIBffdxdy.pointer(),13, -13, 1 );
688  setup_DS6_process_gamma3body( 6, "tau+","tau-",BEreq::dsIBffdxdy.pointer(),15, -15, 1 );
689  setup_DS6_process_gamma3body( 7, "u", "ubar",BEreq::dsIBffdxdy.pointer(), 2, -2, 1 );
690  setup_DS6_process_gamma3body( 8, "d", "dbar",BEreq::dsIBffdxdy.pointer(), 1, -1, 1 );
691  setup_DS6_process_gamma3body( 9, "c", "cbar",BEreq::dsIBffdxdy.pointer(), 4, -4, 1 );
692  setup_DS6_process_gamma3body(10,"s", "sbar",BEreq::dsIBffdxdy.pointer(), 3, -3, 1 );
693  setup_DS6_process_gamma3body(11,"t", "tbar",BEreq::dsIBffdxdy.pointer(), 6, -6, 1 );
694  setup_DS6_process_gamma3body(12,"b", "bbar",BEreq::dsIBffdxdy.pointer(), 5, -5, 1 );
695  };
696 
697 
699  // Import Decay information
701 
702  // Import based on decay table from DecayBit
703  const DecayTable* tbl = &(*Dep::decay_rates);
704 
705  // Set of imported decays - avoids double imports
706  std::set<string> importedDecays;
707 
708  // Option minBranching <double>: Minimum branching fraction of included
709  // processes (default 0.)
710  double minBranching = runOptions->getValueOrDef<double>(0.0,
711  "ProcessCatalog_MinBranching");
712 
713  // Exclude also ttbar final states
714  auto excludeDecays = daFunk::vec<std::string>("Z0", "W+", "W-", "u_3", "ubar_3", "e+_2", "e-_2", "e+_3", "e-_3");
715 
716 
717  // Import relevant decays
719  if(annFinalStates.count("H+") == 1)
720  ImportDecays("H+", catalog, importedDecays, tbl, minBranching, excludeDecays);
721  if(annFinalStates.count("H-") == 1)
722  ImportDecays("H-", catalog, importedDecays, tbl, minBranching, excludeDecays);
723  if(annFinalStates.count("h0_1") == 1)
724  ImportDecays("h0_1", catalog, importedDecays, tbl, minBranching, excludeDecays);
725  if(annFinalStates.count("h0_2") == 1)
726  ImportDecays("h0_2", catalog, importedDecays, tbl, minBranching, excludeDecays);
727  if(annFinalStates.count("A0") == 1)
728  ImportDecays("A0", catalog, importedDecays, tbl, minBranching, excludeDecays);
729 
730 
731  // Add process to process list
732  catalog.processList.push_back(process);
733 
734  // Validate
735  catalog.validate();
736 
737  // Return the finished process catalog
738  result = catalog;
739  } //TH_ProcessCatalog_DS_MSSM
740 
741 
742 
743  void DarkMatter_ID_MSSM(std::string & result)
744  {
745  using namespace Pipes::DarkMatter_ID_MSSM;
746  // TODO: need ask Dep::MSSM_spectrum in future; might have sneutralinos and gravitinos later.
747  result = "~chi0_1";
748  }
749  }
750 }
std::vector< TH_Channel > channelList
List of channels.
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
double DSgamma3bdy(double(*IBfunc)(int &, double &, double &), int IBch, double Eg, double E1, double M_DM, double m_1, double m_2)
Fully initialize DarkSUSY to the current model point.
Definition: MSSM.cpp:61
TH_ParticleProperty getParticleProperty(str) const
Retrieve properties of a given particle involved in one or more processes in this catalog...
void validate()
Validate kinematics and entries.
void ImportDecays(std::string pID, TH_ProcessCatalog &catalog, std::set< std::string > &importedDecays, const DecayTable *tbl, double minBranching, std::vector< std::string > excludeDecays=std::vector< std::string >())
Funk var(std::string arg)
Definition: daFunk.hpp:921
bool isSelfConj
Does the process contain self-conjugate DM? (accounting for correct factors of 1/2 in annihilation sp...
double mass
Particle mass (GeV)
const double sigma
Definition: SM_Z.hpp:43
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
virtual double get(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
std::map< std::string, TH_ParticleProperty > particleProperties
Map from particles involved in the processes of this catalog, to their properties.
A container holding all annihilation and decay initial states relevant for DarkBit.
Funk func(double(*f)(funcargs...), Args... args)
Definition: daFunk.hpp:768
Header file that includes all GAMBIT headers required for a module source file.
A container for the mass and spin of a particle.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
SubSpectrum & get_LE()
Standard getters Return references to internal data members.
Definition: spectrum.cpp:224
All data on a single annihilation or decay channel, e.g.
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
Rollcall header for module DarkBit.
A simple C++ wrapper for the MPI C bindings.
void TH_ProcessCatalog_DS_MSSM(DarkBit::TH_ProcessCatalog &result)
Initialization of Process Catalog based on DarkSUSY 6 calculations.
Definition: MSSM.cpp:417
GAMBIT native decay table class.
Definition: decay_table.hpp:35
Funk cnst(double x, Args... argss)
Definition: daFunk.hpp:606
void DarkMatter_ID_MSSM(std::string &result)
Definition: MSSM.cpp:743
std::string str_flav_to_mass(std::string flav)
invalid_point_exception & invalid_point()
Invalid point exceptions.
A container for a single process.
#define getSMmass(Name, spinX2)
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
std::vector< TH_Process > processList
Vector of all processes in this catalog.
void TH_ProcessCatalog_DS5_MSSM(DarkBit::TH_ProcessCatalog &result)
Initialization of Process Catalog based on DarkSUSY 5 calculations.
Definition: MSSM.cpp:100
#define addParticle(Name, Mass, spinX2)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
SubSpectrum & get_HE()
Definition: spectrum.cpp:225
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110
SMInputs & get_SMInputs()
Definition: spectrum.cpp:226
Utility functions for DarkBit.
Container class for Standard Model input information (defined as in SLHA2)
Definition: sminputs.hpp:29