From 5fa92278ac074ce1cde5335f62dc67d9b677ce97 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Wed, 1 Apr 2020 15:05:30 +0200
Subject: [PATCH] alternate_ as member variable
---
Processes/QGSJetII/Interaction.cc | 149 ++++++++++++++++--------------
Processes/QGSJetII/Interaction.h | 4 +-
2 files changed, 82 insertions(+), 71 deletions(-)
diff --git a/Processes/QGSJetII/Interaction.cc b/Processes/QGSJetII/Interaction.cc
index 92288e76c..17329a2fe 100644
--- a/Processes/QGSJetII/Interaction.cc
+++ b/Processes/QGSJetII/Interaction.cc
@@ -12,9 +12,8 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/FourVector.h>
-#include <corsika/process/qgsjetII/ParticleConversion.h>
+#include <corsika/geometry/QuantityVector.h>
#include <corsika/process/qgsjetII/QGSJetIIFragmentsStack.h>
#include <corsika/process/qgsjetII/QGSJetIIStack.h>
#include <corsika/process/qgsjetII/qgsjet-II-04.h>
@@ -215,13 +214,15 @@ namespace corsika::process::qgsjetII {
if (particles::IsNucleus(corsikaBeamId)) beamA = vP.GetNuclearA();
const HEPEnergyType projectileEnergyLabPerNucleon =
- particles::IsNucleus(corsikaBeamId) ? projectileEnergyLab/beamA : projectileEnergyLab;
-
+ particles::IsNucleus(corsikaBeamId) ? projectileEnergyLab / beamA
+ : projectileEnergyLab;
+
cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << projectileMomentumLab.GetComponents() / 1_GeV
<< endl;
cout << "Interaction: etarget lab: " << targetEnergyLab / 1_GeV << endl
- << "Interaction: ptarget lab: " << targetMomentumLab.GetComponents() / 1_GeV << endl;
+ << "Interaction: ptarget lab: " << targetMomentumLab.GetComponents() / 1_GeV
+ << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
@@ -248,7 +249,7 @@ namespace corsika::process::qgsjetII {
cross_section_of_components[i] = sigProd;
}
- const auto targetCode =
+ const auto targetCode =
mediumComposition.SampleTarget(cross_section_of_components, rng_);
cout << "Interaction: target selected: " << targetCode << endl;
@@ -268,22 +269,21 @@ namespace corsika::process::qgsjetII {
throw std::runtime_error("QgsjetII target outside range.");
// beam id for qgsjetII
- std::uniform_real_distribution<double> select;
-
QgsjetIICode qgsjet_beam_code;
- if (corsikaBeamId == particles::Code::Nucleus) {
- if (select(rng_)>0.5)
- qgsjet_beam_code = QgsjetIICode::Proton;
- else
- qgsjet_beam_code = QgsjetIICode::Neutron;
- } else { // it is a nucleus
+ if (corsikaBeamId == particles::Code::Nucleus) {
+ std::array<QgsjetIICode, 2> constexpr nucleons = {QgsjetIICode::Proton,
+ QgsjetIICode::Neutron};
+ std::uniform_int_distribution select(0, 1);
+ qgsjet_beam_code = nucleons[select(rng_)];
+ } else { // it is an "elementary" particle
qgsjet_beam_code = process::qgsjetII::ConvertToQgsjetII(corsikaBeamId);
// from conex
if (qgsjet_beam_code == QgsjetIICode::Pi0 or
- qgsjet_beam_code == QgsjetIICode::Rho0) { // replace pi0 or rho0 with pi+/pi- in alternating sequence
- static QgsjetIICode alternate = QgsjetIICode::PiPlus;
- qgsjet_beam_code = alternate;
- alternate = (alternate == QgsjetIICode::PiPlus ? QgsjetIICode::PiMinus : QgsjetIICode::PiPlus);
+ qgsjet_beam_code == QgsjetIICode::Rho0) { // replace pi0 or rho0 with pi+/pi-
+ // in alternating sequence
+ qgsjet_beam_code = alternate_;
+ alternate_ = (alternate_ == QgsjetIICode::PiPlus ? QgsjetIICode::PiMinus
+ : QgsjetIICode::PiPlus);
}
// replace lambda by neutron
if (qgsjet_beam_code == QgsjetIICode::Lambda0)
@@ -294,15 +294,17 @@ namespace corsika::process::qgsjetII {
}
int qgsjet_beam_code_int = static_cast<QgsjetIICodeIntType>(qgsjet_beam_code);
-
+
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << endl;
count_++;
- qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode, targetQgsCode);
+ qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode,
+ targetQgsCode);
// this is from CRMC, is this REALLY needed ???
- qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode, targetQgsCode);
+ qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode,
+ targetQgsCode);
qgconf_();
-
+
// bookkeeping
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
@@ -312,8 +314,8 @@ namespace corsika::process::qgsjetII {
// CoM frame definition in QgsjetII projectile: +z
auto const& originalCS = projectileMomentumLab.GetCoordinateSystem();
geometry::CoordinateSystem const zAxisFrame =
- originalCS.RotateToZ(projectileMomentumLab);
-
+ originalCS.RotateToZ(projectileMomentumLab);
+
// nuclear projectile fragments
QGSJetIIFragmentsStack qfs;
for (auto& fragm : qfs) {
@@ -322,30 +324,31 @@ namespace corsika::process::qgsjetII {
int Z = 0;
switch (A) {
case 1: { // proton/neutron
- std::uniform_real_distribution<double> select;
- if (select(rng_)>0.5) {
- idFragm = particles::Code::Proton;
- Z = 1;
- } else {
- idFragm = particles::Code::Neutron;
- Z = 0;
- }
-
- const HEPMassType projectileMass = particles::GetMass(idFragm);
-
- auto momentum = geometry::Vector(
- zAxisFrame,
- corsika::geometry::QuantityVector<hepmomentum_d>{0.0_GeV, 0.0_GeV,
- sqrt((projectileEnergyLab + projectileMass) *
- (projectileEnergyLab - projectileMass))});
-
- auto const energy = sqrt(momentum.squaredNorm() + square(projectileMass));
- momentum.rebase(originalCS); // transform back into standard lab frame
- std::cout << "secondary fragment> id=" << idFragm << " p=" << momentum.GetComponents() << std::endl;
+ std::uniform_real_distribution<double> select;
+ if (select(rng_) > 0.5) {
+ idFragm = particles::Code::Proton;
+ Z = 1;
+ } else {
+ idFragm = particles::Code::Neutron;
+ Z = 0;
+ }
+
+ const HEPMassType projectileMass = particles::GetMass(idFragm);
+
+ auto momentum = geometry::Vector(
+ zAxisFrame, corsika::geometry::QuantityVector<hepmomentum_d>{
+ 0.0_GeV, 0.0_GeV,
+ sqrt((projectileEnergyLab + projectileMass) *
+ (projectileEnergyLab - projectileMass))});
+
+ auto const energy = sqrt(momentum.squaredNorm() + square(projectileMass));
+ momentum.rebase(originalCS); // transform back into standard lab frame
+ std::cout << "secondary fragment> id=" << idFragm
+ << " p=" << momentum.GetComponents() << std::endl;
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- geometry::Point, units::si::TimeType>{
- idFragm, energy, momentum, pOrig, tOrig});
+ geometry::Point, units::si::TimeType>{idFragm, energy, momentum,
+ pOrig, tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
} break;
@@ -365,22 +368,25 @@ namespace corsika::process::qgsjetII {
}
if (idFragm == particles::Code::Nucleus) {
- const HEPMassType nucleusMass = particles::Proton::GetMass() * Z + particles::Neutron::GetMass() * (A-Z);
- auto momentum = geometry::Vector(
- zAxisFrame,
- geometry::QuantityVector<hepmomentum_d>{0.0_GeV, 0.0_GeV,
- sqrt((projectileEnergyLabPerNucleon*A + nucleusMass) *
- (projectileEnergyLabPerNucleon*A - nucleusMass))});
-
- auto const energy = sqrt(momentum.squaredNorm() + square(nucleusMass));
- momentum.rebase(originalCS); // transform back into standard lab frame
- std::cout << "secondary fragment> id=" << idFragm << " p=" << momentum.GetComponents() << " A=" << A << " Z=" << Z << std::endl;
- auto pnew = vP.AddSecondary(
- tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- geometry::Point, units::si::TimeType, unsigned short, unsigned short>{
- idFragm, energy, momentum, pOrig, tOrig, A, Z});
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
+ const HEPMassType nucleusMass =
+ particles::Proton::GetMass() * Z + particles::Neutron::GetMass() * (A - Z);
+ auto momentum = geometry::Vector(
+ zAxisFrame, geometry::QuantityVector<hepmomentum_d>{
+ 0.0_GeV, 0.0_GeV,
+ sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
+ (projectileEnergyLabPerNucleon * A - nucleusMass))});
+
+ auto const energy = sqrt(momentum.squaredNorm() + square(nucleusMass));
+ momentum.rebase(originalCS); // transform back into standard lab frame
+ std::cout << "secondary fragment> id=" << idFragm
+ << " p=" << momentum.GetComponents() << " A=" << A << " Z=" << Z
+ << std::endl;
+ auto pnew = vP.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+ geometry::Point, units::si::TimeType, unsigned short, unsigned short>{
+ idFragm, energy, momentum, pOrig, tOrig, A, Z});
+ Plab_final += pnew.GetMomentum();
+ Elab_final += pnew.GetEnergy();
}
}
@@ -391,14 +397,17 @@ namespace corsika::process::qgsjetII {
auto momentum = psec.GetMomentum(zAxisFrame);
auto const energy = psec.GetEnergy();
- momentum.rebase(originalCS); // transform back into standard lab frame
- std::cout << "secondary fragment> id=" << process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()) << " p=" << momentum.GetComponents() << std::endl;
- auto pnew = vP.AddSecondary(
- tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- geometry::Point, units::si::TimeType>{
- process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()), energy, momentum, pOrig, tOrig});
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
+ momentum.rebase(originalCS); // transform back into standard lab frame
+ std::cout << "secondary fragment> id="
+ << process::qgsjetII::ConvertFromQgsjetII(psec.GetPID())
+ << " p=" << momentum.GetComponents() << std::endl;
+ auto pnew = vP.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{
+ process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()), energy, momentum,
+ pOrig, tOrig});
+ Plab_final += pnew.GetMomentum();
+ Elab_final += pnew.GetEnergy();
}
cout << "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/ << endl
<< "Elab_final=" << Elab_final / 1_GeV
diff --git a/Processes/QGSJetII/Interaction.h b/Processes/QGSJetII/Interaction.h
index 0b2804a23..d06a07843 100644
--- a/Processes/QGSJetII/Interaction.h
+++ b/Processes/QGSJetII/Interaction.h
@@ -13,6 +13,7 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
+#include <corsika/process/qgsjetII/ParticleConversion.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
@@ -25,6 +26,7 @@ namespace corsika::process::qgsjetII {
std::string data_path_;
int count_ = 0;
bool initialized_ = false;
+ QgsjetIICode alternate_ = QgsjetIICode::PiPlus; // for pi0, rho0 projectiles
public:
Interaction(const std::string& dataPath = "");
@@ -58,7 +60,7 @@ namespace corsika::process::qgsjetII {
private:
corsika::random::RNG& rng_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("qgran");
- const int maxMassNumber_ = 208;
+ static constexpr int maxMassNumber_ = 208;
};
} // namespace corsika::process::qgsjetII
--
GitLab