Sub-2GeV mass WIMP energy spectrum calculation

[grischbieter]

There's a problem with the energy spectrum calculation for sub-2GeV mass WIMPs in TestSpectra.cpp. If mass < 2 GeV, the number of sample points becomes large ( > 10,000 ) and the EnergySpec array can only take 10,000 entries. Plus, if the bin size is too small for any WIMP, the calculation can break down and think that the maximum recoil energy is too small, truncating the expected recoil spectrum. I suggest committing replacing the WIMP_prep_spectrum function with these edits:

TestSpectra::WIMP_spectrum_prep TestSpectra::WIMP_prep_spectrum(double mass,
                                                                double eStep) {
  WIMP_spectrum_prep spectrum;
  double divisor, x1, x2;
  vector<double> EnergySpec;
  int numberPoints;
  if (mass < 2.0) {  // GeV/c^2
    divisor = 10. / eStep;
    numberPoints = int(1000. / eStep);
  } else if (mass < 10.) {
    divisor = 10. / eStep;
    numberPoints = int(1000. / eStep);
  } else {
    divisor = 1.0 / eStep;
    numberPoints = int(100. / eStep);
  }
  int nZeros = 0; //keep track of the number of zeros in a row
  for (int i = 0; i < (numberPoints + 1); i++) {
    EnergySpec.push_back( WIMP_dRate(double(i) / divisor, mass) );
    if ( EnergySpec[i] == 0. ) nZeros++;
    else nZeros = 0; //reset the count if EnergySpec[i] != zero
    if ( nZeros == 100 ) break; //quit the for-loop once we're sure we're only getting zeros
  }

  for (long i = 0; i < 1000000; i++) {
    spectrum.integral += WIMP_dRate(double(i) / 1e4, mass) / 1e4;
  }
  for (int i = 0; i < (int) EnergySpec.size() - 2; i++) {
    x1 = double(i) / divisor;
    x2 = double(i + 1) / divisor;
    spectrum.base[i] = EnergySpec[i + 1] *
                       pow(EnergySpec[i + 1] / EnergySpec[i], x2 / (x1 - x2));
    spectrum.exponent[i] = log(EnergySpec[i + 1] / EnergySpec[i]) / (x1 - x2);
    if (spectrum.base[i] > 0. && spectrum.base[i] < DBL_MAX &&
        spectrum.exponent[i] > 0. && spectrum.exponent[i] < DBL_MAX )
      ;  // spectrum.integral+=spectrum.base[i]/spectrum.exponent[i]*(exp(-spectrum.exponent[i]*x1)-exp(-spectrum.exponent[i]*x2));
    else {
      if ( EnergySpec[i+2] > 10. ) { //i.e. the calculation stopped before event rate was low
        cerr << "ERROR: WIMP E_step is too small! Increase it slightly to avoid noise in the calculation." << endl;
        exit(0.);
      }
      spectrum.xMax = double(i - 1) / divisor;
      if (spectrum.xMax <= 0.0) {
        cerr << "ERROR: The maximum possible WIMP recoil is negative, which "
                "usually means your E_step is too small."
             << endl;
        exit(0);
      }
      break;
    }
  }

  spectrum.divisor = divisor;
  return spectrum;
}

where the main changes are the step size for sub-2GeV masses and the warning statement if the energy spectrum calculation breaks a bit early. In addition, for low-mass WIMPs, the current calculation fills EnergySpec with many many zeros. I added a quick if-statement that cuts the code if there have been a bunch of zeros, indicating that the calculation is done.