remove massless particles - for now

carmeloevoli · Nov 17, 2023 · 08f6046 · 08f6046
1 parent 41899cf
commit 08f6046
Show file tree

Hide file tree

Showing 3 changed files with 45 additions and 33 deletions.
diff --git a/examples/printPhotopionRates.cpp b/examples/printPhotopionRates.cpp
@@ -94,51 +94,59 @@ void plot_sampled_epsilon() {
  }
 }
 
+// double sum_pion_energy(const ParticleStack& stack) {
+// double pion_energy = 0.;
+// for (const auto& particle : stack) {
+// if (pidIsPion(particle.getPid())) pion_energy += particle.getGamma() * SI::pionMassC2;
+// }
+
+auto getPionsEnergy = [](double E, const Particle& p) {
+ auto value = (pidIsPion(p.getPid())) ? p.getGamma() * SI::pionMassC2 : 0.;
+ return E + value;
+};
+
 void plot_pion_energies() {
  const auto cmb = std::make_shared<photonfields::CMB>();
  RandomNumberGenerator rng = utils::RNG<double>(1234);
- const auto pppcmb = std::make_shared<interactions::PhotoPionProduction>(cmb);
- utils::OutputFile out("SimProp_photopion_pionenergy_sample.txt");
  const size_t N = 100000;
  const double z = 0;
  std::vector<double> protonEnergies = {1e18 * SI::eV, 1e19 * SI::eV, 1e20 * SI::eV, 1e21 * SI::eV};
- for (size_t i = 0; i < N; ++i) {
- for (auto E : protonEnergies) {
- const auto Gamma = E / SI::protonMassC2;
- auto finalState = pppcmb->finalState({proton, z, Gamma}, z, rng);
- out << finalState[1].getGamma() * SI::pionMassC2 / E << "\t";
+ {
+ const auto pppcmb = std::make_shared<interactions::PhotoPionProduction>(cmb);
+ utils::OutputFile out("SimProp_photopion_pionenergy_sample.txt");
+ for (size_t i = 0; i < N; ++i) {
+ for (auto E : protonEnergies) {
+ const auto Gamma = E / SI::protonMassC2;
+ auto finalState = pppcmb->finalState({proton, z, Gamma}, z, rng);
+ double piEnergy = std::accumulate(finalState.begin(), finalState.end(), 0., getPionsEnergy);
+ out << piEnergy / E << "\t";
+ }
+ out << "\n";
+ }
+ }
+ {
+ const auto pppcmb = std::make_shared<interactions::PhotoPionProductionSophia>(cmb);
+ utils::OutputFile out("SimProp_photopionsophia_pionenergy_sample.txt");
+ for (size_t i = 0; i < N; ++i) {
+ for (auto E : protonEnergies) {
+ const auto Gamma = E / SI::protonMassC2;
+ auto finalState = pppcmb->finalState({proton, z, Gamma}, z, rng);
+ auto piEnergy = std::accumulate(finalState.begin(), finalState.end(), 0., getPionsEnergy);
+ out << piEnergy / E << "\t";
+ }
+ out << "\n";
  }
- out << "\n";
  }
 }
 
-// void plot_inelasticity() {
-// RandomNumberGenerator rng = utils::RNG<double>(234);
-// const auto cmb = std::make_shared<photonfields::CMB>();
-// const auto pppcmb = std::make_shared<interactions::PhotoPionProduction>(cmb);
-// utils::OutputFile out("test_photopion_inelasticity.txt");
-// const auto protonEnergies = utils::LogAxis<double>(1e18 * SI::eV, 1e21 * SI::eV, 3 * 8);
-// const size_t N = 10000;
-// const double z = 0;
-// for (auto& E : protonEnergies) {
-// double Y[N];
-// for (size_t i = 0; i < N; ++i) {
-// auto finalState = pppcmb->finalState({proton, z, E / SI::protonMassC2}, z, rng);
-// Y[i] = finalState[1].getGamma() * SI::pionMassC2 / E;
-// }
-// out << E / SI::eV << " ";
-// out << gsl_stats_mean(Y, 1, N) << " " << gsl_stats_sd(Y, 1, N) << "\n";
-// }
-// }
-
 int main() {
  try {
  utils::startup_information();
- // plot_rates();
- // plot_sampled_s();
- // plot_sampled_epsilon();
+ utils::Timer timer("main timer");
+ plot_rates();
+ plot_sampled_s();
+ plot_sampled_epsilon();
  plot_pion_energies();
- // plot_inelasticity();
  } catch (const std::exception& e) {
  LOGE << "exception caught with message: " << e.what();
  }

diff --git a/src/core/pid.cpp b/src/core/pid.cpp
@@ -75,6 +75,11 @@ double getPidMass(const PID& pid) {
  return SI::electronMassC2;
  else if (pid == pionNeutral || pid == pionMinus || pid == pionPlus)
  return SI::pionMassC2;
+ else if (pid == photon)
+ return 0.;
+ else if (pid == neutrino_e || pid == neutrino_mu || pid == antineutrino_e ||
+ pid == antineutrino_mu)
+ return 0.;
  else
  throw std::invalid_argument("mass not available for this pid");
 }

diff --git a/src/interactions/PhotoPionProductionSophia.cpp b/src/interactions/PhotoPionProductionSophia.cpp
@@ -87,9 +87,8 @@ std::vector<Particle> PhotoPionProductionSophia::finalState(const Particle& inco
  PID pid = ID_sophia_to_SimProp(seo.outPartID[i]);
  auto mass = getPidMass(pid);
  auto E = seo.outPartP[3][i] * SI::GeV;
- outgoingParticle.push_back({pid, zInteractionPoint, E / mass, w});
+ if (mass > 0.) outgoingParticle.push_back({pid, zInteractionPoint, E / mass, w});
  }
-
  return outgoingParticle;
 }