From 4b0f5439b275548be23a0737ec93bb18885e802b Mon Sep 17 00:00:00 2001 From: mszydagis Date: Sat, 11 Jun 2022 12:40:03 -0400 Subject: [PATCH] Piecemeal track method works now at Q_beta_beta, with basement-level customizable Fano-LIKE factor, and 4um step size. -MMS --- examples/rootNEST.cpp | 3 ++- src/NEST.cpp | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/rootNEST.cpp b/examples/rootNEST.cpp index ee3d75f8..27d76bf3 100644 --- a/examples/rootNEST.cpp +++ b/examples/rootNEST.cpp @@ -791,7 +791,8 @@ void GetFile(char* fileName) { S2cor_phd.push_back(-999.); } } - fclose(ifp); if ( E_keV.size() < 100000 ) { skewness = 0; cerr << "WARNING: Not enough stats (at least 10^5 events) for skew fits so doing Gaussian" << endl; } + fclose(ifp); + if ( E_keV.size() < 100000 && numBins > 1 ) { skewness = 0; cerr << "WARNING: Not enough stats (at least 10^5 events) for skew fits so doing Gaussian" << endl; } if (numBins == 1) { TH1F* HistogramArray = new TH1F[3]; diff --git a/src/NEST.cpp b/src/NEST.cpp index 71a9e0ef..ba1dd4a3 100644 --- a/src/NEST.cpp +++ b/src/NEST.cpp @@ -585,8 +585,8 @@ NESTresult NESTcalc::GetYieldERdEOdxBasis(const std::vector &dEOdxParam, 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 +603,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());