IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 70bc86d8 authored by Felix Riehn's avatar Felix Riehn
Browse files

removed explicit call to sibyll from nuclear interaction, added nucleus air cross section

parent 5a21107e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 4
− 2
View file @ 70bc86d8
......@@ -244,9 +244,10 @@ int main() {
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
corsika::process::sibyll::Interaction sibyll(env);
corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
// corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
corsika::process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
corsika::process::sibyll::Decay decay;
ProcessCut cut(200_GeV);
ProcessCut cut(20_GeV);
// corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// corsika::process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env);
......@@ -281,6 +282,7 @@ int main() {
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
......
Processes/Sibyll/Interaction.h
+ 19
− 12
View file @ 70bc86d8
......@@ -31,7 +31,8 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
bool fInitialized = false;
public:
Interaction(corsika::environment::Environment const& env)
: fEnvironment(env) {}
......@@ -47,32 +48,38 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize Sibyll
sibyll_ini_();
if(!fInitialized){
sibyll_ini_();
fInitialized = true;
}
}
std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
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, dummy, dum1, dum2, dum3, dum4;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
const double dEcm = CoMenergy / 1_GeV;
int iTarget = -1;
if (corsika::particles::IsNucleus(TargetId)) {
const int iTarget = corsika::particles::GetNucleusA(TargetId);
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.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
return std::make_tuple(sigProd * 1_mbarn, iTarget);
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == corsika::particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
return std::make_tuple(sigProd * 1_mbarn, 1);
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
iTarget = 1;
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
return std::make_tuple(sigProd * 1_mbarn, 0);
sigEla = std::numeric_limits<double>::infinity();
}
return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
}
template <typename Particle, typename Track>
......@@ -139,7 +146,7 @@ namespace corsika::process::sibyll {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
auto const [productionCrossSection, numberOfNucleons] =
auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
std::cout << "Interaction: "
......@@ -266,7 +273,7 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, 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
......
Processes/Sibyll/NuclearInteraction.h
+ 50
− 116
View file @ 70bc86d8
......@@ -17,9 +17,6 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/process/sibyll/nuclib.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
......@@ -33,10 +30,12 @@ namespace corsika::process::sibyll {
int fNucCount = 0;
public:
NuclearInteraction(corsika::environment::Environment const& env)
: fEnvironment(env) {}
NuclearInteraction(corsika::environment::Environment const& env, corsika::process::sibyll::Interaction& hadint)
: fEnvironment(env),
fHadronicInteraction(hadint)
{}
~NuclearInteraction() {
std::cout << "Sibyll::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
std::cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
<< std::endl;
}
......@@ -47,45 +46,28 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize hadronic interaction module
//sibyll_ini_();
// TODO: save to run multiple initializations?
if(!fHadronicInteraction.WasInitialized())
fHadronicInteraction.Init();
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
}
std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
// 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) {
using namespace corsika::units::si;
double sigProd, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
double sigProd;
std::cout << "NuclearInteraction: GetCrossSection: called with: beamId= "<<BeamId
<< " TargetId= " << TargetId << " CoMenergy= " << CoMenergy / 1_GeV << std::endl;
if(!corsika::particles::IsNucleus(BeamId) ){
// ask sibyll for hadron cross section
// this is needed for sample target, which uses the nucleon cross section
// TODO: generalize to parent hadronic interaction
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
const double dEcm = CoMenergy / 1_GeV;
// check target, hadron or nucleus?
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.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
return std::make_tuple(sigProd * 1_mbarn, iTarget);
} else if (TargetId == corsika::particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
return std::make_tuple(sigProd * 1_mbarn, 1);
} else {
throw std::runtime_error("GetCrossSection: no interaction in sibyll possible");
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
return std::make_tuple(sigProd * 1_mbarn, 0);
}
// ask hadronic interaction model for hadron cross section
return fHadronicInteraction.GetCrossSection( BeamId, TargetId, CoMenergy);
}
// use nuclib to calc. nuclear cross sections
......@@ -101,7 +83,9 @@ namespace corsika::process::sibyll {
const double dEcm = CoMenergy / 1_GeV;
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
if (corsika::particles::IsNucleus(TargetId)) {
const int iTarget = corsika::particles::GetNucleusA(TargetId);
if (iTarget > 18 || iTarget == 0)
......@@ -110,8 +94,9 @@ namespace corsika::process::sibyll {
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, iTarget);
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);
}
template <typename Particle, typename Track>
......@@ -134,7 +119,6 @@ namespace corsika::process::sibyll {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
......@@ -180,12 +164,12 @@ namespace corsika::process::sibyll {
for (auto const targetId : mediumComposition.GetComponents()) {
i++;
cout << "NuclearInteraction: get interaction length for target: " << targetId << endl;
auto const [productionCrossSection, numberOfNucleons] =
// 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: sibyll return (mb): "
<< "IntLength: nuclib return (mb): "
<< productionCrossSection / 1_mbarn << std::endl;
weightedProdCrossSection += w[i] * productionCrossSection;
avgTargetMassNumber += w[i] * numberOfNucleons;
......@@ -208,7 +192,7 @@ namespace corsika::process::sibyll {
}
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&) {
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
......@@ -240,7 +224,7 @@ namespace corsika::process::sibyll {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
// position and time of interaction, not used in NUCLIB
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
......@@ -275,7 +259,6 @@ namespace corsika::process::sibyll {
// define target
// always a nucleon
// for Sibyll is always a single nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
// target is always at rest
......@@ -329,7 +312,7 @@ namespace corsika::process::sibyll {
auto const targetId = compVec[i];
cout << "target component: " << targetId << endl;
cout << "beam id: " << targetId << endl;
const auto [sigProd, nNuc] = GetCrossSection(beamId,targetId,EcmNN);
const auto [sigProd, sigEla, nNuc] = GetCrossSection(beamId,targetId,EcmNN);
cross_section_of_components[i] = sigProd;
}
......@@ -344,7 +327,7 @@ namespace corsika::process::sibyll {
int kATarget = -1;
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
if (targetCode == corsika::particles::Proton::GetCode()) kATarget = 1;
cout << "NuclearInteraction: sibyll target code: " << kATarget << endl;
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 "
......@@ -355,19 +338,12 @@ namespace corsika::process::sibyll {
cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "<< endl;
// get nucleon-nucleon cross section
// (needed to determine number of scatterings)
//
// TODO: this is an explicit dependence on sibyll
// should be changed to a call to GetCrossSection() of the had. int. implementation
// this also means GetCrossSection has to return the elastic cross section as well
double sigProd, sigEla, dum1, dum2, dum3;
double dumdif[3];
const double sqsnuc = EcmNN / 1_GeV;
const int iBeam = 1;
// read hadron-proton cross section table (input: hadron id, CoMenergy, output: cross sections)
sib_sigma_hp_(iBeam, sqsnuc, dum1, sigEla, sigProd, dumdif, dum2, dum3);
const auto protonId = corsika::particles::Proton::GetCode();
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.)
// nuclear multiple scattering according to glauber (r.i.p.)
int_nuc_( kATarget, kAProj, sigProd, sigEla);
cout << "number of nucleons in target : " << kATarget << endl
......@@ -408,10 +384,6 @@ namespace corsika::process::sibyll {
<< " pz=" << fragments_.ppp[j][2]
<< endl;
// bookeeping accross nucleon-nucleon interactions
MomentumVector Plab_all(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_all = 0_GeV;
auto Nucleus = []( int iA ){
// int znew = iA / 2.15 + 0.7;
......@@ -468,12 +440,11 @@ namespace corsika::process::sibyll {
auto const Plab = PprojLab * mass_ratio;
p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
}
// 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();
......@@ -485,70 +456,33 @@ namespace corsika::process::sibyll {
auto const Plab = PprojLab * mass_ratio;
p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
}
// add inelastic interactions
for(int j=0; j<nInelNucleons; ++j){
// create nucleon-nucleus inelastic interaction
// TODO: switch to DoInteraction() of had. interaction implementation
// for now use SIBYLL directly
const int kBeamCode = process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode());
cout << "creating interaction no. "<< j << endl;
sibyll_( kBeamCode, kATarget, sqsnuc );
decsib_();
// print final state
int print_unit = 6;
sib_list_(print_unit);
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// skip particles that have decayed in Sibyll
if (psib.HasDecayed()) continue;
// // transform energy to lab. frame
auto const pCoM = psib.GetMomentum();
HEPEnergyType const eCoM = psib.GetEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
auto pnew = p.AddSecondary(process::sibyll::ConvertFromSibyll(psib.GetPID()), Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
// create dummy stack, to store nucleon projectile
Stack inelasticNucleons;
// add one nucleon per inelastic interaction
for(int j=0; j<nInelNucleons; ++j){
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
}
cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV
<< endl
<< "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
<< endl;
// TODO: sample neutron or proton
auto pCode = corsika::particles::Proton::GetCode();
inelasticNucleons.AddParticle( pCode, PprojNucLab.GetTimeLikeComponent(), PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig );
}
cout << "accross all nucleon interactions: Etot lab: " << Elab_all / 1_GeV << endl
<< "accross all nucleon interactions: Ptot lab: " << (Plab_all / 1_GeV).GetComponents()
<< endl;
// delete current particle
// process stack of inelastic nucleons
for(auto& nucl : inelasticNucleons)
// create inelastic nucleon-nucleon interaction
fHadronicInteraction.DoInteraction( nucl, s);
// delete parent particle
p.Delete();
//throw std::runtime_error(" stop here");
// throw std::runtime_error(" stop here");
return process::EProcessReturn::eOk;
}
private:
corsika::environment::Environment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
};
......
Processes/Sibyll/sibyll2.3c.f
+ 6
− 4
View file @ 70bc86d8
......@@ -5516,7 +5516,7 @@ c external type declarations
 
c local type declarations
DOUBLE PRECISION SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO,
& SIGPROD,SIGBDIF,S_RNDM,S,PF,PB,PD,P0,P1,P2,R
& SIGPROD,SIGBDIF,SIGELA,S_RNDM,S,PF,PB,PD,P0,P1,P2,R
DIMENSION SIGDIF(3)
INTEGER K
SAVE
......@@ -5534,7 +5534,7 @@ c read hadron-nucleon cross section from table
c distinguish between nuclear cross sections..
IF(IAFLG.eq.0)THEN
c if target is nucleus calc. hadron-nucleus cross section (slow)
CALL SIB_SIGMA_HNUC(L,IA,SQS,SIGprod,SIGbdif)
CALL SIB_SIGMA_HNUC(L,IA,SQS,SIGprod,SIGbdif,SIGela)
ELSE
c if target is air read hadron-air cross section from table
CALL SIB_SIGMA_HAIR(L,SQS,SIGprod,SIGbdif)
......@@ -14010,7 +14010,7 @@ c internal
END
C=======================================================================
 
SUBROUTINE SIB_SIGMA_HNUC (L,IAT,SQS,SIGprod,SIGbdif)
SUBROUTINE SIB_SIGMA_HNUC (L,IAT,SQS,SIGprod,SIGbdif,SIGela)
 
C-----------------------------------------------------------------------
C calculate Hadron-nucleus cross sections
......@@ -14061,7 +14061,7 @@ C--------------------------------------------------------------------
+ SIGQSD,SIGPPT,SIGPPEL,SIGPPSD,ITG
 
c external
DOUBLE PRECISION SQS,SIGPROD,SIGBDIF
DOUBLE PRECISION SQS,SIGPROD,SIGBDIF,SIGELA
INTEGER L,IAT
 
c internal
......@@ -14092,6 +14092,8 @@ C particle production cross section
SIGprod = SIGT-SIGQE
C quasi elastic + elastic singl. diff cross section
SIGbdif = SIGQSD
c elastic cross section
SIGela = SIGel
if(ndebug.gt.0)THEN
WRITE(LUN,'(1X,A,3F8.2)')
& 'SIB_SIGMA_HNUC: SIGprod, SIGbdif, ALAM:',
......
Processes/Sibyll/sibyll2.3c.h
+ 1
− 1
View file @ 70bc86d8
......@@ -96,7 +96,7 @@ void decsib_();
// interaction length
// double fpni_(double&, int&);
void sib_sigma_hnuc_(const int&, const int&, const 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&);
......
Processes/Sibyll/testSibyll.cc
+ 2
− 21
View file @ 70bc86d8
......@@ -144,27 +144,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
NuclearInteraction model(env);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
SECTION("DecayInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 10_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);
NuclearInteraction model(env);
Interaction hmodel(env);
NuclearInteraction model(env, hmodel);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum