IAP GITLAB

Commit 9530b4f5 authored by Felix Riehn's avatar Felix Riehn Committed by Ralf Ulrich

Date: Sat May 25 11:11:00 2019 +0100

 On branch 188-sibyll-decay-of-sigmabar0-yields-sigmabar0
 You are currently cherry-picking commit 9c650ba.

 Changes to be committed:
	modified:   Documentation/Examples/boundary_example.cc
	modified:   Documentation/Examples/cascade_example.cc
	modified:   Documentation/Examples/vertical_EAS.cc
	modified:   Processes/Sibyll/Decay.cc
	modified:   Processes/Sibyll/Decay.h
	modified:   Processes/Sibyll/Interaction.cc
	modified:   Processes/Sibyll/Interaction.h
	modified:   Processes/Sibyll/sibyll2.3c.f
	modified:   Processes/Sibyll/testSibyll.cc
parent 61bcbb48
......@@ -121,9 +121,7 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll;
process::sibyll::Decay decay{{particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::KPlus, particles::Code::KMinus,
particles::Code::K0Long, particles::Code::K0Short}};
process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
......
......@@ -28,77 +28,62 @@ using SetupParticle = corsika::setup::Stack::ParticleType;
namespace corsika::process::sibyll {
Decay::Decay(vector<particles::Code> pParticles)
: fTrackedParticles(pParticles) {}
Decay::Decay() {}
Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
void Decay::Init() {
SetHadronsUnstable();
SetParticleListStable(fTrackedParticles);
// switch off decays to avoid internal decay chains
SetAllStable();
}
void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
/*
Sibyll is hadronic generator
only hadrons decay
*/
// set particles unstable
SetHadronsUnstable();
cout << "Interaction: setting tracked hadrons stable.." << endl;
for (auto p : particleList) Decay::SetStable(p);
void Decay::SetStable(const vector<particles::Code> vParticleList) {
for (auto p : vParticleList) Decay::SetStable(p);
}
void Decay::SetUnstable(const particles::Code pCode) {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
void Decay::SetUnstable(const vector<particles::Code> vParticleList) {
for (auto p : vParticleList) Decay::SetUnstable(p);
}
bool Decay::IsStable(const particles::Code vCode) {
return abs(process::sibyll::ConvertToSibyllRaw(vCode)) <= 0 ? true : false;
}
bool Decay::IsUnstable(const particles::Code vCode) {
return abs(process::sibyll::ConvertToSibyllRaw(vCode)) > 0 ? true : false;
}
void Decay::SetDecay(const particles::Code vCode, const bool vMakeUnstable) {
vMakeUnstable ? SetUnstable(vCode) : SetStable(vCode);
}
void Decay::SetUnstable(const particles::Code vCode) {
cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << endl;
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
void Decay::SetStable(const particles::Code pCode) {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
void Decay::SetStable(const particles::Code vCode) {
cout << "Sibyll::Interaction: setting " << vCode << " stable.." << endl;
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
void Decay::SetAllStable() {
// name? also makes EM particles stable
cout << "Decay: setting all particles stable.." << endl;
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. particles::ListOfParticles()
for (auto& p : corsika2sibyll) {
// cout << (int)p << endl;
const int sibCode = static_cast<int>(p);
// skip unknown and antiparticles
if (sibCode < 1) continue;
s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
}
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
}
void Decay::SetHadronsUnstable() {
// name? also makes EM particles stable
void Decay::SetAllUnstable() {
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
}
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. particles::ListOfParticles()
cout << "Sibyll: setting hadrons unstable.." << endl;
// make ALL particles unstable, then set EM stable
for (int sibCode : corsika2sibyll) {
if (sibCode < 1) continue;
s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
}
// set Leptons and Proton and Neutron stable
// use stack to loop over particles
constexpr particles::Code particleList[] = {
particles::Code::Proton, particles::Code::Neutron, particles::Code::Electron,
particles::Code::Positron, particles::Code::NuE, particles::Code::NuEBar,
particles::Code::MuMinus, particles::Code::MuPlus, particles::Code::NuMu,
particles::Code::NuMuBar};
for (auto p : particleList) {
const int sibid = process::sibyll::ConvertToSibyllRaw(p);
s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
}
void Decay::PrintDecayConfig(const particles::Code vCode) {
cout << "Decay: Sibyll decay configuration:" << endl;
const int sibCode = process::sibyll::ConvertToSibyllRaw(vCode);
const int absSibCode = abs(sibCode);
cout << vCode << " is ";
if (s_csydec_.idb[absSibCode - 1] <= 0)
cout << "stable" << endl;
else
cout << "unstable" << endl;
}
template <>
......@@ -148,12 +133,16 @@ namespace corsika::process::sibyll {
// remember position
Point const decayPoint = vP.GetPosition();
TimeType const t0 = vP.GetTime();
// set all particles/hadrons unstable
// setHadronsUnstable();
// remember if particles is unstable
// auto const priorIsUnstable = IsUnstable(pCode);
// switch on decay for this particle
SetUnstable(pCode);
PrintDecayConfig(pCode);
// call sibyll decay
cout << "Decay: calling Sibyll decay routine.." << endl;
decsib_();
// reset to stable
SetStable(pCode);
// print output
......
......@@ -21,23 +21,33 @@ namespace corsika::process {
namespace sibyll {
class Decay : public corsika::process::DecayProcess<Decay> {
std::vector<particles::Code> const fTrackedParticles;
int fCount = 0;
public:
Decay(std::vector<particles::Code>);
Decay();
~Decay();
void Init();
void SetParticleListStable(const std::vector<particles::Code>);
void SetStable(const std::vector<particles::Code>);
void SetUnstable(const std::vector<particles::Code>);
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
void SetAllUnstable();
void SetAllStable();
bool IsStable(const corsika::particles::Code);
bool IsUnstable(const corsika::particles::Code);
void SetDecay(const particles::Code, const bool);
void PrintDecayConfig(const corsika::particles::Code);
void SetHadronsUnstable();
template <typename TParticle>
corsika::units::si::TimeType GetLifetime(TParticle const&) const;
/**
In this function SIBYLL is called to produce to decay the input particle.
*/
template <typename TProjectile>
void DoDecay(TProjectile&);
};
......
......@@ -47,46 +47,38 @@ namespace corsika::process::sibyll {
// initialize Sibyll
if (!fInitialized) {
sibyll_ini_();
// any decays in sibyll? if yes need to define which particles
if (fInternalDecays) {
// define which particles are passed to corsika, i.e. which particles make it into
// history even very shortlived particles like charm or pi0 are of interest here
const std::vector<particles::Code> hadronsWeWantTrackedByCorsika = {
particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::Pi0, particles::Code::KMinus,
particles::Code::KPlus, particles::Code::K0Long,
particles::Code::K0Short, particles::Code::SigmaPlus,
particles::Code::SigmaMinus, particles::Code::Lambda0,
particles::Code::Xi0, particles::Code::XiMinus,
particles::Code::OmegaMinus, particles::Code::DPlus,
particles::Code::DMinus, particles::Code::D0,
particles::Code::D0Bar};
Interaction::SetParticleListStable(hadronsWeWantTrackedByCorsika);
}
fInitialized = true;
}
}
void Interaction::SetParticleListStable(
std::vector<particles::Code> const& particleList) {
for (auto p : particleList) Interaction::SetStable(p);
void Interaction::SetStable(std::vector<particles::Code> const& vParticleList) {
for (auto p : vParticleList) Interaction::SetStable(p);
}
void Interaction::SetUnstable(std::vector<particles::Code> const& vParticleList) {
for (auto p : vParticleList) Interaction::SetUnstable(p);
}
void Interaction::SetUnstable(const particles::Code pCode) {
cout << "Sibyll::Interaction: setting " << pCode << " unstable.." << endl;
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
void Interaction::SetUnstable(const particles::Code vCode) {
cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << endl;
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
void Interaction::SetStable(const particles::Code pCode) {
cout << "Sibyll::Interaction: setting " << pCode << " stable.." << endl;
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
void Interaction::SetStable(const particles::Code vCode) {
cout << "Sibyll::Interaction: setting " << vCode << " stable.." << endl;
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
void Interaction::SetAllStable() {
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
}
void Interaction::SetAllUnstable() {
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
}
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
......@@ -325,8 +317,16 @@ namespace corsika::process::sibyll {
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
// running decays
if (fInternalDecays) decsib_();
if (fInternalDecays) {
// particles that decay internally will never appear on the corsika stack
// switch on all decays except for the particles we want to take part in the
// tracking
SetAllUnstable();
SetStable(fTrackedParticles);
decsib_();
// reset
SetAllStable();
}
// print final state
int print_unit = 6;
sib_list_(print_unit);
......
......@@ -31,9 +31,13 @@ namespace corsika::process::sibyll {
void Init();
void SetParticleListStable(std::vector<particles::Code> const&);
void SetStable(std::vector<particles::Code> const&);
void SetUnstable(std::vector<particles::Code> const&);
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
void SetAllUnstable();
void SetAllStable();
bool WasInitialized() { return fInitialized; }
bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) const {
......@@ -65,7 +69,18 @@ namespace corsika::process::sibyll {
private:
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
// FOR NOW keep trackedParticles private, could be configurable
std::vector<particles::Code> const fTrackedParticles = {
particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::Pi0, particles::Code::KMinus,
particles::Code::KPlus, particles::Code::K0Long,
particles::Code::K0Short, particles::Code::SigmaPlus,
particles::Code::SigmaMinus, particles::Code::Lambda0,
particles::Code::Xi0, particles::Code::XiMinus,
particles::Code::OmegaMinus, particles::Code::DPlus,
particles::Code::DMinus, particles::Code::D0,
particles::Code::MuMinus, particles::Code::MuPlus,
particles::Code::D0Bar};
const bool fInternalDecays = true;
const corsika::units::si::HEPEnergyType fMinEnergyCoM =
10. * 1e9 * corsika::units::si::electronvolt;
......
......@@ -478,7 +478,7 @@ C-----------------------------------------------------------------------
CALL BLOCK_INI
CALL NUC_GEOM_INI
CALL SIG_AIR_INI
CALL DEC_INI
c CALL DEC_INI
c... charm frag. normalisation
CALL ZNORMAL
......
......@@ -69,6 +69,7 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
using namespace corsika::units::si;
using namespace corsika::units;
......@@ -161,14 +162,29 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
const std::vector<particles::Code> particleList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
Decay model(particleList);
Decay model;
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
SECTION("DecayConfiguration") {
Decay model;
const std::vector<particles::Code> particleTestList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::Lambda0Bar, particles::Code::NuE, particles::Code::D0Bar};
for (auto& pCode : particleTestList) {
model.SetUnstable(pCode);
// get state of sibyll internal config
REQUIRE(0 <= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
model.SetStable(pCode);
// get state of sibyll internal config
REQUIRE(0 >= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
}
}
}
Markdown is supported
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