diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index cf8df0febd63f6a144c69343bb66fd01d10086c4..772853db877c6f099370b3b8f35ee949b2129c45 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -250,14 +250,15 @@ int main() {
ProcessCut cut(20_GeV);
// corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
- // corsika::process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env);
+ // corsika::process::HadronicElasticModel::HadronicElasticInteraction
+ // hadronicElastic(env);
corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
- //auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
+ // auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
-
+
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
@@ -265,15 +266,16 @@ int main() {
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Proton;
- const HEPEnergyType E0 = 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
+ const HEPEnergyType E0 =
+ 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.;
double phi = 0.;
-
+
{
- auto elab2plab = []( HEPEnergyType Elab, HEPMassType m){
- return sqrt(Elab * Elab - m * m);
- };
- HEPMomentumType P0 = elab2plab( E0, corsika::particles::GetMass( beamCode ) );
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt(Elab * Elab - m * m);
+ };
+ HEPMomentumType P0 = elab2plab(E0, corsika::particles::GetMass(beamCode));
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
index 752164dc7e5b8cf2ad96bf245e80e40e583ce971..5bf2ece2b1e26af2ccfad2b34e5c3a9534250d86 100644
--- a/Framework/StackInterface/ParticleBase.h
+++ b/Framework/StackInterface/ParticleBase.h
@@ -67,7 +67,7 @@ namespace corsika::stack {
protected:
/** @name Access to underlying stack data
- @{
+ @{
*/
auto& GetStackData() { return GetIterator().GetStackData(); }
const auto& GetStackData() const { return GetIterator().GetStackData(); }
@@ -75,7 +75,7 @@ namespace corsika::stack {
const auto& GetStack() const { return GetIterator().GetStack(); }
///@}
- /**
+ /**
* return the index number of the underlying iterator object
*/
int GetIndex() const { return GetIterator().GetIndex(); }
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index 91d671f5fd83d44194343010c6c677490063ec71..bc9aa2e5c0ee0388df07a6ea4a90b4c3865c4485 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -12,7 +12,7 @@
#ifndef _include_Stack_h__
#define _include_Stack_h__
-#include <corsika/stack/StackIteratorInterface.h>
+#include <corsika/stack/StackIteratorInterface.h>
#include <stdexcept>
@@ -29,7 +29,7 @@ namespace corsika::stack {
Important: ParticleInterface must inherit from ParticleBase !
*/
-
+
template <typename>
class ParticleInterface; // forward decl
diff --git a/Framework/StackInterface/StackIteratorInterface.h b/Framework/StackInterface/StackIteratorInterface.h
index 7ac72c34bca96264fa9fed4e9906cecf882166db..ad00c6c0939a4e99c020164ef4cd85cbc81a8788 100644
--- a/Framework/StackInterface/StackIteratorInterface.h
+++ b/Framework/StackInterface/StackIteratorInterface.h
@@ -112,8 +112,8 @@ namespace corsika::stack {
}
public:
- /** @name Iterator interface
- @{
+ /** @name Iterator interface
+ @{
*/
StackIteratorInterface& operator++() {
++fIndex;
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index a621768a1b328437fd8128f581e6315ff194fd7e..557eade5e6ea1ef9f9f85d84718f825c5cf0d270 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -32,7 +32,7 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
bool fInitialized = false;
-
+
public:
Interaction(corsika::environment::Environment const& env)
: fEnvironment(env) {}
@@ -48,17 +48,19 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize Sibyll
- if(!fInitialized){
- sibyll_ini_();
- fInitialized = true;
+ if (!fInitialized) {
+ sibyll_ini_();
+ fInitialized = true;
}
}
- bool WasInitialized(){ return fInitialized;}
-
- std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType, int> GetCrossSection(
- const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
- const corsika::units::si::HEPEnergyType CoMenergy) {
+ bool WasInitialized() { return fInitialized; }
+
+ std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
+ int>
+ GetCrossSection(const corsika::particles::Code BeamId,
+ const corsika::particles::Code TargetId,
+ const corsika::units::si::HEPEnergyType CoMenergy) {
using namespace corsika::units::si;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
@@ -70,14 +72,14 @@ namespace corsika::process::sibyll {
if (iTarget > 18 || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
- sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
+ sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == corsika::particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
- iTarget = 1;
+ iTarget = 1;
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
- sigEla = std::numeric_limits<double>::infinity();
+ sigEla = std::numeric_limits<double>::infinity();
}
return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
}
@@ -192,12 +194,12 @@ namespace corsika::process::sibyll {
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl;
-
- if(corsika::particles::IsNucleus(corsikaBeamId)){
- // nuclei handled by different process, this should not happen
- throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
+
+ if (corsika::particles::IsNucleus(corsikaBeamId)) {
+ // nuclei handled by different process, this should not happen
+ throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
}
-
+
if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -273,7 +275,8 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- const auto [sigProd, sigEla, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ const auto [sigProd, sigEla, nNuc] =
+ GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
int ideleteme = nNuc; // to avoid not used warning in array binding
ideleteme = ideleteme; // to avoid not used warning in array binding
@@ -306,7 +309,7 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << std::endl;
- throw std::runtime_error("energy too low for SIBYLL");
+ throw std::runtime_error("energy too low for SIBYLL");
// p.Delete(); delete later... different process
} else {
fCount++;
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 35ab127c993c9b0dc16ca4cb71f19bf5f0b5c70e..572f17683e88876535898cb8038fa7229c8d7d4d 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -24,16 +24,17 @@
namespace corsika::process::sibyll {
- class NuclearInteraction : public corsika::process::InteractionProcess<NuclearInteraction> {
+ class NuclearInteraction
+ : public corsika::process::InteractionProcess<NuclearInteraction> {
int fCount = 0;
int fNucCount = 0;
public:
- NuclearInteraction(corsika::environment::Environment const& env, corsika::process::sibyll::Interaction& hadint)
- : fEnvironment(env),
- fHadronicInteraction(hadint)
- {}
+ NuclearInteraction(corsika::environment::Environment const& env,
+ corsika::process::sibyll::Interaction& hadint)
+ : fEnvironment(env)
+ , fHadronicInteraction(hadint) {}
~NuclearInteraction() {
std::cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
<< std::endl;
@@ -47,56 +48,63 @@ namespace corsika::process::sibyll {
// initialize hadronic interaction module
// TODO: save to run multiple initializations?
- if(!fHadronicInteraction.WasInitialized())
- fHadronicInteraction.Init();
-
+ if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
+
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
}
// TODO: remove number of nucleons, avg target mass is available in environment
- std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType, int> GetCrossSection(
- const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
- const corsika::units::si::HEPEnergyType CoMenergy) {
+ std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
+ int>
+ GetCrossSection(const corsika::particles::Code BeamId,
+ const corsika::particles::Code TargetId,
+ const corsika::units::si::HEPEnergyType CoMenergy) {
using namespace corsika::units::si;
double sigProd;
- std::cout << "NuclearInteraction: GetCrossSection: called with: beamId= "<<BeamId
- << " TargetId= " << TargetId << " CoMenergy= " << CoMenergy / 1_GeV << std::endl;
+ std::cout << "NuclearInteraction: GetCrossSection: called with: beamId= " << BeamId
+ << " TargetId= " << TargetId << " CoMenergy= " << CoMenergy / 1_GeV
+ << std::endl;
- if(!corsika::particles::IsNucleus(BeamId) ){
- // ask hadronic interaction model for hadron cross section
- return fHadronicInteraction.GetCrossSection( BeamId, TargetId, CoMenergy);
+ if (!corsika::particles::IsNucleus(BeamId)) {
+ // ask hadronic interaction model for hadron cross section
+ return fHadronicInteraction.GetCrossSection(BeamId, TargetId, CoMenergy);
}
-
+
// use nuclib to calc. nuclear cross sections
// TODO: for now assumes air with hard coded composition
// extend to arbitrary mixtures, requires smarter initialization
-
+
const int iBeam = corsika::particles::GetNucleusA(BeamId);
- if( iBeam > 56 || iBeam < 2){
- //std::cout<<"NuclearInteraction: beam nucleus outside allowed range for NUCLIB!" << std::endl;
- throw std::runtime_error("NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for NUCLIB!");
+ if (iBeam > 56 || iBeam < 2) {
+ // std::cout<<"NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
+ // << std::endl;
+ throw std::runtime_error(
+ "NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
+ "NUCLIB!");
}
-
+
const double dEcm = CoMenergy / 1_GeV;
- if(dEcm<10.)
- throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
+ if (dEcm < 10.)
+ throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
// TODO: these limitations are still sibyll specific.
- // available target nuclei depends on the hadronic interaction model and the initialization
+ // available target nuclei depends on the hadronic interaction model and the
+ // initialization
if (corsika::particles::IsNucleus(TargetId)) {
- const int iTarget = corsika::particles::GetNucleusA(TargetId);
- if (iTarget > 18 || iTarget == 0)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 are allowed.");
- std::cout<< "NuclearInteraction: calling signuc.." << std::endl;
- signuc_(iBeam, dEcm, sigProd);
- std::cout << "cross section: " << sigProd << std::endl;
- return std::make_tuple(sigProd * 1_mbarn, 0_mbarn, iTarget);
+ const int iTarget = corsika::particles::GetNucleusA(TargetId);
+ if (iTarget > 18 || iTarget == 0)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 are allowed.");
+ std::cout << "NuclearInteraction: calling signuc.." << std::endl;
+ signuc_(iBeam, dEcm, sigProd);
+ std::cout << "cross section: " << sigProd << std::endl;
+ return std::make_tuple(sigProd * 1_mbarn, 0_mbarn, iTarget);
}
- return std::make_tuple(std::numeric_limits<double>::infinity()* 1_mbarn, std::numeric_limits<double>::infinity()* 1_mbarn, 0);
+ return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
+ std::numeric_limits<double>::infinity() * 1_mbarn, 0);
}
template <typename Particle, typename Track>
@@ -114,14 +122,14 @@ namespace corsika::process::sibyll {
const particles::Code corsikaBeamId = p.GetPID();
- if(!corsika::particles::IsNucleus(corsikaBeamId)){
- // no nuclear interaction
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ if (!corsika::particles::IsNucleus(corsikaBeamId)) {
+ // no nuclear interaction
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
-
+
// read from cross section code table
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
+ corsika::particles::Neutron::GetMass());
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
@@ -134,71 +142,72 @@ namespace corsika::process::sibyll {
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
- (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
+ (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
std::cout << "NuclearInteraction: LambdaInt: \n"
- << " input energy: " << Elab / 1_GeV << endl
- << " input energy CoM: " << ECoM / 1_GeV << endl
- << " beam pid:" << corsikaBeamId << endl;
-
- if ( Elab >= 8.5_GeV && ECoM >= 10_GeV) {
-
- // get target from environment
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- const auto currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- const auto mediumComposition =
- currentNode->GetModelProperties().GetNuclearComposition();
- // determine average interaction length
- // weighted sum
- int i = -1;
- double avgTargetMassNumber = 0.;
- si::CrossSectionType weightedProdCrossSection = 0_mbarn;
- // get weights of components from environment/medium
- const auto w = mediumComposition.GetFractions();
- // loop over components in medium
- for (auto const targetId : mediumComposition.GetComponents()) {
- i++;
- cout << "NuclearInteraction: get interaction length for target: " << targetId << endl;
- // TODO: remove number of nucleons, avg target mass is available in environment
- auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
- GetCrossSection(corsikaBeamId, targetId, ECoM);
-
- std::cout << "NuclearInteraction: "
- << "IntLength: nuclib return (mb): "
- << productionCrossSection / 1_mbarn << std::endl;
- weightedProdCrossSection += w[i] * productionCrossSection;
- avgTargetMassNumber += w[i] * numberOfNucleons;
- }
- cout << "NuclearInteraction: "
- << "IntLength: weighted CrossSection (mb): "
- << weightedProdCrossSection / 1_mbarn << endl;
-
- // calculate interaction length in medium
- GrammageType const int_length =
- avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
- std::cout << "NuclearInteraction: "
- << "interaction length (g/cm2): "
- << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
-
- return int_length;
+ << " input energy: " << Elab / 1_GeV << endl
+ << " input energy CoM: " << ECoM / 1_GeV << endl
+ << " beam pid:" << corsikaBeamId << endl;
+
+ if (Elab >= 8.5_GeV && ECoM >= 10_GeV) {
+
+ // get target from environment
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+ */
+ const auto currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ // determine average interaction length
+ // weighted sum
+ int i = -1;
+ double avgTargetMassNumber = 0.;
+ si::CrossSectionType weightedProdCrossSection = 0_mbarn;
+ // get weights of components from environment/medium
+ const auto w = mediumComposition.GetFractions();
+ // loop over components in medium
+ for (auto const targetId : mediumComposition.GetComponents()) {
+ i++;
+ cout << "NuclearInteraction: get interaction length for target: " << targetId
+ << endl;
+ // TODO: remove number of nucleons, avg target mass is available in environment
+ auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
+ GetCrossSection(corsikaBeamId, targetId, ECoM);
+
+ std::cout << "NuclearInteraction: "
+ << "IntLength: nuclib return (mb): "
+ << productionCrossSection / 1_mbarn << std::endl;
+ weightedProdCrossSection += w[i] * productionCrossSection;
+ avgTargetMassNumber += w[i] * numberOfNucleons;
+ }
+ cout << "NuclearInteraction: "
+ << "IntLength: weighted CrossSection (mb): "
+ << weightedProdCrossSection / 1_mbarn << endl;
+
+ // calculate interaction length in medium
+ GrammageType const int_length =
+ avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
+ std::cout << "NuclearInteraction: "
+ << "interaction length (g/cm2): "
+ << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
+
+ return int_length;
} else {
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
}
-
+
template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&s) {
+ corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
// this routine superimposes different nucleon-nucleon interactions
- // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
+ // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
// this routine superimposes different nucleon-nucleon interactions
- // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
+ // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
using namespace corsika::units;
using namespace corsika::utl;
@@ -209,34 +218,33 @@ namespace corsika::process::sibyll {
const auto corsikaProjId = p.GetPID();
cout << "NuclearInteraction: DoInteraction: called with:" << corsikaProjId << endl
- << "NuclearInteraction: projectile mass: " << corsika::particles::GetMass(corsikaProjId) / 1_GeV
- << endl;
-
- if(!IsNucleus(corsikaProjId)){
- cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
- // this should not happen
- //throw std::runtime_error("Non nuclear projectile in NUCLIB!");
- return process::EProcessReturn::eOk;
+ << "NuclearInteraction: projectile mass: "
+ << corsika::particles::GetMass(corsikaProjId) / 1_GeV << endl;
+
+ if (!IsNucleus(corsikaProjId)) {
+ cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
+ // this should not happen
+ // throw std::runtime_error("Non nuclear projectile in NUCLIB!");
+ return process::EProcessReturn::eOk;
}
fCount++;
-
+
const CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
// position and time of interaction, not used in NUCLIB
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
- cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
- << endl;
- cout << "Interaction: time: " << tOrig << endl;
-
+ cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
+ cout << "Interaction: time: " << tOrig << endl;
+
// projectile nucleon number
const int kAProj = GetNucleusA(corsikaProjId);
- if(kAProj>56)
- throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
-
+ if (kAProj > 56)
+ throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
+
// kinematics
// define projectile nucleus
HEPEnergyType const eProjectileLab = p.GetEnergy();
@@ -244,31 +252,31 @@ namespace corsika::process::sibyll {
const FourVector PprojLab(eProjectileLab, pProjectileLab);
cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl
- << "NuclearInteraction: pProj lab: " << pProjectileLab.GetComponents() / 1_GeV
- << endl;
+ << "NuclearInteraction: pProj lab: " << pProjectileLab.GetComponents() / 1_GeV
+ << endl;
// define projectile nucleon
HEPEnergyType const eProjectileNucLab = p.GetEnergy() / kAProj;
auto const pProjectileNucLab = p.GetMomentum() / kAProj;
const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab);
- cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV << endl
- << "NuclearInteraction: pProjNucleon lab: " << pProjectileNucLab.GetComponents() / 1_GeV
- << endl;
-
+ cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV
+ << endl
+ << "NuclearInteraction: pProjNucleon lab: "
+ << pProjectileNucLab.GetComponents() / 1_GeV << endl;
// define target
// always a nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
+ corsika::particles::Neutron::GetMass());
// target is always at rest
const auto eTargetNucLab = 0_GeV + nucleon_mass;
const auto pTargetNucLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab);
-
+
cout << "NuclearInteraction: etarget lab: " << eTargetNucLab / 1_GeV << endl
- << "NuclearInteraction: ptarget lab: " << pTargetNucLab.GetComponents() / 1_GeV
- << endl;
+ << "NuclearInteraction: ptarget lab: " << pTargetNucLab.GetComponents() / 1_GeV
+ << endl;
// center-of-mass energy in nucleon-nucleon frame
auto const PtotNN4 = PtargNucLab + PprojNucLab;
@@ -284,13 +292,13 @@ namespace corsika::process::sibyll {
auto const PtargNucCoM = boost.toCoM(PtargNucLab);
cout << "Interaction: ebeam CoM: " << PprojNucCoM.GetTimeLikeComponent() / 1_GeV
- << endl
- << "Interaction: pbeam CoM: "
- << PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ << endl
+ << "Interaction: pbeam CoM: "
+ << PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargNucCoM.GetTimeLikeComponent() / 1_GeV
- << endl
- << "Interaction: ptarget CoM: "
- << PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ << endl
+ << "Interaction: ptarget CoM: "
+ << PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
// sample target nucleon number
//
@@ -298,63 +306,65 @@ namespace corsika::process::sibyll {
const auto beamId = corsika::particles::Proton::GetCode();
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
const auto& mediumComposition =
- currentNode->GetModelProperties().GetNuclearComposition();
+ currentNode->GetModelProperties().GetNuclearComposition();
cout << "get cross sections for target materials.." << endl;
// get cross sections for target materials
/*
- Here we read the cross section from the interaction model again,
- should be passed from GetInteractionLength if possible
+ Here we read the cross section from the interaction model again,
+ should be passed from GetInteractionLength if possible
*/
auto const& compVec = mediumComposition.GetComponents();
std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
-
+
for (size_t i = 0; i < compVec.size(); ++i) {
- auto const targetId = compVec[i];
- cout << "target component: " << targetId << endl;
- cout << "beam id: " << targetId << endl;
- const auto [sigProd, sigEla, nNuc] = GetCrossSection(beamId,targetId,EcmNN);
- cross_section_of_components[i] = sigProd;
+ auto const targetId = compVec[i];
+ cout << "target component: " << targetId << endl;
+ cout << "beam id: " << targetId << endl;
+ const auto [sigProd, sigEla, nNuc] = GetCrossSection(beamId, targetId, EcmNN);
+ cross_section_of_components[i] = sigProd;
}
const auto targetCode = currentNode->GetModelProperties().SampleTarget(
- cross_section_of_components, fRNG);
+ cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
- FOR NOW: allow nuclei with A<18 or protons only.
- when medium composition becomes more complex, approximations will have to be
- allowed air in atmosphere also contains some Argon.
+ FOR NOW: allow nuclei with A<18 or protons only.
+ when medium composition becomes more complex, approximations will have to be
+ allowed air in atmosphere also contains some Argon.
*/
int kATarget = -1;
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
if (targetCode == corsika::particles::Proton::GetCode()) kATarget = 1;
cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
if (kATarget > 18 || kATarget < 1)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 or protons are "
- "allowed.");
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 or protons are "
+ "allowed.");
// end of target sampling
- // superposition
- cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "<< endl;
+ // superposition
+ cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "
+ << endl;
// get nucleon-nucleon cross section
- // (needed to determine number of scatterings)
+ // (needed to determine number of nucleon-nucleon scatterings)
const auto protonId = corsika::particles::Proton::GetCode();
- const auto [prodCrossSection, elaCrossSection, dum] = fHadronicInteraction.GetCrossSection(protonId,protonId,EcmNN);
- const double sigProd = prodCrossSection / 1_mbarn;
+ const auto [prodCrossSection, elaCrossSection, dum] =
+ fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
+ const double sigProd = prodCrossSection / 1_mbarn;
const double sigEla = elaCrossSection / 1_mbarn;
// sample number of interactions (only input variables, output in common cnucms)
// nuclear multiple scattering according to glauber (r.i.p.)
- int_nuc_( kATarget, kAProj, sigProd, sigEla);
+ int_nuc_(kATarget, kAProj, sigProd, sigEla);
cout << "number of nucleons in target : " << kATarget << endl
- << "number of wounded nucleons in target : " << cnucms_.na << endl
- << "number of nucleons in projectile : " << kAProj << endl
- << "number of wounded nucleons in project. : " << cnucms_.nb << endl
- << "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl
- << "number of elastic nucleons in target : " << cnucms_.nael << endl
- << "number of elastic nucleons in project. : " << cnucms_.nbel << endl
- << "impact parameter: " << cnucms_.b << endl;
-
+ << "number of wounded nucleons in target : " << cnucms_.na << endl
+ << "number of nucleons in projectile : " << kAProj << endl
+ << "number of wounded nucleons in project. : " << cnucms_.nb << endl
+ << "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl
+ << "number of elastic nucleons in target : " << cnucms_.nael << endl
+ << "number of elastic nucleons in project. : " << cnucms_.nbel << endl
+ << "impact parameter: " << cnucms_.b << endl;
+
// calculate fragmentation
cout << "calculating nuclear fragments.." << endl;
// number of interactions
@@ -364,117 +374,114 @@ namespace corsika::process::sibyll {
const int nIntProj = nInelNucleons + nElasticNucleons;
const double impactPar = cnucms_.b; // only needed to avoid passing common var.
int nFragments;
- // number of fragments is limited to 60
+ // number of fragments is limited to 60
int AFragments[60];
// call fragmentation routine
- // input: target A, projectile A, number of int. nucleons in projectile, impact parameter (fm)
- // output: nFragments, AFragments
- // in addition the momenta ar stored in pf in common fragments, neglected
+ // input: target A, projectile A, number of int. nucleons in projectile, impact
+ // parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
+ // in pf in common fragments, neglected
fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :)
- if(nFragments>60)
- throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
+ if (nFragments > 60)
+ throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
cout << "number of fragments: " << nFragments << endl;
- for(int j=0; j<nFragments; ++j)
- cout << "fragment: " << j << " A=" << AFragments[j]
- << " px=" << fragments_.ppp[j][0]
- << " py=" << fragments_.ppp[j][1]
- << " pz=" << fragments_.ppp[j][2]
- << endl;
-
-
- auto Nucleus = []( int iA ){
- // int znew = iA / 2.15 + 0.7;
- corsika::particles::Code pCode;
- switch( iA ){
- case 2:
- pCode = corsika::particles::Deuterium::GetCode();
- break;
-
- case 3:
- pCode = corsika::particles::Tritium::GetCode();
- break;
-
- case 4:
- pCode = corsika::particles::Helium::GetCode();
- break;
-
- case 12:
- pCode = corsika::particles::Carbon::GetCode();
- break;
-
- case 14:
- pCode = corsika::particles::Nitrogen::GetCode();
- break;
-
- case 16:
- pCode = corsika::particles::Oxygen::GetCode();
- break;
-
- case 33:
- pCode = corsika::particles::Sulphur::GetCode();
- break;
-
- case 40:
- pCode = corsika::particles::Argon::GetCode();
- break;
-
- default:
- pCode = corsika::particles::Proton::GetCode();
- }
- return pCode;
- };
-
+ for (int j = 0; j < nFragments; ++j)
+ cout << "fragment: " << j << " A=" << AFragments[j]
+ << " px=" << fragments_.ppp[j][0] << " py=" << fragments_.ppp[j][1]
+ << " pz=" << fragments_.ppp[j][2] << endl;
+
+ auto Nucleus = [](int iA) {
+ // int znew = iA / 2.15 + 0.7;
+ corsika::particles::Code pCode;
+ switch (iA) {
+ case 2:
+ pCode = corsika::particles::Deuterium::GetCode();
+ break;
+
+ case 3:
+ pCode = corsika::particles::Tritium::GetCode();
+ break;
+
+ case 4:
+ pCode = corsika::particles::Helium::GetCode();
+ break;
+
+ case 12:
+ pCode = corsika::particles::Carbon::GetCode();
+ break;
+
+ case 14:
+ pCode = corsika::particles::Nitrogen::GetCode();
+ break;
+
+ case 16:
+ pCode = corsika::particles::Oxygen::GetCode();
+ break;
+
+ case 33:
+ pCode = corsika::particles::Sulphur::GetCode();
+ break;
+
+ case 40:
+ pCode = corsika::particles::Argon::GetCode();
+ break;
+
+ default:
+ pCode = corsika::particles::Proton::GetCode();
+ }
+ return pCode;
+ };
+
// put nuclear fragments on corsika stack
- for(int j=0; j<nFragments; ++j){
- // here we need the nucleonNumber to corsika Id conversion
- // A = AFragments[j]
- // pnew.SetPID( corsika::particles::GetCode( corsika::particles::Nucleus(A,Z) ) );
- auto pCode = Nucleus( AFragments[j] );
-
- // CORSIKA 7 way
- // spectators inherit momentum from original projectile
- const double mass_ratio = corsika::particles::GetMass( pCode ) / corsika::particles::GetMass( corsikaProjId );
- auto const Plab = PprojLab * mass_ratio;
-
- p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+ for (int j = 0; j < nFragments; ++j) {
+ // here we need the nucleonNumber to corsika Id conversion
+ auto pCode = Nucleus(AFragments[j]);
+
+ // CORSIKA 7 way
+ // spectators inherit momentum from original projectile
+ const double mass_ratio = corsika::particles::GetMass(pCode) /
+ corsika::particles::GetMass(corsikaProjId);
+ auto const Plab = PprojLab * mass_ratio;
+
+ p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
+ pOrig, tOrig);
}
// add elastic nucleons to corsika stack
// TODO: the elastic interaction could be external like the inelastic interaction,
// e.g. use existing ElasticModel
- for(int j=0; j<nElasticNucleons; ++j){
- // TODO: sample proton or neutron
- auto pCode = corsika::particles::Proton::GetCode();
-
- // CORSIKA 7 way
- // elastic nucleons inherit momentum from original projectile
- // neglecting momentum transfer in interaction
- const double mass_ratio = corsika::particles::GetMass( pCode ) / corsika::particles::GetMass( corsikaProjId );
- auto const Plab = PprojLab * mass_ratio;
+ for (int j = 0; j < nElasticNucleons; ++j) {
+ // TODO: sample proton or neutron
+ auto pCode = corsika::particles::Proton::GetCode();
+
+ // CORSIKA 7 way
+ // elastic nucleons inherit momentum from original projectile
+ // neglecting momentum transfer in interaction
+ const double mass_ratio = corsika::particles::GetMass(pCode) /
+ corsika::particles::GetMass(corsikaProjId);
+ auto const Plab = PprojLab * mass_ratio;
+
+ p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
+ pOrig, tOrig);
+ }
- p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+ // add inelastic interactions
+ for (int j = 0; j < nInelNucleons; ++j) {
+
+ // TODO: sample neutron or proton
+ auto pCode = corsika::particles::Proton::GetCode();
+ // temporarily add to stack, will be removed after interaction in DoInteraction
+ auto inelasticNucleon =
+ p.AddSecondary(pCode, PprojNucLab.GetTimeLikeComponent(),
+ PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig);
+ // create inelastic interaction
+ fHadronicInteraction.DoInteraction(inelasticNucleon, s);
}
- // create dummy stack, to store nucleon projectile
- // Stack inelasticNucleons;
- // add one nucleon per inelastic interaction
- for(int j=0; j<nInelNucleons; ++j){
-
- // TODO: sample neutron or proton
- auto pCode = corsika::particles::Proton::GetCode();
- // temporarily add to stack, will be removed after interaction
- auto inelasticNucleon = p.AddSecondary( pCode, PprojNucLab.GetTimeLikeComponent(),
- PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig );
- // create inelastic interaction
- fHadronicInteraction.DoInteraction( inelasticNucleon, s);
- }
-
// delete parent particle
p.Delete();
- // throw std::runtime_error(" stop here");
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Sibyll/nuclib.h b/Processes/Sibyll/nuclib.h
index 36ab4c0c070889a5a0e6893471751d381ed40c68..8653b2cc2aae811f330cd54f1d058456da4742f3 100644
--- a/Processes/Sibyll/nuclib.h
+++ b/Processes/Sibyll/nuclib.h
@@ -14,42 +14,37 @@
extern "C" {
- // nuclib common, NUClear Multiple Scattering
- /*
- COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
- + ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
- + ,JJAEL(IAMAX), JJBEL(IAMAX)
- */
-
- extern struct {
- double b, bmax;
- int ntry, na, nb, ni, nael, nbel;
- int jja[56], jjb[56], jjint[56][56], jjael[56], jjbel[56];
- } cnucms_;
-
- /*
- nuclib common, nuclear FRAGMENTS
-
- COMMON /FRAGMENTS/ PPP(3,60)
- */
- extern struct {
- double ppp[60][3];
- } fragments_;
-
-
- // NUCLIB
-
- // subroutine to initiate nuclib
- void nuc_nuc_ini_();
-
- // subroutine to sample nuclear interaction structure
- void int_nuc_( const int&, const int&, const double&, const double&);
-
- // subroutine to sample nuclear fragments
- void fragm_(const int&, const int&, const int&, const double&, int&, int*);
-
- void signuc_(const int&, const double&, double&);
-
+// nuclib common, NUClear Multiple Scattering
+/*
+ COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
+ + ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
+ + ,JJAEL(IAMAX), JJBEL(IAMAX)
+ */
+
+extern struct {
+ double b, bmax;
+ int ntry, na, nb, ni, nael, nbel;
+ int jja[56], jjb[56], jjint[56][56], jjael[56], jjbel[56];
+} cnucms_;
+
+/*
+ nuclib common, nuclear FRAGMENTS
+
+ COMMON /FRAGMENTS/ PPP(3,60)
+*/
+extern struct { double ppp[60][3]; } fragments_;
+
+// NUCLIB
+
+// subroutine to initiate nuclib
+void nuc_nuc_ini_();
+
+// subroutine to sample nuclear interaction structure
+void int_nuc_(const int&, const int&, const double&, const double&);
+
+// subroutine to sample nuclear fragments
+void fragm_(const int&, const int&, const int&, const double&, int&, int*);
+
+void signuc_(const int&, const double&, double&);
}
#endif
-
diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
index 61268b43dd0e5486775756dbb726b457968a88ee..756a94bd2f05d6b14eaada201b4d75e4476c0f9b 100644
--- a/Processes/Sibyll/sibyll2.3c.h
+++ b/Processes/Sibyll/sibyll2.3c.h
@@ -70,7 +70,6 @@ extern struct {
int lun;
} s_debug_;
-
// lund random generator setup
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
@@ -80,7 +79,6 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
-
// subroutine to SET DECAYS
void dec_ini_();
@@ -96,7 +94,7 @@ void decsib_();
// interaction length
// double fpni_(double&, int&);
- void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&, double&);
+void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&, double&);
void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
double&);
@@ -107,7 +105,5 @@ double get_sibyll_mass2(int&);
// phojet random generator setup
void pho_rndin_(int&, int&, int&, int&);
-
-
}
#endif
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 067a06e481d733322266f4b5654dc1a33eddd737..6251aab11ed5efc5d8b620824ba2ba45847cfff5 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -136,12 +136,12 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
- const HEPEnergyType E0 = 100_GeV;
+ const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
-
+
auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
Interaction hmodel(env);