diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 8159853875f1155931be45057b15dbd90b6cf153..81c6fc535538d6789ee6170fbe44e6b2fb7fedd5 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -90,7 +90,7 @@ target_link_libraries (vertical_EAS
CORSIKAlogging
CORSIKArandom
ProcessSibyll
- # ProcessPythia
+ ProcessPythia
ProcessUrQMD
ProcessSwitch
CORSIKAcascade
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index d985992eb0707911ac2e9bfadd2fed6409442359..f95c9ae2255ac0512a5457a63840d8950e5fd3fc 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -121,17 +121,13 @@ int main() {
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
- const std::vector<particles::Code> trackedHadrons = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::sibyll::Interaction sibyll(env);
process::pythia::Interaction pythia;
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
- // process::sibyll::Decay decay(trackedHadrons);
- process::pythia::Decay decay(trackedHadrons);
+ // process::sibyll::Decay decay;
+ process::pythia::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 455939072950e977698863c52b7c15f3ef6748a7..8dd6f827df46e28c98c2354c2e122734250d87e1 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -20,9 +20,9 @@
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Plane.h>
@@ -32,6 +32,8 @@
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
+#include <corsika/process/pythia/Decay.h>
+
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/process/particle_cut/ParticleCut.h>
@@ -63,7 +65,7 @@ using namespace corsika::units::si;
void registerRandomStreams() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
- //random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+ random::RNGManager::GetInstance().RegisterRandomStream("pythia");
random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
random::RNGManager::GetInstance().SeedAll();
@@ -139,7 +141,30 @@ int main() {
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
- process::sibyll::Decay decay;
+
+ process::pythia::Decay decayPythia;
+
+ // use sibyll decay routine for decays of particles unknown to pythia
+ process::sibyll::Decay decaySibyll({
+ Code::N1440Plus,
+ Code::N1440MinusBar,
+ Code::N1440_0,
+ Code::N1440_0Bar,
+ Code::N1710Plus,
+ Code::N1710MinusBar,
+ Code::N1710_0,
+ Code::N1710_0Bar,
+
+ Code::Pi1300Plus,
+ Code::Pi1300Minus,
+ Code::Pi1300_0,
+
+ Code::KStar0_1430_0,
+ Code::KStar0_1430_0Bar,
+ Code::KStar0_1430_Plus,
+ Code::KStar0_1430_MinusBar,
+ });
+ decaySibyll.PrintDecayConfig();
process::particle_cut::ParticleCut cut(5_GeV);
@@ -157,7 +182,8 @@ int main() {
auto sibyllSequence = sibyll << sibyllNuc;
process::switch_process::SwitchProcess switchProcess(urqmd, sibyllSequence, 55_GeV);
- auto sequence = switchProcess << decay << eLoss << cut << observationLevel
+ auto decaySequence = decayPythia << decaySibyll;
+ auto sequence = switchProcess << decaySequence << eLoss << cut << observationLevel
<< trackWriter;
// define air shower object, run simulation
diff --git a/Framework/Particles/ParticleData.xml b/Framework/Particles/ParticleData.xml
index 0102f3e261ceee58b0b94003ec940d37231dd10c..5de8af669fcfb347e1aed0c2f23f639111c6924c 100644
--- a/Framework/Particles/ParticleData.xml
+++ b/Framework/Particles/ParticleData.xml
@@ -715,7 +715,15 @@
<channel onMode="1" bRatio="0.9995502" meMode="2" products="211 111"/>
<channel onMode="1" bRatio="0.0004498" products="211 22"/>
</particle>
-
+
+<particle id="100211" name="pi(1300)+" antiName="pi(1300)-" spinType="1" chargeType="3" colType="0"
+ m0="1.300" mWidth="0.200">
+</particle>
+
+<particle id="100111" name="pi(1300)0" spinType="1" chargeType="3" colType="0"
+ m0="1.300" mWidth="0.200">
+</particle>
+
<!--
<particle id="215" name="a_2+" antiName="a_2-" spinType="5" chargeType="3" colType="0"
m0="1.31830" mWidth="0.10700" mMin="1.00000" mMax="1.70000">
@@ -793,7 +801,15 @@
<channel onMode="1" bRatio="0.3326633" meMode="2" products="311 111"/>
<channel onMode="1" bRatio="0.0023900" products="311 22"/>
</particle>
-
+
+<particle id="10311" name="K*0(1430)0" antiName="K*0(1430)0bar" spinType="1" chargeType="0" colType="0"
+ m0="1.425" mWidth="0.270" mMin="1.3750" mMax="1.4750">
+</particle>
+
+<particle id="10321" name="K*0(1430)+" antiName="K*0(1430)bar-" spinType="1" chargeType="0" colType="0"
+ m0="1.425" mWidth="0.270" mMin="1.3750" mMax="1.4750">
+</particle>
+
<!--
<particle id="315" name="K*_2(1430)0" antiName="K*_2(1430)bar0" spinType="5" chargeType="0" colType="0"
m0="1.43240" mWidth="0.10900" mMin="1.10000" mMax="1.80000">
@@ -849,7 +865,6 @@
</particle>
-->
-<!--
<particle id="331" name="eta'" spinType="1" chargeType="0" colType="0"
m0="0.95778" mWidth="0.00020" mMin="0.95578" mMax="0.95978">
<channel onMode="1" bRatio="0.4365815" products="211 -211 221"/>
@@ -860,7 +875,6 @@
<channel onMode="1" bRatio="0.0016900" products="111 111 111"/>
<channel onMode="1" bRatio="0.0001076" products="13 -13 22"/>
</particle>
--->
<particle id="333" name="phi" spinType="3" chargeType="0" colType="0"
m0="1.01946" mWidth="0.00426" mMin="1.00000" mMax="1.04000">
@@ -3360,12 +3374,14 @@
<particle id="1103" name="dd_1" antiName="dd_1bar" spinType="3" chargeType="-2" colType="-1"
m0="0.77133">
</particle>
-
+-->
+
<particle id="1114" name="Delta-" antiName="Deltabar+" spinType="4" chargeType="-3" colType="0"
m0="1.23200" mWidth="0.11700" mMin="1.08000" mMax="1.60000">
<channel onMode="1" bRatio="1.0000000" products="2112 -211"/>
</particle>
-
+
+<!--
<particle id="2101" name="ud_0" antiName="ud_0bar" spinType="1" chargeType="1" colType="-1"
m0="0.57933">
</particle>
@@ -3396,7 +3412,45 @@
<particle id="2212" name="p+" antiName="pbar-" spinType="2" chargeType="3" colType="0"
m0="0.93827">
</particle>
-
+
+<particle id="202212" name="N(1440)+" antiName="N(1440)bar-" spinType="2" chargeType="3" colType="0"
+ m0="1.440" mWidth="0.350" mMin="1.4100" mMax="1.4500">
+ <channel onMode="1" bRatio="0.466" products="2212 211"/>
+ <channel onMode="1" bRatio="0.2340" products="2112 111"/>
+ <channel onMode="1" bRatio="0.1990" products="2212 211 111"/>
+ <channel onMode="1" bRatio="0.1010" products="2112 211 -211"/>
+</particle>
+
+<particle id="202112" name="N(1440)0" antiName="N(1440)bar0" spinType="2" chargeType="3" colType="0"
+ m0="1.440" mWidth="0.350" mMin="1.4100" mMax="1.4500">
+ <channel onMode="1" bRatio="0.4660" products="2112 -211"/>
+ <channel onMode="1" bRatio="0.2340" products="2212 111"/>
+ <channel onMode="1" bRatio="0.1990" products="2112 111 -211"/>
+ <channel onMode="1" bRatio="0.1010" products="2212 211 -211"/>
+</particle>
+
+<particle id="212212" name="N(1710)+" antiName="N(1710)bar-" spinType="2" chargeType="3" colType="0"
+ m0="1.710" mWidth="0.100" mMin="1.6800" mMax="1.7400">
+ <channel onMode="1" bRatio="0.3334" products="2212 211 111"/>
+ <channel onMode="1" bRatio="0.1666" products="2112 211 -211"/>
+ <channel onMode="1" bRatio="0.1334" products="2212 211"/>
+ <channel onMode="1" bRatio="0.1300" products="2112 223"/>
+ <channel onMode="1" bRatio="0.1100" products="3122 321"/>
+ <channel onMode="1" bRatio="0.0660" products="2112 111"/>
+ <channel onMode="1" bRatio="0.0660" products="2112 221"/>
+</particle>
+
+<particle id="212112" name="N(1710)0" antiName="N(1710)bar0" spinType="2" chargeType="3" colType="0"
+ m0="1.710" mWidth="0.100" mMin="1.6800" mMax="1.7400">
+ <channel onMode="1" bRatio="0.3334" products="2112 -211 111"/>
+ <channel onMode="1" bRatio="0.1666" products="2212 211 -211"/>
+ <channel onMode="1" bRatio="0.1334" products="2112 -211"/>
+ <channel onMode="1" bRatio="0.1300" products="2212 223"/>
+ <channel onMode="1" bRatio="0.1100" products="3122 311"/>
+ <channel onMode="1" bRatio="0.0660" products="2212 111"/>
+ <channel onMode="1" bRatio="0.0660" products="2212 221"/>
+</particle>
+
<particle id="2214" name="Delta+" antiName="Deltabar-" spinType="4" chargeType="3" colType="0"
m0="1.23200" mWidth="0.11700" mMin="1.08000" mMax="1.60000">
<channel onMode="1" bRatio="0.6630208" products="2212 111"/>
@@ -3429,14 +3483,14 @@
<channel onMode="1" bRatio="0.0000573" meMode="22" products="-12 11 3122"/>
</particle>
-<!--
+
<particle id="3114" name="Sigma*-" antiName="Sigma*bar+" spinType="4" chargeType="-3" colType="0"
m0="1.38720" mWidth="0.03940" mMin="1.34000" mMax="1.70000">
<channel onMode="1" bRatio="0.8814590" products="3122 -211"/>
<channel onMode="1" bRatio="0.0592705" products="3212 -211"/>
<channel onMode="1" bRatio="0.0592705" products="3112 111"/>
</particle>
--->
+
<particle id="3122" name="Lambda0" antiName="Lambdabar0" spinType="2" chargeType="0" colType="0"
m0="1.11568" tau0="7.89000e+01">
@@ -3483,7 +3537,7 @@
<channel onMode="1" bRatio="0.0049751" meMode="11" products="3122 11 -11"/>
</particle>
-<!--
+
<particle id="3214" name="Sigma*0" antiName="Sigma*bar0" spinType="4" chargeType="0" colType="0"
m0="1.38370" mWidth="0.03600" mMin="1.34000" mMax="1.70000">
<channel onMode="1" bRatio="0.8700000" products="3122 111"/>
@@ -3491,7 +3545,7 @@
<channel onMode="1" bRatio="0.0585000" products="3112 211"/>
<channel onMode="1" bRatio="0.0130000" products="3122 22"/>
</particle>
--->
+
<particle id="3222" name="Sigma+" antiName="Sigmabar-" spinType="2" chargeType="3" colType="0"
m0="1.18937" tau0="2.40400e+01">
@@ -3501,14 +3555,14 @@
<channel onMode="1" bRatio="0.0000200" meMode="22" products="-11 12 3122"/>
</particle>
-<!--
+
<particle id="3224" name="Sigma*+" antiName="Sigma*bar-" spinType="4" chargeType="3" colType="0"
m0="1.38280" mWidth="0.03600" mMin="1.34000" mMax="1.70000">
<channel onMode="1" bRatio="0.8814590" products="3122 211"/>
<channel onMode="1" bRatio="0.0592705" products="3222 111"/>
<channel onMode="1" bRatio="0.0592705" products="3212 211"/>
</particle>
--->
+
<!--
<particle id="3303" name="ss_1" antiName="ss_1bar" spinType="3" chargeType="-2" colType="-1"
@@ -3525,13 +3579,12 @@
<channel onMode="1" bRatio="0.0000870" meMode="22" products="-12 11 3212"/>
</particle>
-<!--
+
<particle id="3314" name="Xi*-" antiName="Xi*bar+" spinType="4" chargeType="-3" colType="0"
m0="1.53500" mWidth="0.00990" mMin="1.46000" mMax="1.63000">
<channel onMode="1" bRatio="0.6670000" products="3322 -211"/>
<channel onMode="1" bRatio="0.3330000" products="3312 111"/>
</particle>
--->
<particle id="3322" name="Xi0" antiName="Xibar0" spinType="2" chargeType="0" colType="0"
m0="1.31486" tau0="8.71000e+01">
@@ -3542,13 +3595,11 @@
<channel onMode="1" bRatio="0.0000046" meMode="22" products="-14 13 3222"/>
</particle>
-<!--
<particle id="3324" name="Xi*0" antiName="Xi*bar0" spinType="4" chargeType="0" colType="0"
m0="1.53180" mWidth="0.00910" mMin="1.46000" mMax="1.63000">
<channel onMode="1" bRatio="0.3330000" products="3322 111"/>
<channel onMode="1" bRatio="0.6670000" products="3312 211"/>
</particle>
--->
<particle id="3334" name="Omega-" antiName="Omegabar+" spinType="4" chargeType="-3" colType="0"
m0="1.67245" tau0="2.46100e+01">
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index eed8d5e9e0e9d345bb572b945721d540f69fa112..7586f0b42c63e88f532ba6bba17f27348a09ae9c 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -28,37 +28,128 @@ using Track = Trajectory;
namespace corsika::process::pythia {
- Decay::Decay(vector<particles::Code> pParticles)
- : fTrackedParticles(pParticles) {}
+ Decay::Decay() {}
+ Decay::Decay(std::set<particles::Code> vHandled)
+ : handleAllDecays_(false)
+ , handledDecays_(vHandled) {}
Decay::~Decay() { cout << "Pythia::Decay n=" << fCount << endl; }
void Decay::Init() {
- Decay::SetParticleListStable(fTrackedParticles);
-
// set random number generator in pythia
Pythia8::RndmEngine* rndm = new corsika::process::pythia::Random();
fPythia.setRndmEnginePtr(rndm);
+ /*
+ issue xyz: definition of particles and decay channels use the same mechanism in
+ corsika and pythia we should force pythia to use the file in corsika.
+ */
+ // bool ParticleData::reInit(string startFile, bool xmlFormat = true)
+ // read in particle data from Corsika 8
+ // fPythia.particleData.reInit("/home/felix/ngcorsika/corsika-build/include/corsika/particles/ParticleData.xml");
+ // fPythia.particleData.checkTable();
+
fPythia.readString("Next:numberShowInfo = 0");
fPythia.readString("Next:numberShowProcess = 0");
fPythia.readString("Next:numberShowEvent = 0");
fPythia.readString("Print:quiet = on");
+ fPythia.readString("Check:particleData = 1");
+
+ /*
+ switching off event check in pythia is needed to allow decays that are off-shell
+ according to the mass definition in pythia.
+ the consistency of particle masses between event generators is an unsolved issues
+ */
+ cout << "Pythia::Init: switching off event checking in pythia.." << endl;
+ fPythia.readString("Check:event = 0");
fPythia.readString("ProcessLevel:all = off");
fPythia.readString("ProcessLevel:resonanceDecays = off");
- fPythia.particleData.readString("59:m0 = 101.00");
+ // making sure
+ SetStable(particles::Code::Pi0);
+
+ // fPythia.particleData.readString("59:m0 = 101.00");
- fPythia.init();
+ if (!fPythia.init())
+ throw std::runtime_error("Pythia::Decay: Initialization failed!");
+ }
+
+ bool Decay::CanHandleDecay(const particles::Code vParticleCode) {
+ using namespace corsika::particles;
+ // if known to pythia and not proton, electron or neutrino it can decay
+ if (vParticleCode == Code::Proton || vParticleCode == Code::AntiProton ||
+ vParticleCode == Code::NuE || vParticleCode == Code::NuMu ||
+ vParticleCode == Code::NuTau || vParticleCode == Code::NuEBar ||
+ vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
+ vParticleCode == Code::Electron || vParticleCode == Code::Positron)
+ return false;
+ else if (CanDecay(vParticleCode)) // non-zero for particles known to sibyll
+ return true;
+ else
+ return false;
+ }
+
+ void Decay::SetHandleDecay(const particles::Code vParticleCode) {
+ handleAllDecays_ = false;
+ cout << "Pythia::Decay: set to handle decay of " << vParticleCode << endl;
+ if (Decay::CanHandleDecay(vParticleCode))
+ handledDecays_.insert(vParticleCode);
+ else
+ throw std::runtime_error("this decay can not be handled by pythia!");
+ }
+
+ void Decay::SetHandleDecay(const vector<particles::Code> vParticleList) {
+ handleAllDecays_ = false;
+ for (auto p : vParticleList) Decay::SetHandleDecay(p);
+ }
+
+ bool Decay::IsDecayHandled(const corsika::particles::Code vParticleCode) {
+ if (handleAllDecays_ && Decay::CanHandleDecay(vParticleCode))
+ return true;
+ else
+ return Decay::handledDecays_.find(vParticleCode) != Decay::handledDecays_.end()
+ ? true
+ : false;
}
- void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
+ bool Decay::IsStable(const particles::Code vCode) {
+ return fPythia.particleData.canDecay(static_cast<int>(particles::GetPDG(vCode)));
+ }
+
+ void Decay::PrintDecayConfig(const particles::Code vCode) {
+ cout << "Decay: Pythia decay configuration:" << endl;
+ cout << vCode << " is ";
+ if (IsStable(vCode))
+ cout << "stable" << endl;
+ else
+ cout << "unstable" << endl;
+ }
+
+ void Decay::PrintDecayConfig() {
+ cout << "Pythia::Decay: decay configuration:" << endl;
+ if (handleAllDecays_)
+ cout << " all particles known to Pythia are handled by Pythia::Decay!" << endl;
+ else
+ for (auto& pCode : handledDecays_)
+ cout << "Decay of " << pCode << " is handled by Pythia!" << endl;
+ }
+
+ void Decay::SetStable(const vector<particles::Code> particleList) {
for (auto p : particleList) Decay::SetStable(p);
}
+ bool Decay::CanDecay(const particles::Code pCode) {
+ std::cout << "Pythia::Decay: checking if particle: " << pCode
+ << " can decay in PYTHIA? ";
+ const bool ans =
+ fPythia.particleData.canDecay(static_cast<int>(particles::GetPDG(pCode)));
+ std::cout << ans << std::endl;
+ return ans;
+ }
+
void Decay::SetUnstable(const particles::Code pCode) {
cout << "Pythia::Decay: setting " << pCode << " unstable.." << endl;
fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
@@ -70,18 +161,29 @@ namespace corsika::process::pythia {
}
template <>
- units::si::TimeType Decay::GetLifetime(Particle const& p) {
+ units::si::TimeType Decay::GetLifetime(Particle const& vP) {
using namespace units::si;
- HEPEnergyType E = p.GetEnergy();
- HEPMassType m = p.GetMass();
-
- const double gamma = E / m;
-
- const TimeType t0 = particles::GetLifetime(p.GetPID());
- auto const lifetime = gamma * t0;
-
- return lifetime;
+ const auto pid = vP.GetPID();
+ if (CanDecay(pid)) {
+ HEPEnergyType E = vP.GetEnergy();
+ HEPMassType m = vP.GetMass();
+
+ const double gamma = E / m;
+
+ const TimeType t0 = particles::GetLifetime(pid);
+ auto const lifetime = gamma * t0;
+ cout << "Pythia::Decay: code: " << vP.GetPID() << endl;
+ cout << "Pythia::Decay: MinStep: t0: " << t0 << endl;
+ cout << "Pythia::Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
+ cout << "Pythia::Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV
+ << " GeV" << endl;
+ cout << "Pythia::Decay: MinStep: gamma: " << gamma << endl;
+ cout << "Pythia::Decay: MinStep: tau: " << lifetime << endl;
+
+ return lifetime;
+ } else
+ return std::numeric_limits<double>::infinity() * 1_s;
}
template <>
@@ -103,24 +205,26 @@ namespace corsika::process::pythia {
Pythia8::Event& event = fPythia.event;
event.reset();
+ auto const particleId = vP.GetPID();
+
// set particle unstable
- Decay::SetUnstable(vP.GetPID());
+ Decay::SetUnstable(particleId);
// input particle PDG
- auto const pdgCode = static_cast<int>(particles::GetPDG(vP.GetPID()));
+ auto const pdgCode = static_cast<int>(particles::GetPDG(particleId));
auto const pcomp = vP.GetMomentum().GetComponents();
double px = pcomp[0] / 1_GeV;
double py = pcomp[1] / 1_GeV;
double pz = pcomp[2] / 1_GeV;
double en = vP.GetEnergy() / 1_GeV;
- double m = particles::GetMass(vP.GetPID()) / 1_GeV;
+ double m = particles::GetMass(particleId) / 1_GeV;
// add particle to pythia stack
event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
if (!fPythia.next())
- cout << "Pythia::Decay: decay failed!" << endl;
+ throw std::runtime_error("Pythia::Decay: decay failed!");
else
cout << "Pythia::Decay: particles after decay: " << event.size() << endl;
@@ -146,7 +250,7 @@ namespace corsika::process::pythia {
}
// set particle stable
- Decay::SetStable(vP.GetPID());
+ Decay::SetStable(particleId);
}
} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Decay.h b/Processes/Pythia/Decay.h
index cc971c6b0f169dbfc3141905d2685e97f539acbb..0bce429c4a630d7e9b6c98155c93ad84765417c2 100644
--- a/Processes/Pythia/Decay.h
+++ b/Processes/Pythia/Decay.h
@@ -22,26 +22,59 @@ namespace corsika::process {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
class Decay : public corsika::process::DecayProcess<Decay> {
- const std::vector<particles::Code> fTrackedParticles;
int fCount = 0;
+ bool handleAllDecays_ = true;
public:
- Decay(std::vector<corsika::particles::Code>);
+ Decay();
+ Decay(std::set<particles::Code>);
~Decay();
void Init();
- void SetParticleListStable(const std::vector<particles::Code>);
- void SetUnstable(const corsika::particles::Code);
- void SetStable(const corsika::particles::Code);
+ // is Pythia::Decay set to handle the decay of this particle?
+ bool IsDecayHandled(const corsika::particles::Code);
+
+ // is decay possible in principle?
+ bool CanHandleDecay(const corsika::particles::Code);
+
+ // set Pythia::Decay to handle the decay of this particle!
+ void SetHandleDecay(const corsika::particles::Code);
+ // set Pythia::Decay to handle the decay of this list of particles!
+ void SetHandleDecay(const std::vector<particles::Code>);
+ // set Pythia::Decay to handle all particle decays
+ void SetHandleAllDecays();
+
+ // print internal configuration for this particle
+ void PrintDecayConfig(const corsika::particles::Code);
+ // print configuration of decays in corsika
+ void PrintDecayConfig();
+
+ bool CanDecay(const corsika::particles::Code);
+
+ /**
+ In this function PYTHIA is asked for the lifetime of the input particle.
+ Unknown particles should return an infinite lifetime so that another decay process
+ can act on the particle later on.
+ */
template <typename TParticle>
corsika::units::si::TimeType GetLifetime(TParticle const&);
+ /**
+ In this function PYTHIA is called to execute the decay of the input particle.
+ */
+
template <typename TProjectile>
void DoDecay(TProjectile&);
private:
+ void SetUnstable(const corsika::particles::Code);
+ void SetStable(const corsika::particles::Code);
+ void SetStable(const std::vector<particles::Code>);
+ bool IsStable(const corsika::particles::Code);
+
Pythia8::Pythia fPythia;
+ std::set<particles::Code> handledDecays_;
};
} // namespace pythia
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index f229b907f92968d7c26ff7bed70b4476ce56282d..5dd0cc8bc6c60e23697d722c8d983ce0f44efc89 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -63,13 +63,9 @@ TEST_CASE("Pythia", "[processes]") {
SECTION("pythia interface") {
using namespace corsika;
- const std::vector<particles::Code> particleList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- process::pythia::Decay model(particleList);
+ process::pythia::Decay model;
model.Init();
}
@@ -126,19 +122,40 @@ TEST_CASE("pythia process") {
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::PiPlus, E0, plab, pos, 0_ns});
- const std::vector<particles::Code> particleList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
- process::pythia::Decay model(particleList);
+ process::pythia::Decay model;
model.Init();
- model.DoDecay(projectile);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
+ model.DoDecay(projectile);
+ REQUIRE(stack.GetSize() == 3);
+ }
+
+ SECTION("pythia decay config") {
+ process::pythia::Decay model({particles::Code::PiPlus, particles::Code::PiMinus});
+ REQUIRE(model.IsDecayHandled(particles::Code::PiPlus));
+ REQUIRE(model.IsDecayHandled(particles::Code::PiMinus));
+ REQUIRE_FALSE(model.IsDecayHandled(particles::Code::KPlus));
+
+ const std::vector<particles::Code> particleTestList = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::Lambda0Bar, particles::Code::D0Bar};
+
+ // setup decays
+ model.SetHandleDecay(particleTestList);
+ for (auto& pCode : particleTestList) REQUIRE(model.IsDecayHandled(pCode));
+
+ // individually
+ model.SetHandleDecay(particles::Code::KMinus);
+
+ // possible decays
+ REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Proton));
+ REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Electron));
+ REQUIRE(model.CanHandleDecay(particles::Code::PiPlus));
+ REQUIRE(model.CanHandleDecay(particles::Code::MuPlus));
}
SECTION("pythia interaction") {
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index ea50d4477715e7d765abb000087fb56ef6f3e722..7750bfd383febc8534b95b07ebffbb54bb263798 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -30,11 +30,59 @@ using SetupParticle = corsika::setup::Stack::ParticleType;
namespace corsika::process::sibyll {
- Decay::Decay() {}
- Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
- void Decay::Init() {
+ Decay::Decay() {
// switch off decays to avoid internal decay chains
SetAllStable();
+ // handle all decays by default
+ handleAllDecays_ = true;
+ }
+
+ Decay::Decay(std::set<particles::Code> vHandled)
+ : handleAllDecays_(false)
+ , handledDecays_(vHandled) {
+ SetAllStable();
+ }
+
+ Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
+ void Decay::Init() {}
+
+ bool Decay::CanHandleDecay(const particles::Code vParticleCode) {
+ using namespace corsika::particles;
+ // if known to sibyll and not proton or neutrino it can decay
+ if (vParticleCode == Code::Proton || vParticleCode == Code::AntiProton ||
+ vParticleCode == Code::NuE || vParticleCode == Code::NuMu ||
+ vParticleCode == Code::NuTau || vParticleCode == Code::NuEBar ||
+ vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
+ vParticleCode == Code::Electron || vParticleCode == Code::Positron)
+ return false;
+ else if (process::sibyll::ConvertToSibyllRaw(
+ vParticleCode)) // non-zero for particles known to sibyll
+ return true;
+ else
+ return false;
+ }
+
+ void Decay::SetHandleDecay(const particles::Code vParticleCode) {
+ handleAllDecays_ = false;
+ cout << "Sibyll::Decay: set to handle decay of " << vParticleCode << endl;
+ if (Decay::CanHandleDecay(vParticleCode))
+ handledDecays_.insert(vParticleCode);
+ else
+ throw std::runtime_error("this decay can not be handled by sibyll!");
+ }
+
+ void Decay::SetHandleDecay(const vector<particles::Code> vParticleList) {
+ handleAllDecays_ = false;
+ for (auto p : vParticleList) Decay::SetHandleDecay(p);
+ }
+
+ bool Decay::IsDecayHandled(const corsika::particles::Code vParticleCode) {
+ if (handleAllDecays_ && Decay::CanHandleDecay(vParticleCode))
+ return true;
+ else
+ return Decay::handledDecays_.find(vParticleCode) != Decay::handledDecays_.end()
+ ? true
+ : false;
}
void Decay::SetStable(const vector<particles::Code> vParticleList) {
@@ -58,13 +106,13 @@ namespace corsika::process::sibyll {
}
void Decay::SetUnstable(const particles::Code vCode) {
- cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << endl;
+ cout << "Sibyll::Decay: 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 vCode) {
- cout << "Sibyll::Interaction: setting " << vCode << " stable.." << endl;
+ cout << "Sibyll::Decay: 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]);
}
@@ -88,33 +136,46 @@ namespace corsika::process::sibyll {
cout << "unstable" << endl;
}
+ void Decay::PrintDecayConfig() {
+ cout << "Sibyll::Decay: decay configuration:" << endl;
+ if (handleAllDecays_)
+ cout << " all particles known to Sibyll are handled by Sibyll::Decay!" << endl;
+ else
+ for (auto& pCode : handledDecays_)
+ cout << "Decay of " << pCode << " is handled by Sibyll!" << endl;
+ }
+
template <>
- units::si::TimeType Decay::GetLifetime(SetupParticle const& vP) const {
+ units::si::TimeType Decay::GetLifetime(SetupParticle const& vP) {
using namespace units::si;
- HEPEnergyType E = vP.GetEnergy();
- HEPMassType m = vP.GetMass();
-
- const double gamma = E / m;
-
- const TimeType t0 = particles::GetLifetime(vP.GetPID());
- auto const lifetime = gamma * t0;
-
- const auto mkin =
- (E * E - vP.GetMomentum().squaredNorm()); // delta_mass(vP.GetMomentum(), E, m);
- cout << "Decay: code: " << vP.GetPID() << endl;
- cout << "Decay: MinStep: t0: " << t0 << endl;
- cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
- cout << "Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV << " GeV"
- << endl;
- cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV << " "
- << m / 1_GeV * m / 1_GeV << endl;
- auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
- cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
- cout << "Decay: MinStep: gamma: " << gamma << endl;
- cout << "Decay: MinStep: tau: " << lifetime << endl;
-
- return lifetime;
+ const particles::Code pid = vP.GetPID();
+ if (Decay::IsDecayHandled(pid)) {
+ HEPEnergyType E = vP.GetEnergy();
+ HEPMassType m = vP.GetMass();
+
+ const double gamma = E / m;
+
+ const TimeType t0 = particles::GetLifetime(vP.GetPID());
+ auto const lifetime = gamma * t0;
+
+ const auto mkin =
+ (E * E - vP.GetMomentum().squaredNorm()); // delta_mass(vP.GetMomentum(), E, m);
+ cout << "Sibyll::Decay: code: " << vP.GetPID() << endl;
+ cout << "Sibyll::Decay: MinStep: t0: " << t0 << endl;
+ cout << "Sibyll::Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
+ cout << "Sibyll::Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV
+ << " GeV" << endl;
+ cout << "Sibyll::Decay: momentum: shell mass-kin. inv. mass "
+ << mkin / 1_GeV / 1_GeV << " " << m / 1_GeV * m / 1_GeV << endl;
+ auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
+ cout << "Sibyll::Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
+ cout << "Sibyll::Decay: MinStep: gamma: " << gamma << endl;
+ cout << "Sibyll::Decay: MinStep: tau (s): " << lifetime / 1_s << endl;
+
+ return lifetime;
+ } else
+ return std::numeric_limits<double>::infinity() * 1_s;
}
template <>
@@ -122,10 +183,15 @@ namespace corsika::process::sibyll {
using geometry::Point;
using namespace units::si;
+ const particles::Code pCode = vP.GetPID();
+ // check if sibyll is configured to handle this decay!
+ if (!IsDecayHandled(pCode))
+ throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
+
fCount++;
SibStack ss;
ss.Clear();
- const particles::Code pCode = vP.GetPID();
+
// copy particle to sibyll stack
ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), vP.GetEnergy(),
vP.GetMomentum(),
@@ -145,12 +211,13 @@ namespace corsika::process::sibyll {
cout << "Decay: calling Sibyll decay routine.." << endl;
decsib_();
- // reset to stable
- SetStable(pCode);
// print output
int print_unit = 6;
sib_list_(print_unit);
+ // reset to stable
+ SetStable(pCode);
+
// copy particles from sibyll stack to corsika
for (auto& psib : ss) {
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 05bc0186ee6b2b4231d540cf8a1ccb209f65994d..96743485308fbfb68e5667cc5b554320749955de 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -15,6 +15,7 @@
#include <corsika/process/DecayProcess.h>
#include <corsika/process/SecondariesProcess.h>
+#include <set>
#include <vector>
namespace corsika::process {
@@ -23,27 +24,33 @@ namespace corsika::process {
class Decay : public corsika::process::DecayProcess<Decay> {
int fCount = 0;
+ bool handleAllDecays_ = true;
public:
Decay();
+ Decay(std::set<particles::Code>);
~Decay();
void Init();
- 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 PrintDecayConfig();
void SetHadronsUnstable();
+ // is Sibyll::Decay set to handle the decay of this particle?
+ bool IsDecayHandled(const corsika::particles::Code);
+
+ // is decay possible in principle?
+ bool CanHandleDecay(const corsika::particles::Code);
+
+ // set Sibyll::Decay to handle the decay of this particle!
+ void SetHandleDecay(const corsika::particles::Code);
+ // set Sibyll::Decay to handle the decay of this list of particles!
+ void SetHandleDecay(const std::vector<particles::Code>);
+ // set Sibyll::Decay to handle all particle decays
+ void SetHandleAllDecay();
+
template <typename TParticle>
- corsika::units::si::TimeType GetLifetime(TParticle const&) const;
+ corsika::units::si::TimeType GetLifetime(TParticle const&);
/**
In this function SIBYLL is called to produce to decay the input particle.
@@ -54,6 +61,28 @@ namespace corsika::process {
template <typename TParticleView>
EProcessReturn DoSecondaries(TParticleView&);
+
+ private:
+ // internal routines to set particles stable and unstable in the COMMON blocks in
+ // sibyll
+ void SetStable(const std::vector<particles::Code>);
+ void SetUnstable(const std::vector<particles::Code>);
+
+ void SetStable(const corsika::particles::Code);
+ void SetUnstable(const corsika::particles::Code);
+
+ // internally set all particles to decay/not to decay
+ void SetAllUnstable();
+ void SetAllStable();
+
+ // will this particle be stable in sibyll ?
+ bool IsStable(const corsika::particles::Code);
+ // will this particle decay in sibyll ?
+ bool IsUnstable(const corsika::particles::Code);
+ // set particle with input code to decay or not
+ void SetDecay(const particles::Code, const bool);
+
+ std::set<particles::Code> handledDecays_;
};
} // namespace sibyll
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 553792e8f945ed6fb414c55f0d898650012da9cf..9f0d210ca8db8a9d4875e7a2b42c647148f8c7e6 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -51,34 +51,10 @@ namespace corsika::process::sibyll {
}
}
- 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 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 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,
@@ -317,16 +293,7 @@ namespace corsika::process::sibyll {
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
- if (internalDecays_) {
- // 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(trackedParticles_);
- decsib_();
- // reset
- SetAllStable();
- }
+
// print final state
int print_unit = 6;
sib_list_(print_unit);
@@ -340,8 +307,9 @@ namespace corsika::process::sibyll {
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
- // skip particles that have decayed in Sibyll
- if (psib.HasDecayed()) continue;
+ // abort on particles that have decayed in Sibyll. Should not happen!
+ if (psib.HasDecayed())
+ throw std::runtime_error("found particle that decayed in SIBYLL!");
// transform energy to lab. frame
auto const pCoM = psib.GetMomentum();
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 62302d09190c1044c65ccf52289de0915645dc9c..ebf4ac22aa465f3bb405e3828d6b647ef0a8bff3 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -31,12 +31,6 @@ namespace corsika::process::sibyll {
void Init();
- 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 initialized_; }
@@ -69,20 +63,7 @@ namespace corsika::process::sibyll {
private:
corsika::random::RNG& RNG_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
- // FOR NOW keep trackedParticles private, could be configurable
- std::vector<particles::Code> const trackedParticles_ = {
- 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::Sigma0, particles::Code::Sigma0Bar,
- 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 internalDecays_ = true;
+
const corsika::units::si::HEPEnergyType minEnergyCoM_ =
10. * 1e9 * corsika::units::si::electronvolt;
const corsika::units::si::HEPEnergyType maxEnergyCoM_ =
diff --git a/Processes/Sibyll/sibyll_codes.dat b/Processes/Sibyll/sibyll_codes.dat
index 1bc87904bf53354323ce514761e22a868d0111b7..e24c95495b88838c02806fe1d7ad6c52bd3a6ecf 100644
--- a/Processes/Sibyll/sibyll_codes.dat
+++ b/Processes/Sibyll/sibyll_codes.dat
@@ -24,6 +24,10 @@ PiMinus 8 1 2
RhoPlus 25 0 0
RhoMinus 26 0 0
Eta 23 0 0
+EtaPrime 24 0 0
+Pi1300Plus 62 0 0
+Pi1300Minus 63 0 0
+Pi1300_0 61 0 0
Omega 32 0 0
K0Short 12 1 3
KStar0 30 0 0
@@ -32,6 +36,10 @@ KPlus 9 1 3
KMinus 10 1 3
KStarPlus 28 0 0
KStarMinus 29 0 0
+KStar0_1430_0 66 0 0
+KStar0_1430_0Bar 67 0 0
+KStar0_1430_Plus 64 0 0
+KStar0_1430_MinusBar 65 0 0
DPlus 59 1 3
DMinus 60 1 3
DStarPlus 78 0 0
@@ -49,14 +57,30 @@ Neutron 14 1 1
AntiNeutron -14 1 1
Delta0 42 0 0
Delta0Bar -42 0 0
+DeltaMinus 43 0 0
+DeltaPlusBar -43 0 0
Proton 13 1 1
AntiProton -13 1 1
+N1440Plus 51 0 0
+N1440MinusBar -51 0 0
+N1440_0 52 0 0
+N1440_0Bar -52 0 0
+N1710Plus 53 0 0
+N1710MinusBar -53 0 0
+N1710_0 54 0 0
+N1710_0Bar -54 0 0
DeltaPlus 41 0 0
DeltaMinusBar -41 0 0
DeltaPlusPlus 40 0 0
DeltaMinusMinusBar -40 0 0
SigmaMinus 36 1 1
SigmaPlusBar -36 1 1
+SigmaStarMinus 46 0 0
+SigmaStarPlusBar -46 0 0
+SigmaStarPlus 44 0 0
+SigmaStarMinusBar -44 0 0
+SigmaStar0 45 0 0
+SigmaStar0Bar -45 0 0
Lambda0 39 1 1
Lambda0Bar -39 1 1
Sigma0 35 1 1
@@ -67,6 +91,10 @@ XiMinus 38 1 1
XiPlusBar -38 1 1
Xi0 37 1 1
Xi0Bar -37 1 1
+XiStarMinus 48 0 0
+XiStarPlusBar -48 0 0
+XiStar0 47 0 0
+XiStar0Bar -47 0 0
OmegaMinus 49 0 0
OmegaPlusBar -49 0 0
SigmaC0 86 0 0
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 492b6bf3ed455c8ee3c7c599ab16b98cccec0030..116a77b26a2bb1266341b2f58281fa4da6bfc7e5 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -167,28 +167,41 @@ TEST_CASE("SibyllInterface", "[processes]") {
Decay model;
+ model.PrintDecayConfig();
+
model.Init();
+
+ [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
+
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
+
// run checks
- [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
+ // lambda decays into proton and pi- or neutron and pi+
+ REQUIRE(stack.GetSize() == 3);
}
SECTION("DecayConfiguration") {
- Decay model;
+ Decay model({particles::Code::PiPlus, particles::Code::PiMinus});
+ REQUIRE(model.IsDecayHandled(particles::Code::PiPlus));
+ REQUIRE(model.IsDecayHandled(particles::Code::PiMinus));
+ REQUIRE_FALSE(model.IsDecayHandled(particles::Code::KPlus));
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]);
- }
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::Lambda0Bar, particles::Code::D0Bar};
+
+ // setup decays
+ model.SetHandleDecay(particleTestList);
+ for (auto& pCode : particleTestList) REQUIRE(model.IsDecayHandled(pCode));
+
+ // individually
+ model.SetHandleDecay(particles::Code::KMinus);
+
+ // possible decays
+ REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Proton));
+ REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Electron));
+ REQUIRE(model.CanHandleDecay(particles::Code::PiPlus));
+ REQUIRE(model.CanHandleDecay(particles::Code::MuPlus));
}
}