From 401e303821d6df6be9263297b1413f6606aff61d Mon Sep 17 00:00:00 2001 From: Greg Rischbieter Date: Fri, 8 Jul 2022 10:50:00 -0400 Subject: [PATCH] Synced NESTv2.3.9 and updated GetQuanta Bindings --- HISTORY.md | 16 +++ src/nestpy/DetectorExample_XENON10.hh | 4 +- src/nestpy/LUX_Run03.hh | 4 +- src/nestpy/NEST.cpp | 176 +++++++++++++------------- src/nestpy/NEST.hh | 16 +-- src/nestpy/RandomGen.cpp | 14 +- src/nestpy/RandomGen.hh | 1 + src/nestpy/TestSpectra.hh | 10 +- src/nestpy/__init__.py | 4 +- src/nestpy/analysis.hh | 4 +- src/nestpy/bindings.cpp | 3 +- src/nestpy/execNEST.cpp | 119 ++++++++--------- 12 files changed, 196 insertions(+), 175 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index 81adce1..29b91e5 100644 --- a/HISTORY.md +++ b/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: diff --git a/src/nestpy/DetectorExample_XENON10.hh b/src/nestpy/DetectorExample_XENON10.hh index 7d7b71b..3c6c2b9 100644 --- a/src/nestpy/DetectorExample_XENON10.hh +++ b/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 diff --git a/src/nestpy/LUX_Run03.hh b/src/nestpy/LUX_Run03.hh index abacd86..b582232 100644 --- a/src/nestpy/LUX_Run03.hh +++ b/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; diff --git a/src/nestpy/NEST.cpp b/src/nestpy/NEST.cpp index bec745a..406475b 100644 --- a/src/nestpy/NEST.cpp +++ b/src/nestpy/NEST.cpp @@ -35,7 +35,7 @@ NESTresult NESTcalc::FullCalculation( if (density < 1.) fdetector->set_inGas(true); NESTresult result; result.yields = GetYields(species, energy, density, dfield, A, Z, NRYieldsParam); - result.quanta = GetQuanta(result.yields, density, NRERWidthsParam); + result.quanta = GetQuanta(result.yields, density, NRERWidthsParam, false, -999.); if (do_times) result.photon_times = GetPhotonTimes( species, result.quanta.photons, result.quanta.excitons, dfield, energy); @@ -65,9 +65,9 @@ double NESTcalc::PhotonTime(INTERACTION_TYPE species, bool exciton, if (ValidityTests::nearlyEqual(ATOM_NUM, 18.)) { // copied from 2013 NEST version for LAr on LBNE tau1 = - RandomGen::rndm()->rand_gauss(6.5, 0.8); // error from weighted average - tau3 = RandomGen::rndm()->rand_gauss(1300, 50); // ibid. - tauR = RandomGen::rndm()->rand_gauss(0.8, 0.2); // Kubota 1979 + RandomGen::rndm()->rand_gauss(6.5, 0.8, true); // error from weighted average + tau3 = RandomGen::rndm()->rand_gauss(1300, 50, true); // ibid. + tauR = RandomGen::rndm()->rand_gauss(0.8, 0.2, true); // Kubota 1979 if (species <= Cf) SingTripRatio = 0.22218 * pow(energy, 0.48211); else if (species == ion) { // really only alphas here @@ -78,9 +78,9 @@ double NESTcalc::PhotonTime(INTERACTION_TYPE species, bool exciton, SingTripRatio = 0.2701 + 0.003379 * LET - 4.7338e-5 * pow(LET, 2.) + 8.1449e-6 * pow(LET, 3.); if (LET < 3. && !exciton) - SingTripRatio = RandomGen::rndm()->rand_gauss(0.5, 0.2); + SingTripRatio = RandomGen::rndm()->rand_gauss(0.5, 0.2, true); if (LET < 3. && exciton) - SingTripRatio = RandomGen::rndm()->rand_gauss(0.36, 0.06); + SingTripRatio = RandomGen::rndm()->rand_gauss(0.36, 0.06, true); } // lastly is ER for LAr } else { if (species <= Cf) // NR @@ -232,10 +232,8 @@ double NESTcalc::FanoER(double density, double Nq_mean, double efield, return Fano; } -QuantaResult NESTcalc::GetQuanta( - const YieldResult &yields, double density, - const std::vector &NRERWidthsParam /*={1.,1.,0.1,0.5,0.19,2.25, 0.0015, 0.0553, 0.205, 0.45, -0.2}*/, - bool oldModelER) { +QuantaResult NESTcalc::GetQuanta(const YieldResult &yields, double density, const std::vector &NRERWidthsParam, + bool oldModelER, double SkewnessER) { QuantaResult result{}; bool HighE; int Nq_actual, Ne, Nph, Ni, Nex; @@ -265,26 +263,26 @@ QuantaResult NESTcalc::GetQuanta( } if (ValidityTests::nearlyEqual(yields.Lindhard, 1.)) { - double Fano = FanoER(density, Nq_mean, yields.ElectricField, NRERWidthsParam); - Nq_actual = int(floor( - RandomGen::rndm()->rand_gauss(Nq_mean, sqrt(Fano * Nq_mean)) + 0.5)); - if (Nq_actual < 0 || ValidityTests::nearlyEqual(Nq_mean, 0.)) Nq_actual = 0; - + if (ValidityTests::nearlyEqual(Nq_mean, 0.)) Nq_actual = 0; + else { + double Fano = FanoER(density, Nq_mean, yields.ElectricField, NRERWidthsParam); + Nq_actual = int(floor( + RandomGen::rndm()->rand_gauss(Nq_mean, sqrt(Fano * Nq_mean), true) + 0.5)); + } + Ni = RandomGen::rndm()->binom_draw(Nq_actual, alf); Nex = Nq_actual - Ni; } else { double Fano = NRERWidthsParam[0]; Ni = int(floor(RandomGen::rndm()->rand_gauss(Nq_mean * alf, - sqrt(Fano * Nq_mean * alf)) + + sqrt(Fano * Nq_mean * alf), true) + 0.5)); - if (Ni < 0) Ni = 0; Fano = NRERWidthsParam[1]; Nex = int(floor(RandomGen::rndm()->rand_gauss( Nq_mean * excitonRatio * alf, - sqrt(Fano * Nq_mean * excitonRatio * alf)) + + sqrt(Fano * Nq_mean * excitonRatio * alf), true) + 0.5)); - if (Nex < 0) Nex = 0; Nq_actual = Nex + Ni; } @@ -351,11 +349,18 @@ QuantaResult NESTcalc::GetQuanta( skewness = 0.; if (ValidityTests::nearlyEqual(yields.Lindhard, 1.)) { - skewness = 1. / (1. + exp((engy - E2) / E3)) * - (alpha0 + - cc0 * exp(-1. * fld / F0) * (1. - exp(-1. * engy / E0))) + - 1. / (1. + exp(-1. * (engy - E2) / E3)) * cc1 * - exp(-1. * engy / E1) * exp(-1. * sqrt(fld) / sqrt(F1)); + + if ( SkewnessER != -999. ) { + skewness = SkewnessER; + } + else { + skewness = 1. / (1. + exp((engy - E2) / E3)) * + (alpha0 + + cc0 * exp(-1. * fld / F0) * (1. - exp(-1. * engy / E0))) + + 1. / (1. + exp(-1. * (engy - E2) / E3)) * cc1 * + exp(-1. * engy / E1) * exp(-1. * sqrt(fld) / sqrt(F1)); + } + // if ( std::abs(skewness) <= DBL_MIN ) skewness = DBL_MIN; } else { skewness = NRERWidthsParam[5]; // 2.25 but ~5-20 also good (for NR). All better @@ -377,7 +382,7 @@ QuantaResult NESTcalc::GetQuanta( 0.5)); else { Ne = int(floor( - RandomGen::rndm()->rand_gauss((1. - recombProb) * Ni, sqrt(Variance)) + + RandomGen::rndm()->rand_gauss((1. - recombProb) * Ni, sqrt(Variance), true) + 0.5)); } @@ -579,14 +584,14 @@ NESTresult NESTcalc::GetYieldERdEOdxBasis(const std::vector &dEOdxParam, Ne += result.quanta.electrons; } else - result.quanta = GetQuanta(result.yields, rho, default_NRERWidthsParam, true); + result.quanta = GetQuanta(result.yields, rho, default_NRERWidthsParam, false, -999.); if (eMin > 0.) Nph += result.quanta.photons * (eStep / refEnergy); else { if ( NRERWidthsParam[0] >= 0. ) { Nq = result.quanta.photons + result.quanta.electrons; - double FanoOverall = 0.0015 * sqrt((-1e3*eMin/Wq_eV)*field); - double FanoScint = 20. * FanoOverall; + double FanoOverall = 1.; //use 6 for ~0.7% energy resolution at 2.6 MeV and 190 V/cm + double FanoScint = 200.; //standin for recombination fluctuations in this approach Nq = int(floor(RandomGen::rndm()->rand_gauss(Nq,sqrt(FanoOverall*Nq),false)+0.5)); result.quanta.photons = int(floor(RandomGen::rndm()->rand_gauss( result.quanta.photons,sqrt(FanoScint*result.quanta.photons),false)+0.5)); @@ -603,7 +608,7 @@ NESTresult NESTcalc::GetYieldERdEOdxBasis(const std::vector &dEOdxParam, driftTime = (fdetector->get_TopDrift() - zz) / vD; if ( zz >= fdetector->get_cathode() && NRERWidthsParam[0] >= 0. ) { if (eMin > 0.) - Ne += result.quanta.electrons * (eStep / refEnergy) * exp(-driftTime / fdetector->get_eLife_us()); + Ne += result.quanta.electrons * exp(-driftTime / fdetector->get_eLife_us()) * (eStep / refEnergy); else Ne += result.quanta.electrons * exp(-driftTime / fdetector->get_eLife_us()); @@ -769,8 +774,8 @@ YieldResult NESTcalc::GetYieldNR( Ni = Nq / (1. + NexONi); Nex = Nq - Ni; } - - if (Nex < 0.) + + if ( Nex < 0. && density >= 1. ) cerr << "\nCAUTION: You are approaching the border of NEST's validity for " "high-energy (OR, for LOW) NR, or are beyond it, at " << energy << " keV." << endl; @@ -986,7 +991,6 @@ YieldResult NESTcalc::GetYieldKr83m(double energy, double density, YieldResult NESTcalc::GetYieldBeta(double energy, double density, double dfield) { // OLD - Wvalue wvalue = WorkFunction(density, fdetector->get_molarMass(), fdetector->get_OldW13eV()); double Qy, Nq; @@ -1154,7 +1158,7 @@ YieldResult NESTcalc::GetYields( break; default: // beta, CH3T, 14C, the pp solar neutrino background, and // Compton/PP spectra of fullGamma - if (ValidityTests::nearlyEqual(ATOM_NUM, 18.) || oldModelER) + if (ValidityTests::nearlyEqual(ATOM_NUM, 18.) || oldModelER || fdetector->get_inGas()) return GetYieldBeta(energy, density, dfield); // OLD else return GetYieldBetaGR(energy, density, dfield, NRYieldsParam); // NEW @@ -1164,8 +1168,10 @@ YieldResult NESTcalc::GetYields( YieldResult NESTcalc::YieldResultValidity(YieldResult &res, const double energy, const double Wq_eV) { - assert(res.ElectronYield != -999 && res.PhotonYield != -999 && - res.ExcitonRatio != -999); + assert(res.ElectronYield != -999 && + res.PhotonYield != -999 && + res.ExcitonRatio != -999 + ); if (res.PhotonYield > energy / W_SCINT) res.PhotonYield = energy / W_SCINT; // yields can never exceed 1 / [ W ~ few eV ] @@ -1201,7 +1207,7 @@ const vector &NESTcalc::GetS1( double truthPosZ, double smearPosX, double smearPosY, double smearPosZ, double driftVelocity, double dV_mid, INTERACTION_TYPE type_num, uint64_t evtNum, double dfield, double energy, S1CalculationMode mode, - bool outputTiming, vector &wf_time, vector &wf_amp) { + int outputTiming, vector &wf_time, vector &wf_amp) { int Nph = quanta.photons; double subtract[2] = {0., 0.}; @@ -1269,16 +1275,9 @@ const vector &NESTcalc::GetS1( current_mode = energy < 100. ? S1CalculationMode::Full : S1CalculationMode::Parametric; } - - // Follow https://en.wikipedia.org/wiki/Truncated_normal_distribution - double TruncGaussAlpha = -1. / fdetector->get_sPEres(); - double LittlePhi_Alpha = - inv_sqrt2_PI * exp(-0.5 * TruncGaussAlpha * TruncGaussAlpha); - double BigPhi_Alpha = 0.5 * (1. + erf(TruncGaussAlpha / sqrt2)); - double NewMean = - 1. + (LittlePhi_Alpha / (1. - BigPhi_Alpha)) * fdetector->get_sPEres(); - switch (current_mode) { - case S1CalculationMode::Full: + + double NewMean = RandomGen::rndm()->FindNewMean(fdetector->get_sPEres()); + switch (current_mode) { case S1CalculationMode::Full: case S1CalculationMode::Waveform: { // Step through the pmt hits for (int i = 0; i < nHits; ++i) { @@ -1289,7 +1288,8 @@ const vector &NESTcalc::GetS1( double phe1 = TruncGauss + RandomGen::rndm()->rand_gauss( fdetector->get_noiseBaseline()[0], - fdetector->get_noiseBaseline()[1]); + fdetector->get_noiseBaseline()[1], + false); ++Nphe; if (phe1 > DBL_MAX) { phe1 = 1.; // for Box-Mueller fuck-ups @@ -1310,7 +1310,8 @@ const vector &NESTcalc::GetS1( phe2 = TruncGauss + RandomGen::rndm()->rand_gauss( fdetector->get_noiseBaseline()[0], - fdetector->get_noiseBaseline()[1]); + fdetector->get_noiseBaseline()[1], + false); ++Nphe; if (phe2 > DBL_MAX) { phe2 = 1.; @@ -1388,7 +1389,7 @@ const vector &NESTcalc::GetS1( Nphe, 1. - (1. - eff) / (1. + fdetector->get_P_dphe()))); // take into account the truncation of the PE distributions pulseArea = RandomGen::rndm()->rand_gauss( - Nphe_det * (1. / NewMean), fdetector->get_sPEres() * sqrt(Nphe_det)); + Nphe_det * (1. / NewMean), fdetector->get_sPEres() * sqrt(Nphe_det/NewMean), true); spike = static_cast(RandomGen::rndm()->binom_draw( nHits, eff * (1. - belowThresh_percentile))); @@ -1416,7 +1417,7 @@ const vector &NESTcalc::GetS1( photonstream photon_times = AddPhotonTransportTime( photon_emission_times, truthPosX, truthPosY, truthPosZ); - if (outputTiming && !pulseFile.is_open()) { + if (outputTiming > 0 && !pulseFile.is_open()) { pulseFile.open("photon_times.txt"); // pulseFile << "Event#\tt [ns]\tA [PE]" << endl; pulseFile << "Event#\tt [ns]\tPEb/bin\tPEt/bin\tin win" << endl; @@ -1471,7 +1472,7 @@ const vector &NESTcalc::GetS1( wf_time.emplace_back((ii - numPts / 2) * SAMPLE_SIZE); wf_amp.emplace_back(AreaTable[0][ii] + AreaTable[1][ii]); - if (outputTiming) { + if (outputTiming > 0) { char line[80]; if (AreaTable[0][ii] > PHE_MAX) { subtract[0] = AreaTable[0][ii] - PHE_MAX; @@ -1516,10 +1517,10 @@ const vector &NESTcalc::GetS1( pulseArea, sqrt( pow(fdetector->get_noiseQuadratic()[0] * pulseArea * pulseArea, 2) + - pow(fdetector->get_noiseLinear()[0] * pulseArea, 2))); + pow(fdetector->get_noiseLinear()[0] * pulseArea, 2)), true); } else { pulseArea = RandomGen::rndm()->rand_gauss( - pulseArea, fdetector->get_noiseLinear()[0] * pulseArea); + pulseArea, fdetector->get_noiseLinear()[0] * pulseArea, true); } pulseArea -= (subtract[0] + subtract[1]); if (pulseArea < fdetector->get_sPEthr()) pulseArea = 0.; @@ -1625,7 +1626,7 @@ const vector &NESTcalc::GetS2( int Ne, double truthPosX, double truthPosY, double truthPosZ, double smearPosX, double smearPosY, double smearPosZ, double dt, double driftVelocity, uint64_t evtNum, double dfield, - S2CalculationMode mode, bool outputTiming, vector &wf_time, + S2CalculationMode mode, int outputTiming, vector &wf_time, vector &wf_amp, const vector &g2_params) { double elYield = g2_params[0]; double ExtEff = g2_params[1]; @@ -1728,8 +1729,8 @@ const vector &NESTcalc::GetS2( double min = 1e100; for (i = 0; i < stopPoint; ++i) { elecTravT = 0.; // resetting for the current electron - DL = RandomGen::rndm()->rand_gauss(0., sigmaDL); - DT = std::abs(RandomGen::rndm()->rand_gauss(0., sigmaDT)); + DL = RandomGen::rndm()->rand_gauss(0., sigmaDL, false); + DT = std::abs(RandomGen::rndm()->rand_gauss(0., sigmaDT, false)); phi = two_PI * RandomGen::rndm()->rand_uniform(); sigX = DT * cos(phi); sigY = DT * sin(phi); @@ -1770,8 +1771,8 @@ const vector &NESTcalc::GetS2( SE += (double)RandomGen::rndm()->binom_draw(uint64_t(SE), fdetector->get_P_dphe()); Nphe += uint64_t(SE); - DL = RandomGen::rndm()->rand_gauss(0., sigmaDLg); - DT = RandomGen::rndm()->rand_gauss(0., sigmaDTg); + DL = RandomGen::rndm()->rand_gauss(0., sigmaDLg, false); + DT = std::abs(RandomGen::rndm()->rand_gauss(0., sigmaDTg, false)); phi = two_PI * RandomGen::rndm()->rand_uniform(); sigX = DT * cos(phi); sigY = DT * sin(phi); @@ -1795,7 +1796,7 @@ const vector &NESTcalc::GetS2( } } for (double j = 0.; j < SE; j += 1.) { - double phe = RandomGen::rndm()->rand_gauss(1., fdetector->get_sPEres()); + double phe = RandomGen::rndm()->rand_zero_trunc_gauss(1./RandomGen::rndm()->FindNewMean(fdetector->get_sPEres()), fdetector->get_sPEres()); if (i < Nee || !eTrain) pulseArea += phe; origin = fdetector->get_TopDrift() + gasGap * RandomGen::rndm()->rand_uniform(); @@ -1855,7 +1856,7 @@ const vector &NESTcalc::GetS2( wf_time.emplace_back(k * SAMPLE_SIZE + int64_t(min + SAMPLE_SIZE / 2.)); wf_amp.emplace_back(AreaTableBot[1][k] + AreaTableTop[1][k]); - if (outputTiming) { + if (outputTiming > 0) { char line[80]; if (AreaTableBot[1][k] > PHE_MAX) subtract[0] = AreaTableBot[1][k] - PHE_MAX; @@ -1883,8 +1884,9 @@ const vector &NESTcalc::GetS2( RandomGen::rndm()->binom_draw(Nph, fdetector->get_g1_gas() * posDep); Nphe = nHits + RandomGen::rndm()->binom_draw(nHits, fdetector->get_P_dphe()); + double NewMean = RandomGen::rndm()->FindNewMean(fdetector->get_sPEres()); pulseArea = RandomGen::rndm()->rand_gauss( - Nphe, fdetector->get_sPEres() * sqrt(Nphe)); + Nphe/NewMean, fdetector->get_sPEres() * sqrt(Nphe/NewMean), true); } if (fdetector->get_noiseLinear()[1] >= 0.1) @@ -1895,10 +1897,10 @@ const vector &NESTcalc::GetS2( pulseArea = RandomGen::rndm()->rand_gauss( pulseArea, sqrt(pow(fdetector->get_noiseQuadratic()[1] * pow(pulseArea, 2), 2) + - pow(fdetector->get_noiseLinear()[1] * pulseArea, 2))); + pow(fdetector->get_noiseLinear()[1] * pulseArea, 2)), true); } else { pulseArea = RandomGen::rndm()->rand_gauss( - pulseArea, fdetector->get_noiseLinear()[1] * pulseArea); + pulseArea, fdetector->get_noiseLinear()[1] * pulseArea, true); } pulseArea -= (subtract[0] + subtract[1]); double pulseAreaC = @@ -1909,7 +1911,7 @@ const vector &NESTcalc::GetS2( double S2b = RandomGen::rndm()->rand_gauss( fdetector->FitTBA(truthPosX, truthPosY, truthPosZ)[1] * pulseArea, sqrt(fdetector->FitTBA(truthPosX, truthPosY, truthPosZ)[1] * pulseArea * - (1. - fdetector->FitTBA(truthPosX, truthPosY, truthPosZ)[1]))); + (1. - fdetector->FitTBA(truthPosX, truthPosY, truthPosZ)[1])), true); double S2bc = S2b / exp(-dt / fdetector->get_eLife_us()) / posDepSm; // for detectors using S2 bottom-only in their analyses @@ -1960,7 +1962,7 @@ const vector &NESTcalc::GetS2( return ionization; } -vector NESTcalc::CalculateG2(bool verbosity) { +vector NESTcalc::CalculateG2(int verbosity) { vector g2_params(5); // Set parameters for calculating EL yield and extraction @@ -2063,7 +2065,7 @@ vector NESTcalc::CalculateG2(bool verbosity) { double g2 = ExtEff * SE; double StdDev = 0., Nphe, pulseArea, pulseAreaC, NphdC, phi, posDep, r, x, y; int Nph, nHits; - if (verbosity) { + if (verbosity > 0) { for (int i = 0; i < 10000; ++i) { // calculate properly the width (1-sigma std dev) in the SE size Nph = int(floor(RandomGen::rndm()->rand_gauss( @@ -2081,16 +2083,16 @@ vector NESTcalc::CalculateG2(bool verbosity) { Nphe = nHits + RandomGen::rndm()->binom_draw(nHits, fdetector->get_P_dphe()); pulseArea = RandomGen::rndm()->rand_gauss( - Nphe, fdetector->get_sPEres() * sqrt(Nphe)); + Nphe, fdetector->get_sPEres() * sqrt(Nphe), true); if (fdetector->get_noiseQuadratic()[1] != 0) { pulseArea = RandomGen::rndm()->rand_gauss( pulseArea, sqrt( pow(fdetector->get_noiseQuadratic()[1] * pow(pulseArea, 2), 2) + - pow(fdetector->get_noiseLinear()[1] * pulseArea, 2))); + pow(fdetector->get_noiseLinear()[1] * pulseArea, 2)), true); } else { pulseArea = RandomGen::rndm()->rand_gauss( - pulseArea, fdetector->get_noiseLinear()[1] * pulseArea); + pulseArea, fdetector->get_noiseLinear()[1] * pulseArea, true); } if (fdetector->get_s2_thr() < 0.) pulseArea = RandomGen::rndm()->rand_gauss( @@ -2100,7 +2102,7 @@ vector NESTcalc::CalculateG2(bool verbosity) { fdetector->FitTBA(0.0, 0.0, fdetector->get_TopDrift() / 2.)[1] * pulseArea * (1. - fdetector->FitTBA(0.0, 0.0, - fdetector->get_TopDrift() / 2.)[1]))); + fdetector->get_TopDrift() / 2.)[1])), true); pulseAreaC = pulseArea / posDep; NphdC = pulseAreaC / (1. + fdetector->get_P_dphe()); StdDev += (SE - NphdC) * (SE - NphdC); @@ -2141,12 +2143,7 @@ const vector &NESTcalc::GetSpike(int Nph, double dx, double dy, newSpike[1] = oldScint[5]; return newSpike; } - newSpike[0] = std::abs(oldScint[6]); - double TruncGauss = 0.; - while (TruncGauss <= 0.) - TruncGauss = RandomGen::rndm()->rand_gauss( - newSpike[0], (fdetector->get_sPEres() / 4.) * sqrt(newSpike[0])); - newSpike[0] = TruncGauss; + newSpike[0] = RandomGen::rndm()->rand_zero_trunc_gauss(oldScint[6], (fdetector->get_sPEres()/4.)*sqrt(std::abs(oldScint[6]))); if (XYcorr == 0 || XYcorr == 2) { dx = 0.; dy = 0.; @@ -2500,13 +2497,13 @@ vector NESTcalc::xyResolution(double xPos_mm, double yPos_mm, 9.); // σ_vertex^2 = ((r_vertex/r_S2)*σ_(x,y))^2 + σ_sys^2 double phi = two_PI * RandomGen::rndm()->rand_uniform(); - sigmaR = std::abs(RandomGen::rndm()->rand_gauss(0.0, sigmaR)); + sigmaR = std::abs(RandomGen::rndm()->rand_gauss(0.0, sigmaR, false)); double sigmaX = sigmaR * cos(phi); double sigmaY = sigmaR * sin(phi); if (sigmaR > 1e2 || std::isnan(sigmaR) || sigmaR <= 0. || std::abs(sigmaX) > 1e2 || std::abs(sigmaY) > 1e2) { - if (A_top > 20.) { + if (A_top > 40.) { // roughly two S2 electrons in most detectors (~1e- in LZ) cerr << "WARNING: your position resolution is worse than 10 cm. Is that " "correct?!" << endl; @@ -2528,24 +2525,24 @@ double NESTcalc::PhotonEnergy(bool s2Flag, bool state, double tempK) { if (ValidityTests::nearlyEqual(ATOM_NUM, 18.)) { // liquid Argon - if (state) return RandomGen::rndm()->rand_gauss(9.7, 0.2); - return RandomGen::rndm()->rand_gauss(9.69, 0.22); + if (state) return RandomGen::rndm()->rand_zero_trunc_gauss(9.7, 0.2); + return RandomGen::rndm()->rand_zero_trunc_gauss(9.69, 0.22); } if (!state) // liquid or solid - wavelength = RandomGen::rndm()->rand_gauss( - 178., 14. / 2.355); // taken from Jortner JchPh 42 '65 + wavelength = RandomGen::rndm()->rand_zero_trunc_gauss( + 178., 14. / 2.355); // taken from Jortner JchPh 42 '65 else { // gas wavelength = - RandomGen::rndm()->rand_gauss(175., 5.); // G4S1Light, probably Doke + RandomGen::rndm()->rand_zero_trunc_gauss(175., 5.); // G4S1Light, probably Doke } if (s2Flag) { // S2 different from ordinary gas (or just measurement error?) if (tempK < 200.) // cold gas - wavelength = RandomGen::rndm()->rand_gauss( - 179., 5.); // source is G4S2Light.cc from the old NEST + wavelength = RandomGen::rndm()->rand_zero_trunc_gauss( + 179., 5.); // source is G4S2Light.cc from the old NEST else { - wavelength = RandomGen::rndm()->rand_gauss(174., 5.); // ditto + wavelength = RandomGen::rndm()->rand_zero_trunc_gauss(174., 5.); // ditto } } @@ -2569,7 +2566,7 @@ double NESTcalc::CalcElectronLET(double E, int Z, bool CSDA) { LET=1.8106-0.45086*log10(E)-.33151*pow(log10(E),2.)+.25916*pow(log10(E),3.)-0.2051*pow(log10(E),4.)+0.15279*pow(log10(E),5.)-.084659*pow(log10(E),6.) +0.030441*pow(log10(E),7.)-.0058953*pow(log10(E),8.)+0.00045633*pow(log10(E),9.); LET = pow(10.,LET); - if ( std::isnan(LET) || LET <= 0. ) LET = 1e2; + if ( std::isnan(LET) || LET <= 1. ) LET = 1e2; return LET; } @@ -2689,10 +2686,7 @@ double NESTcalc::GetDiffLong_Liquid( // Use the standard NEST parametrization if (!highFieldModel) { - output = 345.92 * pow(dfield, -0.47880) * - exp(-81.3230 / dfield); // fit to Aprile & Doke review - // paper and to arXiv:1102.2865; - // plus, LUX Run02+03 + output = 57.381*pow(dfield,-0.22221)+127.27*exp(-dfield/32.821); //fit to Aprile & Doke rev & arXiv:1102.2865 (Peter: XENON10/100); plus, LUX Run03 (181V/cm) & 1911.11580 } // Use the Boyle model, which is drastically different at high (>5kV/cm) // fields. Note here that the Boyle model is only at one temperature. First diff --git a/src/nestpy/NEST.hh b/src/nestpy/NEST.hh index 01a96ca..93e9d17 100644 --- a/src/nestpy/NEST.hh +++ b/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 &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 &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 &wf_time, vector &wf_amp); + int outputTiming, vector &wf_time, vector &wf_amp); // Very comprehensive conversion of the "original" intrinsic scintillation // photons into the many possible definitions of S1 as measured by // photo-sensors @@ -349,7 +349,7 @@ 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 &wf_time, vector &wf_amp, + int outputTiming, vector &wf_time, vector &wf_amp, const vector &g2_params); // Exhaustive conversion of the intrinsic ionization electrons into the many // possible definitions of S2 pulse areas as observed in the photo-tubes @@ -357,7 +357,7 @@ class NESTcalc { // electron mean free path or life time caused by electronegative impurities // (exponential) - std::vector CalculateG2(bool verbosity = true); + std::vector 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. diff --git a/src/nestpy/RandomGen.cpp b/src/nestpy/RandomGen.cpp index e6d21be..372ea14 100644 --- a/src/nestpy/RandomGen.cpp +++ b/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; diff --git a/src/nestpy/RandomGen.hh b/src/nestpy/RandomGen.hh index b604fee..acdee89 100644 --- a/src/nestpy/RandomGen.hh +++ b/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); diff --git a/src/nestpy/TestSpectra.hh b/src/nestpy/TestSpectra.hh index 2909a14..f63978f 100644 --- a/src/nestpy/TestSpectra.hh +++ b/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 diff --git a/src/nestpy/__init__.py b/src/nestpy/__init__.py index 59857ad..da650b7 100644 --- a/src/nestpy/__init__.py +++ b/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 * diff --git a/src/nestpy/analysis.hh b/src/nestpy/analysis.hh index 40eca78..6cee251 100644 --- a/src/nestpy/analysis.hh +++ b/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 diff --git a/src/nestpy/bindings.cpp b/src/nestpy/bindings.cpp index 5b0fa5d..11da8d2 100644 --- a/src/nestpy/bindings.cpp +++ b/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({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) diff --git a/src/nestpy/execNEST.cpp b/src/nestpy/execNEST.cpp index fd650e9..2026388 100644 --- a/src/nestpy/execNEST.cpp +++ b/src/nestpy/execNEST.cpp @@ -39,8 +39,8 @@ int main(int argc, char** argv) { // Instantiate your own VDetector class here, then load into NEST class // constructor auto* detector = new DetectorExample_LUX_RUN03(); - if (verbosity) cerr << "*** Detector definition message ***" << endl; - if (verbosity) + if (verbosity > 0) cerr << "*** Detector definition message ***" << endl; + if (verbosity > 0) cerr << "You are currently using the LUX Run03 template detector." << endl << endl; // Custom parameter modification functions @@ -48,7 +48,7 @@ int main(int argc, char** argv) { if (ValidityTests::nearlyEqual(ATOM_NUM, 18.)) { detector->set_molarMass(39.948); - if (verbosity) + if (verbosity > 0) cerr << "\nWARNING: Argon is currently only in alpha testing mode!! Many " "features copied over from Xenon wholesale still. Use models at " "your own risk.\n" @@ -105,7 +105,7 @@ int main(int argc, char** argv) { type = "MIP"; dEOdxBasis = true; eMin = -TestSpectra::CH3T_spectrum(0.,18.6); - //s1CalculationMode = NEST::S1CalculationMode::Parametric; + s1CalculationMode = NEST::S1CalculationMode::Parametric; eMax = -1.; inField = 180.; position = "-1"; @@ -116,7 +116,7 @@ int main(int argc, char** argv) { NRERWidthsParam.clear(); NRYieldsParam.clear(); ERWeightParam.clear(); - verbosity = false; + verbosity = -1; NRYieldsParam.push_back(0.0); NRYieldsParam.push_back(0.); @@ -156,7 +156,7 @@ int main(int argc, char** argv) { } else { numEvts = (uint64_t)atof(argv[1]); if (numEvts <= 0) { - if (verbosity) + if (verbosity > 0) cerr << "ERROR, you must simulate at least 1 event, or 1 kg*day" << endl; return 1; @@ -252,7 +252,7 @@ NESTObservableArray runNESTvec( particleType, // func suggested by Xin Xiang, PD Brown U. for RG, LZ vector eList, vector> pos3dxyz, double inField, int seed) { - verbosity = false; + verbosity = -1; NESTcalc n(detector); NESTresult result; QuantaResult quanta; @@ -346,11 +346,11 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, // Construct NEST class using detector object NESTcalc n(detector); NRYieldsParam = {11., 1.1, 0.0480, -0.0533, 12.6, 0.3, 2., 0.3, 2., 0.5, 1., 1.}; - NRERWidthsParam = {1.,1.,0.1,0.5,0.19,2.25, 0.0015, 0.0553, 0.205, 0.45, -0.2}; - ERWeightParam = {0.23, 0.77, 2.95, -1.44, 1.0, 1.0, 0., 0.}; + NRERWidthsParam = {1.,1.,0.1,0.5,0.19,2.25, 0.0015, 0.0553, 0.205, 0.45, -0.2}; + ERWeightParam = {0.23, 0.77, 2.95, -1.44, 1.0, 1.0, 0., 0.}; if (detector->get_TopDrift() <= 0. || detector->get_anode() <= 0. || detector->get_gate() <= 0.) { - if (verbosity) + if (verbosity > 0) cerr << "ERROR, unphysical value(s) of position within the detector " "geometry."; // negative or 0 for cathode position is OK (e.g., // LZ) @@ -384,7 +384,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, eMax = hiEregime; // the default energy max if (ValidityTests::nearlyEqual(eMax, 0.) && type != "muon" && type != "MIP" && type != "LIP" && type != "mu" && type != "mu-") { - if (verbosity) + if (verbosity > 0) cerr << "ERROR: The maximum energy (or Kr time sep) cannot be 0 keV (or " "0 ns)!" << endl; @@ -399,7 +399,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, type_num = NR; //-1: default particle type is also NR else if (type == "WIMP") { if (eMin < 0.44) { // here eMin is WIMP mass - if (verbosity) cerr << "WIMP mass too low, you're crazy!" << endl; + if (verbosity > 0) cerr << "WIMP mass too low, you're crazy!" << endl; return 1; } type_num = WIMP; @@ -464,18 +464,18 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, numEvts = RandomGen::rndm()->poisson_draw(1.5764e-7 * double(numEvts)); } else if (type == "fullGamma") { type_num = fullGamma; - if (verbosity) + if (verbosity > 0) cerr << "Please choose gamma source. The allowed sources " "are:\n\"Co57\"\n\"Co60\"\n\"Cs137\"\nSource: "; cin >> gamma_source; - if (gamma_source == "Co60" && verbosity) { + if (gamma_source == "Co60" && verbosity > 0) { cerr << "WARNING: This source is in the pair production range. " "Electron/positron pairs are not accounted for after initial " "interaction, and some " << "photons and electrons may go unaccounted." << endl; } } else { - if (verbosity) { + if (verbosity > 0) { string particleTypes = "UNRECOGNIZED PARTICLE TYPE!! VALID OPTIONS ARE:\n" "NR or neutron,\n" @@ -510,11 +510,11 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, eMin != eMax) { maxTimeSep = eMax; if (eMax <= 0.) { - if (verbosity) cerr << "Max t sep must be +." << endl; + if (verbosity > 0) cerr << "Max t sep must be +." << endl; return 1; } } else { - if (verbosity) + if (verbosity > 0) cerr << "ERROR: For Kr83m, put E_min as 9.4, 32.1, or 41.5 keV and " "E_max as the max time-separation [ns] between the two decays " "please." @@ -523,7 +523,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, } } - if ((eMin < 10. || eMax < 10.) && type_num == gammaRay && verbosity) { + if ((eMin < 10. || eMax < 10.) && type_num == gammaRay && verbosity > 0) { cerr << "WARNING: Typically beta model works better for ER BG at low " "energies as in a WS." << endl; @@ -535,7 +535,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, double rho = n.SetDensity(detector->get_T_Kelvin(), detector->get_p_bar()); if (rho <= 0. || detector->get_T_Kelvin() <= 0. || detector->get_p_bar() <= 0.) { - if (verbosity) cerr << "ERR: Unphysical thermodynamic property!"; + if (verbosity > 0) cerr << "ERR: Unphysical thermodynamic property!"; return 1; } if (rho < 1.75 && ValidityTests::nearlyEqual(ATOM_NUM, 54.)) @@ -584,7 +584,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, double(massNum), double(atomNum), NRYieldsParam); } if ((g1 * yieldsMax.PhotonYield) > (2. * maxS1) && eMin != eMax && - type_num != Kr83m && verbosity && !dEOdxBasis) + type_num != Kr83m && verbosity > 0 && !dEOdxBasis) cerr << "\nWARNING: Your energy maximum may be too high given your maxS1.\n"; @@ -690,8 +690,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if (eMax > eMin) keV = eMin + (eMax - eMin) * RandomGen::rndm()->rand_uniform(); else // polymorphic feature: Gauss E distribution - keV = RandomGen::rndm()->rand_gauss(eMin, eMax); - if (keV < 0.) keV = 1e-3 * Wq_eV; + keV = RandomGen::rndm()->rand_zero_trunc_gauss(eMin, eMax); } else { // negative eMax signals to NEST that you want to use an // exponential energy spectrum profile if (ValidityTests::nearlyEqual(eMin, 0.)) return 1; @@ -711,12 +710,12 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if (!dEOdxBasis) { if (keV < 0) { - if (verbosity) + if (verbosity > 0) cerr << "ERROR: Get ready for time travel or FTL--negative energy!" << endl; return 1; } - if (ValidityTests::nearlyEqual(keV, 0.) && verbosity) { + if (ValidityTests::nearlyEqual(keV, 0.) && verbosity > 0) { cerr << "WARNING: Zero energy has occurred, and this is not normal" << endl; } @@ -749,12 +748,12 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if (pos_x != -999. && pos_y != -999. && (pos_x * pos_x + pos_y * pos_y) > detector->get_radius() * detector->get_radius() && - j == 0 && verbosity) + j == 0 && verbosity > 0) cerr << "WARNING: outside fiducial radius." << endl; if (pos_x != -999. && pos_y != -999. && (pos_x * pos_x + pos_y * pos_y) > detector->get_radmax() * detector->get_radmax()) { - if (verbosity) + if (verbosity > 0) cerr << "\nERROR: outside physical radius!!!" << endl; return EXIT_FAILURE; } @@ -781,17 +780,17 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, field = inField; // no fringing if (field < 0. || detector->get_E_gas() < 0.) { - if (verbosity) + if (verbosity > 0) cerr << "\nERROR: Neg field is not permitted. We don't simulate " "field dir (yet). Put in magnitude.\n"; return 1; } if ((ValidityTests::nearlyEqual(field, 0.) || std::isnan(field)) && - verbosity) + verbosity > 0) cerr << "\nWARNING: A LITERAL ZERO (or undefined) FIELD MAY YIELD WEIRD " "RESULTS. USE A SMALL VALUE INSTEAD.\n"; - if ((field > 12e3 || detector->get_E_gas() > 17e3) && verbosity) + if ((field > 12e3 || detector->get_E_gas() > 17e3) && verbosity > 0) cerr << "\nWARNING: Your field is >12,000 V/cm. No data out here. Are " "you sure about this?\n"; @@ -808,7 +807,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, n.SetDriftVelocity(detector->get_T_Kelvin(), rho, inField); vD = n.SetDriftVelocity(detector->get_T_Kelvin(), rho, field); } - if (verbosity) { + if (verbosity > 0) { cout << "Density = " << rho << " g/mL" << "\t"; cout << "central vDrift = " << vD_middle << " mm/us\n"; @@ -823,7 +822,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if ((ValidityTests::nearlyEqual(eMax, eMin) || type_num == Kr83m) && numBins == 1 && MCtruthE) { MCtruthE = false; - if (verbosity) + if (verbosity > 0) fprintf(stderr, "Simulating a mono-E peak; setting MCtruthE false.\n"); } @@ -863,7 +862,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if (inField != -1. && detector->get_dt_min() > (detector->get_TopDrift() - 0.) / vD && field >= FIELD_MIN) { - if (verbosity) + if (verbosity > 0) cerr << "ERROR: dt_min is too restrictive (too large)" << endl; return 1; } @@ -874,7 +873,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, field >= FIELD_MIN) goto Z_NEW; if (detector->get_dt_max() > (detector->get_TopDrift() - 0.) / vD && !j && - field >= FIELD_MIN && verbosity) { + field >= FIELD_MIN && verbosity > 0) { cerr << "WARNING: dt_max is greater than max possible" << endl; } @@ -882,7 +881,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, // code-block dealing with user error (rounding another possible culprit) if (!dEOdxBasis) { if (pos_z <= 0) { - if (verbosity) + if (verbosity > 0) cerr << "WARNING: unphysically low Z coordinate (vertical axis of " "detector) of " << pos_z << " mm" << endl; // warn user on screen @@ -890,7 +889,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, } if ((pos_z > (detector->get_TopDrift() + z_step) || driftTime < 0.0) && field >= FIELD_MIN) { - if (verbosity) + if (verbosity > 0) cerr << "WARNING: unphysically big Z coordinate (vertical axis of " "detector) of " << pos_z << " mm" << endl; // give the specifics @@ -931,7 +930,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, } else { if (keV > 0.001 * Wq_eV) { if (type == "ER") { - if (verbosity && j == 0) { + if (verbosity > 0 && j == 0) { cerr << "CAUTION: Are you sure you don't want beta model instead " "of ER? This is a weighted average of the beta and gamma " "models" @@ -961,7 +960,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, 0.0015, 0.0553, 0.205, 0.45, -0.2}; // zero out non-binom recomb fluct & skew (NR) } - if (!dEOdxBasis) quanta = n.GetQuanta(yields, rho, NRERWidthsParam); + if (!dEOdxBasis) quanta = n.GetQuanta(yields, rho, NRERWidthsParam, false, -999.); } else { yields.PhotonYield = 0.; yields.ElectronYield = 0.; @@ -980,7 +979,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, detector->get_noiseBaseline()[3] != 0.) quanta.electrons += int(floor( RandomGen::rndm()->rand_gauss(detector->get_noiseBaseline()[2], - detector->get_noiseBaseline()[3]) + + detector->get_noiseBaseline()[3],false) + 0.5)); // If we want the smeared positions (non-MC truth), then implement @@ -1034,7 +1033,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if ((useS2 == 0 && logMax <= log10(maxS2 / maxS1)) || (useS2 == 1 && logMax <= log10(maxS2)) || (useS2 == 2 && logMax <= log10(maxS2 / maxS1))) { - if (j == 0 && verbosity) + if (j == 0 && verbosity > 0) cerr << "err: You may be chopping off the upper half of your (ER?) " "band; increase logMax and/or maxS2" << endl; @@ -1042,7 +1041,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if ((useS2 == 0 && logMin >= log10(lowest2 / lowest1)) || (useS2 == 1 && logMin >= log10(lowest2)) || (useS2 == 2 && logMin >= log10(lowest2 / lowest1))) { - if (j == 0 && verbosity) + if (j == 0 && verbosity > 0) cerr << "err: You may be chopping off the lower half of your (NR?) " "band; decrease logMin and/or minS2" << endl; @@ -1057,10 +1056,10 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if (ValidityTests::nearlyEqual(eMin, eMax) || type_num == Kr83m) { if ((scint[3] > maxS1 || scint[5] > maxS1 || scint[7] > maxS1) && - j < 10 && verbosity) + j < 10 && verbosity > 0) cerr << "WARNING: Some S1 pulse areas are greater than maxS1" << endl; if ((scint2[5] > maxS2 || scint2[7] > maxS2) && j < 10 && - verbosity) // don't repeat too much: only if within first 10 events + verbosity > 0) // don't repeat too much: only if within first 10 events // then show (+above) cerr << "WARNING: Some S2 pulse areas are greater than maxS2" << endl; } @@ -1160,7 +1159,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, signal1.pop_back(); signal2.pop_back(); signalE.pop_back(); - if (verbosity) { + if (verbosity > 0) { cerr << endl << "CAUTION: Efficiency seems to have been zero, so trying " "again with full S1 and S2 ranges." @@ -1216,10 +1215,10 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, // scint2[8] = g2; // g2 = ExtEff * SE, light collection efficiency of EL // in gas gap (from CalculateG2) - if (truthPos[2] < detector->get_cathode() && verbosity && !BeenHere && + if (truthPos[2] < detector->get_cathode() && verbosity > 0 && !BeenHere && !dEOdxBasis) { BeenHere = true; - if (verbosity) + if (verbosity > 0) fprintf( stderr, "gamma-X i.e. MSSI may be happening. This may be why even high-E " @@ -1237,9 +1236,9 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, // other suggestions: minS1, minS2 (or s2_thr) for tighter cuts // depending on analysis.hh settings (think of as analysis v. trigger // thresholds) and using max's too, pinching both ends - if (type_num == Kr83m && verbosity && (eMin < 10. || eMin > 40.)) + if (type_num == Kr83m && verbosity > 0 && (eMin < 10. || eMin > 40.)) printf("%.6f\t", yields.DeltaT_Scint); - if (timeStamp > tZero && verbosity) { + if (timeStamp > tZero && verbosity > 0) { printf("%.6f\t", timeStamp); if (keV > AnnModERange[0] && keV < AnnModERange[1]) SaveTheDates[int(timeStamp) % tMax]++; @@ -1281,7 +1280,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, // place-holder // in case you want to drop all sub-threshold data } catch (exception& e) { - if (verbosity) cerr << e.what() << endl; + if (verbosity > 0) cerr << e.what() << endl; return 1; } } // end of the gigantic primary event number loop @@ -1292,7 +1291,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, } } - if (verbosity) { + if (verbosity >= 0) { if (eMin != eMax && type_num != Kr83m) { if (useS2 == 2) GetBand(signal2, signal1, false, detector->get_coinLevel()); @@ -1310,7 +1309,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, std::isnan(band[j][0]) || std::isnan(band[j][1]) || std::isnan(band[j][2]) || std::isnan(band[j][3]) || std::isnan(band[j][4]) || std::isnan(band[j][5])) { - if (eMax != -999. && verbosity) { + if (eMax != -999. && verbosity > 0) { if (((g1 * yieldsMax.PhotonYield) < maxS1 || (g2 * yieldsMax.ElectronYield) < maxS2) && j != 0) @@ -1331,7 +1330,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, if (type_num == NR || type_num == WIMP || type_num == B8 || type_num == DD || type_num == AmBe || type_num == Cf || type_num == ion) { - fprintf(stderr, + if (verbosity > 0) fprintf(stderr, "S1 Mean\t\tS1 Res [%%]\tS2 Mean\t\tS2 Res [%%]\tEc " "[keVnr]\tEc Res[%%]\tEff[%%>thr]\tEc [keVee]\n"); keVee /= numEvts; @@ -1354,17 +1353,19 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, std::isnan(band[j][1]) || std::isnan(band[j][2]) || std::isnan(band[j][3])) && field >= FIELD_MIN) { - if (numEvts > 1) - cerr << "CAUTION: YOUR S1 and/or S2 MIN and/or MAX may be set to " - "be too restrictive, please check.\n"; - else - cerr << "CAUTION: Poor stats. You must have at least 2 events to " - "calculate S1 and S2 and E resolutions.\n"; + if (verbosity > 0) { + if (numEvts > 1) + cerr << "CAUTION: YOUR S1 and/or S2 MIN and/or MAX may be set to " + "be too restrictive, please check.\n"; + else + cerr << "CAUTION: Poor stats. You must have at least 2 events to " + "calculate S1 and S2 and E resolutions.\n"; + } } else if ((ValidityTests::nearlyEqual(energies[0], eMin) || ValidityTests::nearlyEqual(energies[0], eMax) || energies[1] <= 1E-6) && field >= FIELD_MIN) - if (verbosity) + if (verbosity > 0) cerr << "If your energy resolution is 0% then you probably still " "have MC truth energy on." << endl; @@ -1377,7 +1378,7 @@ int execNEST(VDetector* detector, uint64_t numEvts, const string& type, vector> GetBand(vector S1s, vector S2s, bool resol, int nFold) { if (numBins > NUMBINS_MAX) { - if (verbosity) + if (verbosity > 0) cerr << "ERROR: Too many bins. Decrease numBins (analysis.hh) or " "increase NUMBINS_MAX (TestSpectra.hh)" << endl;