From c647c89a0ae1bd9a4d1b1bb5acaa48f058f8f28d Mon Sep 17 00:00:00 2001 From: sophiaandaloro Date: Mon, 1 Mar 2021 12:28:33 -0600 Subject: [PATCH 1/2] Sync with 2.2.1 --- src/nestpy/GammaHandler.cpp | 314 ++++++++++++++++++------------------ src/nestpy/GammaHandler.hh | 94 +++++------ src/nestpy/LUX_Run03.hh | 9 +- src/nestpy/NEST.cpp | 249 ++++++++++++++-------------- src/nestpy/NEST.hh | 20 +-- src/nestpy/RandomGen.cpp | 4 +- src/nestpy/RandomGen.hh | 2 +- src/nestpy/TestSpectra.cpp | 8 +- src/nestpy/TestSpectra.hh | 2 +- src/nestpy/VDetector.hh | 6 + src/nestpy/ValidityTests.hh | 2 +- src/nestpy/__init__.py | 2 +- src/nestpy/bindings.cpp | 17 +- src/nestpy/execNEST.cpp | 164 +++++++++++-------- src/nestpy/execNEST.hh | 2 +- 15 files changed, 475 insertions(+), 420 deletions(-) diff --git a/src/nestpy/GammaHandler.cpp b/src/nestpy/GammaHandler.cpp index f4983c8..fa3e21d 100644 --- a/src/nestpy/GammaHandler.cpp +++ b/src/nestpy/GammaHandler.cpp @@ -1,157 +1,157 @@ -#include -#include "GammaHandler.hh" -#include -#include "TestSpectra.hh" - -using namespace std; - -double yMax = 1.0; //arbitrary y max, might need to change -double brThresh = 0.1; - -const vector>& GammaHandler::sourceLookupTable(const std::string& source) -{ - // energy container vector orginized as {energy, branch ratio, PE mass attenuation coef, Compton coef, Pair Production - // coef} - typedef vector> LookupTable; - static const LookupTable co57Info { - { 122.0, 0.856, 1.793, 0.1081, 0.00 }, - { 136.0, 0.1068, 0.5651, 0.1019, 0.00 }, - { 14.0, 0.0916, 55.5, 0.0744, 0.00 }, - }; - static const LookupTable co60Info { - { 1332.0, 0.9998, 0.001991, 0.04244, 0.0008853 }, - { 1173.0, 0.9985, 0.004126, 0.05156, 0.00 }, - }; - static const LookupTable cs137Info { - { 662.0, 0.851, 0.01338, 0.06559, 0.00 }, - { 284.0, 0.0006, 0.08009, 0.08495, 0.00 }, - }; - if (source == "Co57") - { - return co57Info; - } - else if (source == "Co60") - { - return co60Info; - } - else if (source == "Cs137") - { - return cs137Info; - } - cerr << source << " Is not a valid option!" << endl; - throw std::invalid_argument { source + " is not a valid option!" }; - return co57Info; -} - -const vector GammaHandler::combineSpectra(double emin, double emax, string source) { - double brSum = 0.0; - double fValue = 0.0; - double pe = 0.0; - double compton = 0.0; - double pp = 0.0; - - vector> sourceInfo = GammaHandler::sourceLookupTable(source); - - vector xyTry = { - emin + (emax - emin) * RandomGen::rndm()->rand_uniform(), - yMax * RandomGen::rndm()->rand_uniform(), 1.}; - - while(xyTry[2] > 0.) { - pe = GammaHandler::photoIonization(sourceInfo, xyTry); - compton = GammaHandler::compton(sourceInfo, xyTry); - pp = GammaHandler::pairProduction(sourceInfo, xyTry); - fValue = pe + compton + pp; - xyTry = RandomGen::rndm()->VonNeumann(emin, emax, 0., yMax, xyTry[0], - xyTry[1], fValue); - } - - vector keV_vec = {-1, -1, -1}; - - if( pe > 0.0) { - keV_vec[0] = xyTry[0]; - } - if(compton > 0.0) { - keV_vec[1] = xyTry[0]; - } - if(pp > 0.0) { - keV_vec[2] = xyTry[0]; - } - return keV_vec; -} - -double GammaHandler::photoIonization(const vector>& sourceInfo, const vector& xyTry) -{ - // implement simple delta function to the spectrum - std::size_t index { 0 }; - bool found = false; - for (int i = 0; i < sourceInfo.size(); i++) - { - double initialEnergy = sourceInfo[i][0]; - if (abs(xyTry[0] - initialEnergy) < brThresh) - { - index = i; - found = true; - } - } - double br = sourceInfo[index][1]; - double pe = sourceInfo[index][2]; - double co = sourceInfo[index][3]; - double pp = sourceInfo[index][4]; - if(found) return yMax * br * (pe / (pe + co + pp)); - return 0.0; -} - -double GammaHandler::compton(const vector>& sourceInfo, const vector& xyTry) { - double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV - double thetaMin = 0.0; - double thetaMax = M_PI; - int simpIterations = 100; - double simpStep = (thetaMax - thetaMin)/simpIterations; - double simpCurrentStep = thetaMin; - double shiftedEnergy, simpResult, kn, B, rY, rPsi, initialEnergy; - bool draw = true; - double a = 1.0/137.04; - double re = pow(0.38616, -12); - - //loop over gamma energies - for(int i = 0; i < sourceInfo.size(); i++) { - double initialEnergy = sourceInfo[i][0]; - double br = sourceInfo[i][1]; - double pe = sourceInfo[i][2]; - double co = sourceInfo[i][3]; - double pp = sourceInfo[i][4]; - //get shifted energy with MC - bool draw = true; - while(draw){ - rPsi = M_PI * RandomGen::rndm()->rand_uniform(); - rY = 10* RandomGen::rndm()->rand_uniform(); - - B = 1.0/(1.0+initialEnergy/energyScaleFactor*(1-cos(rPsi))); - kn = M_PI*pow(B,2)*(B+1.0/B-pow(sin(rPsi),2))*sin(rPsi); //Klein-Nishina - if(rY>& sourceInfo, const vector& xyTry) { - double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV - double initialEnergy, shiftedEnergy; - //loop over allowed gamma energies - for(int i = 0; i < sourceInfo.size(); i++) { - double initialEnergy = sourceInfo[i][0]; - double br = sourceInfo[i][1]; - double pe = sourceInfo[i][2]; - double co = sourceInfo[i][3]; - double pp = sourceInfo[i][4]; - shiftedEnergy = (0.5)*(initialEnergy - 2*energyScaleFactor); //Pair production energy - if(abs(xyTry[0]-shiftedEnergy) < brThresh) { - return yMax*br*(pp/(pe+co+pp)); - } - } - return 0.0; -} +#include +#include "GammaHandler.hh" +#include +#include "TestSpectra.hh" + +using namespace std; + +double yMax = 1.0; //arbitrary y max, might need to change +double brThresh = 0.1; + +const vector>& GammaHandler::sourceLookupTable(const std::string& source) +{ + // energy container vector orginized as {energy, branch ratio, PE mass attenuation coef, Compton coef, Pair Production + // coef} + typedef vector> LookupTable; + static const LookupTable co57Info { + { 122.0, 0.856, 1.793, 0.1081, 0.00 }, + { 136.0, 0.1068, 0.5651, 0.1019, 0.00 }, + { 14.0, 0.0916, 55.5, 0.0744, 0.00 }, + }; + static const LookupTable co60Info { + { 1332.0, 0.9998, 0.001991, 0.04244, 0.0008853 }, + { 1173.0, 0.9985, 0.004126, 0.05156, 0.00 }, + }; + static const LookupTable cs137Info { + { 662.0, 0.851, 0.01338, 0.06559, 0.00 }, + { 284.0, 0.0006, 0.08009, 0.08495, 0.00 }, + }; + if (source == "Co57") + { + return co57Info; + } + else if (source == "Co60") + { + return co60Info; + } + else if (source == "Cs137") + { + return cs137Info; + } + cerr << source << " Is not a valid option!" << endl; + throw std::invalid_argument { source + " is not a valid option!" }; + return co57Info; +} + +const vector GammaHandler::combineSpectra(double emin, double emax, string source) { + double brSum = 0.0; + double fValue = 0.0; + double pe = 0.0; + double compton = 0.0; + double pp = 0.0; + + vector> sourceInfo = GammaHandler::sourceLookupTable(source); + + vector xyTry = { + emin + (emax - emin) * RandomGen::rndm()->rand_uniform(), + yMax * RandomGen::rndm()->rand_uniform(), 1.}; + + while(xyTry[2] > 0.) { + pe = GammaHandler::photoIonization(sourceInfo, xyTry); + compton = GammaHandler::compton(sourceInfo, xyTry); + pp = GammaHandler::pairProduction(sourceInfo, xyTry); + fValue = pe + compton + pp; + xyTry = RandomGen::rndm()->VonNeumann(emin, emax, 0., yMax, xyTry[0], + xyTry[1], fValue); + } + + vector keV_vec = {-1, -1, -1}; + + if( pe > 0.0) { + keV_vec[0] = xyTry[0]; + } + if(compton > 0.0) { + keV_vec[1] = xyTry[0]; + } + if(pp > 0.0) { + keV_vec[2] = xyTry[0]; + } + return keV_vec; +} + +double GammaHandler::photoIonization(const vector>& sourceInfo, const vector& xyTry) +{ + // implement simple delta function to the spectrum + std::size_t index { 0 }; + bool found = false; + for (int i = 0; i < sourceInfo.size(); i++) + { + double initialEnergy = sourceInfo[i][0]; + if (abs(xyTry[0] - initialEnergy) < brThresh) + { + index = i; + found = true; + } + } + double br = sourceInfo[index][1]; + double pe = sourceInfo[index][2]; + double co = sourceInfo[index][3]; + double pp = sourceInfo[index][4]; + if(found) return yMax * br * (pe / (pe + co + pp)); + return 0.0; +} + +double GammaHandler::compton(const vector>& sourceInfo, const vector& xyTry) { + double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV + double thetaMin = 0.0; + double thetaMax = M_PI; + int simpIterations = 100; + double simpStep = (thetaMax - thetaMin)/simpIterations; + double simpCurrentStep = thetaMin; + double shiftedEnergy, simpResult, kn, B, rY, rPsi, initialEnergy; + bool draw = true; + double a = 1.0/137.04; + double re = pow(0.38616, -12); + + //loop over gamma energies + for(int i = 0; i < sourceInfo.size(); i++) { + double initialEnergy = sourceInfo[i][0]; + double br = sourceInfo[i][1]; + double pe = sourceInfo[i][2]; + double co = sourceInfo[i][3]; + double pp = sourceInfo[i][4]; + //get shifted energy with MC + bool draw = true; + while(draw){ + rPsi = M_PI * RandomGen::rndm()->rand_uniform(); + rY = 10* RandomGen::rndm()->rand_uniform(); + + B = 1.0/(1.0+initialEnergy/energyScaleFactor*(1-cos(rPsi))); + kn = M_PI*pow(B,2)*(B+1.0/B-pow(sin(rPsi),2))*sin(rPsi); //Klein-Nishina + if(rY>& sourceInfo, const vector& xyTry) { + double energyScaleFactor = ElectronRestMassEnergy; //mc^2 for electron mass in keV + double initialEnergy, shiftedEnergy; + //loop over allowed gamma energies + for(int i = 0; i < sourceInfo.size(); i++) { + double initialEnergy = sourceInfo[i][0]; + double br = sourceInfo[i][1]; + double pe = sourceInfo[i][2]; + double co = sourceInfo[i][3]; + double pp = sourceInfo[i][4]; + shiftedEnergy = (0.5)*(initialEnergy - 2*energyScaleFactor); //Pair production energy + if(abs(xyTry[0]-shiftedEnergy) < brThresh) { + return yMax*br*(pp/(pe+co+pp)); + } + } + return 0.0; +} diff --git a/src/nestpy/GammaHandler.hh b/src/nestpy/GammaHandler.hh index 653ddae..4067a5a 100644 --- a/src/nestpy/GammaHandler.hh +++ b/src/nestpy/GammaHandler.hh @@ -1,47 +1,47 @@ -#ifndef GAMMAHANDLER_HH -#define GAMMAHANDLER_HH - -#include -#include -#include -#include -#include -#include -#include "RandomGen.hh" -#include -#include -#include -//#include "GammaContainer.hh" -using namespace std; - -class GammaHandler { -public: - GammaHandler() {}; - /* - The main function that combines all spectra from photoionization, compton scattering, pair production, etc... - Takes min and max energies, a vector of monoenergetic gamma energies and a vector of their branching ratios. - The index of the gamma energy must correspond with the index of the branch ratio - */ - const vector combineSpectra(double emin, double emax, string source); - - /* - Get y value for given x value in xyTry. Function is just delta functions at the gammaEnergies with amplitudes given - by yMax*branchRatio at that energy - */ - double photoIonization(const vector>& sourceInfo, const vector& xyTry); - - /*Return compton spectrum from KN formula and shifted energies */ - double compton(const vector>& sourceInfo, const vector& xyTry); - - /*Return pair production spectrum from pair production energy equation */ - double pairProduction(const vector>& sourceInfo, const vector& xyTry); - - /*return energies, branching ratios, and mass attenuation coefficients for a given source*/ - - const vector>& sourceLookupTable(const std::string& source); -}; - - - - -#endif +#ifndef GAMMAHANDLER_HH +#define GAMMAHANDLER_HH + +#include +#include +#include +#include +#include +#include +#include "RandomGen.hh" +#include +#include +#include +//#include "GammaContainer.hh" +using namespace std; + +class GammaHandler { +public: + GammaHandler() {}; + /* + The main function that combines all spectra from photoionization, compton scattering, pair production, etc... + Takes min and max energies, a vector of monoenergetic gamma energies and a vector of their branching ratios. + The index of the gamma energy must correspond with the index of the branch ratio + */ + const vector combineSpectra(double emin, double emax, string source); + + /* + Get y value for given x value in xyTry. Function is just delta functions at the gammaEnergies with amplitudes given + by yMax*branchRatio at that energy + */ + double photoIonization(const vector>& sourceInfo, const vector& xyTry); + + /*Return compton spectrum from KN formula and shifted energies */ + double compton(const vector>& sourceInfo, const vector& xyTry); + + /*Return pair production spectrum from pair production energy equation */ + double pairProduction(const vector>& sourceInfo, const vector& xyTry); + + /*return energies, branching ratios, and mass attenuation coefficients for a given source*/ + + const vector>& sourceLookupTable(const std::string& source); +}; + + + + +#endif diff --git a/src/nestpy/LUX_Run03.hh b/src/nestpy/LUX_Run03.hh index 958d18c..1ec2ead 100644 --- a/src/nestpy/LUX_Run03.hh +++ b/src/nestpy/LUX_Run03.hh @@ -4,6 +4,7 @@ #include "VDetector.hh" using namespace std; + //NOTES: best g1 for DD 0.1193, but for tritium 0.1146; S2 noise 1.9, 7.5%; g1_gas 0.1019, 0.1012 //s2fano 3.6, 0.9; eField in gas 6.25, 6.2; e- life 650, 750 us; fid vol 80-130, 38-305 us; gasGap 4.25, 4.5 mm //DISCLAIMER: Slight differences from official published values due to private LUX algorithms @@ -13,6 +14,8 @@ class DetectorExample_LUX_RUN03: public VDetector { public: DetectorExample_LUX_RUN03() { + if ( verbosity ) cerr << "*** Detector definition message ***" << endl; + if ( verbosity ) cerr << "You are currently using the LUX Run03 template detector." << endl << endl; // Call the initialization of all the parameters Initialization(); @@ -27,7 +30,7 @@ public: sPEres = 0.37; //arXiv:1910.04211 sPEthr = (0.3*1.173)/0.915; //arXiv:1910.04211 sPEeff = 1.00; //arXiv:1910.04211 - noiseB[0] =-0.01; //arXiv:1910.04211 + noiseB[0] = 0.00; //arXiv:1910.04211 says -0.01 noiseB[1] = 0.08; //arXiv:1910.04211 noiseB[2] = 0.; noiseB[3] = 0.; @@ -39,10 +42,10 @@ public: extraPhot =false; //default noiseL[0]=1.4e-2; //1910.04211 p.12, to match 1610.02076 Fig. 8 - noiseL[1]=6.0e-2; //1910.04211 p.12, to match 1610.02076 Fig. 8 + noiseL[1]=5.5e-2; //1910.04211 p.12, to match 1610.02076 Fig. 8 // Ionization and Secondary Scintillation (S2) parameters - g1_gas = 0.1016; //0.1 in 1910.04211 + g1_gas = 0.1033; //0.1 in 1910.04211 s2Fano = 2.2; //3.7 in 1910.04211; this matches 1608.05381 better s2_thr = 165.;//(150.*1.173)/0.915; //65-194 pe in 1608.05381 E_gas = 6.23; //6.55 in 1910.04211 diff --git a/src/nestpy/NEST.cpp b/src/nestpy/NEST.cpp index 6b3daf5..66c618b 100644 --- a/src/nestpy/NEST.cpp +++ b/src/nestpy/NEST.cpp @@ -14,7 +14,7 @@ using namespace NEST; const std::vector NESTcalc::default_NuisParam = {11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}; const std::vector NESTcalc::default_FreeParam = {1.,1.,0.1,0.5,0.19,2.25}; -long NESTcalc::BinomFluct(long N0, double prob) { +int64_t NESTcalc::BinomFluct(int64_t N0, double prob) { double mean = N0 * prob; double sigma = sqrt(N0 * prob * (1. - prob)); int N1 = 0; @@ -40,8 +40,8 @@ long NESTcalc::BinomFluct(long N0, double prob) { NESTresult NESTcalc::FullCalculation(INTERACTION_TYPE species, double energy, double density, double dfield, double A, double Z, - const vector& NuisParam /*={11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/, - const vector& FreeParam /*={1.,1.,0.1,0.5,0.19,2.25}*/, + const std::vector& NuisParam /*={11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/, + const std::vector& FreeParam /*={1.,1.,0.1,0.5,0.19,2.25}*/, bool do_times /*=true*/) { if ( density < 1. ) fdetector->set_inGas(true); @@ -104,7 +104,7 @@ double NESTcalc::PhotonTime(INTERACTION_TYPE species, bool exciton, } if (fdetector->get_inGas() || energy < W_DEFAULT * 0.001) { SingTripRatio = 0.1; - tauR = 0.; + if ( fdetector->get_inGas() && !exciton ) tauR = 28e3; else tauR = 0.; //28 microseconds comes from Henrique: https://doi.org/10.1016/j.astropartphys.2018.04.006 if ( ValidityTests::nearlyEqual(ATOM_NUM, 18.)) { tau3 = 1600.; tau1 = 6.; } // from old G4S2Light } if ( tauR < 0. ) tauR = 0.; //in case varied with Gaussian earlier @@ -146,7 +146,7 @@ photonstream NESTcalc::GetPhotonTimes(INTERACTION_TYPE species, return return_photons; } -double NESTcalc::RecombOmegaNR(double elecFrac,const vector& FreeParam/*={1.,1.,0.1,0.5,0.19,2.25}*/) +double NESTcalc::RecombOmegaNR(double elecFrac,const std::vector& FreeParam/*={1.,1.,0.1,0.5,0.19,2.25}*/) { double omega = FreeParam[2]*exp(-0.5*pow(elecFrac-FreeParam[3],2.)/(FreeParam[4]*FreeParam[4])); if ( omega < 0. ) @@ -154,15 +154,15 @@ double NESTcalc::RecombOmegaNR(double elecFrac,const vector& FreeParam/* return omega; } -double NESTcalc::RecombOmegaER(double efield, double elecFrac) +double NESTcalc::RecombOmegaER(double efield, double elecFrac, const std::vector& FreeParam) { - double ampl = 0.14+(0.043-0.14)/(1.+pow(efield/1210.,1.25)); //0.14+(0.05-0.14)/(1.+pow(efield/500.,6.)); //pair with GregR mean yields model + double ampl = 0.14+(0.043-0.14)/(1.+pow(efield/1210.,1.25)); //0.086036+(0.0553-0.086036)/pow(1.+pow(efield/295.2,251.6),0.0069114); //pair with GregR mean yields model if ( ampl < 0. ) ampl = 0.; - double wide = 0.205; - double cntr = 0.5; //0.41 agrees better with Dahl thesis. Odd! Reduces fluctuations for high e-Frac (high EF,low E). Also works with GregR LUX Run04 model + double wide = 0.205; //or: FreeParam #2, like amplitude (#1) + double cntr = 0.5; //0.41-45 agrees better with Dahl thesis. Odd! Reduces fluctuations for high e-Frac (high EF,low E). Also works with GregR LUX Run04 model. FreeParam #3 //for gamma-rays larger than 100 keV at least in XENON10 use 0.43 as the best fit. 0.62-0.37 for LUX Run03 - double skew = -0.2; + double skew = -0.2; //FreeParam #4 double mode = cntr + sqrt(2./M_PI)*skew*wide/sqrt(1.+skew*skew); double norm = 1./(exp(-0.5*pow(mode-cntr,2.)/(wide*wide))*(1.+erf(skew*(mode-cntr)/(wide*sqrt(2.))))); //makes sure omega never exceeds ampl double omega = norm*ampl*exp(-0.5*pow(elecFrac-cntr,2.)/(wide*wide))*(1.+erf(skew*(elecFrac-cntr)/(wide*sqrt(2.)))); @@ -189,15 +189,15 @@ double NESTcalc::FanoER(double density, double Nq_mean,double efield) QuantaResult NESTcalc::GetQuanta(const YieldResult& yields, double density, - const vector& FreeParam/*={1.,1.,0.1,0.5,0.19,2.25}*/) { + const std::vector& FreeParam/*={1.,1.,0.1,0.5,0.19,2.25}*/) { QuantaResult result{}; bool HighE; int Nq_actual, Ne, Nph, Ni, Nex; - + if ( FreeParam.size() < 6 ) { throw std::runtime_error("ERROR: You need a minimum of 6 free parameters for the resolution model."); } - + double excitonRatio = yields.ExcitonRatio; double Nq_mean = yields.PhotonYield + yields.ElectronYield; @@ -211,7 +211,7 @@ QuantaResult NESTcalc::GetQuanta(const YieldResult& yields, double density, } else{ HighE = false; } - + double alf = 1. / (1. + excitonRatio); double recombProb = 1. - (excitonRatio + 1.) * elecFrac; if (recombProb < 0.){ @@ -286,11 +286,11 @@ QuantaResult NESTcalc::GetQuanta(const YieldResult& yields, double density, } //set omega (non-binomial recombination fluctuations parameter) according to whether the Lindhard <1, i.e. this is NR. - double omega = yields.Lindhard <1 ? RecombOmegaNR(elecFrac, FreeParam) : RecombOmegaER(yields.ElectricField, elecFrac); + double omega = yields.Lindhard <1 ? RecombOmegaNR(elecFrac, FreeParam) : RecombOmegaER(yields.ElectricField, elecFrac, FreeParam); if ( ValidityTests::nearlyEqual(ATOM_NUM, 18.) ) omega = 0.0; // Ar has no non-binom sauce double Variance = recombProb * (1. - recombProb) * Ni + omega * omega * Ni * Ni; - + double skewness; if ( (yields.PhotonYield+yields.ElectronYield) > 1e4 || yields.ElectricField > 4e3 || yields.ElectricField < 50. ) { skewness = 0.00; //make it a constant 0 when outside the range of Vetri Velan's Run04 models. @@ -300,14 +300,14 @@ QuantaResult NESTcalc::GetQuanta(const YieldResult& yields, double density, double Wq_eV = wvalue.Wq_eV; double engy = 1e-3 * Wq_eV * (yields.PhotonYield + yields.ElectronYield); double fld = yields.ElectricField; - + double alpha0 = 1.39; double cc0 = 4.0, cc1 = 22.1; double E0 = 7.7, E1 = 54., E2 = 26.7, E3 = 6.4; double F0 = 225., F1 = 71.; - + 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)); @@ -317,14 +317,14 @@ QuantaResult NESTcalc::GetQuanta(const YieldResult& yields, double density, skewness = FreeParam[5]; //2.25 but ~5-20 also good (for NR). All better than zero, but 0 is OK too } //note to self: find way to make 0 for ion (wall BG) incl. alphas? } - + double widthCorrection = sqrt( 1. - (2./M_PI) * skewness*skewness/(1. + skewness*skewness)); double muCorrection = (sqrt(Variance)/widthCorrection)*(skewness/sqrt(1.+skewness*skewness))*sqrt(2./M_PI); if ( fabs(skewness) > DBL_MIN && ValidityTests::nearlyEqual(ATOM_NUM, 54.) ) //skewness model only for Xenon! Ne = int(floor(RandomGen::rndm()->rand_skewGauss((1.-recombProb)*Ni-muCorrection,sqrt(Variance)/widthCorrection,skewness)+0.5)); else Ne = int(floor(RandomGen::rndm()->rand_gauss((1.-recombProb)*Ni,sqrt(Variance))+0.5)); - + if (Ne < 0) Ne = 0; if (Ne > Ni) Ne = Ni; @@ -335,7 +335,7 @@ QuantaResult NESTcalc::GetQuanta(const YieldResult& yields, double density, if ((Nph + Ne) != (Nex + Ni)) { throw std::runtime_error("ERROR: Quanta not conserved. Tell Matthew Immediately!"); } - + if ( fdetector->get_extraPhot() ) { if ( yields.Lindhard != 1. && density >= 3.10 ) { Nph = int(floor(double(Nph)*InfraredNR+0.5)); //IR photons for NR (in SXe) but happens for alphas/ions too @@ -358,7 +358,7 @@ YieldResult NESTcalc::GetYieldGamma(double energy, double density, double dfield double Wq_eV = wvalue.Wq_eV; double alpha = wvalue.alpha; constexpr double m3 = 2., m4 = 2., m6 = 0.; - + const double m1 = 33.951 + (3.3284 - 33.951) / (1. + pow(dfield / 165.34, .72665)); double m2 = 1000 / Wq_eV; @@ -367,7 +367,7 @@ YieldResult NESTcalc::GetYieldGamma(double energy, double density, double dfield double densCorr = 240720. / pow(density, 8.2076); double m7 = 66.825 + (829.25 - 66.825) / (1. + pow(dfield / densCorr, .83344)); - + double Nq = energy * 1000. / Wq_eV; double m8 = 2.; if (fdetector->get_inGas()) m8 = -2.; @@ -386,9 +386,9 @@ YieldResult NESTcalc::GetYieldGamma(double energy, double density, double dfield } YieldResult NESTcalc::GetYieldNROld ( double energy, int option ) { // possible anti-correlation in NR ignored totally - + double Ne, Nph, FakeField; - + if ( option == 0 ) { // with old-school L_eff*S_nr conversion, and more explicit Thomas-Imel box model formulae but approximation where Qy can sail off as energy to 0 FakeField = 1.00; Nph = 0.95*64.*energy * (0.050295*pow(energy,1.3483)*(log(1.+(0.84074*pow(energy,1.3875)))/(0.84074*pow(energy,1.3875)))); // NESTv0.97beta (2011) ~0 V/cm arXiv:1106.1613 @@ -409,7 +409,7 @@ YieldResult NESTcalc::GetYieldNROld ( double energy, int option ) { // possible Ne = ( -3.8780 + 12.372 * pow ( energy, 0.73502 ) ) * exp ( -0.0034329 * energy ); // slide 68 of Matthew's private Google Slides mega-deck Nph = 0.81704 + 3.8584 * pow ( energy, 1.30180 ); // ditto } - + YieldResult result{}; if(Nph<0.)Nph=0.; if(Ne<0.)Ne=0.; result.PhotonYield = Nph; result.ElectronYield = Ne; result.ExcitonRatio = 1.; @@ -417,28 +417,15 @@ YieldResult NESTcalc::GetYieldNROld ( double energy, int option ) { // possible result.ElectricField = FakeField; result.DeltaT_Scint = -999.; return YieldResultValidity ( result, energy, W_DEFAULT ); - + } YieldResult NESTcalc::GetYieldNR(double energy, double density, double dfield, double massNum, - vector NuisParam/*{11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/) + const std::vector &NuisParam/*{11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/) { - - if ( ValidityTests::nearlyEqual(ATOM_NUM, 18.) ) { // liquid Ar - NuisParam[0] = 11.1025; // +/-1.10 Everything from https://docs.google.com/document/d/1vLg8vvY5bcdl4Ah4fzyE182DGWt0Wr7_FJ12_B10ujU - NuisParam[1] = 1.087399; // +/-0.025 - NuisParam[2] = 0.1; // +/-0.005 - NuisParam[3] = -0.0932; // +/-0.0095 - NuisParam[4] = 2.998; // +/-1.026 - NuisParam[5] = 0.3; // Fixed - NuisParam[6] = 2.94; // +/-0.12 - NuisParam[7] = W_DEFAULT / 1000.; - NuisParam[8] = DBL_MAX; - NuisParam[9] = 0.5; // square root - NuisParam[10] = 1.0; - NuisParam[11] = 1.0; - massNum = fdetector->get_molarMass(); - } + + if ( ValidityTests::nearlyEqual ( ATOM_NUM, 18. ) ) massNum = fdetector->get_molarMass(); + if ( NuisParam.size() < 12 ) { throw std::runtime_error("ERROR: You need a minimum of 12 nuisance parameters for the mean yields."); @@ -474,11 +461,11 @@ YieldResult NESTcalc::GetYieldNR(double energy, double density, double dfield, d throw std::runtime_error("ERROR: Quanta not conserved. Tell Matthew Immediately!"); } double NexONi = Nex / Ni; - + Wvalue wvalue = WorkFunction(density,fdetector->get_molarMass()); double Wq_eV = wvalue.Wq_eV; double L = (Nq / energy) * Wq_eV * 1e-3; - + YieldResult result{}; result.PhotonYield = Nph; result.ElectronYield = Ne; @@ -490,7 +477,7 @@ YieldResult NESTcalc::GetYieldNR(double energy, double density, double dfield, d } YieldResult NESTcalc::GetYieldIon(double energy, double density, double dfield, double massNum, double atomNum, - const vector& NuisParam/*{11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/) + const std::vector& NuisParam/*{11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/) { //simply uses the original Lindhard model, but as cited by Hitachi in: //https://indico.cern.ch/event/573069/sessions/230063/attachments/1439101/2214448/Hitachi_XeSAT2017_DM.pdf @@ -536,7 +523,7 @@ YieldResult NESTcalc::GetYieldIon(double energy, double density, double dfield, double Ne = Nq - Nph; if ( ValidityTests::nearlyEqual(A1, 206.) && ValidityTests::nearlyEqual(Z1, 82.) ) Ne = RandomGen::rndm()->rand_gauss(Ne/ChargeLoss,2.*sqrt(Ne/(ChargeLoss*ChargeLoss))); // to compensate for accidentally including Q-loss in fits to Xed data - + if ( ValidityTests::nearlyEqual(Z2, 18.) && ValidityTests::nearlyEqual(Z1, 2.) && ValidityTests::nearlyEqual(A1, 4.) ) { //alphas in argon double factorE = pow((4.71598+pow(dfield,7.72848)),0.109802); double Qy = (1./6200.) * @@ -549,7 +536,7 @@ YieldResult NESTcalc::GetYieldIon(double energy, double density, double dfield, Nph = Ly * energy; L = 0.0; NexONi = 0.21; Wq_eV = 19.5; } - + YieldResult result{}; result.PhotonYield = Nph; result.ElectronYield = Ne; @@ -570,15 +557,21 @@ YieldResult NESTcalc::GetYieldKr83m(double energy, double density, double dfield double alpha = wvalue.alpha; constexpr double deltaT_ns_halflife = 154.4; if (ValidityTests::nearlyEqual(energy, 9.4)) - { + { if (deltaT_ns < 0) { while (deltaT_ns > maxTimeSeparation || deltaT_ns < 0) { deltaT_ns = RandomGen::rndm()->rand_exponential(deltaT_ns_halflife); } } + if (deltaT_ns < 100 && energy < 41.5) { + cerr << "\tWARNING! Past Kr83m model fit validity region. Details: " + << " deltaT_ns is <100 ns and your input energy is either 32.1 or 9.4 keV. " + << " Data for separated Kr83m decays does not yet exist for deltaT_ns <100 ns. " + << " 9.4 & 32.1 keV yields are still summed to physically accurate result, but individually will be nonsensical." << endl; + } Nq = energy * 1e3 / Wq_eV; double medTlevel = 57.462 + (69.201 - 57.462 ) / pow(1. + pow(dfield / 250.13, 0.9), 1.); - double lowTdrop = 35. + (75. - 35.) / pow(1. + pow(dfield/60, 1), 1); + double lowTdrop = 35. + (75. - 35.) / pow(1. + pow(dfield/60, 1), 1); double lowTpeak = 6.2831e4 - (6.2831e4 - 5.949e4 ) / pow(1. + pow(dfield/60, 1.), 1); Nph = energy * (lowTpeak * pow(2. * deltaT_ns + 10., -1.5) + medTlevel) / (1. + pow(deltaT_ns / lowTdrop, -1.*lowTdrop/5.)); Ne = Nq - Nph; @@ -588,7 +581,7 @@ YieldResult NESTcalc::GetYieldKr83m(double energy, double density, double dfield } else { if (ValidityTests::nearlyEqual(energy, 32.1) ) { - Nq = energy * 1e3 / Wq_eV; + Nq = energy * 1e3 / Wq_eV; Nph = energy * (6. + (66.742 - 6.) / pow(1. + pow(dfield / 115.037, 0.6409),0.3215)); Ne = Nq - Nph; if (Ne < 0.) @@ -615,7 +608,7 @@ YieldResult NESTcalc::GetYieldKr83m(double energy, double density, double dfield if ( Ne < 0. ) Ne = 0.; } } - + YieldResult result{}; result.PhotonYield = Nph; result.ElectronYield = Ne; @@ -632,7 +625,7 @@ YieldResult NESTcalc::GetYieldBeta(double energy, double density, double dfield) double Qy, Nq; double Wq_eV = wvalue.Wq_eV; //double alpha = wvalue.alpha; // duplicate definition below. We don't even need this here (it is Nex/Ni) - + if ( ValidityTests::nearlyEqual(ATOM_NUM, 18.) ) { // Liquid Argon double alpha = 32.988 - 552.988/(15.5578+pow(dfield/(-4.7+0.025115*exp(1.3954/0.265360653)), 0.208889)); double beta = 2.01952 + 20.9/pow(1.105 + pow(dfield/0.4, 4.55), 7.502); @@ -662,11 +655,11 @@ YieldResult NESTcalc::GetYieldBeta(double energy, double density, double dfield) Qy = QyLvlmedE + (QyLvllowE - QyLvlmedE) / pow(1. + 1.304 * pow(energy, 2.1393), 0.35535) + QyLvlhighE / (1. + DokeBirks * pow(energy, LET_power)); if (Qy > QyLvllowE && energy > 1. && dfield > 1e4) Qy = QyLvllowE; } - + double Ly = Nq / energy - Qy; double Ne = Qy * energy; double Nph = Ly * energy; - + YieldResult result{}; result.PhotonYield = Nph; result.ElectronYield = Ne; @@ -677,25 +670,25 @@ YieldResult NESTcalc::GetYieldBeta(double energy, double density, double dfield) return YieldResultValidity(result,energy,Wq_eV); // everything needed to calculate fluctuations; } -YieldResult NESTcalc::GetYieldBetaGR ( double energy, double density, double dfield ) { - - if ( RecombOmegaER ( 0.0, 0.5 ) > 0.042 ) +YieldResult NESTcalc::GetYieldBetaGR ( double energy, double density, double dfield, const std::vector& NuisParam ) { + + if ( RecombOmegaER ( 0.0, 0.5, NuisParam ) < 0.05 ) cerr << "WARNING! You need to change RecombOmegaER to go along with GetYieldBetaGR" << endl; - + Wvalue wvalue = WorkFunction(density,fdetector->get_molarMass()); double Wq_eV = wvalue.Wq_eV; double alpha = wvalue.alpha; - + double Nq = energy * 1e3 / Wq_eV; - double m1 = 35.*(1.-1./(1.+(dfield/160.)));//(14.10181492*log10(dfield)-13.1164354516); if ( m1 > 30.66 ) { m1 = 30.66; } + double m1 = 30.66+(6.1978-30.66)/pow(1.+pow(dfield/73.855,2.0318),0.41883); //NuisParam[0]; double m5 = Nq/energy/(1 + alpha*erf(0.05 * energy))-m1; double m10 = (0.0508273937+(0.1166087199-0.0508273937)/(1+pow(dfield/1.39260460e+02,-0.65763592))); - + double Qy = m1 + (77.2931084-m1)/ pow((1.+pow(energy/(log10(dfield)*0.13946236+0.52561312),1.82217496+(2.82528809-1.82217496)/(1+pow(dfield/144.65029656,-2.80532006)))),0.3344049589) + m5 +(0.- m5)/pow((1.+pow(energy/(7.02921301+(98.27936794-7.02921301)/ (1.+pow(dfield/256.48156448,1.29119251))),4.285781736)), m10); - + double coeff_TI = pow(1./DENSITY, 0.3); double coeff_Ni = pow(1./DENSITY, 1.4); double coeff_OL = pow(1./DENSITY, -1.7)/log(1.+coeff_TI*coeff_Ni*pow(DENSITY,1.7)); @@ -704,7 +697,7 @@ YieldResult NESTcalc::GetYieldBetaGR ( double energy, double density, double dfi double Ne = Qy * energy; double Nph = Ly * energy; double lux_NexONi = alpha * erf(0.05 * energy); - + YieldResult result{}; result.PhotonYield = Nph; result.ElectronYield = Ne; @@ -713,11 +706,11 @@ YieldResult NESTcalc::GetYieldBetaGR ( double energy, double density, double dfi result.ElectricField = dfield; result.DeltaT_Scint = -999; return YieldResultValidity(result,energy,Wq_eV); // everything needed to calculate fluctuations; - + } YieldResult NESTcalc::GetYields(INTERACTION_TYPE species, double energy, double density, double dfield, double massNum, - double atomNum, const vector& NuisParam + double atomNum, const std::vector& NuisParam /*={11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/) { switch (species) { case NR: @@ -749,7 +742,7 @@ YieldResult NESTcalc::GetYields(INTERACTION_TYPE species, double energy, double break; default: // beta, CH3T, 14C, the pp solar neutrino background, and Compton/PP spectra of fullGamma return GetYieldBeta(energy,density,dfield); - //return GetYieldBetaGR(energy,density,dfield); + //return GetYieldBetaGR(energy,density,dfield,NuisParam); break; } @@ -787,13 +780,13 @@ NESTcalc::~NESTcalc() { vector NESTcalc::GetS1(const QuantaResult &quanta, double truthPosX, double truthPosY, double truthPosZ, double smearPosX, double smearPosY, double smearPosZ, double driftVelocity, double dV_mid, INTERACTION_TYPE type_num, - long evtNum, double dfield, double energy, + uint64_t evtNum, double dfield, double energy, int useTiming, bool outputTiming, - vector& wf_time, + vector& wf_time, vector& wf_amp) { double truthPos[3] = { truthPosX, truthPosY, truthPosZ }; double smearPos[3] = { smearPosX, smearPosY, smearPosZ }; int Nph = quanta.photons; double subtract[2] = { 0., 0. }; - + wf_time.clear(); wf_amp.clear(); @@ -823,20 +816,20 @@ vector NESTcalc::GetS1(const QuantaResult &quanta, double truthPosX, dou // generate a number of PMT hits drawn from a binomial distribution. // Initialize number of photo-electrons int nHits = BinomFluct ( Nph, g1_XYZ ), Nphe = 0; - + double eff = fdetector->get_sPEeff(); if ( eff < 1. ) eff += ((1.-eff)/(2.*double(fdetector->get_numPMTs())))*double(nHits); if ( eff > 1. ) eff = 1.; if ( eff < 0. ) eff = 0.; - + // Initialize the pulse area and spike count variables double pulseArea = 0., spike = 0., prob; // If single photo-electron efficiency is under 1 and the threshold is above 0 // (some phe will be below threshold) if ( useTiming != -1 ) { // digital nHits eventually becomes spikes (spike++) based upon threshold - + // Follow https://en.wikipedia.org/wiki/Truncated_normal_distribution double TruncGaussAlpha = -1. / fdetector->get_sPEres(); double TruncGaussBeta = DBL_MAX; @@ -846,7 +839,7 @@ vector NESTcalc::GetS1(const QuantaResult &quanta, double truthPosX, dou double BigPhi_Beta = 0.5*(1.+erf(TruncGaussBeta/sqrt(2.))); double CapitalZ = BigPhi_Beta - BigPhi_Alpha; double NewMean = 1. + ( ( LittlePhi_Alpha - LittlePhi_Beta ) / CapitalZ ) * fdetector->get_sPEres(); - + // Step through the pmt hits for (int i = 0; i < nHits; ++i) { // generate photo electron, integer count and area @@ -973,14 +966,14 @@ vector NESTcalc::GetS1(const QuantaResult &quanta, double truthPosX, dou wf_time.push_back((ii - numPts / 2) * SAMPLE_SIZE); wf_amp.push_back(AreaTable[0][ii] + AreaTable[1][ii]); - + if (outputTiming) { char line[80]; if ( AreaTable[0][ii] > PHE_MAX ) subtract[0] = AreaTable[0][ii] - PHE_MAX; else subtract[0] = 0.0; if ( AreaTable[1][ii] > PHE_MAX ) subtract[1] = AreaTable[1][ii] - PHE_MAX; else subtract[1] = 0.0; - sprintf(line, "%lu\t%ld\t%.3f\t%.3f", evtNum, wf_time.back() + (long)tRandOffset, + sprintf(line, "%llu\t%lld\t%.3f\t%.3f", evtNum, wf_time.back() + (int64_t)tRandOffset, AreaTable[0][ii]-subtract[0], AreaTable[1][ii]-subtract[1]); pulseFile << line << flush; } @@ -1001,11 +994,17 @@ vector NESTcalc::GetS1(const QuantaResult &quanta, double truthPosX, dou // TimeTable[0][ii]+24., TimeTable[1][ii] ); pulseFile << line << endl; } } - + if ( fdetector->get_noiseL()[0] >= 0.1 ) cerr << " !!WARNING!! S1 linear noise term is greater than or equal to 10% (i.e. 0.1) Did you mistake fraction for percent??" << endl; - pulseArea = RandomGen::rndm()->rand_gauss( + + if(fdetector->get_noiseQ()[0] != 0) { + pulseArea = RandomGen::rndm()->rand_gauss( + pulseArea, sqrt(pow(fdetector->get_noiseQ()[0] * pow(pulseArea, 2), 2) + pow(fdetector->get_noiseL()[0] * pulseArea, 2))); + }else { + pulseArea = RandomGen::rndm()->rand_gauss( pulseArea, fdetector->get_noiseL()[0] * pulseArea); + } pulseArea -= ( subtract[0] + subtract[1] ); if (pulseArea < fdetector->get_sPEthr()) pulseArea = 0.; if (spike < 0) spike = 0; @@ -1101,9 +1100,9 @@ vector NESTcalc::GetS1(const QuantaResult &quanta, double truthPosX, dou } vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, double truthPosZ, double smearPosX, double smearPosY, double smearPosZ, - double dt, double driftVelocity, long evtNum, + double dt, double driftVelocity, uint64_t evtNum, double dfield, int useTiming, bool outputTiming, - vector& wf_time, + vector& wf_time, vector& wf_amp, const vector& g2_params) { double truthPos[3] = { truthPosX, truthPosY, truthPosZ }; double smearPos[3] = { smearPosX, smearPosY, smearPosZ }; @@ -1136,18 +1135,18 @@ vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, doubl posDep /= fdetector->FitS2(0., 0., VDetector::fold); posDepSm /= fdetector->FitS2(0., 0., VDetector::unfold); double dz = fdetector->get_TopDrift() - dt * driftVelocity; - + if ( ValidityTests::nearlyEqual(fdetector->get_g1_gas(), 1. )) { posDep = 1.; posDepSm = 1.; } if ( ValidityTests::nearlyEqual(fdetector->get_g1_gas(), 0. )) { posDep = 0.; posDepSm = 0.; } - + int Nee = BinomFluct(Ne, ExtEff * exp(-dt / fdetector->get_eLife_us())); //MAKE this 1 for SINGLE e- DEBUG - - long Nph = 0, nHits = 0, Nphe = 0; + + uint64_t Nph = 0, nHits = 0, Nphe = 0; double pulseArea = 0.; - + if ( useTiming >= 1 ) { - long k; + uint64_t k; int stopPoint; double tau1, tau2, E_liq, amp2; vector electronstream, AreaTableBot[2], AreaTableTop[2], TimeTable; @@ -1209,9 +1208,9 @@ vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, doubl SE = floor(RandomGen::rndm()->rand_gauss( elYield, sqrt(fdetector->get_s2Fano() * elYield)) + 0.5); - Nph += long(SE); - SE = (double)BinomFluct(long(SE), fdetector->get_g1_gas() * posDep); - nHits += long(SE); + Nph += uint64_t(SE); + SE = (double)BinomFluct(uint64_t(SE), fdetector->get_g1_gas() * posDep); + nHits += uint64_t(SE); double KE = 0.5 * 9.109e-31 * driftVelocity_gas * driftVelocity_gas * 1e6 / 1.602e-16; double origin = fdetector->get_TopDrift() + gasGap / 2.; @@ -1220,12 +1219,12 @@ vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, doubl quanta.electrons = 0; quanta.ions = 0; quanta.excitons = int(floor(0.0566 * SE + 0.5)); - photonstream photon_emission_times = GetPhotonTimes(NEST::beta, quanta.photons, + photonstream photon_emission_times = GetPhotonTimes(NEST::beta, quanta.photons, quanta.excitons, dfield, KE); photonstream photon_times = AddPhotonTransportTime(photon_emission_times, newX, newY, origin); - SE += (double)BinomFluct(long(SE), fdetector->get_P_dphe()); - Nphe += long(SE); + SE += (double)BinomFluct(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); phi = 2. * M_PI * RandomGen::rndm()->rand_uniform(); @@ -1254,7 +1253,7 @@ vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, doubl if (i < Nee || !eTrain) pulseArea += phe; origin = fdetector->get_TopDrift() + gasGap * RandomGen::rndm()->rand_uniform(); - k = long(j); + k = uint64_t(j); if (k >= photon_times.size()) k -= photon_times.size(); double offset = ((fdetector->get_anode() - origin) / driftVelocity_gas + electronstream[i]) * @@ -1306,7 +1305,7 @@ vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, doubl } for (k = 0; k < numPts; ++k) { if ((AreaTableBot[1][k] + AreaTableTop[1][k]) <= PULSEHEIGHT) continue; - wf_time.push_back(k * SAMPLE_SIZE + long(min + SAMPLE_SIZE / 2.)); + wf_time.push_back(k * SAMPLE_SIZE + int64_t(min + SAMPLE_SIZE / 2.)); wf_amp.push_back(AreaTableBot[1][k] + AreaTableTop[1][k]); if (outputTiming) { @@ -1315,13 +1314,13 @@ vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, doubl else subtract[0] = 0.0; if ( AreaTableTop[1][k] > PHE_MAX ) subtract[1] = AreaTableTop[1][k] - PHE_MAX; else subtract[1] = 0.0; - sprintf(line, "%lu\t%ld\t%.3f\t%.3f", evtNum, wf_time.back(), + sprintf(line, "%llu\t%lld\t%.3f\t%.3f", evtNum, wf_time.back(), AreaTableBot[1][k]-subtract[0], AreaTableTop[1][k]-subtract[1]); pulseFile << line << endl; } } } else { - Nph = long(floor(RandomGen::rndm()->rand_gauss( + Nph = uint64_t(floor(RandomGen::rndm()->rand_gauss( elYield * double(Nee), sqrt(fdetector->get_s2Fano() * elYield * double(Nee))) + 0.5)); @@ -1330,11 +1329,16 @@ vector NESTcalc::GetS2(int Ne, double truthPosX, double truthPosY, doubl pulseArea = RandomGen::rndm()->rand_gauss( Nphe, fdetector->get_sPEres() * sqrt(Nphe)); } - + if ( fdetector->get_noiseL()[1] >= 0.1 ) cerr << " !!WARNING!! S2 linear noise term is greater than or equal to 10% (i.e. 0.1) Did you mistake fraction for percent??" << endl; - pulseArea = RandomGen::rndm()->rand_gauss( + if(fdetector->get_noiseQ()[1] != 0) { + pulseArea = RandomGen::rndm()->rand_gauss( + pulseArea, sqrt(pow(fdetector->get_noiseQ()[1] * pow(pulseArea, 2), 2) + pow(fdetector->get_noiseL()[1] * pulseArea, 2))); + }else { + pulseArea = RandomGen::rndm()->rand_gauss( pulseArea, fdetector->get_noiseL()[1] * pulseArea); + } pulseArea -= (subtract[0]+subtract[1]); double pulseAreaC = pulseArea / exp(-dt / fdetector->get_eLife_us()) / posDepSm; @@ -1424,7 +1428,7 @@ vector NESTcalc::CalculateG2(bool verbosity) { } if (ExtEff > 1. || fdetector->get_inGas()) ExtEff = 1.; if (ExtEff < 0. || E_liq <= 0.) ExtEff = 0.; - + double gasGap = fdetector->get_anode() - fdetector @@ -1472,7 +1476,13 @@ vector NESTcalc::CalculateG2(bool verbosity) { nHits = BinomFluct ( Nph, fdetector->get_g1_gas() * posDep ); Nphe = nHits+BinomFluct(nHits,fdetector->get_P_dphe()); pulseArea = RandomGen::rndm()->rand_gauss(Nphe,fdetector->get_sPEres()*sqrt(Nphe)); - pulseArea = RandomGen::rndm()->rand_gauss(pulseArea,fdetector->get_noiseL()[1]*pulseArea); + if(fdetector->get_noiseQ()[1] != 0) { + pulseArea = RandomGen::rndm()->rand_gauss( + pulseArea, sqrt(pow(fdetector->get_noiseQ()[1] * pow(pulseArea, 2), 2) + pow(fdetector->get_noiseL()[1] * pulseArea, 2))); + }else { + pulseArea = RandomGen::rndm()->rand_gauss( + pulseArea, fdetector->get_noiseL()[1] * pulseArea); + } if ( fdetector->get_s2_thr() < 0. ) pulseArea = RandomGen::rndm()->rand_gauss(fdetector->FitTBA(0.0,0.0,fdetector->get_TopDrift()/2.)[1]*pulseArea,sqrt (fdetector->FitTBA(0.0,0.0,fdetector->get_TopDrift()/2.)[1]* @@ -1481,14 +1491,14 @@ vector NESTcalc::CalculateG2(bool verbosity) { NphdC = pulseAreaC/(1.+fdetector->get_P_dphe()); StdDev += (SE-NphdC)*(SE-NphdC); } StdDev = sqrt(StdDev)/sqrt(9999.); // N-1 from above (10,000) - + cout << endl << "g1 = " << fdetector->get_g1() << " phd per photon\tg2 = " << g2 << " phd per electron (e-EE = "; cout << ExtEff * 100. << "%, SE_mean,width = " << SE << "," << StdDev << ")\t"; } - + // Store the g2 parameters in a vector for later (calculated once per // detector) g2_params[0] = elYield; @@ -1547,15 +1557,15 @@ double NESTcalc::GetDensity(double Kelvin, else return DENSITY; } - + //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 && ValidityTests::nearlyEqual(ATOM_NUM, 54. )) { // solid Xenon cerr << "\nWARNING: SOLID PHASE. IS THAT WHAT YOU WANTED?\n"; return 3.41; // from Yoo at 157K // other sources say 3.1 (Wikipedia, 'minimum') and 3.640g/mL at an unknown temperature } - + double VaporP_bar; // we will calculate using NIST if (Kelvin < 289.7) { if ( ValidityTests::nearlyEqual(ATOM_NUM, 54. ) ) VaporP_bar = pow(10., 4.0519 - 667.16 / Kelvin); @@ -1602,7 +1612,7 @@ double NESTcalc::GetDriftVelocity_Liquid(double Kelvin, double Density, 0.0; // returns drift speed in mm/usec. based on Fig. 14 arXiv:1712.08607 int i, j; double vi, vf, slope, Ti, Tf, offset; - + if ( ValidityTests::nearlyEqual(ATOM_NUM, 18.) ) { //replace eventually with Kate's work, once T's splined as done for LXe and SXe /*x_nest is field in kV/cm: Liquid: @@ -1619,7 +1629,7 @@ double NESTcalc::GetDriftVelocity_Liquid(double Kelvin, double Density, speed = 0.097384*pow(log10(eField),3.0622)-0.018614*sqrt(eField); if ( speed < 0. ) speed = 0.; return speed; } - + double polyExp[11][7] = { {-3.1046, 27.037, -2.1668, 193.27, -4.8024, 646.04, 9.2471}, // 100K {-2.7394, 22.760, -1.7775, 222.72, -5.0836, 724.98, 8.7189}, // 120 @@ -1637,14 +1647,14 @@ double NESTcalc::GetDriftVelocity_Liquid(double Kelvin, double Density, double Temperatures[11] = {100., 120., 140., 155., 157., 163., 165., 167., 184., 200., 230.}; - + if ( Kelvin < 100. || Kelvin > 230. ) { cerr << "\nWARNING: TEMPERATURE OUT OF RANGE (100-230 K) for vD\n"; if ( Kelvin < 100. ) Kelvin = 100.; if ( Kelvin > 230. ) Kelvin = 230.; cerr << "Using value at closest temp for a drift speed estimate\n"; } - + if (Kelvin >= Temperatures[0] && Kelvin < Temperatures[1]) i = 0; else if (Kelvin >= Temperatures[1] && Kelvin < Temperatures[2]) @@ -1694,7 +1704,7 @@ double NESTcalc::GetDriftVelocity_Liquid(double Kelvin, double Density, slope = (vf - vi) / (Tf - Ti); speed = slope * (Kelvin - Ti) + vi; } - + if ( speed <= 0. ) { cerr << "\nWARNING: DRIFT SPEED NON-POSITIVE. Setting to 0.1 mm/us\t" << "Line Number ~1700 of NEST.cpp, in function NESTcalc::GetDriftVelocity_Liquid\t" << @@ -1792,7 +1802,7 @@ vector NESTcalc::xyResolution(double xPos_mm, double yPos_mm, sigmaR = fabs(RandomGen::rndm()->rand_gauss(0.0, sigmaR)); double sigmaX = sigmaR * cos(phi); double sigmaY = sigmaR * sin(phi); - + if ( sigmaR > 1e2 || std::isnan(sigmaR) || sigmaR <= 0. || fabs(sigmaX) > 1e2 || fabs(sigmaY) > 1e2 ) { if ( A_top > 20. ) { @@ -1800,7 +1810,7 @@ vector NESTcalc::xyResolution(double xPos_mm, double yPos_mm, cerr << "Setting resolution to perfect." << endl; sigmaX=0.; sigmaY=0.; } // this is only a problem if the area is large and the res is still bad } - + xySmeared[0] = xPos_mm + sigmaX; xySmeared[1] = yPos_mm + sigmaY; @@ -1810,12 +1820,12 @@ vector NESTcalc::xyResolution(double xPos_mm, double yPos_mm, double NESTcalc::PhotonEnergy(bool s2Flag, bool state, double tempK) { double wavelength, E_keV; // wavelength is in nanometers - + 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) // liquid or solid wavelength = RandomGen::rndm()->rand_gauss( 178., 14. / 2.355); // taken from Jortner JchPh 42 '65 @@ -1842,7 +1852,7 @@ double NESTcalc::PhotonEnergy(bool s2Flag, bool state, double tempK) { double NESTcalc::CalcElectronLET ( double E, int Z ) { double LET; - + // use a spline fit to online ESTAR data if ( ValidityTests::nearlyEqual(ATOM_NUM, 54.)) { if ( E >= 1. ) @@ -1868,7 +1878,7 @@ double NESTcalc::CalcElectronLET ( double E, int Z ) { else LET = 0.; } - + return LET; } @@ -1886,8 +1896,7 @@ NESTcalc::Wvalue NESTcalc::WorkFunction(double density, double MolarMass) { +xi_se/(1.+alpha);*/ double eDensity = ( density / MolarMass ) * NEST_AVO * ATOM_NUM; Wq_eV = 20.7 - 1.01e-23 * eDensity; - - return Wvalue{.Wq_eV=Wq_eV,.alpha=alpha}; //W and Nex/Ni together + return Wvalue{.Wq_eV=Wq_eV,.alpha=alpha}; //W and Nex/Ni together } double NESTcalc::NexONi(double energy, double density) diff --git a/src/nestpy/NEST.hh b/src/nestpy/NEST.hh index e50d8bf..2384812 100644 --- a/src/nestpy/NEST.hh +++ b/src/nestpy/NEST.hh @@ -170,7 +170,7 @@ class NESTcalc { explicit NESTcalc(VDetector* detector); virtual ~NESTcalc(); - static long BinomFluct(long, double); + static int64_t BinomFluct(int64_t, double); static const std::vector default_NuisParam; /* = {11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}*/ static const std::vector default_FreeParam; /* = {1.,1.,0.1,0.5,0.19,2.25} */ @@ -206,17 +206,17 @@ class NESTcalc { virtual YieldResult GetYieldGamma(double energy, double density, double dfield); // Called by GetYields in the Gamma/x-ray/Photoabsorption Case virtual YieldResult GetYieldNR(double energy, double density, double dfield, double massNum, - std::vector NuisParam={11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}); + const std::vector& NuisParam={11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}); // Called by GetYields in the NR (and related) cases virtual YieldResult GetYieldNROld ( double energy, int alt ); // Quick and dirty simple analytical approximations saved for earlier NEST versions that were first principles: power laws, ln, sigmoid, exponentials - virtual YieldResult GetYieldIon(double energy, double density, double dfield, double massNum, double atomNum, const vector& NuisParam={11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}); + virtual YieldResult GetYieldIon(double energy, double density, double dfield, double massNum, double atomNum, const std::vector& NuisParam={11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1.,1.}); // Called by GetYields in the ion case virtual YieldResult GetYieldKr83m(double energy, double density, double dfield, double maxTimeSeparation, double deltaT_ns); // Called by GetYields in the K383m case virtual YieldResult GetYieldBeta(double energy, double density, double dfield); // Called by GetYields in the Beta/Compton/etc.(IC,Auger,EC) Case - virtual YieldResult GetYieldBetaGR(double energy, double density, double dfield); + virtual YieldResult GetYieldBetaGR(double energy, double density, double dfield, const std::vector& NuisParam); // Greg R. version: arXiv:1910.04211 virtual YieldResult YieldResultValidity(YieldResult& res, const double energy, const double Wq_eV); // Confirms and sometimes adjusts YieldResult to make physical sense @@ -225,18 +225,18 @@ class NESTcalc { // quanta (photons+electrons) with a Fano-like factor, and the "slosh" between // photons and electrons // Namely, the recombination fluctuations - virtual double RecombOmegaNR(double elecFrac,const vector& FreeParam/*={1.,1.,0.1,0.5,0.19,2.25}*/); + virtual double RecombOmegaNR(double elecFrac,const std::vector& FreeParam/*={1.,1.,0.1,0.5,0.19,2.25}*/); //Calculates the Omega parameter governing non-binomial recombination fluctuations for nuclear recoils and ions (Lindhard<1) - virtual double RecombOmegaER(double efield, double elecFrac); + virtual double RecombOmegaER(double efield, double elecFrac, const std::vector& FreeParam); //Calculates the Omega parameter governing non-binomial recombination fluctuations for gammas and betas (Lindhard==1) virtual double FanoER(double density, double Nq_mean,double efield); //Fano-factor (and Fano-like additional energy resolution model) for gammas and betas (Lindhard==1) std::vector GetS1(const QuantaResult& quanta, double truthPosX, double truthPosY, double truthPosZ, double smearPosX, double smearPosY, double smearPosZ, double driftSpeed, double dS_mid, INTERACTION_TYPE species, - long evtNum, double dfield, double energy, + uint64_t evtNum, double dfield, double energy, int useTiming, bool outputTiming, - vector& wf_time, vector& wf_amp); + 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 @@ -246,9 +246,9 @@ class NESTcalc { // GetSpike takes the extremely basic digital/integer number of spike counts // provided by GetS1 and does more realistic smearing std::vector GetS2(int Ne, double truthPosX, double truthPosY, double truthPosZ, double smearPosX, double smearPosY, double smearPosZ, - double dt, double driftSpeed, long evtNum, + double dt, double driftSpeed, uint64_t evtNum, double dfield, int useTiming, bool outputTiming, - vector& wf_time, vector& wf_amp, + 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 diff --git a/src/nestpy/RandomGen.cpp b/src/nestpy/RandomGen.cpp index 676d81e..98c5f43 100644 --- a/src/nestpy/RandomGen.cpp +++ b/src/nestpy/RandomGen.cpp @@ -19,8 +19,8 @@ std::uint64_t splitmix64(std::uint64_t z) { return z ^ (z >> 31); } -void RandomGen::SetSeed(unsigned long int s) { - unsigned long int s1 = splitmix64(s); +void RandomGen::SetSeed(uint64_t s) { + uint64_t s1 = splitmix64(s); rng = xoroshiro128plus64(s1, splitmix64(s1)); } diff --git a/src/nestpy/RandomGen.hh b/src/nestpy/RandomGen.hh index c7834de..7e8190a 100644 --- a/src/nestpy/RandomGen.hh +++ b/src/nestpy/RandomGen.hh @@ -16,7 +16,7 @@ using namespace std; class RandomGen { public: static RandomGen* rndm(); - void SetSeed(unsigned long int s); + void SetSeed(uint64_t s); double rand_uniform(); double rand_gauss(double mean, double sigma); double rand_exponential(double half_life); diff --git a/src/nestpy/TestSpectra.cpp b/src/nestpy/TestSpectra.cpp index eca27e9..81fd84e 100644 --- a/src/nestpy/TestSpectra.cpp +++ b/src/nestpy/TestSpectra.cpp @@ -186,7 +186,7 @@ double TestSpectra::DD_spectrum( } double TestSpectra::ppSolar_spectrum ( double xMin, double xMax ) { - + if ( xMax > 250. ) xMax = 250.; if ( xMin < 0.00 ) xMin = 0.00; double yMax = 0.000594; @@ -202,7 +202,7 @@ double TestSpectra::ppSolar_spectrum ( double xMin, double xMax ) { } double TestSpectra::atmNu_spectrum ( double xMin, double xMax ) { - + if ( xMax > 85. ) xMax = 85.; if ( xMin < 0.0 ) xMin = 0.0; vector xyTry = { @@ -240,7 +240,7 @@ double TestSpectra::WIMP_dRate(double ER, double mWimp, double dayNum) { double v_esc = V_ESCAPE * cmPerkm; // escape velocity double v_e = ( V_SUN + ( 0.49 * 29.8 * cos ( ( dayNum * 2. * M_PI / 365.24 ) - ( 0.415 * 2. * M_PI ) ) ) ) * cmPerkm; // the Earth's velocity // used Eq. 18 for SHM w/ June 1 as reference date (MAX!) from arXiv 0607121 [Savage, Freese, Gondolo 2006] - Juergen Reichenbacher 09/17/2020 - + // Define the detector Z and A and the mass of the target nucleus double Z = ATOM_NUM; double A = (double)RandomGen::rndm()->SelectRanXeAtom(); @@ -378,7 +378,7 @@ TestSpectra::WIMP_spectrum_prep TestSpectra::WIMP_prep_spectrum(double mass, dou if ( nZeros == 100 ) break; //quit the for-loop once we're sure we're only getting zeros } - for (long i = 0; i < 1000000; ++i) { + for (uint64_t i = 0; i < 1000000; ++i) { spectrum.integral += WIMP_dRate(double(i) / 1e4, mass, dayNum) / 1e4; } spectrum.xMax = ( (double) EnergySpec.size() - 1. )/divisor; diff --git a/src/nestpy/TestSpectra.hh b/src/nestpy/TestSpectra.hh index 5166758..abcd51f 100644 --- a/src/nestpy/TestSpectra.hh +++ b/src/nestpy/TestSpectra.hh @@ -60,7 +60,7 @@ class TestSpectra { static double WIMP_spectrum(WIMP_spectrum_prep wprep, double mass, double day); static const vector Gamma_spectrum(double xMin, double xMax, string source); - double ZeplinBackground(); // an example of how to do a better (non-flat) ER + static double ZeplinBackground(); // an example of how to do a better (non-flat) ER // BG spectrum for a WS, from Henrique Araujo }; diff --git a/src/nestpy/VDetector.hh b/src/nestpy/VDetector.hh index e980740..9cbd216 100644 --- a/src/nestpy/VDetector.hh +++ b/src/nestpy/VDetector.hh @@ -27,6 +27,7 @@ class VDetector { double get_sPEeff() const { return sPEeff; } const double* get_noiseB() { return &noiseB[0]; } const double* get_noiseL() { return &noiseL[0]; } + const double* get_noiseQ() { return &noiseQ[0]; } double get_P_dphe() const { return P_dphe; } bool get_extraPhot() const{ return extraPhot; } @@ -80,6 +81,10 @@ class VDetector { noiseL[0] = p1; noiseL[1] = p2; } + void set_noiseQ(double p1, double p2) { + noiseQ[0] = p1; + noiseQ[1] = p2; + } void set_P_dphe(double param) { P_dphe = param; } void set_extraPhot(bool param){ extraPhot = param;} @@ -171,6 +176,7 @@ protected: bool extraPhot=false; // for matching EXO-200's W measurement //"Linear noise" terms as defined in Dahl thesis and by D. McK double noiseL[2] = {3e-2,3e-2}; // S1->S1 Gaussian-smeared w/ noiseL[0]*S1. Ditto S2 + double noiseQ[2] = {0.0, 0.0}; //(n)EXO quadratic noise term // Ionization and Secondary Scintillation (S2) parameters double g1_gas = 0.06; // phd per S2 photon in gas, used to get SE size diff --git a/src/nestpy/ValidityTests.hh b/src/nestpy/ValidityTests.hh index cf62f98..27faded 100644 --- a/src/nestpy/ValidityTests.hh +++ b/src/nestpy/ValidityTests.hh @@ -2,7 +2,7 @@ #define VALIDITYTESTS_HH 1 #include -#include +#include #include using namespace std; diff --git a/src/nestpy/__init__.py b/src/nestpy/__init__.py index 79f73da..5d4e4fc 100644 --- a/src/nestpy/__init__.py +++ b/src/nestpy/__init__.py @@ -1,5 +1,5 @@ __version__ = '1.4.4' -__nest_version__ = '2.2.0' +__nest_version__ = '2.2.1' from .nestpy import * diff --git a/src/nestpy/bindings.cpp b/src/nestpy/bindings.cpp index 6d70850..433e775 100644 --- a/src/nestpy/bindings.cpp +++ b/src/nestpy/bindings.cpp @@ -3,7 +3,6 @@ #include "VDetector.hh" #include "execNEST.hh" #include "DetectorExample_XENON10.hh" -#include "LUX_Run03.hh" #include "RandomGen.hh" #include #include @@ -152,14 +151,14 @@ PYBIND11_MODULE(nestpy, m) { .def("FitTBA", &DetectorExample_XENON10::FitTBA) .def("OptTrans", &DetectorExample_XENON10::OptTrans) .def("SinglePEWaveForm", &DetectorExample_XENON10::SinglePEWaveForm); - - // Binding for example LUX_Run03 - py::class_>(m, "LUX_RUN3") - .def(py::init<>()) - .def("Initialization", &DetectorExample_LUX_RUN03::Initialization) - .def("FitTBA", &DetectorExample_LUX_RUN03::FitTBA) - .def("OptTrans", &DetectorExample_LUX_RUN03::OptTrans) - .def("SinglePEWaveForm", &DetectorExample_LUX_RUN03::SinglePEWaveForm); +// +// // Binding for example LUX_Run03 +// py::class_>(m, "LUX_RUN3") +// .def(py::init<>()) +// .def("Initialization", &DetectorExample_LUX_RUN03::Initialization) +// .def("FitTBA", &DetectorExample_LUX_RUN03::FitTBA) +// .def("OptTrans", &DetectorExample_LUX_RUN03::OptTrans) +// .def("SinglePEWaveForm", &DetectorExample_LUX_RUN03::SinglePEWaveForm); // Binding for the NESTcalc class py::class_>(m, "NESTcalc") diff --git a/src/nestpy/execNEST.cpp b/src/nestpy/execNEST.cpp index c19b260..92bc9c4 100644 --- a/src/nestpy/execNEST.cpp +++ b/src/nestpy/execNEST.cpp @@ -36,17 +36,17 @@ int main(int argc, char** argv) { auto* detector = new DetectorExample_LUX_RUN03(); // Custom parameter modification functions // detector->ExampleFunction(); - + if ( ValidityTests::nearlyEqual(ATOM_NUM, 18.) ) { detector->set_molarMass(39.948); 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" << endl; } - + /* vector eList = { 1., 2., 3. }; // fast example--for PLR, ML train vector> pos3dxyz = { {0.,-1.,60.},{-1.,0.,70.},{1.,0.,80.} }; runNESTvec ( detector, NEST::beta, eList, pos3dxyz ); return EXIT_SUCCESS; */ - + if (argc < 7) { cout << "This program takes 6 (or 7) inputs, with Z position in mm from " "bottom of detector:" << endl; @@ -55,7 +55,7 @@ int main(int argc, char** argv) { << endl; cout << "For Kr83m time-dependent 9.4, 32.1, or 41.5 keV yields: " << endl; cout << "\t ./execNEST numEvents Kr83m Energy[keV] maxTimeDiff[ns] " - "field_drift[V/cm] x,y,z-position[mm] {optional:seed}" << endl + "field_drift[V/cm] x,y,z-position[mm] {optional:seed}" << endl << endl; cout << "For 8B or pp or atmNu, numEvts is kg-days of exposure with everything else same. " "For WIMPs:" << endl; @@ -70,13 +70,13 @@ int main(int argc, char** argv) { << endl; return 1; } - - unsigned long int numEvts; + + uint64_t numEvts; string type, position, posiMuon; double eMin, eMax, inField, fPos; int seed; bool no_seed = false; if ( loopNEST ) { - + numEvts = 100000; //10,000 faster but of course less precise if ( loopNEST == 1 ) type = "ER"; @@ -93,9 +93,14 @@ int main(int argc, char** argv) { FreeParam.clear(); NuisParam.clear(); verbosity = false; useTiming = 0; //-1 for faster (less accurate) - + if ( type == "ER" ) { - + + /*detector->set_g1(atof(argv[1])); //an alternate loop approach + detector->set_g1_gas(atof(argv[2])); + inField = atof(argv[3]); + detector->set_noiseL(atof(argv[4]),atof(argv[5]));*/ + FreeParam.push_back(atof(argv[1])); //-0.1 for LUX C-14 ~200V/cm FreeParam.push_back(atof(argv[2])); //0.5 FreeParam.push_back(atof(argv[3])); //0.06 @@ -104,7 +109,7 @@ int main(int argc, char** argv) { FreeParam.push_back(atof(argv[6])); //0.95 FreeParam.push_back(atof(argv[7])); //8e-2 FreeParam.push_back(atof(argv[7])); //repeat - + NuisParam.push_back(11.); NuisParam.push_back(1.1); NuisParam.push_back(0.0480); @@ -117,11 +122,11 @@ int main(int argc, char** argv) { NuisParam.push_back(0.5); NuisParam.push_back(1.0); NuisParam.push_back(1.0); - + } - + else { - + NuisParam.push_back(atof(argv[1])); //11.0 XENON10 NuisParam.push_back(atof(argv[2])); //1.09 NuisParam.push_back(0.0480); @@ -143,15 +148,15 @@ int main(int argc, char** argv) { FreeParam.push_back(0.50); FreeParam.push_back(0.19); FreeParam.push_back(2.25); - + } - + } else { - - numEvts = (unsigned long)atof(argv[1]); + + numEvts = (uint64_t)atof(argv[1]); if ( numEvts <= 0 ) { - cerr << "ERROR, you must simulate at least 1 event" << endl; + cerr << "ERROR, you must simulate at least 1 event, or 1 kg*day" << endl; return 1; } type = argv[2]; @@ -161,7 +166,7 @@ int main(int argc, char** argv) { position = argv[6]; posiMuon = argv[4]; fPos = atof(argv[6]); - + seed = 0; //if not given make 0 if (argc == 8) { seed = atoi(argv[7]); @@ -170,7 +175,7 @@ int main(int argc, char** argv) { RandomGen::rndm()->SetSeed(0); no_seed = true; } - + FreeParam.clear(); NuisParam.clear(); if ( type == "ER" ) { @@ -191,42 +196,57 @@ int main(int argc, char** argv) { FreeParam.push_back(0.19); // width parameter (Gaussian 1-sigma) FreeParam.push_back(2.25); // raw skewness, for NR } - NuisParam.push_back(11.); //alpha, for NR model. See http://nest.physics.ucdavis.edu - NuisParam.push_back(1.1); //beta - NuisParam.push_back(0.0480); //gamma - NuisParam.push_back(-0.0533); //delta - NuisParam.push_back(12.6); //epsilon - NuisParam.push_back(0.3); //zeta - NuisParam.push_back(2.); //eta - NuisParam.push_back(0.3); //theta - NuisParam.push_back(2.); //iota - // last 3 are the secret extra parameters for additional flexibility - NuisParam.push_back(0.5); //changes sqrt in Qy equation - NuisParam.push_back(1.0); //makes low-E sigmoid an asymmetric one, for charge - NuisParam.push_back(1.0); //makes low-E sigmoid an asymmetric one, for light - + if ( ValidityTests::nearlyEqual ( ATOM_NUM, 18. ) ) { // liquid Ar + NuisParam[0] = 11.1025; // +/-1.10 Everything from https://docs.google.com/document/d/1vLg8vvY5bcdl4Ah4fzyE182DGWt0Wr7_FJ12_B10ujU + NuisParam[1] = 1.087399; // +/-0.025 + NuisParam[2] = 0.1; // +/-0.005 + NuisParam[3] = -0.0932; // +/-0.0095 + NuisParam[4] = 2.998; // +/-1.026 + NuisParam[5] = 0.3; // Fixed + NuisParam[6] = 2.94; // +/-0.12 + NuisParam[7] = W_DEFAULT / 1000.; + NuisParam[8] = DBL_MAX; + NuisParam[9] = 0.5; // square root + NuisParam[10] = 1.0; + NuisParam[11] = 1.0; + } + else { + NuisParam.push_back(11.); //alpha, for NR model. See http://nest.physics.ucdavis.edu + NuisParam.push_back(1.1); //beta + NuisParam.push_back(0.0480); //gamma + NuisParam.push_back(-0.0533); //delta + NuisParam.push_back(12.6); //epsilon + NuisParam.push_back(0.3); //zeta + NuisParam.push_back(2.); //eta + NuisParam.push_back(0.3); //theta + NuisParam.push_back(2.); //iota + // last 3 are the secret extra parameters for additional flexibility + NuisParam.push_back(0.5); //changes sqrt in Qy equation + NuisParam.push_back(1.0); //makes low-E sigmoid an asymmetric one, for charge + NuisParam.push_back(1.0); //makes low-E sigmoid an asymmetric one, for light + } } - + auto exec = execNEST ( detector, numEvts, type, eMin, eMax, inField, position, posiMuon, fPos, seed, no_seed, tZero ); delete detector; return exec; - + } NESTObservableArray runNESTvec ( VDetector* detector, INTERACTION_TYPE particleType, //func suggested by Xin Xiang, PD Brown U. for RG, LZ vector eList, vector> pos3dxyz, double inField, int seed ) { - + verbosity = false; NESTcalc n(detector); NESTresult result; QuantaResult quanta; double x, y, z, driftTime, vD; RandomGen::rndm()->SetSeed(seed); NuisParam = {11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1., 1.}; FreeParam = {1.,1.,0.10,0.5,0.19,2.25}; - vector scint, scint2, wf_amp; vector wf_time; + vector scint, scint2, wf_amp; vector wf_time; NESTObservableArray OutputResults; double useField; vector g2_params = n.CalculateG2(verbosity); double rho = n.SetDensity(detector->get_T_Kelvin(),detector->get_p_bar()); - - for ( long i = 0; i < eList.size(); ++i ) { + + for ( uint64_t i = 0; i < eList.size(); ++i ) { x = pos3dxyz[i][0]; y = pos3dxyz[i][1]; z = pos3dxyz[i][2]; @@ -281,18 +301,18 @@ NESTObservableArray runNESTvec ( VDetector* detector, INTERACTION_TYPE particleT OutputResults.s2c_phd.push_back(0.); } } - + delete detector; return OutputResults; } -int execNEST(VDetector* detector, unsigned long int numEvts, const string& type, +int execNEST(VDetector* detector, uint64_t numEvts, const string& type, double eMin, double eMax, double inField, string position, const string& posiMuon, double fPos, int seed, bool no_seed, double dayNumber ) { // Construct NEST class using detector object NESTcalc n(detector); - //Needed for python runability. THese are only valid for NR + // Here to get nestpy bindings to run... be careful of these! NuisParam = {11.,1.1,0.0480,-0.0533,12.6,0.3,2.,0.3,2.,0.5,1., 1.}; FreeParam = {1.,1.,0.10,0.5,0.19,2.25}; @@ -396,7 +416,7 @@ vector signal1, signal2, signalE, vTable; cerr << "Please choose gamma source. The allowed sources are:\n\"Co57\"\n\"Co60\"\n\"Cs137\"\nSource: "; cin >> gamma_source; if ( gamma_source == "Co60" ) { - cerr << "WARNING: This source is in the pair production range. Electron/positron pairs are not accounted for after initial interaction, and some" + 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 { @@ -422,10 +442,10 @@ vector signal1, signal2, signalE, vTable; return 1; } - + double maxTimeSep = DBL_MAX; if (type_num == Kr83m) { - if ( ValidityTests::nearlyEqual(eMin, 9.4) || ValidityTests::nearlyEqual(eMin, 32.1) || ValidityTests::nearlyEqual(eMin, 41.5) || ValidityTests::nearlyEqual(eMin, 41.55) ||ValidityTests::nearlyEqual(eMin, 41.6) && eMin != eMax ) { + if ( (ValidityTests::nearlyEqual(eMin, 9.4) || ValidityTests::nearlyEqual(eMin, 32.1) || ValidityTests::nearlyEqual(eMin, 41.5) || ValidityTests::nearlyEqual(eMin, 41.55) || ValidityTests::nearlyEqual(eMin, 41.6)) && eMin != eMax ) { maxTimeSep = eMax; if ( eMax <= 0. ) { cerr << "Max t sep must be +." << endl; return 1; } } else { @@ -455,12 +475,12 @@ vector signal1, signal2, signalE, vTable; //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) vector g2_params = n.CalculateG2(verbosity); g2 = fabs(g2_params[3]); double g1 = detector->get_g1(); - + double centralZ = (detector->get_gate() * 0.8 + detector->get_cathode() * 1.03) / 2.; // fid vol def usually shave more off the top, because of gas @@ -488,13 +508,13 @@ vector signal1, signal2, signalE, vTable; } if ( ( g1 * yieldsMax.PhotonYield ) > ( 2. * maxS1 ) && eMin != eMax && type_num != Kr83m ) cerr << "\nWARNING: Your energy maximum may be too high given your maxS1.\n"; - + if ( type_num < 6 ) massNum = 0; - if ( type_num == Kr83m ) massNum = maxTimeSep; + if ( type_num == Kr83m ) massNum = maxTimeSep; //use massNum to input maxTimeSep into GetYields(...) double keV = -999.; double timeStamp = dayNumber; vector keV_vec; - for (unsigned long int j = 0; j < numEvts; ++j) { + for (uint64_t j = 0; j < numEvts; ++j) { try { //timeStamp += tStep; //detector->set_eLife_us(5e1+1e3*(timeStamp/3e2)); //for E-recon when you've changed g1,g2-related stuff, redo line 341+ @@ -573,6 +593,8 @@ vector signal1, signal2, signalE, vTable; } break; default: + //keV = TestSpectra::ZeplinBackground(); //example of continuous ER/beta/gamma BG spec + //break; if ( eMin < 0. ) return 1; if ( eMax > 0. ) { if ( eMax > eMin ) @@ -702,11 +724,10 @@ vector signal1, signal2, signalE, vTable; if(MCtruthE) fprintf(stdout, "E_truth [keV]"); else fprintf(stdout, "E_recon [keV]"); } - fprintf(stdout, - "\tfield [V/cm]\ttDrift [us]\tX,Y,Z " - "[mm]\tNph\tNe-\tS1 [PE or phe]\tS1_3Dcor " - "[phd]\tspikeC(NON-INT)\tNe-Extr\tS2_rawArea [PE]\tS2_3Dcorr " - "[phd]\n"); + if ( seed == -2 ) + printf("\tfield [V/cm]\ttDrift [us]\tX,Y,Z [mm]\tNph\t\tNe-\t\tS1 [PE or phe]\tS1_3Dcor [phd]\tspikeC(NON-INT)\tNe-Extr\tS2_rawArea [PE]\tS2_3Dcorr [phd]\n"); + else + printf("\tfield [V/cm]\ttDrift [us]\tX,Y,Z [mm]\tNph\tNe-\tS1 [PE or phe]\tS1_3Dcor [phd]\tspikeC(NON-INT)\tNe-Extr\tS2_rawArea [PE]\tS2_3Dcorr [phd]\n"); } } if(ValidityTests::nearlyEqual(inField, -1.)) { @@ -884,7 +905,7 @@ vector signal1, signal2, signalE, vTable; smearPos[1] = xySmeared[1]; } - vector wf_time; + vector wf_time; vector wf_amp; vector scint = n.GetS1(quanta, truthPos[0], truthPos[1], truthPos[2], smearPos[0], smearPos[1], smearPos[2], @@ -908,7 +929,24 @@ vector signal1, signal2, signalE, vTable; signal1.push_back(scint[7]); else signal1.push_back(-999.); - + + double lowest1 = (double)detector->get_coinLevel(); + double lowest2 = minS2; + if ( lowest1 <= 0. ) lowest1 = 1.; + if ( lowest2 <= 0. ) lowest2 = 1.; + if ( (useS2 == 0 && logMax <= log10(maxS2/maxS1)) || + (useS2 == 1 && logMax <= log10(maxS2)) || + (useS2 == 2 && logMax <= log10(maxS2/maxS1)) ) { + if ( j == 0 ) + cerr << "err: You may be chopping off the upper half of your (ER?) band; increase logMax and/or maxS2" << endl; + } + if ( (useS2 == 0 && logMin >= log10(lowest2/lowest1)) || + (useS2 == 1 && logMin >= log10(lowest2)) || + (useS2 == 2 && logMin >= log10(lowest2/lowest1)) ) { + if ( j == 0 ) + cerr << "err: You may be chopping off the lower half of your (NR?) band; decrease logMin and/or minS2" << endl; + } + if(usePD == 0 && fabs(scint2[5]) > minS2 && scint2[5] < maxS2) signal2.push_back(scint2[5]); else if(usePD >= 1 && fabs(scint2[7]) > minS2 && scint2[7] < maxS2) @@ -1053,12 +1091,12 @@ vector signal1, signal2, signalE, vTable; // adjusted for 2-PE effect (LUX phd units) // scint2[8] = g2; // g2 = ExtEff * SE, light collection efficiency of EL in // gas gap (from CalculateG2) - + if ( truthPos[2] < detector->get_cathode() && verbosity && !BeenHere ) { BeenHere = true; fprintf ( stderr, "gamma-X i.e. MSSI may be happening. This may be why even high-E eff is <100%%. Check your cathode position definition.\n\n" ); } - + if(PrintSubThr || ( scint[0] > PHE_MIN && scint[1] > PHE_MIN && scint[2] > PHE_MIN && scint[3] > PHE_MIN && scint[4] > PHE_MIN && @@ -1086,7 +1124,7 @@ vector signal1, signal2, signalE, vTable; type == "muon" || type == "MIP" || type == "LIP" || type == "mu" || type == "mu-") { printf("%e\t%e\t%e\t", scint[2], scint[5], scint[7]); - printf("%li\t%e\t%e\n", (long) scint2[0], scint2[4], scint2[7]); + printf("%lli\t%e\t%e\n", (int64_t) scint2[0], scint2[4], scint2[7]); } else { printf("%.6f\t%.6f\t%.6f\t", scint[2], scint[5], scint[7]); // see GetS1 inside of NEST.cpp for full explanation @@ -1105,7 +1143,7 @@ vector signal1, signal2, signalE, vTable; return 1; } } - + if (verbosity) { if (eMin != eMax && type_num != Kr83m) { if (useS2 == 2) @@ -1204,7 +1242,7 @@ vector> GetBand(vector S1s, vector S2s, } int i = 0, j = 0; double s1c, numPts; - unsigned long reject[NUMBINS_MAX] = {0}; + uint64_t reject[NUMBINS_MAX] = {0}; if (resol) { numBins = 1; diff --git a/src/nestpy/execNEST.hh b/src/nestpy/execNEST.hh index 249b2ba..67f5130 100644 --- a/src/nestpy/execNEST.hh +++ b/src/nestpy/execNEST.hh @@ -31,7 +31,7 @@ vector> GetBand(vector S1s, vector S2s, void GetEnergyRes(vector Es); -int execNEST(VDetector* detector, unsigned long int numEvts, const string& type, +int execNEST(VDetector* detector, uint64_t numEvts, const string& type, double eMin, double eMax, double inField, string position, const string& posiMuon, double fPos, int seed, bool no_seed, double dayNumber); NESTObservableArray runNESTvec(VDetector* detector, INTERACTION_TYPE scatterType, From e4558356edaad8d4bb705b4d6237f0144bab2cc6 Mon Sep 17 00:00:00 2001 From: sophiaandaloro Date: Mon, 1 Mar 2021 12:41:52 -0600 Subject: [PATCH 2/2] remove unneeded includes --- MANIFEST.in | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/MANIFEST.in b/MANIFEST.in index 143f18e..0ac0666 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,6 +1,6 @@ recursive-include src *.cpp recursive-include src *.hh recursive-include include *.h +recursive-include src *.py include *.md -include *.txt -global-include * \ No newline at end of file +include *.txt \ No newline at end of file