Added softsusy-specific options to MssmSoftsusy object. They now get …

…printed out in SLHA output

README.md

@@ -11,7 +11,7 @@ can be used in conjunction with other programs for [many different particle
 physics calculations](https://arxiv.org/abs/0805.2088). SOFTSUSY now has a mode with 3 loop RGEs and some 2-loop threshold correction and 2-loop SUSY QCD corrections to gluino and squark pole masses. SOFTSUSY *now computes decay branching ratios for the MSSM and NMSSM*. It also ships with (and links to) Himalaya-1.0 for three-loop corrections to mh0. 
 
 ## Download
-* <a id="raw-url" href="https://raw.githubusercontent.com/BAllanach/softsusy/master/tags/softsusy-4.1.20.tar.gz">SOFTSUSY-4.1.20</a> (??/8/23): Fixed SLHA output which was giving the wrong RPV inputs. Also switched MSUSY scale RPV BC on in rpvHouchesInput: this gives lighter neutrinos.
+* <a id="raw-url" href="https://raw.githubusercontent.com/BAllanach/softsusy/master/tags/softsusy-4.1.20.tar.gz">SOFTSUSY-4.1.20</a> (??/8/23): Fixed SLHA output which was giving the wrong RPV inputs. SOFTSUSY-specific SLHA input options given automatically printed out now in SLHA output.
 * <a id="raw-url" href="https://raw.githubusercontent.com/BAllanach/softsusy/master/tags/softsusy-4.1.19.tar.gz">SOFTSUSY-4.1.19</a> (16/8/23): improved 2/3 body neutralino/chargino decays. Previously, there were occasionally numerical problems yielding junk for a (very heavy) RH sneutrino contribution. The RH sneutrino parts have now been explicitly decoupled. Thanks to S Kraml and T Pascal for reporting the bug.
 * <a id="raw-url" href="https://raw.githubusercontent.com/BAllanach/softsusy/master/tags/softsusy-4.1.18.tar.gz">SOFTSUSY-4.1.18</a> (25/7/23): improved 2/3 body neutralino decays. Occasionally, SOFTSUSY would find neutralino2 to be stable when in fact it really had 3 body decays. This occured when the decay of the running mass of the Z would allow 2 body decays but not the pole mass. Thanks to S Kraml and T Pascal for reporting the bug.
 * <a id="raw-url" href="https://raw.githubusercontent.com/BAllanach/softsusy/master/tags/softsusy-4.1.17.tar.gz">SOFTSUSY-4.1.17</a> (2/2/23): improved documentation for github

inOutFiles/lesHouchesInput

@@ -31,7 +31,7 @@ Block SOFTSUSY               # Optional SOFTSUSY-specific parameters
    13	0.000000000e+00	     # Set a tachyonic A^0 to zero mass
    19   1.000000000e+00      # Include 3-loop SUSY RGEs
    20   0.000000000e+00      # Include 2-loop g/Yuk corrections: 31 for all
-   21   1.000000000e+00	      # Calculate uncertainties in Higgs mass prediction
+   21   1.000000000e+00	     # Calculate uncertainties in Higgs mass prediction
 #   22   1.000000000e+00      # Include 2-loop sparticle mass thresholds
 #   23   0.000000000e+00      # No expansion of 2-loop gluino terms
 #   24   1.000000000e-06      # If decay BR is below this number, don't output

inOutFiles/lesHouchesOutput

@@ -4,6 +4,21 @@
 Block SPINFO          # Program information
      1    SOFTSUSY    # spectrum calculator
      2    4.1.20       # version number
+Block SOFTSUSY               # SOFTSUSY specific SLHA options
+    0   1.000000000e+00      # Calculate decays in output (only for RPC (N)MSSM)
+    1   1.000000000e-03      # Numerical precision: suggested range 10^(-3...-6)
+    2   0.000000000e+00	     # Quark mixing parameter: see manual
+    3   0.000000000e+00	     # Additional verbose output?
+    4   1.000000000e+00	     # Change electroweak symmetry breaking scale?
+    6   1.000000000e-04	     # Numerical precision
+    7	2.000000000e+00	     # Number of loops in Higgs mass computation
+   10	0.000000000e+00	     # Force it to SLHA***1*** output?
+   11	1.000000000e+19	     # Gravitino mass
+   12   0.000000000e+00	     # Print spectrum even when point disallowed
+   13	0.000000000e+00	     # Set a tachyonic A^0 to zero mass
+   19   1.000000000e+00      # Include 3-loop SUSY RGEs
+   20   0.000000000e+00      # Include 2-loop g/Yuk corrections: 31 for all
+   21   1.000000000e+00	     # Calculate uncertainties in Higgs mass prediction
 Block MODSEL  # Select model
      1    1   # sugra
 Block SMINPUTS             # Standard Model inputs

inOutFiles/neutOutputTest

@@ -7,6 +7,7 @@
 Block SPINFO          # Program information
      1    SOFTSUSY    # spectrum calculator
      2    4.1.20       # version number
+Block SOFTSUSY               # SOFTSUSY specific SLHA options
 Block MODSEL  # Select model
      1    1   # sugra
      4    1   # R-parity violating

inOutFiles/nmssmSLHAZ3Output

@@ -5,6 +5,18 @@
 Block SPINFO          # Program information
      1    SOFTSUSY    # spectrum calculator
      2    4.1.20       # version number
+Block SOFTSUSY               # SOFTSUSY specific SLHA options
+    0   1.000000000e+00      # Calculate decays in output (only for RPC (N)MSSM)
+    1   1.000000000e-04      # tolerance
+    2   2.000000000e+00      # up-quark mixing (=1) or down (=2)
+    5   1.000000000E+00      # 2-loop running
+    3   0.000000000E+00      # printout
+   15   1.000000000E+00      # NMSSMTools compatible output (default: 0)
+   16   0.000000000E+00      # Call micrOmegas (default: 0 = no,
+   17   1.000000000E+00      # Sparticle decays via NMSDECAY (default: 0)
+   18   0.000000000E+00      # use soft Higgs masses as EWSB output
+   24   1.000000000e-06      # If decay BR is below this number, don't output
+   25   1.000000000e+00	     # If set to 0, don't calculate 3-body decays (1=default)
 Block MODSEL  # Select model
      1    1   # sugra
      3    1   # NMSSM

inOutFiles/nmssmSLHAnoZ3Output

@@ -5,6 +5,18 @@
 Block SPINFO          # Program information
      1    SOFTSUSY    # spectrum calculator
      2    4.1.20       # version number
+Block SOFTSUSY               # SOFTSUSY specific SLHA options
+    0   1.000000000e+00      # Calculate decays in output (only for RPC (N)MSSM)
+    1   1.000000000e-04      # tolerance
+    2   2.000000000e+00      # up-quark mixing (=1) or down (=2)
+    5   1.000000000E+00      # 2-loop running
+    3   0.000000000E+00      # printout
+   15   1.000000000E+00      # NMSSMTools compatible output (default: 0)
+   16   0.000000000E+00      # Call micrOmegas (default: 0 = no,
+   17   0.000000000E+00      # Sparticle decays via NMSDECAY (default: 0)
+   18   0.000000000E+00      # use soft Higgs masses as EWSB output
+   24   1.000000000e-06      # If decay BR is below this number, don't output
+   25   1.000000000e+00	     # If set to 0, don't calculate 3-body decays (1=default)
 Block MODSEL  # Select model
      1    1   # sugra
      3    1   # NMSSM

inOutFiles/rpvHouchesOutput

@@ -7,6 +7,11 @@
 Block SPINFO          # Program information
      1    SOFTSUSY    # spectrum calculator
      2    4.1.20       # version number
+Block SOFTSUSY               # SOFTSUSY specific SLHA options
+    1   1.000000000e-03
+    2   2.000000000e+00      # quark mixing: none (=0), up (=1) or down (=2)
+    8   0.000000000e+00	     # Switch MSUSY scale RPV boundary condition on?
+    9   0.000000000e+00	     # Use normal hierarchy(=0.0) or inverted (=1.0)
 Block MODSEL  # Select model
      1    1   # sugra
      4    1   # R-parity violating

inOutFiles/slha2Output

@@ -3,6 +3,11 @@
 Block SPINFO          # Program information
      1    SOFTSUSY    # spectrum calculator
      2    4.1.20       # version number
+Block SOFTSUSY               # SOFTSUSY specific SLHA options
+    1   1.000000000e-03      # tolerance
+    2   2.000000000e+00      # up-quark mixing (=1) or down (=2)
+    5   1.000000000E+00      # 2-loop running
+    3   0.000000000E+00	     # printout
 Block MODSEL  # Select model
      1    1   # sugra
      6    1   # flavour violating MSSM

src/nmssmsoftsusy.cpp

@@ -9512,6 +9512,7 @@ namespace {
       bool GUTsVev = displayGUTsVev();
       int MICROMEGAS = displayMICROMEGAS();
       int NMSDECAY   = displayNMSDECAY();
+      string sopts = displaySoftsusyOpts();
 
       setSoftsusy(empty); /// Always starts from an empty object
       /// These are things that are re-written by the new initialisation
@@ -9527,6 +9528,7 @@ namespace {
       setGUTsVev(GUTsVev);
       setMICROMEGAS(MICROMEGAS);
       setNMSDECAY(NMSDECAY);
+      setSoftsusyOpts(sopts);
       double mz = displayMz();
 
       if (mxGuess > 0.0)

src/rpvneut.cpp

@@ -2793,6 +2793,8 @@ void RpvNeutrino::spinfoSLHA(ostream & out) {
 	<< displayFracDiff() << endl;
   if (displayProblem().testSeriousProblem()) 
     out << "     4   Point invalid: " << displayProblem() << endl;
+  out << "Block SOFTSUSY               # SOFTSUSY specific SLHA options" << endl;
+  out << displaySoftsusyOpts();
 }
 
 void RpvNeutrino::hmixSLHA(ostream & out) {

src/softpoint.cpp

@@ -1132,6 +1132,7 @@ int main(int argc, char *argv[]) {
 		  }
 		}
 		else if (block == "SOFTSUSY") {
+		  r->setSoftsusyOpts(r->displaySoftsusyOpts() + line + "\n");
 		  int i; double d; kk >> i >> d;
 		  switch(i) {
 		  case 0: if (int(d) > 0) calcDecays = true;
@@ -1418,10 +1419,13 @@ int main(int argc, char *argv[]) {
          }
       }
 
-      if (RPVflag) {
+      if (RPVflag) {	
 	kw.setNumRpcBcs(pars.displayEnd());
 	kw.rpvDisplay(pars);
+	string sopts = r->displaySoftsusyOpts();
+	/// Bug in next line - temporarily fixed by following line
 	kw.setFlavourSoftsusy(k);
+	kw.setSoftsusyOpts(sopts);
 	r = &kw;
 
 	if (boundaryCondition == &sugraBcs) 
@@ -1464,6 +1468,8 @@ int main(int argc, char *argv[]) {
       if (useThreeLoopRge) {
 	r->setLoops(3); k.setLoops(3); nmssm.setLoops(3);
       }
+      k.setSoftsusyOpts(r->displaySoftsusyOpts());
+      nmssm.setSoftsusyOpts(r->displaySoftsusyOpts());
 
       /// Run to scale at which MSUSY and QEDxQCD are matched: by default it's
       /// mt for MSSM, MZ for NMSSM

src/softsusy.cpp

@@ -58,6 +58,7 @@ const MssmSoftsusy& MssmSoftsusy::operator=(const MssmSoftsusy& s) {
   t1OV1Ms1loop = s.displayTadpole1Ms1loop(); 
   t2OV2Ms1loop = s.displayTadpole2Ms1loop(); 
   qewsb = s.displayQewsb();
+  softsusyOpts = s.softsusyOpts;
   mxBC = s.displayMxBC();
   included_thresholds = s.included_thresholds;
 
@@ -71,7 +72,8 @@ const MssmSoftsusy& MssmSoftsusy::operator=(const MssmSoftsusy& s) {
       mw(0.0), dataSet(), fracDiff(1.), setTbAtMX(false), altEwsb(false),
       altMt(false), altAlphaS(false), altAlphaEm(false),
       predMzSq(0.), t1OV1Ms(0.), t2OV2Ms(0.),
-      t1OV1Ms1loop(0.), t2OV2Ms1loop(0.), qewsb(1.), mxBC(mxDefault)
+      t1OV1Ms1loop(0.), t2OV2Ms1loop(0.), qewsb(1.), softsusyOpts(),
+      mxBC(mxDefault)
   {
     setPars(110);
     setMu(0.0);
@@ -101,6 +103,7 @@ const MssmSoftsusy& MssmSoftsusy::operator=(const MssmSoftsusy& s) {
       t1OV1Ms(s.displayTadpole1Ms()), t2OV2Ms(s.displayTadpole2Ms()), 
       t1OV1Ms1loop(s.displayTadpole1Ms1loop()), 
       t2OV2Ms1loop(s.displayTadpole2Ms1loop()), qewsb(s.displayQewsb()),
+      softsusyOpts(s.displaySoftsusyOpts()), 
       mxBC(s.displayMxBC()) {
 
     setPars(110);
@@ -116,7 +119,7 @@ MssmSoftsusy::MssmSoftsusy(const MssmSusyRGE &s)
 	setTbAtMX(false), altEwsb(false), altMt(false), altAlphaS(false),
         altAlphaEm(false), predMzSq(0.),
         t1OV1Ms(0.), t2OV2Ms(0.), t1OV1Ms1loop(0.), t2OV2Ms1loop(0.),
-	qewsb(1.0), mxBC(mxDefault) { 
+	  qewsb(1.0), softsusyOpts(), mxBC(mxDefault) { 
     setPars(110);
     setMu(s.displayMu()); 
 
@@ -141,7 +144,7 @@ MssmSoftsusy::MssmSoftsusy(const MssmSusyRGE &s)
       altEwsb(false), altMt(false), altAlphaS(false), altAlphaEm(false),
       predMzSq(0.), t1OV1Ms(0.),
       t2OV2Ms(0.), t1OV1Ms1loop(0.), t2OV2Ms1loop(0.),
-	qewsb(1.0), mxBC(mxDefault){
+      qewsb(1.0), softsusyOpts(), mxBC(mxDefault){
     setHvev(hv);
     setPars(110);
     setMu(mu);
@@ -7383,6 +7386,7 @@ void MssmSoftsusy::fixedPointIteration
     double maCondFirst = displayMaCond();
     double qqewsb      = displayQewsb();
     int lpnum = displayLoops();
+    string sopts = displaySoftsusyOpts();
 
     int enabled_thresholds = included_thresholds;    
     setSoftsusy(empty); /// Always starts from an empty object
@@ -7398,6 +7402,7 @@ void MssmSoftsusy::fixedPointIteration
     setMuCond(muCondFirst);
     setMaCond(maCondFirst);
     setQewsb(qqewsb);
+    setSoftsusyOpts(sopts);
 
     double mz = displayMz();
 
@@ -10012,6 +10017,8 @@ void MssmSoftsusy::spinfoSLHA(ostream & out) {
   }
   if (displayProblem().testSeriousProblem()) 
     out << "     4    " << displayProblem() << " # Point invalid" << endl;
+  out << "Block SOFTSUSY               # SOFTSUSY specific SLHA options" << endl;
+  out << displaySoftsusyOpts();
 }
 
 void MssmSoftsusy::softsusySLHA(ostream & out) {

src/softsusy.h

@@ -103,7 +103,9 @@ namespace softsusy {
     /// there are two possible conventions: if QEWSB > MZ, its value is assumed
     /// in GeV and used as a constant MSUSY. Otherwise, it MULTIPLIES the usual 
     /// MSUSY value, of root(mstop1 mstop2)
-    double qewsb;       
+    double qewsb;
+    /// SOFTSUSY 
+    string softsusyOpts;
 
   protected:
     void setT1OV1Ms(double t1) { t1OV1Ms = t1; } 
@@ -169,8 +171,9 @@ namespace softsusy {
     /// Return a trilinear element in "SUGRA style"
     double displaySoftA(trilinears, int, int) const;
 
-    /// Displays iteration accuracy attained 
-     double displayFracDiff() const { return fracDiff; }; 
+    /// Displays iteration accuracy attained
+    string displaySoftsusyOpts() const { return softsusyOpts; }
+    double displayFracDiff() const { return fracDiff; }; 
     double displayMinpot() const;    ///< Returns minimum of Higgs potential
     double displayMsusy() const; ///< Returns Higgs minimisation scale
     double displayMw() const; ///< Returns predicted pole MW
@@ -242,7 +245,9 @@ namespace softsusy {
     void flagProblemThrown(bool a) { problem.problemThrown = a; }
     /// Sets all problem flags 
     void setProblem(const sProblem a) { problem = a; }
-
+
+    /// Sets SOFTSUSY-specific option string
+    void setSoftsusyOpts(const string & s) { softsusyOpts = s; }
     /// Sets whole object equal to another  
     void setSoftsusy(const MssmSoftsusy & s) { *this = s; };
     /// Sets low energy Standard Model fermion mass and gauge coupling data