Better conditions for approximating binom w/ Gauss, and explicit extr…

…aPhot option to match EXO-200. -MMS

include/Detectors/DetectorExample_XENON10.hh

@@ -47,6 +47,7 @@ class DetectorExample_XENON10 : public VDetector {
     numPMTs = 89;    // For coincidence calculation
 
     //"Linear noise" terms as defined in Dahl thesis and by D. McK
+    extraPhot=false;  // for matching EXO-200's W measurement
     noise[2] = 3e-2;  // S1 -> S1 Gaussian-smeared with noise[2]*S1
     noise[3] = 3e-2;  // S2 -> S2 Gaussian-smeared with noise[3]*S2
 

include/Detectors/VDetector.hh

@@ -27,7 +27,8 @@ class VDetector {
   double get_sPEeff() { return sPEeff; }
   double* get_noise() { return &noise[0]; }
   double get_P_dphe() { return P_dphe; }
-
+
+  bool get_extraPhot(){ return extraPhot; }
   double get_coinWind() { return coinWind; }
   int get_coinLevel() { return coinLevel; }
   int get_numPMTs() { return numPMTs; }
@@ -73,6 +74,7 @@ class VDetector {
   }
   void set_P_dphe(double param) { P_dphe = param; }
 
+  void set_extraPhot(bool param){ extraPhot = param;}
   void set_coinWind(double param) { coinWind = param; }
   void set_coinLevel(int param) { coinLevel = param; }
   void set_numPMTs(int param) { numPMTs = param; }
@@ -144,7 +146,7 @@ class VDetector {
   double g1_gas, s2Fano, s2_thr, E_gas, eLife_us;
 
   // Thermodynamic Properties
-  bool inGas;
+  bool inGas, extraPhot;
   double T_Kelvin, p_bar;
 
   // Data Analysis Parameters and Geometry

include/NEST/analysis.hh

@@ -3,7 +3,7 @@
 bool verbosity = true;
 
 // General parameters of importance changing global behavior
-bool MCtruthE = true;    // false means reconstructed energy
+bool MCtruthE = false;    // false means reconstructed energy
 bool MCtruthPos = true;  // false means reconstructed position
 
 int useTiming = 0;  // photon arrival times + pulse shapes (2=eTrains)

src/NEST.cpp

@@ -15,7 +15,8 @@ long NESTcalc::BinomFluct(long N0, double prob) {
   if (prob <= 0.00) return N1;
   if (prob >= 1.00) return N0;
 
-  if (N0 < 10) {
+  if ( N0 < 5 || fabs(1.-2.*prob)/sqrt(N0*prob*(1.-prob)) > (1./3.) ) {
+    //https://en.wikipedia.org/wiki/Binomial_distribution#Normal_approximation
     for (int i = 0; i < N0; i++) {
       if (RandomGen::rndm()->rand_uniform() < prob) N1++;
     }
@@ -246,11 +247,11 @@ QuantaResult NESTcalc::GetQuanta(YieldResult yields, double density,
     exit(EXIT_FAILURE);
   }
 
-  if ( density > 3. ) { //solid OR enriched. Units of g/mL
+  if ( fdetector->get_extraPhot() ) {
     if ( yields.Lindhard != 1. )
       Nph = int(floor(double(Nph)*7.00+0.5)); //IR photons for NR
     else
-      Nph = int(floor(double(Nph)*1.35+0.5));
+      Nph = int(floor(double(Nph)*1.35+0.5)); //EXO
   }
   result.photons = Nph;
   result.electrons = Ne;
@@ -1238,7 +1239,7 @@ double NESTcalc::GetDensity(double Kelvin,
                                             // (saturated vapor pressure); will
                                             // add pressure dependence later
 
-  if (MOLAR_MASS == 136.) //enriched Xe for 0vBB experiment (136)
+  if (MOLAR_MASS > 134.5) //enrichment for 0vBB expt (~0.8 Xe-136)
     return 3.0305; // ±0.0077 g/cm^3, EXO-200 @167K: arXiv:1908.04128
 
   if (Kelvin < 161.40) {  // solid Xenon

src/VDetector.cpp

@@ -25,11 +25,16 @@ void VDetector::Initialization() {
   noise[0] = 0.0;  // baseline noise mean and width in PE (Gaussian)
   noise[1] = 0.0;  // baseline noise mean and width in PE (Gaussian)
   P_dphe = 0.2;    // chance 1 photon makes 2 phe instead of 1 in Hamamatsu PMT
-
+  
   coinWind = 100;  // S1 coincidence window in ns
   coinLevel = 2;   // how many PMTs have to fire for an S1 to count
   numPMTs = 89;    // For coincidence calculation
 
+  //"Linear noise" terms as defined in Dahl thesis and by D. McK
+  extraPhot=false;  // for matching EXO-200's W measurement
+  noise[2] = 3e-2;  // S1 -> S1 Gaussian-smeared with noise[2]*S1
+  noise[3] = 3e-2;  // S2 -> S2 Gaussian-smeared with noise[3]*S2
+
   // Ionization and Secondary Scintillation (S2) parameters
   g1_gas = 0.06;  // phd per S2 photon in gas, used to get SE size
   s2Fano = 3.61;  // Fano-like fudge factor for SE width

src/testNEST.cpp

@@ -213,7 +213,9 @@ int testNEST(VDetector* detector, unsigned long int numEvts, string type,
   if (rho < 1.75) detector->set_inGas(true);
 
   double Wq_eV = n.WorkFunction(rho).Wq_eV;
-  if ( MOLAR_MASS == 136. ) Wq_eV = 11.5; //11.5±0.5(syst.)±0.1(stat.) from EXO
+  //if ( rho > 3. ) detector->set_extraPhot(true); //solid OR enriched. Units of g/mL
+  if ( detector->get_extraPhot() )
+    Wq_eV = 11.5; //11.5±0.5(syst.)±0.1(stat.) from EXO
 
 
   // Calculate and print g1, g2 parameters (once per detector)
@@ -604,7 +606,7 @@ int testNEST(VDetector* detector, unsigned long int numEvts, string type,
       if (signal1.back() <= 0.) Nph = 0.;
       if (signal2.back() <= 0.) Ne = 0.;
       if (yields.Lindhard > DBL_MIN && Nph > 0. && Ne > 0.) {
-	if ( rho > 3. ) yields.Lindhard = 1.;
+	if ( detector->get_extraPhot() ) yields.Lindhard = 1.;
         keV = (Nph/FudgeFactor[0] + Ne/FudgeFactor[1]) *
 	  Wq_eV * 1e-3 / yields.Lindhard;
         keVee += (Nph + Ne) * Wq_eV * 1e-3;  // as alternative, use W_DEFAULT in
@@ -658,9 +660,9 @@ int testNEST(VDetector* detector, unsigned long int numEvts, string type,
     // including bottom PMTs, NO XYZ correction
     // scint2[5] = S2bc; // floating real# smeared pulse areas in phe ONLY
     // including bottom PMTs, WITH XYZ correction
-    // scint2[6] = S2b / (1.+fdetector->get_P_dphe()); // same as S2b, but
+    // scint2[6] = S2b / (1.+detector->get_P_dphe()); // same as S2b, but
     // adjusted for 2-PE effect (LUX phd units)
-    // scint2[7] = S2bc / (1.+fdetector->get_P_dphe()); // same as S2bc, but
+    // scint2[7] = S2bc / (1.+detector->get_P_dphe()); // same as S2bc, but
     // adjusted for 2-PE effect (LUX phd units)
     // scint2[8] = g2; // g2 = ExtEff * SE, light collection efficiency of EL in
     // gas gap (from CalculateG2)