IAP GITLAB

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

Merge branch '128-add-process-nuclearinteraction' into 'master' 

Resolve "add sibyll process NuclearInteraction"

Closes #128

See merge request !65
parents d3c55725 93aec428
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 4474 additions and 3579 deletions
Documentation/Examples/cascade_example.cc
+ 25
− 6
View file @ 09a4bb14
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
#include <corsika/process/sibyll/Decay.h> #include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h> #include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/track_writer/TrackWriter.h> #include <corsika/process/track_writer/TrackWriter.h>
...@@ -124,6 +125,14 @@ public: ...@@ -124,6 +125,14 @@ public:
is_inv = true; is_inv = true;
break; break;
case Code::Neutron:
is_inv = true;
break;
case Code::AntiNeutron:
is_inv = true;
break;
default: default:
break; break;
} }
...@@ -235,16 +244,20 @@ int main() { ...@@ -235,16 +244,20 @@ int main() {
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
corsika::process::sibyll::Interaction sibyll(env); corsika::process::sibyll::Interaction sibyll(env);
// corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
corsika::process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
corsika::process::sibyll::Decay decay; corsika::process::sibyll::Decay decay;
ProcessCut cut(8_GeV); ProcessCut cut(20_GeV);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); // 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"); corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list // 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() // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n"; // << "\n";
...@@ -252,12 +265,17 @@ int main() { ...@@ -252,12 +265,17 @@ int main() {
// setup particle stack, and add primary particle // setup particle stack, and add primary particle
setup::Stack stack; setup::Stack stack;
stack.Clear(); stack.Clear();
const Code beamCode = Code::Proton;
const HEPEnergyType E0 = const HEPEnergyType E0 =
100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later) 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.; double theta = 0.;
double phi = 0.; double phi = 0.;
{ {
HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass()); 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) { auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi), return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta)); -ptot * cos(theta));
...@@ -265,10 +283,11 @@ int main() { ...@@ -265,10 +283,11 @@ int main() {
auto const [px, py, pz] = auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0); momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz}); 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 angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl; cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m); Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(Code::Proton, E0, plab, pos, 0_ns); stack.AddParticle(beamCode, E0, plab, pos, 0_ns);
} }
// define air shower object, run simulation // define air shower object, run simulation
......
Framework/Particles/NuclearData.xml
+ 3
− 0
View file @ 09a4bb14
...@@ -45,6 +45,9 @@ ...@@ -45,6 +45,9 @@
<particle id="1000100220" name="neon" A="22" Z="10" > <particle id="1000100220" name="neon" A="22" Z="10" >
</particle> </particle>
<particle id="1000160330" name="sulphur" A="33" Z="16" >
</particle>
<particle id="1000180400" name="argon" A="40" Z="18" > <particle id="1000180400" name="argon" A="40" Z="18" >
</particle> </particle>
......
Processes/Sibyll/CMakeLists.txt
+ 4
− 0
View file @ 09a4bb14
...@@ -16,6 +16,8 @@ set ( ...@@ -16,6 +16,8 @@ set (
MODEL_SOURCES MODEL_SOURCES
ParticleConversion.cc ParticleConversion.cc
sibyll2.3c.f sibyll2.3c.f
nuclib.f
signuc.f
sibyll2.3c.cc sibyll2.3c.cc
gasdev.f gasdev.f
) )
...@@ -24,9 +26,11 @@ set ( ...@@ -24,9 +26,11 @@ set (
MODEL_HEADERS MODEL_HEADERS
ParticleConversion.h ParticleConversion.h
sibyll2.3c.h sibyll2.3c.h
nuclib.h
SibStack.h SibStack.h
Decay.h Decay.h
Interaction.h Interaction.h
NuclearInteraction.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
) )
......
Processes/Sibyll/Interaction.h
+ 30
− 13
View file @ 09a4bb14
...@@ -31,6 +31,7 @@ namespace corsika::process::sibyll { ...@@ -31,6 +31,7 @@ namespace corsika::process::sibyll {
int fCount = 0; int fCount = 0;
int fNucCount = 0; int fNucCount = 0;
bool fInitialized = false;
public: public:
Interaction(corsika::environment::Environment const& env) Interaction(corsika::environment::Environment const& env)
...@@ -47,32 +48,40 @@ namespace corsika::process::sibyll { ...@@ -47,32 +48,40 @@ namespace corsika::process::sibyll {
using std::endl; using std::endl;
// initialize Sibyll // initialize Sibyll
sibyll_ini_(); if (!fInitialized) {
sibyll_ini_();
fInitialized = true;
}
} }
std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection( bool WasInitialized() { return fInitialized; }
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; using namespace corsika::units::si;
double sigProd, dummy, dum1, dum2, dum3, dum4; double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3]; double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId); const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
const double dEcm = CoMenergy / 1_GeV; const double dEcm = CoMenergy / 1_GeV;
int iTarget = -1;
if (corsika::particles::IsNucleus(TargetId)) { if (corsika::particles::IsNucleus(TargetId)) {
const int iTarget = corsika::particles::GetNucleusA(TargetId); iTarget = corsika::particles::GetNucleusA(TargetId);
if (iTarget > 18 || iTarget == 0) if (iTarget > 18 || iTarget == 0)
throw std::runtime_error( throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed."); "Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy); sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
return std::make_tuple(sigProd * 1_mbarn, iTarget);
} else if (TargetId == corsika::particles::Proton::GetCode()) { } else if (TargetId == corsika::particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4); sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
return std::make_tuple(sigProd * 1_mbarn, 1); iTarget = 1;
} else { } else {
// no interaction in sibyll possible, return infinite cross section? or throw? // no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity(); 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> template <typename Particle, typename Track>
...@@ -139,7 +148,7 @@ namespace corsika::process::sibyll { ...@@ -139,7 +148,7 @@ namespace corsika::process::sibyll {
i++; i++;
cout << "Interaction: get interaction length for target: " << targetId << endl; cout << "Interaction: get interaction length for target: " << targetId << endl;
auto const [productionCrossSection, numberOfNucleons] = auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
GetCrossSection(corsikaBeamId, targetId, ECoM); GetCrossSection(corsikaBeamId, targetId, ECoM);
std::cout << "Interaction: " std::cout << "Interaction: "
...@@ -185,6 +194,12 @@ namespace corsika::process::sibyll { ...@@ -185,6 +194,12 @@ namespace corsika::process::sibyll {
cout << "ProcessSibyll: " cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? " << "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl; << 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 (process::sibyll::CanInteract(corsikaBeamId)) { if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS = const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
...@@ -260,7 +275,8 @@ namespace corsika::process::sibyll { ...@@ -260,7 +275,8 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) { for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[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; cross_section_of_components[i] = sigProd;
int ideleteme = nNuc; // to avoid not used warning in array binding int ideleteme = nNuc; // to avoid not used warning in array binding
ideleteme = ideleteme; // to avoid not used warning in array binding ideleteme = ideleteme; // to avoid not used warning in array binding
...@@ -293,6 +309,7 @@ namespace corsika::process::sibyll { ...@@ -293,6 +309,7 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: " std::cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. " << " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << std::endl; << "THIS IS AN ERROR" << std::endl;
throw std::runtime_error("energy too low for SIBYLL");
// p.Delete(); delete later... different process // p.Delete(); delete later... different process
} else { } else {
fCount++; fCount++;
......
Processes/Sibyll/NuclearInteraction.h 0 → 100644
+ 498
− 0
View file @ 09a4bb14
This diff is collapsed.
Processes/Sibyll/nuclib.f 0 → 100644
+ 3554
− 0
View file @ 09a4bb14
This diff is collapsed.
Processes/Sibyll/nuclib.h 0 → 100644
+ 50
− 0
View file @ 09a4bb14
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_nuclib_interface_h_
#define _include_nuclib_interface_h_
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&);
}
#endif
Processes/Sibyll/sibyll2.3c.f
+ 6
− 3557
View file @ 09a4bb14
This diff is collapsed.
Processes/Sibyll/sibyll2.3c.h
+ 5
− 2
View file @ 09a4bb14
...@@ -70,6 +70,7 @@ extern struct { ...@@ -70,6 +70,7 @@ extern struct {
int lun; int lun;
} s_debug_; } s_debug_;
// lund random generator setup // lund random generator setup
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_; // extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
...@@ -78,7 +79,7 @@ void sibyll_(const int&, const int&, const double&); ...@@ -78,7 +79,7 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll // subroutine to initiate sibyll
void sibyll_ini_(); void sibyll_ini_();
// subroutine to SET DECAYS // subroutine to SET DECAYS
void dec_ini_(); void dec_ini_();
...@@ -94,7 +95,7 @@ void decsib_(); ...@@ -94,7 +95,7 @@ void decsib_();
// interaction length // interaction length
// double fpni_(double&, int&); // 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&, void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
double&); double&);
...@@ -105,5 +106,7 @@ double get_sibyll_mass2(int&); ...@@ -105,5 +106,7 @@ double get_sibyll_mass2(int&);
// phojet random generator setup // phojet random generator setup
void pho_rndin_(int&, int&, int&, int&); void pho_rndin_(int&, int&, int&, int&);
} }
#endif #endif
Processes/Sibyll/signuc.f 0 → 100644
+ 276
− 0
View file @ 09a4bb14
*-- Author : D. HECK IK FZK KARLSRUHE 06/12/1996
C=======================================================================
SUBROUTINE SIGNUC( IA,E0,SSIGNUC )
C-----------------------------------------------------------------------
C SIG(MA) NUC(LEUS) INI(TIALIZATION) 2
C
C IN ANALOGY WITH SUBROUT. SIGNUC_INI OF SIBYLL PACKAGE
C THIS SUBROUT. RECEIVES IN INPUT E0 (GEV) ENERGY PER NUCLEON AND
C INTERPOLATES THE CROSS-SECTIONS FOR NUCLEI WITH A < 56.
C UPDATED FOR SIBYLL-2.3 Sept. 2015 by D. Heck
C THIS SUBROUTINE IS CALLED FROM SIBSIG.
C ARGUMENTS:
C IA = MASS NUMBER OF PROJECTILE NUCLEUS
C E0 = LAB.ENERGY/NUCLEON OF PROJECTILE (GEV)
C SSIGNUC= CROSS-SECTION OF NUCLEUS IA WITH AIR
C
C. ATTENTION: THE TABULATED CROSS-SECTIONS ARE OBTAINED WITH
C. NEW P-P CROSS-SECTIONS AS USED IN SIBYLL 2.1,
C. IN ADDITION FIELD DIMENSIONS CHANGED (RE 04/2000)
C-----------------------------------------------------------------------
IMPLICIT NONE
DIMENSION SIGMA(6,56), SIGQE(6,56)
DOUBLE PRECISION SIGMA,SIGQE
DIMENSION AA(6)
DOUBLE PRECISION AA,DA,AMIN,ABEAM,S1,S2,ASQS
DOUBLE PRECISION E0,SSIGNUC
INTEGER IA,J,JE,NE
DOUBLE PRECISION QUAD_INT
EXTERNAL QUAD_INT
SAVE
DATA NE /6/, AMIN /1.D0/, DA /1.D0/
DATA AA /1.D0,2.D0,3.D0,4.D0,5.D0,6.D0/
* DATA AVOG /6.0221367D-04/
* DATA ATARGET /14.514D0/ ! EFFECTIVE MASSS OF AIR
C DATA ON 'INELASTIC-PRODUCTION' NUCLEUS-AIR CROSS-SECTION
DATA (SIGMA(J, 2),J=1,6) /
&3.842D+02,4.287D+02,4.940D+02,5.887D+02,6.922D+02,7.767D+02/
DATA (SIGMA(J, 3),J=1,6) /
&4.601D+02,5.149D+02,5.595D+02,6.663D+02,7.641D+02,8.446D+02/
DATA (SIGMA(J, 4),J=1,6) /
&4.881D+02,5.373D+02,6.005D+02,6.895D+02,7.716D+02,8.967D+02/
DATA (SIGMA(J, 5),J=1,6) /
&5.874D+02,6.176D+02,7.181D+02,7.993D+02,9.089D+02,1.031D+03/
DATA (SIGMA(J, 6),J=1,6) /
&7.054D+02,7.399D+02,8.388D+02,9.463D+02,1.080D+03,1.197D+03/
DATA (SIGMA(J, 7),J=1,6) /
&7.192D+02,7.611D+02,8.449D+02,9.539D+02,1.061D+03,1.176D+03/
DATA (SIGMA(J, 8),J=1,6) /
&7.550D+02,7.975D+02,9.153D+02,9.944D+02,1.126D+03,1.236D+03/
DATA (SIGMA(J, 9),J=1,6) /
&7.929D+02,8.392D+02,9.265D+02,1.059D+03,1.167D+03,1.262D+03/
DATA (SIGMA(J, 10),J=1,6) /
&8.157D+02,8.644D+02,9.512D+02,1.058D+03,1.182D+03,1.298D+03/
DATA (SIGMA(J, 11),J=1,6) /
&8.039D+02,8.587D+02,9.534D+02,1.055D+03,1.182D+03,1.298D+03/
DATA (SIGMA(J, 12),J=1,6) /
&8.515D+02,8.957D+02,9.869D+02,1.122D+03,1.253D+03,1.366D+03/
DATA (SIGMA(J, 13),J=1,6) /
&8.769D+02,9.100D+02,1.018D+03,1.119D+03,1.252D+03,1.341D+03/
DATA (SIGMA(J, 14),J=1,6) /
&9.058D+02,9.532D+02,1.057D+03,1.171D+03,1.302D+03,1.391D+03/
DATA (SIGMA(J, 15),J=1,6) /
&9.555D+02,9.799D+02,1.098D+03,1.201D+03,1.342D+03,1.444D+03/
DATA (SIGMA(J, 16),J=1,6) /
&1.009D+03,1.058D+03,1.149D+03,1.290D+03,1.414D+03,1.520D+03/
DATA (SIGMA(J, 17),J=1,6) /
&9.907D+02,1.045D+03,1.166D+03,1.290D+03,1.384D+03,1.516D+03/
DATA (SIGMA(J, 18),J=1,6) /
&1.036D+03,1.121D+03,1.198D+03,1.328D+03,1.470D+03,1.592D+03/
DATA (SIGMA(J, 19),J=1,6) /
&1.083D+03,1.162D+03,1.250D+03,1.371D+03,1.516D+03,1.661D+03/
DATA (SIGMA(J, 20),J=1,6) /
&1.146D+03,1.215D+03,1.295D+03,1.443D+03,1.544D+03,1.744D+03/
DATA (SIGMA(J, 21),J=1,6) /
&1.158D+03,1.234D+03,1.292D+03,1.467D+03,1.618D+03,1.750D+03/
DATA (SIGMA(J, 22),J=1,6) /
&1.153D+03,1.205D+03,1.329D+03,1.451D+03,1.596D+03,1.734D+03/
DATA (SIGMA(J, 23),J=1,6) /
&1.210D+03,1.274D+03,1.356D+03,1.493D+03,1.655D+03,1.803D+03/
DATA (SIGMA(J, 24),J=1,6) /
&1.212D+03,1.273D+03,1.398D+03,1.489D+03,1.641D+03,1.800D+03/
DATA (SIGMA(J, 25),J=1,6) /
&1.236D+03,1.315D+03,1.423D+03,1.561D+03,1.669D+03,1.855D+03/
DATA (SIGMA(J, 26),J=1,6) /
&1.279D+03,1.345D+03,1.431D+03,1.595D+03,1.734D+03,1.889D+03/
DATA (SIGMA(J, 27),J=1,6) /
&1.228D+03,1.304D+03,1.438D+03,1.546D+03,1.714D+03,1.836D+03/
DATA (SIGMA(J, 28),J=1,6) /
&1.289D+03,1.370D+03,1.451D+03,1.597D+03,1.754D+03,1.913D+03/
DATA (SIGMA(J, 29),J=1,6) /
&1.411D+03,1.469D+03,1.613D+03,1.777D+03,1.910D+03,2.075D+03/
DATA (SIGMA(J, 30),J=1,6) /
&1.347D+03,1.401D+03,1.498D+03,1.642D+03,1.816D+03,1.975D+03/
DATA (SIGMA(J, 31),J=1,6) /
&1.359D+03,1.448D+03,1.551D+03,1.694D+03,1.858D+03,2.007D+03/
DATA (SIGMA(J, 32),J=1,6) /
&1.358D+03,1.460D+03,1.559D+03,1.698D+03,1.842D+03,1.974D+03/
DATA (SIGMA(J, 33),J=1,6) /
&1.418D+03,1.448D+03,1.578D+03,1.727D+03,1.872D+03,2.047D+03/
DATA (SIGMA(J, 34),J=1,6) /
&1.433D+03,1.466D+03,1.605D+03,1.738D+03,1.892D+03,2.019E+03/
DATA (SIGMA(J, 35),J=1,6) /
&1.430D+03,1.511D+03,1.602D+03,1.752D+03,1.935D+03,2.060D+03/
DATA (SIGMA(J, 36),J=1,6) /
&1.462D+03,1.499D+03,1.653D+03,1.805D+03,1.920D+03,2.057D+03/
DATA (SIGMA(J, 37),J=1,6) /
&1.470D+03,1.520D+03,1.656D+03,1.818D+03,1.946D+03,2.131D+03/
DATA (SIGMA(J, 38),J=1,6) /
&1.470D+03,1.542D+03,1.691D+03,1.800D+03,1.968D+03,2.133D+03/
DATA (SIGMA(J, 39),J=1,6) /
&1.495D+03,1.588D+03,1.676D+03,1.834D+03,1.969D+03,2.163D+03/
DATA (SIGMA(J, 40),J=1,6) /
&1.525D+03,1.551D+03,1.722D+03,1.833D+03,2.020D+03,2.192D+03/
DATA (SIGMA(J, 41),J=1,6) /
&1.526D+03,1.615D+03,1.709D+03,1.899D+03,2.040D+03,2.181D+03/
DATA (SIGMA(J, 42),J=1,6) /
&1.510D+03,1.567D+03,1.716D+03,1.892D+03,2.056D+03,2.197D+03/
DATA (SIGMA(J, 43),J=1,6) /
&1.557D+03,1.658D+03,1.776D+03,1.898D+03,2.092D+03,2.200D+03/
DATA (SIGMA(J, 44),J=1,6) /
&1.556D+03,1.645D+03,1.752D+03,1.920D+03,2.091D+03,2.243D+03/
DATA (SIGMA(J, 45),J=1,6) /
&1.583D+03,1.663D+03,1.798D+03,1.940D+03,2.051D+03,2.263D+03/
DATA (SIGMA(J, 46),J=1,6) /
&1.599D+03,1.642D+03,1.799D+03,1.941D+03,2.107D+03,2.268D+03/
DATA (SIGMA(J, 47),J=1,6) /
&1.611D+03,1.692D+03,1.811D+03,1.956D+03,2.107D+03,2.264D+03/
DATA (SIGMA(J, 48),J=1,6) /
&1.625D+03,1.706D+03,1.819D+03,1.986D+03,2.139D+03,2.354D+03/
DATA (SIGMA(J, 49),J=1,6) /
&1.666D+03,1.737D+03,1.854D+03,1.971D+03,2.160D+03,2.318D+03/
DATA (SIGMA(J, 50),J=1,6) /
&1.648D+03,1.747D+03,1.856D+03,2.023D+03,2.181D+03,2.352D+03/
DATA (SIGMA(J, 51),J=1,6) /
&1.653D+03,1.763D+03,1.868D+03,2.015D+03,2.203D+03,2.386D+03/
DATA (SIGMA(J, 52),J=1,6) /
&1.690D+03,1.720D+03,1.902D+03,2.027D+03,2.189D+03,2.357D+03/
DATA (SIGMA(J, 53),J=1,6) /
&1.690D+03,1.750D+03,1.921D+03,2.059D+03,2.208D+03,2.417D+03/
DATA (SIGMA(J, 54),J=1,6) /
&1.705D+03,1.781D+03,1.911D+03,2.073D+03,2.242D+03,2.411D+03/
DATA (SIGMA(J, 55),J=1,6) /
&1.714D+03,1.806D+03,1.896D+03,2.100D+03,2.253D+03,2.411D+03/
DATA (SIGMA(J, 56),J=1,6) /
&1.774D+03,1.813D+03,1.954D+03,2.098D+03,2.280D+03,2.482D+03/
C DATA ON 'QUASI-ELASTIC' NUCLEUS-AIR CROSS-SECTION
DATA (SIGQE(J, 2),J=1,6) /
&4.141D+01,3.708D+01,5.428D+01,8.696D+01,1.403D+02,1.885D+02/
DATA (SIGQE(J, 3),J=1,6) /
&4.357D+01,3.894D+01,5.177D+01,9.675D+01,1.447D+02,2.029D+02/
DATA (SIGQE(J, 4),J=1,6) /
&4.123D+01,3.933D+01,6.070D+01,9.482D+01,1.474D+02,2.023D+02/
DATA (SIGQE(J, 5),J=1,6) /
&4.681D+01,4.287D+01,6.381D+01,1.050D+02,1.519D+02,2.198D+02/
DATA (SIGQE(J, 6),J=1,6) /
&5.407D+01,5.195D+01,6.723D+01,1.108D+02,1.750D+02,2.368D+02/
DATA (SIGQE(J, 7),J=1,6) /
&4.975D+01,4.936D+01,6.880D+01,1.162D+02,1.689D+02,2.329D+02/
DATA (SIGQE(J, 8),J=1,6) /
&5.361D+01,5.027D+01,6.858D+01,1.177D+02,1.759D+02,2.412D+02/
DATA (SIGQE(J, 9),J=1,6) /
&4.980D+01,5.063D+01,7.210D+01,1.196D+02,1.806D+02,2.299D+02/
DATA (SIGQE(J, 10),J=1,6) /
&5.170D+01,5.070D+01,7.105D+01,1.182D+02,1.679D+02,2.411D+02/
DATA (SIGQE(J, 11),J=1,6) /
&4.950D+01,4.950D+01,7.286D+01,1.137D+02,1.769D+02,2.477D+02/
DATA (SIGQE(J, 12),J=1,6) /
&5.262D+01,5.133D+01,7.110D+01,1.204D+02,1.789D+02,2.501D+02/
DATA (SIGQE(J, 13),J=1,6) /
&5.320D+01,5.378D+01,6.847D+01,1.200D+02,1.805D+02,2.442D+02/
DATA (SIGQE(J, 14),J=1,6) /
&5.638D+01,5.271D+01,6.985D+01,1.209D+02,1.867D+02,2.610D+02/
DATA (SIGQE(J, 15),J=1,6) /
&5.294D+01,5.353D+01,7.435D+01,1.211D+02,1.899D+02,2.612D+02/
DATA (SIGQE(J, 16),J=1,6) /
&5.668D+01,5.254D+01,7.557D+01,1.269D+02,1.917D+02,2.707D+02/
DATA (SIGQE(J, 17),J=1,6) /
&5.456D+01,5.721D+01,7.481D+01,1.208D+02,1.859D+02,2.658D+02/
DATA (SIGQE(J, 18),J=1,6) /
&5.901D+01,5.382D+01,7.591D+01,1.246D+02,1.872D+02,2.874D+02/
DATA (SIGQE(J, 19),J=1,6) /
&6.328D+01,6.116D+01,8.451D+01,1.318D+02,2.088D+02,2.749D+02/
DATA (SIGQE(J, 20),J=1,6) /
&5.779D+01,5.924D+01,8.382D+01,1.370D+02,2.062D+02,2.837D+02/
DATA (SIGQE(J, 21),J=1,6) /
&7.155D+01,5.732D+01,8.231D+01,1.363D+02,2.047D+02,2.820D+02/
DATA (SIGQE(J, 22),J=1,6) /
&6.699D+01,5.651D+01,8.511D+01,1.477D+02,2.031D+02,2.921D+02/
DATA (SIGQE(J, 23),J=1,6) /
&6.179D+01,6.269D+01,9.395D+01,1.437D+02,2.195D+02,2.964D+02/
DATA (SIGQE(J, 24),J=1,6) /
&6.784D+01,6.028D+01,8.622D+01,1.279D+02,2.214D+02,2.867D+02/
DATA (SIGQE(J, 25),J=1,6) /
&6.589D+01,5.795D+01,8.890D+01,1.385D+02,2.055D+02,2.988D+02/
DATA (SIGQE(J, 26),J=1,6) /
&6.364D+01,6.325D+01,8.942D+01,1.421D+02,2.128D+02,3.083D+02/
DATA (SIGQE(J, 27),J=1,6) /
&6.449D+01,6.664D+01,8.986D+01,1.453D+02,2.140D+02,2.932D+02/
DATA (SIGQE(J, 28),J=1,6) /
&7.284D+01,6.139D+01,8.867D+01,1.425D+02,2.179D+02,2.978D+02/
DATA (SIGQE(J, 29),J=1,6) /
&7.221D+01,7.085D+01,9.079D+01,1.482D+02,2.277D+02,2.913D+02/
DATA (SIGQE(J, 30),J=1,6) /
&6.928D+01,6.294D+01,8.935D+01,1.463D+02,2.265D+02,2.834D+02/
DATA (SIGQE(J, 31),J=1,6) /
&6.611D+01,6.586D+01,9.133D+01,1.461D+02,2.201D+02,2.959D+02/
DATA (SIGQE(J, 32),J=1,6) /
&6.401D+01,6.177D+01,8.971D+01,1.480D+02,2.155D+02,3.152D+02/
DATA (SIGQE(J, 33),J=1,6) /
&7.057D+01,6.918D+01,8.410D+01,1.465D+02,2.288D+02,3.088D+02/
DATA (SIGQE(J, 34),J=1,6) /
&6.453D+01,7.020D+01,9.272D+01,1.517D+02,2.189D+02,2.999D+02/
DATA (SIGQE(J, 35),J=1,6) /
&6.741D+01,6.295D+01,9.323D+01,1.536D+02,2.190D+02,2.930D+02/
DATA (SIGQE(J, 36),J=1,6) /
&6.807D+01,7.046D+01,1.025D+02,1.565D+02,2.315D+02,3.090D+02/
DATA (SIGQE(J, 37),J=1,6) /
&8.082D+01,6.565D+01,9.160D+01,1.572D+02,2.229D+02,3.125D+02/
DATA (SIGQE(J, 38),J=1,6) /
&6.494D+01,6.964D+01,9.089D+01,1.653D+02,2.336D+02,3.120D+02/
DATA (SIGQE(J, 39),J=1,6) /
&6.833D+01,6.860D+01,8.933D+01,1.601D+02,2.261D+02,3.167D+02/
DATA (SIGQE(J, 40),J=1,6) /
&7.021D+01,6.866D+01,8.437D+01,1.588D+02,2.249D+02,2.941D+02/
DATA (SIGQE(J, 41),J=1,6) /
&7.122D+01,6.205D+01,9.545D+01,1.582D+02,2.335D+02,3.395D+02/
DATA (SIGQE(J, 42),J=1,6) /
&7.265D+01,6.936D+01,9.486D+01,1.505D+02,2.379D+02,3.248D+02/
DATA (SIGQE(J, 43),J=1,6) /
&7.048D+01,7.539D+01,9.192D+01,1.566D+02,2.532D+02,3.182D+02/
DATA (SIGQE(J, 44),J=1,6) /
&6.650D+01,7.139D+01,9.862D+01,1.602D+02,2.289D+02,3.077D+02/
DATA (SIGQE(J, 45),J=1,6) /
&7.511D+01,6.893D+01,9.245D+01,1.641D+02,2.519D+02,3.381D+02/
DATA (SIGQE(J, 46),J=1,6) /
&6.437D+01,6.894D+01,8.697D+01,1.544D+02,2.391D+02,3.213D+02/
DATA (SIGQE(J, 47),J=1,6) /
&7.980D+01,6.958D+01,1.022D+02,1.609D+02,2.408D+02,3.246D+02/
DATA (SIGQE(J, 48),J=1,6) /
&7.265D+01,7.313D+01,8.989D+01,1.578D+02,2.387D+02,3.235D+02/
DATA (SIGQE(J, 49),J=1,6) /
&6.959D+01,6.337D+01,9.084D+01,1.656D+02,2.331D+02,3.226D+02/
DATA (SIGQE(J, 50),J=1,6) /
&7.371D+01,6.807D+01,9.726D+01,1.535D+02,2.445D+02,3.189D+02/
DATA (SIGQE(J, 51),J=1,6) /
&7.882D+01,6.680D+01,9.377D+01,1.629D+02,2.448D+02,3.297D+02/
DATA (SIGQE(J, 52),J=1,6) /
&7.223D+01,6.794D+01,9.925D+01,1.738D+02,2.446D+02,3.162D+02/
DATA (SIGQE(J, 53),J=1,6) /
&7.703D+01,6.971D+01,9.601D+01,1.595D+02,2.484D+02,3.265D+02/
DATA (SIGQE(J, 54),J=1,6) /
&7.549D+01,7.459D+01,8.984D+01,1.645D+02,2.348D+02,3.201D+02/
DATA (SIGQE(J, 55),J=1,6) /
&7.891D+01,6.840D+01,1.017D+02,1.698D+02,2.501D+02,3.429D+02/
DATA (SIGQE(J, 56),J=1,6) /
&7.545D+01,6.673D+01,1.057D+02,1.684D+02,2.424D+02,3.181D+02/
C-----------------------------------------------------------------------
C ENERGY E0 IN GEV
ASQS = 0.5D0*LOG10(1.876D0*E0)
JE = MIN( INT( (ASQS-AMIN)/DA )+1, NE-2 )
ABEAM = DBLE(IA)
C INELASTIC CROSS-SECTION
S1 = QUAD_INT( ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ SIGMA(JE,IA),SIGMA(JE+1,IA),SIGMA(JE+2,IA) )
C QUASI ELASTIC CROSS-SECTION
S2 = QUAD_INT( ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ SIGQE(JE,IA),SIGQE(JE+1,IA),SIGQE(JE+2,IA) )
C ADD INELASTIC AND QUASI-ELASTIC CROSS-SECTIONS
SSIGNUC = S1 + S2
RETURN
END
Processes/Sibyll/testSibyll.cc
+ 23
− 1
View file @ 09a4bb14
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
#include <corsika/process/sibyll/Decay.h> #include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h> #include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/sibyll/ParticleConversion.h> #include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/random/RNGManager.h> #include <corsika/random/RNGManager.h>
...@@ -116,7 +117,7 @@ TEST_CASE("SibyllInterface", "[processes]") { ...@@ -116,7 +117,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("InteractionInterface") { SECTION("InteractionInterface") {
setup::Stack stack; setup::Stack stack;
const HEPEnergyType E0 = 10_GeV; const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 = HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass()); sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0}); auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
...@@ -132,6 +133,27 @@ TEST_CASE("SibyllInterface", "[processes]") { ...@@ -132,6 +133,27 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.GetInteractionLength(particle, track); model.GetInteractionLength(particle, track);
} }
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
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);
NuclearInteraction model(env, hmodel);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
SECTION("DecayInterface") { SECTION("DecayInterface") {
setup::Stack stack; setup::Stack stack;
......
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