Piecemeal track method works now at Q_beta_beta, with basement-level …

…customizable Fano-LIKE factor, and 4um step size. -MMS
NESTCollaboration · Jun 11, 2022 · 4b0f543 · 4b0f543
1 parent e626956
commit 4b0f543
Show file tree

Hide file tree

Showing 2 changed files with 5 additions and 4 deletions.
diff --git 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
@@ -585,8 +585,8 @@ NESTresult NESTcalc::GetYieldERdEOdxBasis(const std::vector<double> &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<double> &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());