IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 5fa92278 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

alternate_ as member variable

parent c607f023
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Processes/QGSJetII/Interaction.cc
+ 79
− 70
View file @ 5fa92278
......@@ -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
......
Processes/QGSJetII/Interaction.h
+ 3
− 1
View file @ 5fa92278
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum