Synced NESTv2.3.9 and updated GetQuanta Bindings

HISTORY.md

@@ -5,6 +5,22 @@ History
 Patch releases mean (the Z number in X.Y.Z version) that the underlying physics has not changed.  Changes to the NEST version will always trigger a minor or major release.  If this library changes such that end users have to change their code, this may also trigger a minor or major release.
 
 
+1.5.5 (2022-07-08)
+-----------------
+Synced with NEST v2.3.9
+  * New Physics Modeling:
+  ** Skewness can be turned off and on now for ER just like for NR. For on -> old model or fixed (Quentin Riffard, LZ/LBNL)
+  ** Older beta model is default for gaseous Xenon, a better fit to old world data at the keV scale (Eric Church, DUNE/PNNL)
+  ** New dark matter halo model defaults, bringing NEST up to date on WIMP and Sun v (Baxter et al., arXiv:2105.00599)
+
+  * Miscellaneous Bug Fixes:
+  ** Fluctuations adjust for difference in width from truncated Gaussians for PE not just mean
+  ** Complaint that position resolution too poor does not activate until above S2 (top) of 40 PE
+  ** In the dE/dx-based model the minimum LET is now 1.0 MeV/cm not 0 to avoid weirdness
+
+  * Updated binding for GetQuanta to allow for nestpy control over ER skewness. 
+
+
 1.5.4 (2022-04-16)
 -----------------
 nestpy specific code quality:

src/nestpy/DetectorExample_XENON10.hh

@@ -154,9 +154,9 @@ class DetectorExample_XENON10 : public VDetector {
     }
 
     double sig = RandomGen::rndm()->rand_gauss(
-        3.84, .09);  // includes stat unc but not syst
+					       3.84, .09, true);  // includes stat unc but not syst
     phoTravT += RandomGen::rndm()->rand_gauss(
-        0.00, sig);  // the overall width added to photon time spectra by the
+					      0.00, sig, false);  // the overall width added to photon time spectra by the
                      // effects in the electronics and the data reduction
                      // pipeline
 

src/nestpy/LUX_Run03.hh

@@ -196,10 +196,10 @@ class DetectorExample_LUX_RUN03 : public VDetector {
     }
 
     double sig = RandomGen::rndm()->rand_gauss(
-        3.84, .09);  // includes stat unc but not syst
+					       3.84, .09, true);  // includes stat unc but not syst
     phoTravT += RandomGen::rndm()->rand_gauss(
         0.00,
-        sig);  // the overall width added to photon time spectra by the effects
+        sig, false);  // the overall width added to photon time spectra by the effects
                // in the electronics and the data reduction pipeline
 
     if (phoTravT > DBL_MAX) phoTravT = tau_a;

src/nestpy/NEST.hh

@@ -57,7 +57,7 @@
    calculations with other doubles as a mixture
    The interaction type ("species") is a new type defined for use in NEST that
    attempts to cover all possible kinds of particle/scatter
-   For waveforms (with output timing set to TRUE) the unit is PE per bin
+   For waveforms (with output timing set to >=1) the unit is PE per bin
    (sample, usually 10 ns) but has option of total area in PE
 
 */
@@ -303,10 +303,10 @@ class NESTcalc {
                                           const double Wq_eV);
   // Confirms and sometimes adjusts YieldResult to make physical sense
 
-  virtual QuantaResult GetQuanta(
-      const YieldResult &yields, double density,
-      const std::vector<double> &NRERWidthsParam = {1.,1.,0.1,0.5,0.19,2.25, 0.0015, 0.0553, 0.205, 0.45, -0.2},
-      bool oldModelER = false);
+  virtual QuantaResult GetQuanta(const YieldResult &yields, double density,
+                                 const std::vector<double> &NRERWidthsParam = {1.,1.,0.1,0.5,0.19,2.25, 0.0015, 0.0553, 0.205, 0.45, -0.2},
+                                 bool oldModelER = false,
+                                 double SkewnessER = -999.);
   // GetQuanta takes the yields from above and fluctuates them, both the total
   // quanta (photons+electrons) with a Fano-like factor, and the "slosh" between
   // photons and electrons
@@ -334,7 +334,7 @@ class NESTcalc {
       double truthPosZ, double smearPosX, double smearPosY, double smearPosZ,
       double driftSpeed, double dS_mid, INTERACTION_TYPE species,
       uint64_t evtNum, double dfield, double energy, S1CalculationMode mode,
-      bool outputTiming, vector<int64_t> &wf_time, vector<double> &wf_amp);
+      int outputTiming, vector<int64_t> &wf_time, vector<double> &wf_amp);
   // Very comprehensive conversion of the "original" intrinsic scintillation
   // photons into the many possible definitions of S1 as measured by
   // photo-sensors
@@ -349,15 +349,15 @@ class NESTcalc {
       int Ne, double truthPosX, double truthPosY, double truthPosZ,
       double smearPosX, double smearPosY, double smearPosZ, double dt,
       double driftSpeed, uint64_t evtNum, double dfield, S2CalculationMode mode,
-      bool outputTiming, vector<int64_t> &wf_time, vector<double> &wf_amp,
+      int outputTiming, vector<int64_t> &wf_time, vector<double> &wf_amp,
       const vector<double> &g2_params);
   // Exhaustive conversion of the intrinsic ionization electrons into the many
   // possible definitions of S2 pulse areas as observed in the photo-tubes
   // This function also applies the extraction efficiency (binomial) and finite
   // electron mean free path or life time caused by electronegative impurities
   // (exponential)
 
-  std::vector<double> CalculateG2(bool verbosity = true);
+  std::vector<double> CalculateG2(int verbosity = 1);
   // Calculates "g2" by combining the single electron size with the extraction
   // efficiency. Called by GetS2 above. Includes helper variables like gas gap
   // and SE width.

src/nestpy/RandomGen.cpp

@@ -41,13 +41,23 @@ double RandomGen::rand_gauss(double mean, double sigma, bool zero_min) {
 }
 
 double RandomGen::rand_zero_trunc_gauss(double mean, double sigma) {
-  double r = rand_gauss(mean, sigma);
+  double r = rand_gauss(mean, sigma, false);
   while (r <= 0) {
-    r = rand_gauss(mean, sigma);
+    r = rand_gauss(mean, sigma, false);
   }
   return r;
 }
 
+double RandomGen::FindNewMean(double sigma) {
+  // Follow https://en.wikipedia.org/wiki/Truncated_normal_distribution
+  double TruncGaussAlpha = -1. / sigma;
+  double LittlePhi_Alpha =
+    exp(-0.5 * TruncGaussAlpha * TruncGaussAlpha) / sqrt(2. * M_PI);
+  double BigPhi_Alpha = 0.5 * (1. + erf(TruncGaussAlpha / sqrt2));
+  return
+    1. + (LittlePhi_Alpha / (1. - BigPhi_Alpha)) * sigma;
+}
+
 double RandomGen::rand_exponential(double half_life, double t_min, double t_max) {
   double r = rand_uniform();
   double r_min = 0;

src/nestpy/RandomGen.hh

@@ -21,6 +21,7 @@ class RandomGen {
   double rand_uniform();
   double rand_gauss(double mean, double sigma, bool zero_min=false);
   double rand_zero_trunc_gauss(double mean, double sigma);
+  double FindNewMean(double sigma); //shift Gaussian of mean 1 for 0 truncation
   double rand_exponential(double half_life, double t_min=-1, double t_max=-1);
   double rand_skewGauss(double xi, double omega, double alpha);
   int poisson_draw(double mean);

src/nestpy/TestSpectra.hh

@@ -24,12 +24,10 @@
   6.0221409e+23       // good to keep in sync w/ NEST.hh, can't define twice
 #define ATOM_NUM 54.  // ibid.
 
-#define RHO_NAUGHT 0.3  // local DM halo density in [GeV/cm^3]
-#define V_SUN 233.
-// works out to mean V_EARTH of 245 for LUX Run03; Run04, 230 km/s
-// (arXiv:1705.03380)
-#define V_WIMP 220.
-#define V_ESCAPE 544.
+#define RHO_NAUGHT 0.3  // local DM halo density in [GeV/cm^3]. Lewin & Smith
+#define V_SUN 250.2  // +/- 1.4: arXiv:2105.00599, page 12 (V_EARTH 29.8km/s)
+#define V_WIMP 238.  // +/- 1.5: page 12 and Table 1
+#define V_ESCAPE 544.  // M.C. Smith et al. RAVE Survey
 
 #define NUMBINS_MAX 1000
 

src/nestpy/__init__.py

@@ -1,5 +1,5 @@
-__version__ = '1.5.4'
-__nest_version__ = '2.3.6'
+__version__ = '1.5.5'
+__nest_version__ = '2.3.9'
 
 from .nestpy import *
 

src/nestpy/analysis.hh

@@ -1,7 +1,7 @@
 #include "NEST.hh"
 
 // Verbosity flag (for limiting output to yields; no timing)
-bool verbosity = true;
+int verbosity = 1;  // use -1 for off, but 0 for efficiency
 // Loop for execNEST and rootNEST to find the best-fit model parameters
 unsigned loopNEST = 0;
 // 0 for no or off, 1 for ER, 2 for NR
@@ -33,7 +33,7 @@ int usePD = 2;
 // band style: log(S2) with 1, while 0 means log(S2/S1)
 int useS2 = 0;  // xtra feature: 2 means S2 x-axis energy scale
 
-double minS1 = 1.5;  // units are controlled by the usePE flag
+double minS1 = 1.5;  // units are controlled by the usePD flag
 // this is separate from S1 thresholds controlled by detector
 double maxS1 = 99.5;
 int numBins = 98;  // for LUXRun03 DD, change these to 1.7,110.6,99

src/nestpy/bindings.cpp

@@ -289,7 +289,8 @@ PYBIND11_MODULE(nestpy, m) {
 	    py::arg("yields"),
 	    py::arg("density") = 2.9,
 	    py::arg("free_parameters") = std::vector<double>({1.,1.,0.1,0.5,0.19,2.25, 0.0015, 0.0553, 0.205, 0.45, -0.2}), 
-        py::arg("oldModelER") = false
+            py::arg("oldModelER") = false,
+        py::arg("SkewnessER") = -999.
      )   
 	.def("GetS1", &NEST::NESTcalc::GetS1)
 	.def("GetSpike", &NEST::NESTcalc::GetSpike)