CHEIDI: Heidi.m

(* Spectrum calculator for compact HEIDI. For more information, please take    *
 * a look at                                                                   *
 *                                                                             *
 *   "Exploring the golden channel for HEIDI models using an interface between *
 *      WHIZARD and FeynRules", N. D. Christensen, C. Duhr, B. Fuks, J. Reuter,*
 *    C. Speckner, arXiv:1010.3251                                             *
 *                                                                             *
 * Please cite above reference if you use this program.                        *)

BeginPackage["Heidi`"];
Begin["Heidi`p`"];
Needs["Hdecay`"];

(* Initialize the package and clear all caches. *)
Heidi`InitHeidi[mv_, mmh_, mcs_, mmb_, mg2_, mprec_] := Block[{},
   v = mv;
   mh = mmh;
   cs = mcs;
   mb = mmb;
   g2 = mg2;
   g = N[Sqrt[g2]];
   alpha = N[g2 / mb * Coth[mb * Pi / cs]];
   lambda = N[mh^2 / v^2 / 2];
   prec = mprec;
   Clear[MassCache, WFCache, WidthCache];
   If[lambda <= alpha,
      mh = 1.1 * Sqrt[2*alpha] * v;
      Print["WARNING: choosen Higgs mass is to small; using mh "
         <> ToString [mh] <> "GeV instead!"];
      InitHeidi[v, mh, cs, mb, g2, prec];
   ];
];

Heidi`InitHeidi[v_, mh_, cs_, mb_, g2_] := InitHeidi[v, mh, cs, mb, g2, 0.000001];

Heidi`InitHeidi::usage = "usage: InitHeidi [v, mh, cs, mb, g2, prec]
   v    : higgs VeV
   mh   : higgs mass
   cs   : compactification scale
   mb   : bulk mass
   g2   : 5D mixing squared
   prec : precision for the numeric determination of the masses
          (optional; default: 0.000001)";


(* Calculate the nth mass eigenvalue *)
Heidi`HeidiMass[n_] := Block[{fun, l, r, x, oldx},
   If[NumericQ[MassCache[n]],
      MassCache[n]
   ,
      If[n > 0,
         fun = (mh^2 - mb^2 + 2*g2*v^2/# * Cot[Pi/cs*#] - #^2)&;
         l = (n - 1) * cs; r = n * cs;
      ,
         fun = (mb^2 - mh^2 + 2*g2*v^2/# * Coth[Pi/cs*#] - #^2)&;
         l = 0; r = mb;
      ];
      oldx = l;
      x = (l + r) / 2;
      While[Abs[(x - oldx) / x] > prec,
         oldx = x;
         If[fun[x] > 0, l = x, r = x];
         x = (l + r) / 2;
      ];
      MassCache[n] = If[n == 0, Sqrt[mb^2 - N[x]^2], Sqrt[N[x]^2 + mb^2]]
   ]
]/;(n >= 0);

Heidi`HeidiMass::usage = "usage: HeidiMass [n]
   n    : mode index";


(* Calculate the mth component of the nth wavefunction *)
Heidi`HeidiWavefunction[n_, m_] := Block[{},
   If[NumericQ[WFCache[n, m]],
      WFCache[n, m]
   ,
      If[m == 0,
         WFCache[n, m] = N[1/Sqrt[
            1 + g2*v^2*Pi/cs/(HeidiMass[n]^2 - mb^2) +
               (HeidiMass[n]^2 - mh^2)/2/(HeidiMass[n]^2 - mb^2) +
               (HeidiMass[n]^2 - mh^2)^2 * Pi/4/g2/v^2/cs
            ]]
      ,
         WFCache[n, m] = N[HeidiWavefunction[n, 0] *
            2*g*v*Sqrt[cs/Pi/If[m == 1, 2, 1]] /
            (HeidiMass[n]^2 - mb^2 - (m - 1)^2*cs^2)]
      ]
   ]
]/;((n >= 0) && (m >= 0));

Heidi`HeidiWavefunction::usage = "usage: HeidiWavefunction [n, m]
   n : mode index
   m : wavefunction component index (0 = pre-mixing higgs)";


Heidi`HeidiScalarCoupling[i__] := Block[{},
   -I*6 * (Times @@ (Heidi`HeidiWavefunction[#, 0]& /@ {i}))
]/;((Length[{i}] == 4) && (And @@ ((# >= 0)& /@ {i})));

Heidi`HeidiScalarCoupling[i__] := Block [{thing},
   -2*I*v * (Times @@ (Heidi`HeidiWavefunction[#, 0]& /@ {i})) *
      (3*lambda - ((Plus @@ ((HeidiMass[#]^2)& /@ {i})) - 3 * mh^2)/2/v^2)
]/;((Length[{i}] == 3) && (And @@ ((# >= 0)& /@ {i})));

Heidi`HeidiScalarCoupling::usage = "usage: HeidiScalarCoupling [i1, ..., in]
   i1, ..., in : modes meeting at the vertex (n = 3 or n = 4)";


Heidi`HeidiScalarWidth[i_] := Block[{calc},
   calc[{x_, j_, k_}] := Which[
      j == i, {x, j, k},
      (HeidiMass[i] < HeidiMass[j] + HeidiMass[k]) || (k > j), {x, j+1, 0},
      True, {
         Sqrt[
            (HeidiMass[i]^2 - (HeidiMass[j] + HeidiMass[k])^2) *
            (HeidiMass[i]^2 - (HeidiMass[j] - HeidiMass[k])^2)]
         /16/Pi/HeidiMass[i]^3 * Abs[HeidiScalarCoupling[i, j, k]]^2
         /If[j == k, 2, 1] + x, j, k+1}
   ];
   If[NumericQ[WidthCache[i]], WidthCache[i], WidthCache[i] = FixedPoint[calc, {0, 0, 0}][[1]]]
]/;(i >= 0);

Heidi`HeidiScalarWidth::usage = "usage: HeidiScalarWidth [n]
   n : mode index";


Heidi`HeidiSMWidth[i_] := Block[{},
   Hdecay[HeidiMass[i]] * HeidiWavefunction[i, 0]^2
]/;(i >= 0);

Heidi`HeidiSMWidth::usage = "usage: HeidiSMWidth [n]
   n : mode index";


Heidi`HeidiWidth[i_] := HeidiSMWidth[i] + HeidiScalarWidth[i];

Heidi`HeidiWidth::usage = "usage: HeidiWidth [n]
   n : mode index";

End[];
Protect[InitHeidi, HeidiMass, HeidiWavefunction, HeidiScalarCoupling, HeidiScalarWidth,
   HeidiSMWidth, HeidiWidth];
EndPackage[];