diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index e221b28fbc3b5d977cd6ae8ac193708e6e149b7a..cf6cd03b5b59b6d130388422b73385c7b2476a5a 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -119,8 +119,10 @@ public:
template <typename Particle>
LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
+ cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
+ cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
const Code pid = p.GetPID();
- if (isEmParticle(pid) || isInvisible(pid)) {
+ if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
return 0_m;
} else {
@@ -134,40 +136,26 @@ public:
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
const Code pid = p.GetPID();
+ EProcessReturn ret = EProcessReturn::eOk;
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += p.GetEnergy();
fEmCount += 1;
p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += p.GetEnergy();
fInvCount += 1;
p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy();
p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
}
- // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
- // cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
- // cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
- // const Code pid = p.GetPID();
- // if( isEmParticle( pid ) ){
- // cout << "removing em. particle..." << endl;
- // fEmEnergy += p.GetEnergy();
- // fEmCount += 1;
- // p.Delete();
- // return EProcessReturn::eParticleAbsorbed;
- // }
- // if ( isInvisible( pid ) ){
- // cout << "removing inv. particle..." << endl;
- // fInvEnergy += p.GetEnergy();
- // fInvCount += 1;
- // p.Delete();
- // return EProcessReturn::eParticleAbsorbed;
- // }
- return EProcessReturn::eOk;
+ return ret;
}
void Init() {
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 23954fef70ffe62140074010a1207624b1cfcaa2..24fd807d2723343be6a827f0783fa8fd5fdd0383 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -42,13 +42,23 @@ namespace corsika::process {
}
}
+ void setUnstable(const corsika::particles::Code pCode) const {
+ int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
+ }
+
+ void setStable(const corsika::particles::Code pCode) const {
+ int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ }
+
void setAllStable() const {
+ // name? also makes EM particles stable
- // using namespace corsika::io;
using std::cout;
using std::endl;
- // name? also makes EM particles stable
+ std::cout << "Decay: setting all particles stable.." << std::endl;
// loop over all particles in sibyll
// should be changed to loop over human readable list
@@ -58,11 +68,7 @@ namespace corsika::process {
const int sibCode = (int)p;
// skip unknown and antiparticles
if (sibCode < 1) continue;
- const corsika::particles::Code pid =
- ConvertFromSibyll(static_cast<SibyllCode>(sibCode));
- std::cout << "Sibyll: Decay: setting " << pid << " stable" << std::endl;
s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
- std::cout << "decay table value: " << s_csydec_.idb[sibCode - 1] << std::endl;
}
}
@@ -124,6 +130,7 @@ namespace corsika::process {
corsika::particles::GetLifetime(p.GetPID());
cout << "Decay: code: " << p.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
+ cout << "Decay: MinStep: energy: " << E / 1_GeV << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
// cout << "Decay: MinStep: density: " << density << endl;
// return as column density
@@ -139,20 +146,29 @@ namespace corsika::process {
template <typename Particle, typename Stack>
void DoDecay(Particle& p, Stack& s) const {
+ using corsika::geometry::Point;
+ using namespace corsika::units::si;
SibStack ss;
ss.Clear();
// copy particle to sibyll stack
auto pin = ss.NewParticle();
- pin.SetPID(process::sibyll::ConvertToSibyllRaw(p.GetPID()));
+ const corsika::particles::Code pCode = p.GetPID();
+ pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
pin.SetEnergy(p.GetEnergy());
pin.SetMomentum(p.GetMomentum());
+ // remember position
+ Point decayPoint = p.GetPosition();
+ TimeType t0 = p.GetTime();
// remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable
- setHadronsUnstable();
+ // setHadronsUnstable();
+ setUnstable(pCode);
// call sibyll decay
std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
decsib_();
+ // reset to stable
+ setStable(pCode);
// print output
int print_unit = 6;
sib_list_(print_unit);
@@ -169,6 +185,8 @@ namespace corsika::process {
pnew.SetEnergy(psib.GetEnergy());
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
pnew.SetMomentum(psib.GetMomentum());
+ pnew.SetPosition(decayPoint);
+ pnew.SetTime(t0);
}
// empty sibyll stack
ss.Clear();
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index a82dac38034c26a6373fe6ff07a4a278dcdb400e..10e45a0c3c1036acf0c91dd7d25800ab1ae334ed 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -56,9 +56,8 @@ namespace corsika::process::sibyll {
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
- int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
-
- bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
+ const int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
+ const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/*
the target should be defined by the Environment,
@@ -67,7 +66,7 @@ namespace corsika::process::sibyll {
*/
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
- int kTarget = 16;
+ const int kTarget = corsika::particles::Oxygen::GetNucleusA();
hep::EnergyType Etot =
p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
@@ -77,8 +76,8 @@ namespace corsika::process::sibyll {
Ptot += p.GetMomentum();
Ptot += pTarget;
// calculate cm. energy
- hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
- double Ecm = sqs / 1_GeV;
+ const hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
+ const double Ecm = sqs / 1_GeV;
std::cout << "Interaction: "
<< "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
@@ -149,8 +148,8 @@ namespace corsika::process::sibyll {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
- Point pOrig(rootCS, coordinates);
+ Point pOrig = p.GetPosition();
+ TimeType tOrig = p.GetTime();
/*
the target should be defined by the Environment,
@@ -160,8 +159,9 @@ namespace corsika::process::sibyll {
here we need: GetTargetMassNumber() or GetTargetPID()??
GetTargetMomentum() (zero in EAS)
*/
- // FOR NOW: set target to proton
- int kTarget = 1; // env.GetTargetParticle().GetPID();
+ // FOR NOW: set target to oxygen
+ const int kTarget = corsika::particles::Oxygen::
+ GetNucleusA(); // env.GetTargetParticle().GetPID();
cout << "defining target momentum.." << endl;
// FOR NOW: target is always at rest
@@ -171,6 +171,8 @@ namespace corsika::process::sibyll {
<< pTarget.GetComponents() / 1_GeV * constants::c << endl;
cout << "beam momentum (GeV/c): "
<< p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
+ cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
+ cout << "time: " << tOrig << endl;
// get energy of particle from stack
/*
@@ -209,7 +211,7 @@ namespace corsika::process::sibyll {
if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
<< " DoInteraction: dropping particle.." << std::endl;
- p.Delete();
+ // p.Delete(); delete later... different process
} else {
// Sibyll does not know about units..
double sqs = Ecm / 1_GeV;
@@ -265,6 +267,8 @@ namespace corsika::process::sibyll {
corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]};
pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
+ pnew.SetPosition(pOrig);
+ pnew.SetTime(tOrig);
Ptot_final += pnew.GetMomentum();
}
// cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
index 5d2d1e69644e812246b731aa4bf276d090d79910..0a0ecc1133ce2fce469e93747368ab9224568449 100644
--- a/Processes/Sibyll/sibyll2.3c.h
+++ b/Processes/Sibyll/sibyll2.3c.h
@@ -59,7 +59,7 @@ extern struct {
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
// sibyll main subroutine
-void sibyll_(int&, int&, double&);
+void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
@@ -79,8 +79,9 @@ void decsib_();
// interaction length
// double fpni_(double&, int&);
-void sib_sigma_hnuc_(int&, int&, double&, double&, double&);
-void sib_sigma_hp_(int&, double&, double&, double&, double&, double*, double&, double&);
+void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&);
+void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
+ double&);
double s_rndm_(int&);