gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
lep_mssm_xsecs.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
24 
25 
26 #include <iostream>
27 #include <fstream>
28 #include <math.h>
29 
34 
35 
36 #define pow2(a) ((a)*(a)) // Get speedy
37 
38 namespace Gambit
39 {
40 
41  namespace ColliderBit
42  {
43 
49  void get_sigma_ee_ll(triplet<double>& result, const double sqrts, const int generation, const int l_chirality,
50  const int lbar_chirality, const double gtol, const double ftol, const bool gpt_error,
51  const bool fpt_error, const Spectrum& spec, const double gammaZ, const bool l_are_gauge_es)
52  {
53  static const str genmap[3][2] =
54  {
55  {"~e_L", "~e_R" },
56  {"~mu_L", "~mu_R" },
57  {"~tau_L", "~tau_R"}
58  };
59 
60  // Subspectrum
61  const SubSpectrum& mssm = spec.get_HE();
62  // PDG codes
63  const int id1 = 1000000*l_chirality + 11 +2*(generation-1);
64  const int id2 = -(1000000*lbar_chirality + 11 +2*(generation-1));
65 
66  // SM parameters
67  const double mZ = spec.safeget(Par::Pole_Mass,23,0);
68  const double g2 = mssm.safeget(Par::dimensionless,"g2");
69  const double sinW2 = mssm.safeget(Par::dimensionless,"sinW2");
70  const double alpha = 0.25*sinW2*g2*g2/pi;
71  // MSSM parameters
72  const double tanb = mssm.safeget(Par::dimensionless,"tanbeta");
73 
74  // Get the mass eigenstate strings and 2x2 slepton generation mass mixing matrix
75  str mass_es1, mass_es2;
76  MixMatrix sleptonmix(2,std::vector<double>(2));
77 
78  if (l_are_gauge_es)
79  {
80  // Requested final states are gauge eigenstates. Pass diagonal mixing matrix to low-level routine.
81  sleptonmix[0][0] = 1.0;
82  sleptonmix[0][1] = 0.0;
83  sleptonmix[1][0] = 0.0;
84  sleptonmix[1][1] = 1.0;
85  mass_es1 = slhahelp::mass_es_from_gauge_es(genmap[generation-1][l_chirality-1], mssm, gtol, LOCAL_INFO, gpt_error);
86  mass_es2 = slhahelp::mass_es_from_gauge_es(genmap[generation-1][lbar_chirality-1], mssm, gtol, LOCAL_INFO, gpt_error);
87  }
88  else
89  {
90  // Requested final states are family mass eigenstates. Pass 2x2 family mass mixing matrix to low-level routine.
91  str m_light, m_heavy;
92  std::vector<double> slepton4vec = slhahelp::family_state_mix_matrix("~e-", generation, m_light, m_heavy, mssm, ftol, LOCAL_INFO, fpt_error);
93  mass_es1 = (l_chirality == 1) ? m_light : m_heavy;
94  mass_es2 = (lbar_chirality == 1) ? m_light : m_heavy;
95  sleptonmix[0][0] = slepton4vec[0];
96  sleptonmix[0][1] = slepton4vec[1];
97  sleptonmix[1][0] = slepton4vec[2];
98  sleptonmix[1][1] = slepton4vec[3];
99  }
100  const double m1 = spec.safeget(Par::Pole_Mass,mass_es1);
101  // FIXME when spectrum object has separate pole mass getters for antiparticles
102  //const double m2 = spec.safeget(Par::Pole_Mass,Models::ParticleDB().get_antiparticle(mass_es2));
103  // until then
104  const double m2 = spec.safeget(Par::Pole_Mass,mass_es2);
105  std::pair<double,double> m1_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, mass_es1),
106  mssm.safeget(Par::Pole_Mass_1srd_low, mass_es1));
107  std::pair<double,double> m2_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, mass_es2),
108  mssm.safeget(Par::Pole_Mass_1srd_low, mass_es2));
109 
110  // If the final state is kinematically inaccessible *even* if both masses
111  // are 2simga lower than their central values, then return zero.
112  if (m1*(1.0-2.0*m1_uncerts.second) + m2*(1.0-2.0*m2_uncerts.second) > sqrts)
113  {
114  result.central = 0.0;
115  result.upper = 0.0;
116  result.lower = 0.0;
117  return;
118  }
119 
120  // Get the neutralino masses
121  const double neutmass[4] = { spec.safeget(Par::Pole_Mass,1000022,0), spec.safeget(Par::Pole_Mass,1000023,0),
122  spec.safeget(Par::Pole_Mass,1000025,0), spec.safeget(Par::Pole_Mass,1000035,0) };
123  // Get the 4x4 neutralino mixing matrix
124  MixMatrix neutmix(4,std::vector<double>(4));
125  for (int i=0; i<4; i++) for (int j=0; j<4; j++) neutmix[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi0",i+1,j+1);
126 
127  // Convert neutralino mixing matrix to BFM convention
128  SLHA2BFM_NN(neutmix, tanb, sinW2);
129 
130  // Calculate the cross-section
131  result.central = xsec_sleislej(id1, id2, sqrts, m1, m2, sleptonmix, neutmix, neutmass, alpha, mZ, gammaZ, sinW2);
132 
133  // Calculate the uncertainty on the cross-section due to final state masses varying by +/- 1 sigma
134  std::vector<double> xsecs;
135  xsecs.push_back(result.central);
136  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.+m2_uncerts.first), sleptonmix, neutmix,
137  neutmass, alpha, mZ, gammaZ, sinW2, false));
138  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.+m2_uncerts.first), sleptonmix, neutmix,
139  neutmass, alpha, mZ, gammaZ, sinW2, false));
140  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.-m2_uncerts.second), sleptonmix, neutmix,
141  neutmass, alpha, mZ, gammaZ, sinW2, false));
142  xsecs.push_back(xsec_sleislej(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.-m2_uncerts.second), sleptonmix, neutmix,
143  neutmass, alpha, mZ, gammaZ, sinW2, false));
144  result.upper = *std::max_element(xsecs.begin(), xsecs.end());
145  result.lower = *std::min_element(xsecs.begin(), xsecs.end());
146 
147  }
148 
149 
151  void get_sigma_ee_chi00(triplet<double>& result, const double sqrts, const int chi_first, const int chi_second,
152  const double tol, const bool pt_error, const Spectrum& spec, const double gammaZ)
153  {
154  // Subspectrum
155  const SubSpectrum& mssm = spec.get_HE();
156 
157  // PDG codes
158  const int id1 = 1000021 + chi_first + (chi_first > 2 ? 1 + (chi_first -3)*9 : 0);
159  const int id2 = 1000021 + chi_second + (chi_second > 2 ? 1 + (chi_second-3)*9 : 0);
160 
161  // SM parameters
162  const double mZ = spec.safeget(Par::Pole_Mass,23,0);
163  const double g2 = mssm.safeget(Par::dimensionless,"g2");
164  const double sinW2 = mssm.safeget(Par::dimensionless,"sinW2");
165  const double alpha = 0.25*sinW2*g2*g2/pi;
166 
167  // MSSM parameters
168  const double tanb = mssm.safeget(Par::dimensionless,"tanbeta");
169  // Get the mass eigenstates best corresponding to ~eL and ~eR.
170  const str mass_esL = slhahelp::mass_es_from_gauge_es("~e_L", mssm, tol, LOCAL_INFO, pt_error);
171  const str mass_esR = slhahelp::mass_es_from_gauge_es("~e_R", mssm, tol, LOCAL_INFO, pt_error);
172  // Get the slepton masses
173  const double mS[2] = {spec.safeget(Par::Pole_Mass,mass_esL), spec.safeget(Par::Pole_Mass,mass_esR)};
174  // Get the neutralino masses
175  const double m1 = spec.safeget(Par::Pole_Mass,id1,0);
176  const double m2 = spec.safeget(Par::Pole_Mass,id2,0);
177  std::pair<double,double> m1_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id1, 0),
178  mssm.safeget(Par::Pole_Mass_1srd_low, id1, 0));
179  std::pair<double,double> m2_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id2, 0),
180  mssm.safeget(Par::Pole_Mass_1srd_low, id2, 0));
181 
182  // Just return zero if the final state is kinematically inaccessible
183  // *even* if both masses are 2simga lower than their central values
184  if (std::abs(m1)*(1.0-2.0*m1_uncerts.second) + std::abs(m2)*(1.0-2.0*m2_uncerts.second) > sqrts)
185  {
186  result.central = 0.0;
187  result.upper = 0.0;
188  result.lower = 0.0;
189  return;
190  }
191 
192  // Get the 4x4 neutralino mixing matrix
193  MixMatrix neutmix(4,std::vector<double>(4));
194  for (int i=0; i<4; i++) for (int j=0; j<4; j++) neutmix[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi0",i+1,j+1);
195 
196  // Convert neutralino mixing matrix to BFM convention
197  SLHA2BFM_NN(neutmix, tanb, sinW2);
198 
199  // Calculate the cross-section
200  result.central = xsec_neuineuj(id1, id2, sqrts, m1, m2, neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2);
201 
202  // Calculate the uncertainty on the cross-section due to final state masses varying by +/- 1 sigma
203  std::vector<double> xsecs;
204  xsecs.push_back(result.central);
205  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.+m2_uncerts.first),
206  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
207  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.-m2_uncerts.second),
208  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
209  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.+m2_uncerts.first),
210  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
211  xsecs.push_back(xsec_neuineuj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.-m2_uncerts.second),
212  neutmix, mS, 1./tanb, alpha, mZ, gammaZ, sinW2));
213  result.upper = *std::max_element(xsecs.begin(), xsecs.end());
214  result.lower = *std::min_element(xsecs.begin(), xsecs.end());
215 
216  }
217 
219  void get_sigma_ee_chipm(triplet<double>& result, const double sqrts, const int chi_plus, const int chi_minus,
220  const double tol, const bool pt_error, const Spectrum& spec, const double gammaZ)
221  {
222  // Subspectrum
223  const SubSpectrum& mssm = spec.get_HE();
224 
225  // PDG codes
226  const int id1 = 1000023 + chi_plus + (chi_plus - 1)*12;
227  const int id2 = -(1000023 + chi_minus + (chi_minus - 1)*12);
228 
229  // SM parameters
230  const double mZ = spec.safeget(Par::Pole_Mass,23,0);
231  const double g2 = mssm.safeget(Par::dimensionless,"g2");
232  const double sinW2 = mssm.safeget(Par::dimensionless,"sinW2");
233  const double alpha = 0.25*sinW2*g2*g2/pi;
234 
235  // MSSM parameters
236  // Get the mass eigenstate best corresponding to ~nu_e_L.
237  const str mass_snue = slhahelp::mass_es_from_gauge_es("~nu_e_L", mssm, tol, LOCAL_INFO, pt_error);
238  // Get the electron sneutrino mass
239  const double msn = spec.safeget(Par::Pole_Mass,mass_snue);
240  // Get the chargino masses
241  const double m1 = spec.safeget(Par::Pole_Mass,id1,0);
242  const double m2 = spec.safeget(Par::Pole_Mass,id2,0);
243  std::pair<double,double> m1_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id1, 0),
244  mssm.safeget(Par::Pole_Mass_1srd_low, id1, 0));
245  std::pair<double,double> m2_uncerts(mssm.safeget(Par::Pole_Mass_1srd_high, id2, 0),
246  mssm.safeget(Par::Pole_Mass_1srd_low, id2, 0));
247 
248  // Just return zero if the final state is kinematically inaccessible
249  // *even* if both masses are 2simga lower than their central values
250  if (std::abs(m1)*(1.0-2.0*m1_uncerts.second) + std::abs(m2)*(1.0-2.0*m2_uncerts.second) > sqrts)
251  {
252  result.central = 0.0;
253  result.upper = 0.0;
254  result.lower = 0.0;
255  return;
256  }
257 
258  // Get the 2x2 chargino mixing matrices
259  MixMatrix charginomixV(2,std::vector<double>(2));
260  MixMatrix charginomixU(2,std::vector<double>(2));
261  for (int i=0; i<2; i++) for (int j=0; j<2; j++)
262  {
263  charginomixV[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi+",i+1,j+1);
264  charginomixU[i][j] = mssm.safeget(Par::Pole_Mixing,"~chi-",i+1,j+1);
265  }
266 
267  // Convert chargino mixing matrices to BFM convention
268  SLHA2BFM_VV(charginomixV);
269  SLHA2BFM_VV(charginomixU);
270 
271  // Calculate the cross-section
272  result.central = xsec_chaichaj(id1, id2, sqrts, m1, m2, charginomixV, charginomixU,
273  msn, alpha, mZ, gammaZ, sinW2);
274 
275  // Calculate the uncertainty on the cross-section due to final state masses varying by +/- 1 sigma
276  std::vector<double> xsecs;
277  xsecs.push_back(result.central);
278  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.+m2_uncerts.first), charginomixV, charginomixU,
279  msn, alpha, mZ, gammaZ, sinW2));
280  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.+m1_uncerts.first), m2*(1.-m2_uncerts.second), charginomixV, charginomixU,
281  msn, alpha, mZ, gammaZ, sinW2));
282  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.+m2_uncerts.first), charginomixV, charginomixU,
283  msn, alpha, mZ, gammaZ, sinW2));
284  xsecs.push_back(xsec_chaichaj(id1, id2, sqrts, m1*(1.-m1_uncerts.second), m2*(1.-m2_uncerts.second), charginomixV, charginomixU,
285  msn, alpha, mZ, gammaZ, sinW2));
286  result.upper = *std::max_element(xsecs.begin(), xsecs.end());
287  result.lower = *std::min_element(xsecs.begin(), xsecs.end());
288 
289  }
290 
295  double I1(double s, double m1, double m2, double mk, double ml)
296  {
297  double S = sqrt(s-pow2(m1+m2))*sqrt(s-pow2(m1-m2));
298  double m1sq = pow2(m1);
299  double m2sq = pow2(m2);
300  double mksq = pow2(mk);
301  double mlsq = pow2(ml);
302 
303  double I1 = 0;
304  // Careful with degenerate masses!
305  if( fabs(mksq-mlsq) < 0.1 ){
306  I1 = (m1sq+m2sq-2.*mksq-s)*log((m1sq+m2sq-2.*mksq-s+S)/(m1sq+m2sq-2.*mksq-s-S))-4*S*((m1sq-mksq)*(m2sq-mksq)+mksq*s)/(m1sq+m2sq-2.*mksq-s-S)/(m1sq+m2sq-2.*mksq-s+S)-S;
307  }
308  else{
309  I1 = (m1sq*(m2sq-mksq)+mksq*(mksq-m2sq+s))/(mksq-mlsq)*log((m1sq+m2sq-2.*mksq-s-S)/(m1sq+m2sq-2.*mksq-s+S));
310  I1 += (m1sq*(m2sq-mlsq)+mlsq*(mlsq-m2sq+s))/(mlsq-mksq)*log((m1sq+m2sq-2.*mlsq-s-S)/(m1sq+m2sq-2.*mlsq-s+S));
311  I1 -= S;
312  }
313  return I1;
314  }
315  double I2(double s, double m1, double m2, double mk, double ml)
316  {
317  double S = sqrt(s-pow2(m1+m2))*sqrt(s-pow2(m1-m2));
318  double m1sq = pow2(m1);
319  double m2sq = pow2(m2);
320  double mksq = pow2(mk);
321  double mlsq = pow2(ml);
322 
323  double I2 = 0;
324  // Careful with degenerate masses!
325  if( fabs(mksq-mlsq) < 0.1 )
326  {
327  I2 = S/(m1sq*(m2sq-mksq)+mksq*(-m2sq+mksq+s));
328  }
329  else
330  {
331  I2 = log((m1sq+m2sq-2.*mksq-(s-S))/(m1sq+m2sq-2.*mksq-(s+S)));
332  I2 += log((m1sq+m2sq-2.*mlsq-(s+S))/(m1sq+m2sq-2.*mlsq-(s-S)));
333  I2 *= 1./(mksq-mlsq);
334  }
335  return I2;
336  }
337  double I3(double s, double m1, double m2, double mk)
338  {
339  double S = sqrt(s-pow2(m1+m2))*sqrt(s-pow2(m1-m2));
340  double m1sq = pow2(m1);
341  double m2sq = pow2(m2);
342  double mksq = pow2(mk);
343 
344  double I3 = 0;
345  I3 = log((m1sq+m2sq-2.*mksq-(s+S))/(m1sq+m2sq-2.*mksq-(s-S)));
346  I3 *= m1sq*(m2sq-mksq)+mksq*(-m2sq+mksq+s);
347  I3 += (m1sq+m2sq-2.*mksq-s)*S/2.;
348  return I3;
349  }
351 
352 
355  double xsec_sleislej(int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix F,
356  MixMatrix N, const double mN[4], double alpha, double mZ, double gZ,
357  double sin2thetaW, bool CP_lock)
358  {
359 
360  // Just return zero if the final state isn't kinematically accessible
361  if (m1+m2 > sqrts) return 0.0;
362 
363  // Slepton mixing
364  double cosphi = F[0][0];
365  double sinphi = F[0][1];
366 
367  // Figure out what we are calculating
368  double tempphi;
369  bool bMixed = false;
370  bool bSelectron = false;
371  // ~e_L ~e_L^*
372  if((pid1 == 1000011 && pid2 == -1000011) || (pid1 == -1000011 && pid2 == 1000011)){
373  bSelectron = true;
374  if(m1 != m2 and CP_lock) ColliderBit_warning().raise(LOCAL_INFO, "You are using a different mass for antiparticle!");
375  }
376  // ~e_L ~e_R^*
377  else if((pid1 == 1000011 && pid2 == -2000011) || (pid1 == -2000011 && pid2 == 1000011)){
378  bSelectron = true;
379  bMixed = true;
380  }
381  // ~e_R ~e_L^*
382  else if((pid1 == 2000011 && pid2 == -1000011) || (pid1 == -1000011 && pid2 == 2000011)){
383  bSelectron = true;
384  bMixed = true;
385  }
386  // ~e_R ~e_R^*
387  else if((pid1 == 2000011 && pid2 == -2000011) || (pid1 == -2000011 && pid2 == 2000011)){
388  bSelectron = true;
389  tempphi = cosphi; cosphi = sinphi; sinphi = tempphi;
390  if(m1 != m2 and CP_lock) ColliderBit_warning().raise(LOCAL_INFO, "You are using a different mass for antiparticle!");
391  }
392  // ~mu_L ~mu_L^*
393  else if((pid1 == 1000013 && pid2 == -1000013) || (pid1 == -1000013 && pid2 == 1000013)){
394  if(m1 != m2 and CP_lock) ColliderBit_warning().raise(LOCAL_INFO, "You are using a different mass for antiparticle!");
395  }
396  // ~mu_L ~mu_R^*
397  else if((pid1 == 1000013 && pid2 == -2000013) || (pid1 == -2000013 && pid2 == 1000013)){
398  bMixed = true;
399  ColliderBit_warning().raise(LOCAL_INFO, "Will give zero cross section unless there is left-right smuon mixing!");
400  }
401  // ~mu_R ~mu_L^*
402  else if((pid1 == 2000013 && pid2 == -1000013) || (pid1 == -1000013 && pid2 == 2000013)){
403  bMixed = true;
404  ColliderBit_warning().raise(LOCAL_INFO, "Will give zero cross section unless there is left-right smuon mixing!");
405  }
406  // ~mu_R ~mu_R^*
407  else if((pid1 == 2000013 && pid2 == -2000013) || (pid1 == -2000013 && pid2 == 2000013)){
408  tempphi = cosphi; cosphi = sinphi; sinphi = tempphi;
409  if(m1 != m2 and CP_lock) ColliderBit_warning().raise(LOCAL_INFO, "You are using a different mass for antiparticle!");
410  }
411  // ~tau_1 ~tau_1^*
412  else if((pid1 == 1000015 && pid2 == -1000015) || (pid1 == -1000015 && pid2 == 1000015)){
413  if(m1 != m2 and CP_lock) ColliderBit_warning().raise(LOCAL_INFO, "You are using a different mass for antiparticle!");
414  }
415  // ~tau_1 ~tau_2^*
416  else if((pid1 == 1000015 && pid2 == -2000015) || (pid1 == -2000015 && pid2 == 1000015)){
417  bMixed = true;
418  }
419  // ~tau_2 ~tau_1^*
420  else if((pid1 == 2000015 && pid2 == -1000015) || (pid1 == -1000015 && pid2 == 2000015)){
421  bMixed = true;
422  }
423  // ~tau_2 ~tau_2^*
424  else if((pid1 == 2000015 && pid2 == -2000015) || (pid1 == -2000015 && pid2 == 2000015)){
425  tempphi = cosphi; cosphi = sinphi; sinphi = tempphi;
426  if(m1 != m2 and CP_lock) ColliderBit_warning().raise(LOCAL_INFO, "You are using a different mass for antiparticle!");
427  }
428  else
429  {
430  std::stringstream ss;
431  ss << "I don't know that process!" << endl
432  << "You asked me to calculate slepton cross section with final states"
433  << "PID1 " << pid1 << " PID2 " << pid2;
434  ColliderBit_warning().raise(LOCAL_INFO, ss.str());
435  return -1;
436  }
437 
438  // Couplings
439  double T3l = -0.5;
440  double Le = T3l+sin2thetaW;
441  double Re = sin2thetaW;
442  // Left-right mixing
443  double cos2phi = pow2(cosphi);
444  double sin2phi = pow2(sinphi);
445 
446  double fL[4], fR[4];
447  for(int k = 0; k < 4; k++){
448  fL[k] = -sqrt(2.) * (1./sqrt(1.-sin2thetaW)*(T3l+sin2thetaW)*N[k][1]-sqrt(sin2thetaW)*N[k][0]);
449  fR[k] = sqrt(2.) * sqrt(sin2thetaW) * (sqrt(sin2thetaW/(1.-sin2thetaW))*N[k][1]-N[k][0]);
450  }
451 
452  // Kinematics
453  double s, S, DZ2, ReDZ;
454  s = pow2(sqrts);
455  S = sqrt(s-pow2(m1+m2))*sqrt(s-pow2(m1-m2));
456  DZ2 = 1./(pow2(s-pow2(mZ))+pow2(mZ*gZ)); // Breit-Wigner for Z
457  ReDZ = (s-pow2(mZ))*DZ2;
458 
459  // Cross sections per diagram and interference terms
460  double sigma, sigma_Z, sigma_Z_mix, sigma_g, sigma_gZ, sigma_N, sigma_N_mix, sigma_gN, sigma_ZN, sigma_ZN_mix;
461  // gamma
462  sigma_g = 2.*pi*pow2(alpha)/pow(s,4) * pow(S,3)/6.;
463  // Z
464  sigma_Z = pi*pow2(alpha)/pow2(s)/pow2(sin2thetaW)/pow2(1.-sin2thetaW) * DZ2 * pow(S,3)/6.;
465  sigma_Z *= (pow2(Le)+pow2(Re))*pow2(Le*cos2phi+Re*sin2phi);
466  sigma_Z_mix = sigma_Z/pow2(Le*cos2phi+Re*sin2phi)*pow2(Le-Re)*cos2phi*sin2phi;
467  // Interference
468  sigma_gZ = 2*pi*pow2(alpha)/pow(s,3)/sin2thetaW/(1.-sin2thetaW) * ReDZ;
469  sigma_gZ *= (Le+Re)*(Le*cos2phi+Re*sin2phi) * pow(S,3)/6.;
470  // Neutralino
471  // Loop over neutralinos
472  sigma_N = 0;
473  for(int k = 0; k < 4; k++){
474  for(int l = 0; l < 4; l++){
475  sigma_N += pow2(cos2phi)*I1(s,m1,m2,mN[k],mN[l])*pow2(fL[k]*fL[l]);
476  sigma_N += pow2(sin2phi)*I1(s,m1,m2,mN[k],mN[l])*pow2(fR[k]*fR[l]);
477  sigma_N += 2.*cos2phi*sin2phi*s*mN[k]*mN[l]*I2(s,m1,m2,mN[k],mN[l])*fL[k]*fL[l]*fR[k]*fR[l];
478  }
479  }
480  sigma_N *= pi*pow2(alpha)/4./pow2(sin2thetaW)/pow2(s);
481  sigma_N_mix = 0;
482  for(int k = 0; k < 4; k++)
483  {
484  for(int l = 0; l < 4; l++)
485  {
486  sigma_N_mix += cos2phi*sin2phi*I1(s,m1,m2,mN[k],mN[l])*pow2(fL[k]*fL[l]);
487  sigma_N_mix += cos2phi*sin2phi*I1(s,m1,m2,mN[k],mN[l])*pow2(fR[k]*fR[l]);
488  sigma_N_mix += (pow2(cos2phi)+pow2(sin2phi))*s*mN[k]*mN[l]*I2(s,m1,m2,mN[k],mN[l])*fL[k]*fL[l]*fR[k]*fR[l];
489  }
490  }
491  sigma_N_mix *= pi*pow2(alpha)/4./pow2(sin2thetaW)/pow2(s);
492  // Neutralino interference terms
493  sigma_gN = 0;
494  for(int k = 0; k < 4; k++){
495  sigma_gN += I3(s,m1,m2,mN[k])*(cos2phi*pow2(fL[k])+sin2phi*pow2(fR[k]));
496  }
497  sigma_gN *= pi*pow2(alpha)/sin2thetaW/pow(s,3);
498  sigma_ZN = 0;
499  for(int k = 0; k < 4; k++){
500  sigma_ZN += I3(s,m1,m2,mN[k])*(Le*cos2phi*pow2(fL[k])+Re*sin2phi*pow2(fR[k]));
501  }
502  sigma_ZN *= pi*pow2(alpha)/pow2(sin2thetaW)/(1.-sin2thetaW)/pow2(s)*(Le*cos2phi+Re*sin2phi)*ReDZ;
503  sigma_ZN_mix = 0;
504  for(int k = 0; k < 4; k++){
505  sigma_ZN_mix += I3(s,m1,m2,mN[k])*(Le*pow2(fL[k])-Re*pow2(fR[k]));
506  }
507  sigma_ZN_mix *= pi*pow2(alpha)/pow2(sin2thetaW)/(1.-sin2thetaW)/pow2(s)*sin2phi*cos2phi*(Le-Re)*ReDZ;
508 
509  // Total cross section
510  if( bMixed ) { sigma = sigma_Z_mix; }
511  else { sigma = sigma_g + sigma_Z + sigma_gZ; }
512  if( bSelectron && !bMixed ) { sigma += sigma_N + sigma_gN + sigma_ZN; }
513  else if( bSelectron && bMixed ) { sigma += sigma_N_mix + sigma_ZN_mix; }
514 
515  // Fix units
516  sigma *= gev2pb;
517 
518  // Return zero in corner cases where numerical roundoff has sent sigma negative.
519  return std::max(sigma, 0.0);
520 
521  }
522 
526  double xsec_neuineuj(int pid1, int pid2, double sqrts, double mi, double mj, MixMatrix N,
527  const double mS[2], double tanb, double alpha, double mZ, double gZ, double sin2thetaW)
528  {
529 
530  // Just return zero if the final state isn't kinematically accessible
531  if (std::abs(mi)+std::abs(mj) > sqrts) return 0.0;
532 
533  // Translate from PDG codes to neutralino indices (starting at zero)
534  int i, j;
535  if(pid1 == 1000022) i = 0;
536  else if(pid1 == 1000023) i = 1;
537  else if(pid1 == 1000025) i = 2;
538  else if(pid1 == 1000035) i = 3;
539  else
540  {
541  std::stringstream ss;
542  ss << "Invalid final state neutralino PDG code " << pid1;
543  ColliderBit_error().raise(LOCAL_INFO, ss.str());
544  return -1;
545  }
546  if(pid2 == 1000022) j = 0;
547  else if(pid2 == 1000023) j = 1;
548  else if(pid2 == 1000025) j = 2;
549  else if(pid2 == 1000035) j = 3;
550  else
551  {
552  std::stringstream ss;
553  ss << "Invalid final state neutralino PDG code " << pid2;
554  ColliderBit_error().raise(LOCAL_INFO, ss.str());
555  return -1;
556  }
557 
558  // Set slepton masses
559  double msL = mS[0];
560  double msR = mS[1];
561 
562  // Couplings
563  // e = g \sin\theta_W = g' \cos\theta_W
564  // alpha = e^2 / 4\pi
565  int deltaij = 0;
566  if (i == j) deltaij = 1;
567  double cos2b = (1.-pow2(tanb))/(1.+pow2(tanb));
568  double sin2b = 2.*tanb/(1.+pow2(tanb));
569  double T3l = -0.5;
570  double Le = T3l+sin2thetaW;
571  double Re = sin2thetaW;
572  double OL[4][4];
573  for(int k = 0; k < 4; k++){
574  for(int l = 0; l < 4; l++){
575  OL[k][l] = 0.5*(N[k][2]*N[l][2]-N[k][3]*N[l][3])*cos2b-0.5*(N[k][2]*N[l][3]+N[k][3]*N[l][2])*sin2b;
576  }
577  }
578  double fL[4], fR[4];
579  for(int k = 0; k < 4; k++){
580  fL[k] = -sqrt(2.) * (Le*N[k][1]/sqrt(1.-sin2thetaW)+sqrt(sin2thetaW)*N[k][0]);
581  fR[k] = sqrt(2.) * sqrt(sin2thetaW) * (sqrt(sin2thetaW/(1.-sin2thetaW))*N[k][1]-N[k][0]);
582  }
583 
584  // Kinematics
585  double s, q, Ei, Ej, DZ2, ReDZ;
586  s = pow2(sqrts);
587  DZ2 = 1./(pow2(s-pow2(mZ))+pow2(mZ*gZ)); // Breit-Wigner for Z
588  ReDZ = (s-pow2(mZ))*DZ2;
589  Ei = (s+pow2(mi)-pow2(mj))/2./sqrts; // Energy of \tilde\chi^0_i in e+e- CoM system
590  q = sqrt(pow2(Ei)-pow2(mi)); // Momentum of \tilde\chi^0_i in e+e- CoM system
591  Ej = sqrt(pow2(q)+pow2(mj));
592 
593  double dL, dR;
594  dL = 0.5/s * (s + 2*pow2(msL) - pow2(mi) - pow2(mj));
595  dR = 0.5/s * (s + 2*pow2(msR) - pow2(mi) - pow2(mj));
596 
597  // Cross sections per diagram and interference terms
598  double sigma, sigma_Z, sigma_s, sigma_Zs;
599  // Z
600  sigma_Z = 4.*pi*pow2(alpha)/pow2(sin2thetaW)/pow2(1.-sin2thetaW) * DZ2 * q/sqrts * pow2(OL[i][j]) * (pow2(Le)+pow2(Re));
601  sigma_Z *= Ei*Ej + 1/3.*pow2(q)-mi*mj;
602  // selectrons
603  sigma_s = pow2(fL[i]*fL[j]) * ((Ei*Ej-s*dL+pow2(q))/(s*pow2(dL)-pow2(q)) + 2. + 0.5*sqrts/q*(1.-2.*dL-mi*mj/s/dL)*log(fabs((dL+q/sqrts)/(dL-q/sqrts))));
604  sigma_s += pow2(fR[i]*fR[j]) * ((Ei*Ej-s*dR+pow2(q))/(s*pow2(dR)-pow2(q)) + 2. + 0.5*sqrts/q*(1.-2.*dR-mi*mj/s/dR)*log(fabs((dR+q/sqrts)/(dR-q/sqrts))));
605  sigma_s *= pi*pow2(alpha)/pow2(sin2thetaW) * q/s/sqrts;
606  // Interference
607  sigma_Zs = Le*fL[i]*fL[j] * (1./q/sqrts*(Ei*Ej-s*dL*(1.-dL)-mi*mj)*log(fabs((dL+q/sqrts)/(dL-q/sqrts)))+2.*(1-dL));
608  sigma_Zs += -Re*fR[i]*fR[j] * (1./q/sqrts*(Ei*Ej-s*dR*(1.-dR)-mi*mj)*log(fabs((dR+q/sqrts)/(dR-q/sqrts)))+2.*(1-dR));
609  sigma_Zs *= -2.*pi*pow2(alpha)/pow2(sin2thetaW)/(1.-sin2thetaW) * q/sqrts * ReDZ * OL[i][j];
610 
611  // Total cross section
612  sigma = 0.5*(sigma_Z + sigma_s + sigma_Zs)*(2.-deltaij);
613 
614  // Fix units
615  sigma *= gev2pb;
616 
617  // Return zero in corner cases where numerical roundoff has sent sigma negative.
618  return std::max(sigma, 0.0);
619  }
620 
621 
624  double xsec_chaichaj(int pid1, int pid2, double sqrts, double mi, double mj, MixMatrix V,
625  MixMatrix U, double ms, double alpha, double mZ, double gZ, double sin2thetaW)
626  {
627  // Just return zero if the final state isn't kinematically accessible
628  if (std::abs(mi)+std::abs(mj) > sqrts) return 0.0;
629 
630  // Translate from PDG codes to chargino indices (silly paper convention that i=2 lighter than i=1!)
631  int i, j;
632  pid1 = abs(pid1); pid2 = abs(pid2);
633  if(pid1 == 1000024) i = 1;
634  else if(pid1 == 1000037) i = 0;
635  else
636  {
637  std::stringstream ss;
638  ss << "Invalid final state chargino PDG code " << pid1;
639  ColliderBit_error().raise(LOCAL_INFO, ss.str());
640  return -1;
641  }
642  if(pid2 == 1000024) j = 1;
643  else if(pid2 == 1000037) j = 0;
644  else
645  {
646  std::stringstream ss;
647  ss << "Invalid final state chargino PDG code " << pid2;
648  ColliderBit_error().raise(LOCAL_INFO, ss.str());
649  return -1;
650  }
651 
652  // Couplings
653  int deltaij = 0;
654  if (i == j) deltaij = 1;
655  // e = g \sin\theta_W = g' \cos\theta_W
656  double T3l = -0.5;
657  double Le = T3l+sin2thetaW;
658  double Re = sin2thetaW;
659  double OL[2][2], OR[2][2];
660  for(int k = 0; k < 2; k++){
661  for(int l = 0; l < 2; l++){
662  OL[k][l] = -V[k][0]*V[l][0]-0.5*V[k][1]*V[l][1]+deltaij*sin2thetaW;
663  OR[k][l] = -U[k][0]*U[l][0]-0.5*U[k][1]*U[l][1]+deltaij*sin2thetaW;
664  }
665  }
666 
667  // Kinematics
668  double s, q, Ei, Ej, DZ2, ReDZ;
669  s = pow2(sqrts);
670  Ei = (s+pow2(mi)-pow2(mj))/2./sqrts; // Energy of \tilde\chi^+_i in e+e- CoM system
671  q = sqrt(pow2(Ei)-pow2(mi)); // Momentum of \tilde\chi^+_i in e+e- CoM system
672  Ej = sqrt(pow2(q)+pow2(mj));
673  DZ2 = 1./(pow2(s-pow2(mZ))+pow2(mZ*gZ)); // Breit-Wigner for Z
674  ReDZ = (s-pow2(mZ))*DZ2;
675 
676  double aL, bL, h;
677  aL = 0.5/pow2(ms)*(2*pow2(ms)+s-pow2(mi)-pow2(mj));
678  bL = q*sqrts/pow2(ms);
679  h = 2.*q*sqrts-2.*pow2(q)*aL/bL+(Ei*Ej+pow2(q*aL/bL)-q*sqrts*aL/bL)*log(fabs((aL+bL)/(aL-bL)));
680 
681  // Cross sections per diagram and interference terms
682  double sigma, sigma_g, sigma_Z, sigma_s, sigma_gZ, sigma_gs, sigma_Zs;
683  // Gamma
684  sigma_g = 8*pi*pow2(alpha) * q*sqrts/pow(s,3) * deltaij * (Ei*Ej+pow2(q)/3.+fabs(mi*mj));
685  // Z
686  sigma_Z = 2.*pi*pow2(alpha)/pow2(sin2thetaW)/pow2(1.-sin2thetaW) * q/sqrts * DZ2;
687  sigma_Z *= (pow2(OL[i][j])+pow2(OR[i][j]))*(pow2(Le)+pow2(Re))*(Ei*Ej+pow2(q)/3.)+2.*(pow2(Le)+pow2(Re))*OL[i][j]*OR[i][j]*mi*mj;
688  // Sneutrino
689  sigma_s = pi*pow2(alpha)/2./pow2(sin2thetaW)*pow2(V[i][0]*V[j][0])/pow(ms,4) * q/sqrts;
690  sigma_s *= (Ei*Ej+pow2(q)-q*sqrts*aL/bL)/(pow2(aL)-pow2(bL)) + 2.*pow2(q/bL) + 0.5/pow2(bL)*(q*sqrts-2.*pow2(q)*aL/bL)*log(fabs((aL+bL)/(aL-bL)));
691  // Interference
692  sigma_gZ = 4*pi*pow2(alpha)/(1.-sin2thetaW)/sin2thetaW * q*sqrts/pow2(s)*ReDZ*deltaij*(Le+Re);
693  sigma_gZ *= (OL[i][j]+OR[i][j])*(Ei*Ej+pow2(q)/3.+fabs(mi*mj));
694  sigma_gs = -pi*pow2(alpha)/sin2thetaW*pow2(V[i][0]) * deltaij/pow2(s);
695  sigma_gs *= h + fabs(mi*mj)*log(fabs((aL+bL)/(aL-bL)));
696  sigma_Zs = -pi*pow2(alpha)/pow2(sin2thetaW)/(1.-sin2thetaW)*V[i][0]*V[j][0] * ReDZ/s * Le;
697  sigma_Zs *= OL[i][j]*h + OR[i][j]*mi*mj*log(fabs((aL+bL)/(aL-bL)));
698 
699  // Total cross section with interference terms
700  sigma = sigma_g + sigma_Z + sigma_s+ sigma_gZ + sigma_gs + sigma_Zs;
701 
702  // Units
703  sigma *= gev2pb;
704 
705  // Return zero in corner cases where numerical roundoff has sent sigma negative.
706  return std::max(sigma, 0.0);
707 
708  }
709 
710 
715 
717  void SLHA2BFM_NN(MixMatrix &NN, double tanb, double sin2thetaW)
718  {
719  // Define conversion matrix
720  double sinthetaW = sqrt(sin2thetaW);
721  double costhetaW = sqrt(1.-sin2thetaW);
722  double tanv = 1./tanb; // Needed because of convention difference
723  double sinv = sin(atan(tanv));
724  double cosv = cos(atan(tanv));
725  MixMatrix T(4,std::vector<double>(4));
726  T[0][0] = costhetaW; T[0][1] = -sinthetaW;
727  T[1][0] = sinthetaW; T[1][1] = costhetaW;
728  T[2][2] = sinv; T[2][3] = cosv;
729  T[3][2] = -cosv; T[3][3] = sinv;
730  // Multiply N_{BFM} = N_{SLHA} T
731  NN = multiply(NN,T);
732  }
733 
736  {
737  // Define conversion matrix (\sigma_3)
738  MixMatrix T(2,std::vector<double>(2));
739  T[0][0] = 1; T[0][1] = 0;
740  T[1][0] = 0; T[1][1] = -1;
741  // Multiply V_{BFM} = \sigma_3 V_{SLHA}
742  VV = multiply(T,VV);
743  }
744 
746  void BFM2SLHA_NN(MixMatrix &NN, double tanb, double sin2thetaW)
747  {
748  // Define conversion matrix
749  double sinthetaW = sqrt(sin2thetaW);
750  double costhetaW = sqrt(1.-sin2thetaW);
751  double tanv = 1./tanb; // Needed because of convention difference
752  double sinv = sin(atan(tanv));
753  double cosv = cos(atan(tanv));
754  MixMatrix T(4,std::vector<double>(4));
755  T[0][0] = costhetaW; T[0][1] = -sinthetaW;
756  T[1][0] = sinthetaW; T[1][1] = costhetaW;
757  T[2][2] = sinv; T[2][3] = cosv;
758  T[3][2] = -cosv; T[3][3] = sinv;
759  // Multiply N_{SLHA} = N_{BFM} T^T
760  NN = multiply(NN,transpose(T));
761  }
762 
765  {
766  SLHA2BFM_VV(VV);
767  }
768 
771  {
772  int dim = A.size();
773  MixMatrix C(dim,std::vector<double>(dim));
774  for(int i = 0; i < dim; i++){
775  for(int j = 0; j < dim; j++){
776  for(int k = 0; k < dim; k++){
777  C[i][j] += A[i][k] * B[k][j];
778  }
779  }
780  }
781  return C;
782  }
783 
786  {
787  double temp;
788  int dim = A.size();
789  for(int i = 0; i < dim; i++){
790  for(int j = i+1; j < dim; j++){
791  temp = A[i][j];
792  A[i][j] = A[j][i];
793  A[j][i] = temp;
794  }
795  }
796  return A;
797  }
798 
800  void print(MixMatrix A)
801  {
802  int dim = A.size();
803  cout << "Matrix dimension: " << dim << endl;
804  for(int i = 0; i < dim; i++){
805  for(int j = 0; j < dim; j++){
806  cout << A[i][j] << " ";
807  if(j == dim-1) cout << endl;
808  }
809  }
810  }
811 
813 
814  }
815 }
void print(MixMatrix A)
Helper function to print a matrix.
Rollcall header for ColliderBit module.
std::vector< double > family_state_mix_matrix(str type, int generation, str &mass_es1, str &mass_es2, const SubSpectrum &mssm, double tol, str context, bool pterror)
Get the family mixing matrix and corresponding mass eigenstates, then check for interfamily mixing...
void SLHA2BFM_NN(MixMatrix &NN, double tanb, double sin2thetaW)
Conversion between SLHA and BFM conventions.
double xsec_chaichaj(int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix V, MixMatrix U, double msn, double alpha, double mZ, double gZ, double sin2thetaW)
Cross section [pb] for Masses mi and mj for the charginos are signed.
MixMatrix transpose(MixMatrix A)
Helper function to find matrix transpose.
void BFM2SLHA_NN(MixMatrix &NN, double tanb, double sin2thetaW)
Converts a neutralino mixing matrix in BFM conventions to SLHA conventions, is as defined in SLHA...
void BFM2SLHA_VV(MixMatrix &VV)
Converts the chargino mixing matrix V in BFM conventions to SLHA conventions.
#define LOCAL_INFO
Definition: local_info.hpp:34
MixMatrix multiply(MixMatrix A, MixMatrix B)
Helper function to multiply matrices.
Sparticle production cross-section calculators for LEP.
void get_sigma_ee_chipm(triplet< double > &result, const double sqrts, const int chi_plus, const int chi_minus, const double tol, const bool pt_error, const Spectrum &spec, const double gammaZ)
Retrieve the production cross-section at an e+e- collider for chargino pairs.
Routines to help translate between SLHA2 sfermions and SLHA1 (or similar) sfermions.
std::chrono::milliseconds ms
void get_sigma_ee_ll(triplet< double > &result, const double sqrts, const int generation, const int l_chirality, const int lbar_chirality, const double gtol, const double ftol, const bool gpt_error, const bool fpt_error, const Spectrum &spec, const double gammaZ, const bool l_are_gauge_es)
High-level cross section routines.
START_MODEL alpha
double I2(double s, double m1, double m2, double mk, double ml)
double xsec_sleislej(int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix F, MixMatrix N, const double mN[4], double alpha, double mZ, double gZ, double sin2thetaW, bool warn_on_CP_violating_masses=true)
Low-level cross section routines.
const double sigma
Definition: SM_Z.hpp:43
const double pi
double safeget(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const
safeget functions, by Abram
void SLHA2BFM_VV(MixMatrix &VV)
Converts the chargino mixing matrix V in SLHA conventions to BFM conventions.
double safeget(const Par::Tags partype, const std::string &mass) const
Getters which first check the sanity of the thing they are returning.
Definition: spectrum.cpp:357
Header file that includes all GAMBIT headers required for a module source file.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
double xsec_neuineuj(int pid1, int pid2, double sqrts, double m1, double m2, MixMatrix N, const double mS[2], double tanb, double alpha, double mZ, double gZ, double sin2thetaW)
Cross section [pb] for Masses mi and mj for the neutralinos are signed.
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
double pow(const double &a)
Outputs a^i.
const double gev2pb
double I1(double s, double m1, double m2, double mk, double ml)
Integrals for t-channel neutralino diagrams m1 and m2 are masses of final state sleptons mk and ml ar...
void get_sigma_ee_chi00(triplet< double > &result, const double sqrts, const int chi_first, const int chi_second, const double tol, const bool pt_error, const Spectrum &spec, const double gammaZ)
Retrieve the production cross-section at an e+e- collider for neutralino pairs.
str mass_es_from_gauge_es(str gauge_es, double &max_mixing, std::vector< double > &gauge_composition, const SubSpectrum &mssm)
indentifies the state with largest gauge_es content also fills largest max_mixing and full gauge_comp...
std::vector< std::vector< double > > MixMatrix
TODO: see if we can use this one:
Definition: Analysis.hpp:33
SubSpectrum & get_HE()
Definition: spectrum.cpp:225
#define pow2(a)
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110
double I3(double s, double m1, double m2, double mk)