diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f439cf2ab3fd895457b49610720c46da1ddd4bd8..abd95beef9c1c5e7f816cb53bf03ef3161757b76 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -255,16 +255,15 @@ int main() {
stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll(env);
process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
process::sibyll::Decay decay(trackedHadrons);
- //random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- //process::pythia::Decay decay(trackedHadrons);
+ // random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+ // process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
@@ -275,9 +274,10 @@ int main() {
process::EnergyLoss::EnergyLoss eLoss;
// assemble all processes into an ordered process list
- // auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut << trackWriter;
- // auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut << trackWriter;
- //auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
+ // auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut <<
+ // trackWriter; auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss <<
+ // cut << trackWriter;
+ // auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
@@ -290,9 +290,7 @@ int main() {
const int nuclA = 4;
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
- const HEPEnergyType E0 =
- nuclA *
- 100_GeV;
+ const HEPEnergyType E0 = nuclA * 100_GeV;
double theta = 0.;
double phi = 0.;
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index 41b485adee0bdc1fac6dcfea79391d25d4c6b4a9..472baf40dfff9707fb52302043adcd10a4449c3b 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -78,7 +78,8 @@ namespace corsika::particles {
*/
int16_t constexpr GetChargeNumber(Code const p) {
if (p == Code::Nucleus)
- throw std::runtime_error("Cannot GetChargeNumber() of particle::Nucleus -> unspecified");
+ throw std::runtime_error(
+ "Cannot GetChargeNumber() of particle::Nucleus -> unspecified");
// electric_charges stores charges in units of (e/3), e.g. 3 for a proton
return detail::electric_charges[static_cast<CodeIntType>(p)] / 3;
}
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 94c686ff2d4186582afa2120d73549965df5c340..1e23f72062d449c27f0b9ac175e3dacfbb2958b0 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -37,15 +37,14 @@ namespace corsika::process::EnergyLoss {
EnergyLoss::EnergyLoss()
: fEnergyLossTot(0_GeV) {}
-
/**
* PDG2018, passage of particles through matter
*
- * Note, that \f$I_{\mathrm{eff}}\f$ of composite media a determined from \f$ \ln I = \sum_i
- * a_i \ln(I_i) \f$ where \f$ a_i \f$ is the fraction of the electron population
+ * Note, that \f$I_{\mathrm{eff}}\f$ of composite media a determined from \f$ \ln I =
+ * \sum_i a_i \ln(I_i) \f$ where \f$ a_i \f$ is the fraction of the electron population
* (\f$\sim Z_i\f$) of the \f$i\f$-th element. This can also be used for shell
* corrections or density effects.
- *
+ *
* The \f$I_{\mathrm{eff}}\f$ of compounds is not better than a few percent, if not
* measured explicitly.
*
@@ -57,9 +56,9 @@ namespace corsika::process::EnergyLoss {
// all these are material constants and have to come through Environment
// right now: values for nitrogen_D
// 7 nitrogen_gas 82.0 0.49976 D E 0.0011653 0.0 1.7378 4.1323 0.15349 3.2125 10.54
- auto Ieff = 82.0_eV;
- [[maybe_unused]] auto Zmat = 7;
- auto ZoverA = 0.49976_mol/1_g;
+ auto Ieff = 82.0_eV;
+ [[maybe_unused]] auto Zmat = 7;
+ auto ZoverA = 0.49976_mol / 1_g;
const double x0 = 1.7378;
const double x1 = 4.1323;
const double Cbar = 10.54;
@@ -68,9 +67,8 @@ namespace corsika::process::EnergyLoss {
const double sk = 3.2125;
// end of material constants
-
// this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
- auto const K = 0.307075_MeV / 1_mol * square(1_cm);
+ auto const K = 0.307075_MeV / 1_mol * square(1_cm);
HEPEnergyType const E = p.GetEnergy();
HEPMassType const m = p.GetMass();
double const gamma = E / m;
@@ -80,30 +78,34 @@ namespace corsika::process::EnergyLoss {
auto const m2 = m * m;
auto const me2 = me * me;
double const gamma2 = pow(gamma, 2);
-
- double const beta2 = (gamma2-1)/gamma2; // 1-1/gamma2 (1-1/gamma)*(1+1/gamma); (gamma_2-1)/gamma_2 = (1-1/gamma2);
- double const c2 = 1; // HEP convention here c=c2=1
+
+ double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2 (1-1/gamma)*(1+1/gamma);
+ // (gamma_2-1)/gamma_2 = (1-1/gamma2);
+ double const c2 = 1; // HEP convention here c=c2=1
cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
- [[maybe_unused]] double const eta2 = beta2/(1 - beta2);
- HEPMassType const Wmax = 2*me*c2*beta2*gamma2 / (1 + 2*gamma*me/m + me2/m2);
- // approx, but <<1% HEPMassType const Wmax = 2*me*c2*beta2*gamma2; for HEAVY PARTICLES
- // Wmax ~ 2me v2 for non-relativistic particles
+ [[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
+ HEPMassType const Wmax =
+ 2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
+ // approx, but <<1% HEPMassType const Wmax = 2*me*c2*beta2*gamma2; for HEAVY
+ // PARTICLES Wmax ~ 2me v2 for non-relativistic particles
cout << "BetheBloch Wmax=" << Wmax << endl;
-
+
// Sternheimer parameterization, density corrections towards high energies
- // NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization -> MISSING
- cout << "BetheBloch p.GetMomentum().GetNorm()/m=" << p.GetMomentum().GetNorm()/m << endl;
- double const x = log10(p.GetMomentum().GetNorm()/m);
+ // NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
+ // MISSING
+ cout << "BetheBloch p.GetMomentum().GetNorm()/m=" << p.GetMomentum().GetNorm() / m
+ << endl;
+ double const x = log10(p.GetMomentum().GetNorm() / m);
double delta = 0;
- if (x>=x1) {
- delta = 2*(log(10))*x - Cbar;
- } else if (x<x1 && x>=x0) {
- delta = 2*(log(10))*x - Cbar + aa*pow((x1-x), sk);
- } else if (x<x0) { // and IF conductor (otherwise, this is 0)
- delta = delta0*pow(100,2*(x-x0));
+ if (x >= x1) {
+ delta = 2 * (log(10)) * x - Cbar;
+ } else if (x < x1 && x >= x0) {
+ delta = 2 * (log(10)) * x - Cbar + aa * pow((x1 - x), sk);
+ } else if (x < x0) { // and IF conductor (otherwise, this is 0)
+ delta = delta0 * pow(100, 2 * (x - x0));
}
cout << "BetheBloch delta=" << delta << endl;
-
+
// with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
// shell correction, <~100MeV
@@ -114,36 +116,38 @@ namespace corsika::process::EnergyLoss {
HEPEnergyType Iadj = 12_eV * Z + 7_eV; // Iadj<163eV
if (Iadj>=163_eV)
Iadj = 9.76_eV * Z + 58.8_eV * pow(Z, -0.19); // Iadj>=163eV
- double const Cadj = (0.422377/eta2 + 0.0304043/(eta2*eta2) - 0.00038106/(eta2*eta2*eta2)) * 1e-6 * Iadj*Iadj +
- (3.858019/eta2 - 0.1667989/(eta2*eta2) + 0.00157955/(eta2*eta2*eta2)) * 1e-9 * Iadj*Iadj*Iadj;
+ double const Cadj = (0.422377/eta2 + 0.0304043/(eta2*eta2) -
+ 0.00038106/(eta2*eta2*eta2)) * 1e-6 * Iadj*Iadj + (3.858019/eta2 -
+ 0.1667989/(eta2*eta2) + 0.00157955/(eta2*eta2*eta2)) * 1e-9 * Iadj*Iadj*Iadj;
*/
-
+
// Barkas correction O(Z3) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
double A = 1;
- if (p.GetPID() == particles::Code::Nucleus)
- A = p.GetNuclearA();
- //double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
- //double const Llow = 0.01 * Erel;
- //double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
- //double const barkas = Z * Llow*Lhigh/(Llow+Lhigh); // RU, I think the Z was missing...
+ if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
+ // double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
+ // double const Llow = 0.01 * Erel;
+ // double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
+ // double const barkas = Z * Llow*Lhigh/(Llow+Lhigh); // RU, I think the Z was
+ // missing...
double const barkas = 1; // does not work yet
// Bloch correction for O(Z4) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
- const double alpha = 1./137.035999173;
- double const y2 = Z*Z * alpha*alpha/beta2;
- double const bloch = -y2 * (1.202 - y2*(1.042-0.855*y2+0.343*y2*y2) );
-
- // cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch << endl;
-
- double const aux = 2*me*c2*beta2*gamma2*Wmax / (Ieff*Ieff);
- return -K * Z2 * ZoverA / beta2 * (0.5 * log(aux) - beta2 - Cadj/Z - delta/2 + barkas + bloch) * dX;
+ const double alpha = 1. / 137.035999173;
+ double const y2 = Z * Z * alpha * alpha / beta2;
+ double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
+
+ // cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch <<
+ // endl;
+
+ double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
+ return -K * Z2 * ZoverA / beta2 *
+ (0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
}
-
+
process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) {
- if (p.GetChargeNumber()==0)
- return process::EProcessReturn::eOk;
+ if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index c9fe4aa61d617d8a42efa6000017564f8d61ee68..75971bc8c0762b53a254f401da55ecf137407aef 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -44,8 +44,10 @@ namespace corsika::process::EnergyLoss {
void SaveSave();
private:
- corsika::units::si::HEPEnergyType BetheBloch(corsika::setup::Stack::ParticleType& p, const corsika::units::si::GrammageType dX);
-
+ corsika::units::si::HEPEnergyType BetheBloch(
+ corsika::setup::Stack::ParticleType& p,
+ const corsika::units::si::GrammageType dX);
+
int GetXbin(corsika::setup::Stack::ParticleType& p,
const corsika::units::si::HEPEnergyType dE);
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 8bdfb554fb44538c816e76440acb8e0cf0845609..77241b106dff5d6f674f60ec9a2acfe763ca09be 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -8,9 +8,9 @@
* the license.
*/
+#include <Pythia8/Pythia.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/pythia/Random.h>
-#include <Pythia8/Pythia.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -28,45 +28,44 @@ using Track = Trajectory;
namespace corsika::process::pythia {
Decay::Decay(vector<particles::Code> pParticles)
- : fTrackedParticles(pParticles) {}
-
+ : fTrackedParticles(pParticles) {}
+
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 );
+ fPythia.setRndmEnginePtr(rndm);
fPythia.readString("Next:numberShowInfo = 0");
fPythia.readString("Next:numberShowProcess = 0");
fPythia.readString("Next:numberShowEvent = 0");
fPythia.readString("Print:quiet = on");
-
+
fPythia.readString("ProcessLevel:all = off");
fPythia.readString("ProcessLevel:resonanceDecays = off");
fPythia.particleData.readString("59:m0 = 101.00");
-
+
fPythia.init();
}
void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
- for (auto p : particleList)
- Decay::SetStable( p );
+ for (auto p : particleList) Decay::SetStable(p);
}
void Decay::SetUnstable(const particles::Code pCode) {
cout << "Pythia::Decay: setting " << pCode << " unstable.." << endl;
- fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , true);
+ fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
}
void Decay::SetStable(const particles::Code pCode) {
cout << "Pythia::Decay: setting " << pCode << " stable.." << endl;
- fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , false);
+ fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), false);
}
template <>
@@ -95,29 +94,29 @@ namespace corsika::process::pythia {
// coordinate system, get global frame of reference
geometry::CoordinateSystem& rootCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
fCount++;
// pythia stack
Pythia8::Event& event = fPythia.event;
event.reset();
-
+
// set particle unstable
- Decay::SetUnstable( p.GetPID() );
+ Decay::SetUnstable(p.GetPID());
// input particle PDG
- auto const pdgCode = static_cast<int>( particles::GetPDG( p.GetPID() ));
+ auto const pdgCode = static_cast<int>(particles::GetPDG(p.GetPID()));
- auto const pcomp = p.GetMomentum().GetComponents();
- double px = pcomp[0] / 1_GeV ;
- double py = pcomp[1] / 1_GeV ;
- double pz = pcomp[2] / 1_GeV ;
+ auto const pcomp = p.GetMomentum().GetComponents();
+ double px = pcomp[0] / 1_GeV;
+ double py = pcomp[1] / 1_GeV;
+ double pz = pcomp[2] / 1_GeV;
double en = p.GetEnergy() / 1_GeV;
- double m = particles::GetMass( p.GetPID() ) / 1_GeV;
+ double m = particles::GetMass(p.GetPID()) / 1_GeV;
// add particle to pythia stack
- event.append( pdgCode, 1, 0, 0, px, py, pz, en, m);
+ event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
if (!fPythia.next())
cout << "Pythia::Decay: decay failed!" << endl;
@@ -130,27 +129,27 @@ namespace corsika::process::pythia {
// loop over final state
for (int i = 0; i < event.size(); ++i)
if (event[i].isFinal()) {
- auto const pyId = particles::ConvertFromPDG(static_cast<particles::PDGCode>(event[i].id()));
- HEPEnergyType pyEn = event[i].e() * 1_GeV;
- MomentumVector pyP(rootCS, { event[i].px() * 1_GeV,
- event[i].py() * 1_GeV,
- event[i].pz() * 1_GeV});
-
- cout << "particle: id=" << pyId << " momentum=" << pyP.GetComponents() / 1_GeV
- << " energy=" << pyEn << endl;
-
- p.AddSecondary( tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- pyId, pyEn, pyP, decayPoint, t0} );
+ auto const pyId =
+ particles::ConvertFromPDG(static_cast<particles::PDGCode>(event[i].id()));
+ HEPEnergyType pyEn = event[i].e() * 1_GeV;
+ MomentumVector pyP(rootCS, {event[i].px() * 1_GeV, event[i].py() * 1_GeV,
+ event[i].pz() * 1_GeV});
+
+ cout << "particle: id=" << pyId << " momentum=" << pyP.GetComponents() / 1_GeV
+ << " energy=" << pyEn << endl;
+
+ p.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ pyId, pyEn, pyP, decayPoint, t0});
}
-
+
// set particle stable
- Decay::SetStable( p.GetPID() );
+ Decay::SetStable(p.GetPID());
// remove original particle from corsika stack
p.Delete();
// if (fCount>10) throw std::runtime_error("stop here");
}
-
} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Decay.h b/Processes/Pythia/Decay.h
index 06f1838867a9be760128eefadaf3e653e8c89598..fc76137160177afabe0d6624b2bc1b708f545ab4 100644
--- a/Processes/Pythia/Decay.h
+++ b/Processes/Pythia/Decay.h
@@ -12,28 +12,28 @@
#ifndef _include_corsika_process_pythia_decay_h_
#define _include_corsika_process_pythia_decay_h_
+#include <Pythia8/Pythia.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/DecayProcess.h>
-#include <Pythia8/Pythia.h>
namespace corsika::process {
namespace pythia {
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;
public:
- Decay(std::vector<corsika::particles::Code> );
+ Decay(std::vector<corsika::particles::Code>);
~Decay();
void Init();
void SetParticleListStable(const std::vector<particles::Code>);
- void SetUnstable(const corsika::particles::Code );
- void SetStable(const corsika::particles::Code );
+ void SetUnstable(const corsika::particles::Code);
+ void SetStable(const corsika::particles::Code);
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle const& p);
@@ -44,7 +44,7 @@ namespace corsika::process {
private:
Pythia8::Pythia fPythia;
};
-
+
} // namespace pythia
} // namespace corsika::process
diff --git a/Processes/Pythia/Random.cc b/Processes/Pythia/Random.cc
index 1432cc3cb4267a641118baf8c22e1e37f45ba051..8fcbcf4ee7274ba36252f7969306073f3cffddd5 100644
--- a/Processes/Pythia/Random.cc
+++ b/Processes/Pythia/Random.cc
@@ -12,10 +12,9 @@
namespace corsika::process::pythia {
- double Random::flat(){
+ double Random::flat() {
std::uniform_real_distribution<double> dist;
return dist(fRNG);
}
-
} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Random.h b/Processes/Pythia/Random.h
index 885bb598a983079b9b0fe438a78762f109a88067..276fc532ca53ac84b708e8ad12cb16fbaed01c9a 100644
--- a/Processes/Pythia/Random.h
+++ b/Processes/Pythia/Random.h
@@ -12,8 +12,8 @@
#ifndef _include_corsika_process_pythia_random_h_
#define _include_corsika_process_pythia_random_h_
-#include <corsika/random/RNGManager.h>
#include <Pythia8/Pythia.h>
+#include <corsika/random/RNGManager.h>
namespace corsika::process {
@@ -21,11 +21,12 @@ namespace corsika::process {
class Random : public Pythia8::RndmEngine {
double flat();
+
private:
corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
+ corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
};
-
+
} // namespace pythia
} // namespace corsika::process
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index 05895edaca71bc5a54ad2dccebd19fc7cdd883dc..9e2d5dd39cd7a02297c22f298802cc79d5f040a9 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -9,8 +9,8 @@
* the license.
*/
-#include <corsika/process/pythia/Decay.h>
#include <Pythia8/Pythia.h>
+#include <corsika/process/pythia/Decay.h>
#include <corsika/random/RNGManager.h>
@@ -58,26 +58,23 @@ TEST_CASE("Pythia", "[processes]") {
// loop over final state
for (int i = 0; i < pythia.event.size(); ++i)
if (pythia.event[i].isFinal()) {
- cout << "particle: id=" << pythia.event[i].id() << endl;
+ cout << "particle: id=" << pythia.event[i].id() << endl;
}
}
-
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);
-
- model.Init();
+ model.Init();
}
-
}
#include <corsika/geometry/Point.h>
@@ -97,7 +94,7 @@ TEST_CASE("Pythia", "[processes]") {
using namespace corsika;
using namespace corsika::units::si;
-TEST_CASE("pytia decay"){
+TEST_CASE("pytia decay") {
// setup environment, geometry
environment::Environment env;
@@ -126,7 +123,7 @@ TEST_CASE("pytia decay"){
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
SECTION("pythia decay") {
-
+
setup::Stack stack;
const HEPEnergyType E0 = 10_GeV;
HEPMomentumType P0 =
@@ -138,20 +135,17 @@ TEST_CASE("pytia decay"){
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");
-
+
process::pythia::Decay model(particleList);
-
+
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
stack);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
-
}
-
}
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 6816a350bb611fa0f5e17dc6bbd02fdc03e57bd8..3cd009c580d10845ea2219b59322dc1855a88295 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -29,12 +29,12 @@ using Track = Trajectory;
namespace corsika::process::sibyll {
- Decay::Decay(vector<particles::Code>pParticles)
- : fTrackedParticles(pParticles) {}
+ Decay::Decay(vector<particles::Code> pParticles)
+ : fTrackedParticles(pParticles) {}
Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
void Decay::Init() {
SetHadronsUnstable();
- SetParticleListStable( fTrackedParticles );
+ SetParticleListStable(fTrackedParticles);
}
void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
@@ -45,8 +45,7 @@ namespace corsika::process::sibyll {
// set particles unstable
SetHadronsUnstable();
cout << "Interaction: setting tracked hadrons stable.." << endl;
- for (auto p : particleList)
- Decay::SetStable( p );
+ for (auto p : particleList) Decay::SetStable(p);
}
void Decay::SetUnstable(const particles::Code pCode) {
@@ -177,7 +176,6 @@ namespace corsika::process::sibyll {
ss.Clear();
// remove original particle from corsika stack
p.Delete();
-
}
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 5a1f6323b0cb0dde9d73395d9ecc77f460860c53..76d61f7b5a44928a4a43d6d97aac7cfd231514eb 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -29,8 +29,8 @@ namespace corsika::process {
void Init();
void SetParticleListStable(const std::vector<particles::Code>);
- void SetUnstable(const corsika::particles::Code );
- void SetStable(const corsika::particles::Code );
+ void SetUnstable(const corsika::particles::Code);
+ void SetStable(const corsika::particles::Code);
void SetAllStable();
void SetHadronsUnstable();
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 5594ec39de31270b0e1014f5aa1e4b64c917c99f..06d301e3a9f4d264e3b6b56c969fd8d3742269cc 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -35,7 +35,7 @@ using Track = Trajectory;
namespace corsika::process::sibyll {
Interaction::Interaction(environment::Environment const& env)
- : fEnvironment(env) {}
+ : fEnvironment(env) {}
Interaction::~Interaction() {
cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl;
@@ -50,29 +50,30 @@ namespace corsika::process::sibyll {
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);
+ 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(const std::vector<particles::Code> particleList) {
- for (auto p : particleList)
- Interaction::SetStable( p );
+ void Interaction::SetParticleListStable(
+ const std::vector<particles::Code> particleList) {
+ for (auto p : particleList) Interaction::SetStable(p);
}
void Interaction::SetUnstable(const particles::Code pCode) {
@@ -87,8 +88,6 @@ namespace corsika::process::sibyll {
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
-
-
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
@@ -97,8 +96,9 @@ namespace corsika::process::sibyll {
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
- if( !IsValidCoMEnergy(CoMenergy) ){
- throw std::runtime_error("Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
+ if (!IsValidCoMEnergy(CoMenergy)) {
+ throw std::runtime_error(
+ "Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
}
const double dEcm = CoMenergy / 1_GeV;
if (particles::IsNucleus(TargetId)) {
@@ -177,7 +177,7 @@ namespace corsika::process::sibyll {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
- auto const [productionCrossSection, elaCrossSection ] =
+ auto const [productionCrossSection, elaCrossSection] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
[[maybe_unused]] auto elaCrossSectionCopy =
elaCrossSection; // ONLY TO AVOID COMPILER WARNING
@@ -193,8 +193,8 @@ namespace corsika::process::sibyll {
// calculate interaction length in medium
//#warning check interaction length units
- GrammageType const int_length =
- mediumComposition.GetAverageMassNumber() * units::constants::u / weightedProdCrossSection;
+ GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
+ units::constants::u / weightedProdCrossSection;
cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< endl;
@@ -299,8 +299,7 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- const auto [sigProd, sigEla ] =
- GetCrossSection(corsikaBeamId, targetId, Ecm);
+ const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy =
sigEla; // to avoid not used warning in array binding
@@ -329,10 +328,10 @@ namespace corsika::process::sibyll {
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl;
- if( Ecm > GetMaxEnergyCoM() )
- throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
+ if (Ecm > GetMaxEnergyCoM())
+ throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
// FR: removed eProjectileLab < 8.5_GeV ||
- if ( Ecm < GetMinEnergyCoM() ) {
+ if (Ecm < GetMinEnergyCoM()) {
cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 0f2fc58efc2d27aa337501c9b7d5174c6ef3f460..21d4675b98f6d1733646aaf7184ff615243556ba 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -37,8 +37,8 @@ namespace corsika::process::sibyll {
void Init();
void SetParticleListStable(const std::vector<particles::Code>);
- void SetUnstable(const corsika::particles::Code );
- void SetStable(const corsika::particles::Code );
+ void SetUnstable(const corsika::particles::Code);
+ void SetStable(const corsika::particles::Code);
bool WasInitialized() { return fInitialized; }
bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
@@ -47,13 +47,11 @@ namespace corsika::process::sibyll {
int GetMaxTargetMassNumber() { return fMaxTargetMassNumber; }
corsika::units::si::HEPEnergyType GetMinEnergyCoM() { return fMinEnergyCoM; }
corsika::units::si::HEPEnergyType GetMaxEnergyCoM() { return fMaxEnergyCoM; }
- bool IsValidTarget(corsika::particles::Code TargetId)
- {
- return ( corsika::particles::GetNucleusA(TargetId) < fMaxTargetMassNumber )
- && corsika::particles::IsNucleus( TargetId );
- }
-
-
+ bool IsValidTarget(corsika::particles::Code TargetId) {
+ return (corsika::particles::GetNucleusA(TargetId) < fMaxTargetMassNumber) &&
+ corsika::particles::IsNucleus(TargetId);
+ }
+
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
@@ -76,10 +74,10 @@ namespace corsika::process::sibyll {
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
const bool fInternalDecays = true;
- const corsika::units::si::HEPEnergyType fMinEnergyCoM =
- 10. * 1e9 * corsika::units::si::electronvolt;
- const corsika::units::si::HEPEnergyType fMaxEnergyCoM =
- 1.e6 * 1e9 * corsika::units::si::electronvolt;
+ const corsika::units::si::HEPEnergyType fMinEnergyCoM =
+ 10. * 1e9 * corsika::units::si::electronvolt;
+ const corsika::units::si::HEPEnergyType fMaxEnergyCoM =
+ 1.e6 * 1e9 * corsika::units::si::electronvolt;
const int fMaxTargetMassNumber = 18;
};
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 6c543c3581538976d0c6e84189dba272d773538b..8b37868b6e879d5614cf88379b4222615f7ce140 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -55,10 +55,11 @@ namespace corsika::process::sibyll {
if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
// check compatibility of energy ranges, someone could try to use low-energy model..
- if(!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM())||
- !fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
- throw std::runtime_error("NuclearInteraction: hadronic interaction model incompatible!");
-
+ if (!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
+ !fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
+ throw std::runtime_error(
+ "NuclearInteraction: hadronic interaction model incompatible!");
+
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
@@ -67,117 +68,116 @@ namespace corsika::process::sibyll {
InitializeNuclearCrossSections();
}
- void NuclearInteraction::InitializeNuclearCrossSections()
- {
+ void NuclearInteraction::InitializeNuclearCrossSections() {
using namespace corsika::particles;
using namespace units::si;
auto& universe = *(fEnvironment.GetUniverse());
auto const allElementsInUniverse = std::invoke([&]() {
- std::set<particles::Code> allElementsInUniverse;
- auto collectElements = [&](auto& vtn) {
- if (auto const mp = vtn.GetModelPropertiesPtr();
- mp != nullptr) { // do not query Universe it self, it has no ModelProperties
- auto const& comp = mp->GetNuclearComposition().GetComponents();
- for (auto const c : comp)
- allElementsInUniverse.insert(c);
- // std::for_each(comp.cbegin(), comp.cend(),
- // [&](particles::Code c) { allElementsInUniverse.insert(c); });
- }
- };
- universe.walk(collectElements);
- return allElementsInUniverse;
- });
+ std::set<particles::Code> allElementsInUniverse;
+ auto collectElements = [&](auto& vtn) {
+ if (auto const mp = vtn.GetModelPropertiesPtr();
+ mp != nullptr) { // do not query Universe it self, it has no ModelProperties
+ auto const& comp = mp->GetNuclearComposition().GetComponents();
+ for (auto const c : comp) allElementsInUniverse.insert(c);
+ // std::for_each(comp.cbegin(), comp.cend(),
+ // [&](particles::Code c) { allElementsInUniverse.insert(c); });
+ }
+ };
+ universe.walk(collectElements);
+ return allElementsInUniverse;
+ });
-
cout << "NuclearInteraction: initializing nuclear cross sections..." << endl;
-
+
// loop over target components, at most 4!!
- int k =-1;
- for(auto &ptarg: allElementsInUniverse){
+ int k = -1;
+ for (auto& ptarg : allElementsInUniverse) {
++k;
cout << "NuclearInteraction: init target component: " << ptarg << endl;
- const int ib = GetNucleusA( ptarg );
- if(!fHadronicInteraction.IsValidTarget( ptarg )){
- cout << "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id="
- << ptarg << endl;
- throw std::runtime_error(" target can not be handled by hadronic interaction model! ");
+ const int ib = GetNucleusA(ptarg);
+ if (!fHadronicInteraction.IsValidTarget(ptarg)) {
+ cout << "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id="
+ << ptarg << endl;
+ throw std::runtime_error(
+ " target can not be handled by hadronic interaction model! ");
}
- fTargetComponentsIndex.insert( std::pair<Code,int>(ptarg, k) );
+ fTargetComponentsIndex.insert(std::pair<Code, int>(ptarg, k));
// loop over energies, fNEnBins log. energy bins
- for(int i=0; i<GetNEnergyBins(); ++i){
- // hard coded energy grid, has to be aligned to definition in signuc2!!, no comment..
- const units::si::HEPEnergyType Ecm = pow(10., 1. + 1.*i ) * 1_GeV;
- // get p-p cross sections
- auto const protonId = Code::Proton;
- auto const [siginel, sigela ] =
- fHadronicInteraction.GetCrossSection(protonId, protonId, Ecm);
- const double dsig = siginel / 1_mbarn;
- const double dsigela = sigela / 1_mbarn;
- // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
- for(int j=1; j<fMaxNucleusAProjectile; ++j){
- const int jj = j+1;
- double sig_out, dsig_out, sigqe_out, dsigqe_out;
- sigma_mc_(jj,ib,dsig,dsigela,fNSample,sig_out,dsig_out,sigqe_out,dsigqe_out);
- // write to table
- cnucsignuc_.sigma[j][k][i] = sig_out;
- cnucsignuc_.sigqe[j][k][i] = sigqe_out;
- }
+ for (int i = 0; i < GetNEnergyBins(); ++i) {
+ // hard coded energy grid, has to be aligned to definition in signuc2!!, no
+ // comment..
+ const units::si::HEPEnergyType Ecm = pow(10., 1. + 1. * i) * 1_GeV;
+ // get p-p cross sections
+ auto const protonId = Code::Proton;
+ auto const [siginel, sigela] =
+ fHadronicInteraction.GetCrossSection(protonId, protonId, Ecm);
+ const double dsig = siginel / 1_mbarn;
+ const double dsigela = sigela / 1_mbarn;
+ // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
+ for (int j = 1; j < fMaxNucleusAProjectile; ++j) {
+ const int jj = j + 1;
+ double sig_out, dsig_out, sigqe_out, dsigqe_out;
+ sigma_mc_(jj, ib, dsig, dsigela, fNSample, sig_out, dsig_out, sigqe_out,
+ dsigqe_out);
+ // write to table
+ cnucsignuc_.sigma[j][k][i] = sig_out;
+ cnucsignuc_.sigqe[j][k][i] = sigqe_out;
+ }
}
}
- cout << "NuclearInteraction: cross sections for " << fTargetComponentsIndex.size() << " components initialized!" << endl;
- for(auto &ptarg: allElementsInUniverse){
+ cout << "NuclearInteraction: cross sections for " << fTargetComponentsIndex.size()
+ << " components initialized!" << endl;
+ for (auto& ptarg : allElementsInUniverse) {
cout << "cross section table: " << ptarg << endl;
- PrintCrossSectionTable( ptarg );
+ PrintCrossSectionTable(ptarg);
}
-
}
- void NuclearInteraction::PrintCrossSectionTable( corsika::particles::Code pCode)
- {
+ void NuclearInteraction::PrintCrossSectionTable(corsika::particles::Code pCode) {
using namespace corsika::particles;
- const int k = fTargetComponentsIndex.at( pCode );
- Code pNuclei [] = {Code::Helium, Code::Lithium7, Code::Oxygen,
- Code::Neon, Code::Argon, Code::Iron};
+ const int k = fTargetComponentsIndex.at(pCode);
+ Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
+ Code::Neon, Code::Argon, Code::Iron};
cout << "en/A ";
- for(auto &j: pNuclei)
- cout << std::setw(9) << j;
+ for (auto& j : pNuclei) cout << std::setw(9) << j;
cout << endl;
// loop over energy bins
- for(int i=0; i<GetNEnergyBins(); ++i){
+ for (int i = 0; i < GetNEnergyBins(); ++i) {
cout << " " << i << " ";
- for(auto &n: pNuclei){
- auto const j= GetNucleusA( n );
- cout << " " << std::setprecision(5) << std::setw(8) << cnucsignuc_.sigma[j-1][k][i];
+ for (auto& n : pNuclei) {
+ auto const j = GetNucleusA(n);
+ cout << " " << std::setprecision(5) << std::setw(8)
+ << cnucsignuc_.sigma[j - 1][k][i];
}
- cout << endl;
+ cout << endl;
}
}
-
- units::si::CrossSectionType NuclearInteraction::ReadCrossSectionTable(const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc)
- {
+
+ units::si::CrossSectionType NuclearInteraction::ReadCrossSectionTable(
+ const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc) {
using namespace corsika::particles;
using namespace units::si;
- const int ib = fTargetComponentsIndex.at( pTarget )+1; // table index in fortran
- auto const ECoMNuc = sqrt( 2. * corsika::units::constants::nucleonMass * elabnuc );
- if( ECoMNuc < GetMinEnergyPerNucleonCoM() || ECoMNuc > GetMaxEnergyPerNucleonCoM() )
+ const int ib = fTargetComponentsIndex.at(pTarget) + 1; // table index in fortran
+ auto const ECoMNuc = sqrt(2. * corsika::units::constants::nucleonMass * elabnuc);
+ if (ECoMNuc < GetMinEnergyPerNucleonCoM() || ECoMNuc > GetMaxEnergyPerNucleonCoM())
throw std::runtime_error("NuclearInteraction: energy outside tabulated range!");
const double e0 = elabnuc / 1_GeV;
double sig;
cout << "ReadCrossSectionTable: " << ia << " " << ib << " " << e0 << endl;
- signuc2_(ia,ib,e0,sig);
+ signuc2_(ia, ib, e0, sig);
cout << "ReadCrossSectionTable: sig=" << sig << endl;
return sig * 1_mbarn;
}
-
+
// TODO: remove elastic cross section?
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
NuclearInteraction::GetCrossSection(Particle& p, const particles::Code TargetId) {
using namespace units::si;
- if ( p.GetPID() != particles::Code::Nucleus)
+ if (p.GetPID() != particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: particle not a nucleus!");
@@ -200,11 +200,11 @@ namespace corsika::process::sibyll {
}
if (fHadronicInteraction.IsValidTarget(TargetId)) {
- auto const sigProd = ReadCrossSectionTable( iBeamA, TargetId, LabEnergyPerNuc );
+ auto const sigProd = ReadCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc);
cout << "cross section (mb): " << sigProd / 1_mbarn << endl;
return std::make_tuple(sigProd, 0_mbarn);
} else {
- throw std::runtime_error( "target outside range.");
+ throw std::runtime_error("target outside range.");
}
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
std::numeric_limits<double>::infinity() * 1_mbarn);
@@ -261,7 +261,8 @@ namespace corsika::process::sibyll {
// energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific
- if (ElabNuc >= 8.5_GeV && ECoMNN >= fMinEnergyPerNucleonCoM && ECoMNN < fMaxEnergyPerNucleonCoM) {
+ if (ElabNuc >= 8.5_GeV && ECoMNN >= fMinEnergyPerNucleonCoM &&
+ ECoMNN < fMaxEnergyPerNucleonCoM) {
// get target from environment
/*
@@ -358,7 +359,7 @@ namespace corsika::process::sibyll {
// projectile nucleon number
const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
- if (kAProj > GetMaxNucleusAProjectile() )
+ if (kAProj > GetMaxNucleusAProjectile())
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics
@@ -460,7 +461,7 @@ namespace corsika::process::sibyll {
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) kATarget = 1;
cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
- if(!fHadronicInteraction.IsValidTarget(targetCode))
+ if (!fHadronicInteraction.IsValidTarget(targetCode))
throw std::runtime_error("target outside range. ");
// end of target sampling
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index ad1d0bc7e06e064d859f0817d0128d6b1d73170e..c62249bbefb0193d603714cc960743e921ee1a5d 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -42,10 +42,15 @@ namespace corsika::process::sibyll {
void Init();
void InitializeNuclearCrossSections();
void PrintCrossSectionTable(corsika::particles::Code);
- corsika::units::si::CrossSectionType ReadCrossSectionTable(const int, corsika::particles::Code, corsika::units::si::HEPEnergyType);
- corsika::units::si::HEPEnergyType GetMinEnergyPerNucleonCoM() { return fMinEnergyPerNucleonCoM; }
- corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() { return fMaxEnergyPerNucleonCoM; }
- int GetMaxNucleusAProjectile() {return fMaxNucleusAProjectile; }
+ corsika::units::si::CrossSectionType ReadCrossSectionTable(
+ const int, corsika::particles::Code, corsika::units::si::HEPEnergyType);
+ corsika::units::si::HEPEnergyType GetMinEnergyPerNucleonCoM() {
+ return fMinEnergyPerNucleonCoM;
+ }
+ corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() {
+ return fMaxEnergyPerNucleonCoM;
+ }
+ int GetMaxNucleusAProjectile() { return fMaxNucleusAProjectile; }
int GetMaxNFragments() { return fMaxNFragments; }
int GetNEnergyBins() { return fNEnBins; }
template <typename Particle>
@@ -61,7 +66,7 @@ namespace corsika::process::sibyll {
private:
corsika::environment::Environment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
- std::map<corsika::particles::Code,int> fTargetComponentsIndex;
+ std::map<corsika::particles::Code, int> fTargetComponentsIndex;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
const int fNSample = 500; // number of samples in MC estimation of cross section
@@ -70,10 +75,10 @@ namespace corsika::process::sibyll {
const int fMaxNFragments = 60;
// energy limits defined by table used for cross section in signuc.f
// 10**1 GeV to 10**6 GeV
- const corsika::units::si::HEPEnergyType fMinEnergyPerNucleonCoM =
- 10. * 1e9 * corsika::units::si::electronvolt;
- const corsika::units::si::HEPEnergyType fMaxEnergyPerNucleonCoM =
- 1.e6 * 1e9 * corsika::units::si::electronvolt;
+ const corsika::units::si::HEPEnergyType fMinEnergyPerNucleonCoM =
+ 10. * 1e9 * corsika::units::si::electronvolt;
+ const corsika::units::si::HEPEnergyType fMaxEnergyPerNucleonCoM =
+ 1.e6 * 1e9 * corsika::units::si::electronvolt;
};
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/nuclib.h b/Processes/Sibyll/nuclib.h
index ddd1856cd2eea972e257ba2958f6705ae15bcf96..44e86f50ef825bc8bfd7d98763d86be1ab903fa4 100644
--- a/Processes/Sibyll/nuclib.h
+++ b/Processes/Sibyll/nuclib.h
@@ -34,13 +34,12 @@ extern struct {
*/
extern struct { double ppp[60][3]; } fragments_;
- // COMMON /cnucsignuc/SIGMA(6,4,56), SIGQE(6,4,56)
- extern struct
- {
- double sigma[56][4][6];
- double sigqe[56][4][6];
- } cnucsignuc_;
-
+// COMMON /cnucsignuc/SIGMA(6,4,56), SIGQE(6,4,56)
+extern struct {
+ double sigma[56][4][6];
+ double sigqe[56][4][6];
+} cnucsignuc_;
+
// NUCLIB
// subroutine to initiate nuclib
@@ -56,6 +55,7 @@ void signuc_(const int&, const double&, double&);
void signuc2_(const int&, const int&, const double&, double&);
-void sigma_mc_(const int&, const int&, const double&, const double&, const int&, double&, double&, double&, double&);
+void sigma_mc_(const int&, const int&, const double&, const double&, const int&, double&,
+ double&, double&, double&);
}
#endif