IAP GITLAB
Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in / Register
Toggle navigation
Menu
Open sidebar
Air Shower Physics
corsika
Commits
322427c4
Commit
322427c4
authored
Oct 09, 2022
by
Felix Riehn
Browse files
clang
parent
c6c77fcc
Pipeline
#9424
passed with stages
in 12 minutes and 15 seconds
Changes
7
Pipelines
2
Hide whitespace changes
Inline
Side-by-side
corsika/detail/modules/sophia/InteractionModel.inl
View file @
322427c4
...
...
@@ -83,8 +83,8 @@ namespace corsika::sophia {
double
theta
=
0.0
;
// set nucleon at rest in collision
int
Imode
;
CORSIKA_LOGGER_DEBUG
(
logger_
,
"calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}"
,
nucleonSophiaCode
,
Enucleon
,
Ephoton
,
theta
);
"calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}"
,
nucleonSophiaCode
,
Enucleon
,
Ephoton
,
theta
);
count_
++
;
// call sophia
eventgen_
(
nucleonSophiaCode
,
Enucleon
,
Ephoton
,
theta
,
Imode
);
...
...
@@ -124,10 +124,10 @@ namespace corsika::sophia {
calculate_kinetic_energy
(
momentum
.
getNorm
(),
get_mass
(
pid
));
CORSIKA_LOGGER_TRACE
(
logger_
,
"SOPHIA: pid={}, p={} GeV"
,
pidSophia
,
momentumSophia
.
getComponents
()
/
1
_GeV
);
momentumSophia
.
getComponents
()
/
1
_GeV
);
CORSIKA_LOGGER_TRACE
(
logger_
,
"CORSIKA: pid={}, p={} GeV"
,
pid
,
momentum
.
getComponents
()
/
1
_GeV
);
momentum
.
getComponents
()
/
1
_GeV
);
auto
pnew
=
secondaries
.
addSecondary
(
std
::
make_tuple
(
pid
,
Ekin
,
momentum
.
normalized
()));
...
...
@@ -136,7 +136,7 @@ namespace corsika::sophia {
E_final
+=
pnew
.
getEnergy
();
}
CORSIKA_LOGGER_TRACE
(
logger_
,
"Efinal={} GeV,Pfinal={} GeV"
,
E_final
/
1
_GeV
,
P_final
.
getComponents
()
/
1
_GeV
);
P_final
.
getComponents
()
/
1
_GeV
);
}
}
// namespace corsika::sophia
\ No newline at end of file
corsika/modules/Sophia.hpp
View file @
322427c4
...
...
@@ -28,7 +28,8 @@ namespace corsika::sophia {
* The sophia::InteractionModel is wrapped as an InteractionProcess here in order
* to provide all the functions for ProcessSequence.
*/
// struct Interaction : public InteractionModel, public InteractionProcess<Interaction> {
// struct Interaction : public InteractionModel, public InteractionProcess<Interaction>
// {
// template <typename TEnvironment>
// Interaction(TEnvironment const& env)
// : InteractionModel{env} {}
...
...
corsika/modules/sophia/ParticleConversion.hpp
View file @
322427c4
...
...
@@ -29,7 +29,7 @@ namespace corsika::sophia {
// Pion = 2,
// Kaon = 3,
// };
//using SophiaXSClassIntType = std::underlying_type<SophiaXSClass>::type;
//
using SophiaXSClassIntType = std::underlying_type<SophiaXSClass>::type;
#include
<corsika/modules/sophia/Generated.inc>
...
...
examples/corsika.cpp
View file @
322427c4
...
...
@@ -82,8 +82,7 @@
using
namespace
corsika
;
using
namespace
std
;
using
EnvironmentInterface
=
IMediumPropertyModel
<
IMagneticFieldModel
<
IMediumModel
>>
;
using
EnvironmentInterface
=
IMediumPropertyModel
<
IMagneticFieldModel
<
IMediumModel
>>
;
using
EnvType
=
Environment
<
EnvironmentInterface
>
;
using
Particle
=
setup
::
Stack
<
EnvType
>::
particle_type
;
...
...
@@ -110,7 +109,7 @@ void registerRandomStreams(int seed) {
template
<
typename
T
>
using
MyExtraEnv
=
MediumPropertyModel
<
UniformMagneticField
<
T
>>
;
int
main
(
int
argc
,
char
**
argv
)
{
int
main
(
int
argc
,
char
**
argv
)
{
// the main command line description
CLI
::
App
app
{
"Simulate standard (downgoing) showers with CORSIKA 8."
};
...
...
@@ -171,8 +170,7 @@ int main(int argc, char **argv) {
CLI11_PARSE
(
app
,
argc
,
argv
);
if
(
app
.
count
(
"--verbosity"
))
{
std
::
string_view
const
loglevel
=
app
[
"--verbosity"
]
->
as
<
std
::
string_view
>
();
std
::
string_view
const
loglevel
=
app
[
"--verbosity"
]
->
as
<
std
::
string_view
>
();
if
(
loglevel
==
"warn"
)
{
logging
::
set_level
(
logging
::
level
::
warn
);
}
else
if
(
loglevel
==
"info"
)
{
...
...
@@ -193,8 +191,7 @@ int main(int argc, char **argv) {
// gets all messed up
if
(
app
.
count
(
"--pdg"
)
==
0
)
{
if
((
app
.
count
(
"-A"
)
==
0
)
||
(
app
.
count
(
"-Z"
)
==
0
))
{
CORSIKA_LOG_ERROR
(
"If --pdg is not provided, then both -A and -Z are required."
);
CORSIKA_LOG_ERROR
(
"If --pdg is not provided, then both -A and -Z are required."
);
return
1
;
}
}
...
...
@@ -204,7 +201,7 @@ int main(int argc, char **argv) {
/* === START: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
EnvType
env
;
CoordinateSystemPtr
const
&
rootCS
=
env
.
getCoordinateSystem
();
CoordinateSystemPtr
const
&
rootCS
=
env
.
getCoordinateSystem
();
Point
const
center
{
rootCS
,
0
_m
,
0
_m
,
0
_m
};
GeomagneticModel
wmm
(
center
,
corsika_data
(
"GeoMag/WMM.COF"
));
...
...
@@ -218,10 +215,9 @@ int main(int argc, char **argv) {
ofstream
atmout
(
"earth.dat"
);
for
(
LengthType
h
=
0
_m
;
h
<
110
_km
;
h
+=
100
_m
)
{
Point
const
ptest
{
rootCS
,
0
_m
,
0
_m
,
constants
::
EarthRadius
::
Mean
+
h
};
auto
rho
=
env
.
getUniverse
()
->
getContainingNode
(
ptest
)
->
getModelProperties
()
.
getMassDensity
(
ptest
);
auto
rho
=
env
.
getUniverse
()
->
getContainingNode
(
ptest
)
->
getModelProperties
().
getMassDensity
(
ptest
);
atmout
<<
h
/
1
_m
<<
" "
<<
rho
/
1
_kg
*
cube
(
1
_m
)
<<
"
\n
"
;
}
atmout
.
close
();
...
...
@@ -257,8 +253,8 @@ int main(int argc, char **argv) {
// convert the momentum to the zenith and azimuth angle of the primary
auto
const
[
px
,
py
,
pz
]
=
std
::
make_tuple
(
P0
*
sin
(
thetaRad
)
*
cos
(
phiRad
),
P0
*
sin
(
thetaRad
)
*
sin
(
phiRad
),
-
P0
*
cos
(
thetaRad
));
std
::
make_tuple
(
P0
*
sin
(
thetaRad
)
*
cos
(
phiRad
),
P0
*
sin
(
thetaRad
)
*
sin
(
phiRad
),
-
P0
*
cos
(
thetaRad
));
auto
plab
=
MomentumVector
(
rootCS
,
{
px
,
py
,
pz
});
/* === END: CONSTRUCT PRIMARY PARTICLE === */
...
...
@@ -270,16 +266,14 @@ int main(int argc, char **argv) {
static_pow
<
2
>
(
injectionHeight
));
Point
const
showerCore
{
rootCS
,
0
_m
,
0
_m
,
observationHeight
};
Point
const
injectionPos
=
showerCore
+
DirectionVector
{
rootCS
,
{
-
sin
(
thetaRad
)
*
cos
(
phiRad
),
-
sin
(
thetaRad
)
*
sin
(
phiRad
),
cos
(
thetaRad
)}}
*
t
;
showerCore
+
DirectionVector
{
rootCS
,
{
-
sin
(
thetaRad
)
*
cos
(
phiRad
),
-
sin
(
thetaRad
)
*
sin
(
phiRad
),
cos
(
thetaRad
)}}
*
t
;
// we make the axis much longer than the inj-core distance since the
// profile will go beyond the core, depending on zenith angle
ShowerAxis
const
showerAxis
{
injectionPos
,
(
showerCore
-
injectionPos
)
*
1.2
,
env
};
ShowerAxis
const
showerAxis
{
injectionPos
,
(
showerCore
-
injectionPos
)
*
1.2
,
env
};
/* === END: CONSTRUCT GEOMETRY === */
// create the output manager that we then register outputs with
...
...
@@ -326,8 +320,7 @@ int main(int argc, char **argv) {
HEPEnergyType
const
emcut
=
50
_GeV
;
HEPEnergyType
const
hadcut
=
50
_GeV
;
ParticleCut
<
SubWriter
<
decltype
(
dEdX
)
>>
cut
(
emcut
,
emcut
,
hadcut
,
hadcut
,
true
,
dEdX
);
ParticleCut
<
SubWriter
<
decltype
(
dEdX
)
>>
cut
(
emcut
,
emcut
,
hadcut
,
hadcut
,
true
,
dEdX
);
// energy threshold for high energy hadronic model. Affects LE/HE switch for
// hadron interactions and the hadronic photon model in proposal
...
...
@@ -351,13 +344,12 @@ int main(int argc, char **argv) {
// assemble all processes into an ordered process list
struct
EnergySwitch
{
HEPEnergyType
cutE_
;
EnergySwitch
(
HEPEnergyType
cutE
)
:
cutE_
(
cutE
)
{}
bool
operator
()(
const
Particle
&
p
)
const
{
return
(
p
.
getKineticEnergy
()
<
cutE_
);
}
EnergySwitch
(
HEPEnergyType
cutE
)
:
cutE_
(
cutE
)
{}
bool
operator
()(
const
Particle
&
p
)
const
{
return
(
p
.
getKineticEnergy
()
<
cutE_
);
}
};
auto
hadronSequence
=
make_select
(
EnergySwitch
(
heHadronModelThreshold
),
urqmdCounted
,
sibyllCounted
);
auto
hadronSequence
=
make_select
(
EnergySwitch
(
heHadronModelThreshold
),
urqmdCounted
,
sibyllCounted
);
auto
decaySequence
=
make_sequence
(
decayPythia
,
decaySibyll
);
TrackWriter
trackWriter
{
tracks
};
...
...
@@ -370,8 +362,8 @@ int main(int argc, char **argv) {
output
.
add
(
"particles"
,
observationLevel
);
// assemble the final process sequence
auto
sequence
=
make_sequence
(
stackInspect
,
hadronSequence
,
decaySequence
,
cut
,
emCascade
,
emContinuous
,
// trackWriter,
auto
sequence
=
make_sequence
(
stackInspect
,
hadronSequence
,
decaySequence
,
cut
,
emCascade
,
emContinuous
,
// trackWriter,
observationLevel
,
longprof
);
/* === END: SETUP PROCESS LIST === */
...
...
@@ -390,8 +382,7 @@ int main(int argc, char **argv) {
CORSIKA_LOG_INFO
(
"Primary Energy: {}"
,
E0
);
CORSIKA_LOG_INFO
(
"Primary Momentum: {}"
,
P0
);
CORSIKA_LOG_INFO
(
"Point of Injection: {}"
,
injectionPos
.
getCoordinates
());
CORSIKA_LOG_INFO
(
"Shower Axis Length: {}"
,
(
showerCore
-
injectionPos
).
getNorm
()
*
1.2
);
CORSIKA_LOG_INFO
(
"Shower Axis Length: {}"
,
(
showerCore
-
injectionPos
).
getNorm
()
*
1.2
);
// trigger the output manager to open the library for writing
output
.
startOfLibrary
();
...
...
@@ -403,10 +394,8 @@ int main(int argc, char **argv) {
// directory for outputs
string
const
outdir
(
app
[
"--filename"
]
->
as
<
std
::
string
>
());
string
const
labHist_file
=
outdir
+
"/inthist_lab_"
+
to_string
(
i_shower
)
+
".npz"
;
string
const
cMSHist_file
=
outdir
+
"/inthist_cms_"
+
to_string
(
i_shower
)
+
".npz"
;
string
const
labHist_file
=
outdir
+
"/inthist_lab_"
+
to_string
(
i_shower
)
+
".npz"
;
string
const
cMSHist_file
=
outdir
+
"/inthist_cms_"
+
to_string
(
i_shower
)
+
".npz"
;
// setup particle stack, and add primary particle
stack
.
clear
();
...
...
@@ -428,16 +417,15 @@ int main(int argc, char **argv) {
HEPEnergyType
const
Efinal
=
dEdX
.
getEnergyLost
()
+
observationLevel
.
getEnergyGround
();
CORSIKA_LOG_INFO
(
"total energy budget (GeV): {} (dEdX={} ground={}), "
"relative difference (%): {}"
,
Efinal
/
1
_GeV
,
dEdX
.
getEnergyLost
()
/
1
_GeV
,
observationLevel
.
getEnergy
Ground
()
/
1
_GeV
,
(
Efinal
/
E0
-
1
)
*
100
);
CORSIKA_LOG_INFO
(
"total energy budget (GeV): {} (dEdX={} ground={}), "
"relative difference (%): {}"
,
Efinal
/
1
_GeV
,
dEdX
.
getEnergy
Lost
()
/
1
_GeV
,
observationLevel
.
getEnergyGround
()
/
1
_GeV
,
(
Efinal
/
E0
-
1
)
*
100
);
// auto const hists = heModelCounted.getHistogram() +
// urqmdCounted.getHistogram();
auto
const
hists
=
sibyllCounted
.
getHistogram
()
+
urqmdCounted
.
getHistogram
();
auto
const
hists
=
sibyllCounted
.
getHistogram
()
+
urqmdCounted
.
getHistogram
();
save_hist
(
hists
.
labHist
(),
labHist_file
,
true
);
save_hist
(
hists
.
CMSHist
(),
cMSHist_file
,
true
);
...
...
examples/em_shower.cpp
View file @
322427c4
...
...
@@ -80,7 +80,7 @@ void registerRandomStreams(int seed) {
template
<
typename
T
>
using
MyExtraEnv
=
MediumPropertyModel
<
UniformMagneticField
<
T
>>
;
int
main
(
int
argc
,
char
**
argv
)
{
int
main
(
int
argc
,
char
**
argv
)
{
logging
::
set_level
(
logging
::
level
::
warn
);
...
...
@@ -92,18 +92,15 @@ int main(int argc, char **argv) {
feenableexcept
(
FE_INVALID
);
int
seed
=
0
;
if
(
argc
>=
3
)
{
seed
=
std
::
stoi
(
std
::
string
(
argv
[
2
]));
}
if
(
argc
>=
3
)
{
seed
=
std
::
stoi
(
std
::
string
(
argv
[
2
]));
}
// initialize random number sequence(s)
registerRandomStreams
(
seed
);
// setup environment, geometry
using
EnvironmentInterface
=
IMediumPropertyModel
<
IMagneticFieldModel
<
IMediumModel
>>
;
using
EnvironmentInterface
=
IMediumPropertyModel
<
IMagneticFieldModel
<
IMediumModel
>>
;
using
EnvType
=
Environment
<
EnvironmentInterface
>
;
EnvType
env
;
CoordinateSystemPtr
const
&
rootCS
=
env
.
getCoordinateSystem
();
CoordinateSystemPtr
const
&
rootCS
=
env
.
getCoordinateSystem
();
Point
const
center
{
rootCS
,
0
_m
,
0
_m
,
0
_m
};
// build a Linsley US Standard atmosphere into `env`
...
...
@@ -148,11 +145,9 @@ int main(int argc, char **argv) {
static_pow
<
2
>
(
injectionHeight
));
Point
const
showerCore
{
rootCS
,
0
_m
,
0
_m
,
observationHeight
};
Point
const
injectionPos
=
showerCore
+
DirectionVector
{
rootCS
,
{
-
sin
(
thetaRad
),
0
,
cos
(
thetaRad
)}}
*
t
;
showerCore
+
DirectionVector
{
rootCS
,
{
-
sin
(
thetaRad
),
0
,
cos
(
thetaRad
)}}
*
t
;
std
::
cout
<<
"point of injection: "
<<
injectionPos
.
getCoordinates
()
<<
std
::
endl
;
std
::
cout
<<
"point of injection: "
<<
injectionPos
.
getCoordinates
()
<<
std
::
endl
;
stack
.
addParticle
(
std
::
make_tuple
(
beamCode
,
calculate_kinetic_energy
(
plab
.
getNorm
(),
get_mass
(
beamCode
)),
...
...
@@ -161,8 +156,8 @@ int main(int argc, char **argv) {
CORSIKA_LOG_INFO
(
"shower axis length: {} "
,
(
showerCore
-
injectionPos
).
getNorm
()
*
1.02
);
ShowerAxis
const
showerAxis
{
injectionPos
,
(
showerCore
-
injectionPos
)
*
1.02
,
env
,
false
,
1000
};
ShowerAxis
const
showerAxis
{
injectionPos
,
(
showerCore
-
injectionPos
)
*
1.02
,
env
,
false
,
1000
};
OutputManager
output
(
"em_shower_outputs"
);
...
...
@@ -172,15 +167,13 @@ int main(int argc, char **argv) {
// setup processes, decays and interactions
ParticleCut
<
SubWriter
<
decltype
(
dEdX
)
>>
cut
(
2
_MeV
,
2
_MeV
,
100
_GeV
,
100
_GeV
,
true
,
dEdX
);
ParticleCut
<
SubWriter
<
decltype
(
dEdX
)
>>
cut
(
2
_MeV
,
2
_MeV
,
100
_GeV
,
100
_GeV
,
true
,
dEdX
);
corsika
::
sibyll
::
Interaction
sibyll
{
env
};
corsika
::
sophia
::
InteractionModel
sophia
;
HEPEnergyType
heThresholdNN
=
60
_GeV
;
corsika
::
proposal
::
Interaction
emCascade
(
env
,
sophia
,
sibyll
.
getHadronInteractionModel
(),
heThresholdNN
);
corsika
::
proposal
::
ContinuousProcess
<
SubWriter
<
decltype
(
dEdX
)
>>
emContinuous
(
env
,
dEdX
);
corsika
::
proposal
::
ContinuousProcess
<
SubWriter
<
decltype
(
dEdX
)
>>
emContinuous
(
env
,
dEdX
);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// NOT possible right now, due to interface differenc in PROPOSAL
...
...
@@ -199,8 +192,7 @@ int main(int argc, char **argv) {
obsPlane
,
DirectionVector
(
rootCS
,
{
1.
,
0.
,
0.
})};
output
.
add
(
"particles"
,
observationLevel
);
auto
sequence
=
make_sequence
(
emCascade
,
emContinuous
,
longprof
,
cut
,
observationLevel
);
auto
sequence
=
make_sequence
(
emCascade
,
emContinuous
,
longprof
,
cut
,
observationLevel
);
// define air shower object, run simulation
setup
::
Tracking
tracking
;
...
...
@@ -212,12 +204,12 @@ int main(int argc, char **argv) {
EAS
.
run
();
HEPEnergyType
const
Efinal
=
dEdX
.
getEnergyLost
()
+
observationLevel
.
getEnergyGround
();
HEPEnergyType
const
Efinal
=
dEdX
.
getEnergyLost
()
+
observationLevel
.
getEnergyGround
();
CORSIKA_LOG_INFO
(
"total energy budget (GeV): {}, "
"relative difference (%): {}"
,
Efinal
/
1
_GeV
,
(
Efinal
/
E0
-
1
)
*
100
);
CORSIKA_LOG_INFO
(
"total energy budget (GeV): {}, "
"relative difference (%): {}"
,
Efinal
/
1
_GeV
,
(
Efinal
/
E0
-
1
)
*
100
);
output
.
endOfLibrary
();
}
examples/mars.cpp
View file @
322427c4
...
...
@@ -12,56 +12,57 @@
// to include it first...
#include
<corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include
<corsika/framework/geometry/PhysicalGeometry.hpp>
#include
<corsika/framework/geometry/Plane.hpp>
#include
<corsika/framework/geometry/Sphere.hpp>
#include
<corsika/framework/geometry/PhysicalGeometry.hpp>
#include
<corsika/framework/core/
Logging
.hpp>
#include
<corsika/framework/core/
Cascade
.hpp>
#include
<corsika/framework/core/EnergyMomentumOperations.hpp>
#include
<corsika/framework/core/Logging.hpp>
#include
<corsika/framework/core/PhysicalUnits.hpp>
#include
<corsika/framework/core/Cascade.hpp>
#include
<corsika/framework/utility/SaveBoostHistogram.hpp>
#include
<corsika/framework/utility/CorsikaFenv.hpp>
#include
<corsika/framework/process/InteractionCounter.hpp>
#include
<corsika/framework/process/ProcessSequence.hpp>
#include
<corsika/framework/process/SwitchProcessSequence.hpp>
#include
<corsika/framework/process/InteractionCounter.hpp>
#include
<corsika/framework/random/RNGManager.hpp>
#include
<corsika/framework/utility/CorsikaFenv.hpp>
#include
<corsika/framework/utility/SaveBoostHistogram.hpp>
#include
<corsika/output/OutputManager.hpp>
#include
<corsika/modules/writers/SubWriter.hpp>
#include
<corsika/modules/writers/EnergyLossWriter.hpp>
#include
<corsika/modules/writers/LongitudinalWriter.hpp>
#include
<corsika/modules/writers/SubWriter.hpp>
#include
<corsika/output/OutputManager.hpp>
#include
<corsika/media/Environment.hpp>
#include
<corsika/media/FlatExponential.hpp>
#include
<corsika/media/HomogeneousMedium.hpp>
#include
<corsika/media/IMagneticFieldModel.hpp>
#include
<corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include
<corsika/media/NuclearComposition.hpp>
#include
<corsika/media/MediumPropertyModel.hpp>
#include
<corsika/media/
UniformMagneticField
.hpp>
#include
<corsika/media/
NuclearComposition
.hpp>
#include
<corsika/media/ShowerAxis.hpp>
#include
<corsika/media/SlidingPlanarExponential.hpp>
#include
<corsika/media/UniformMagneticField.hpp>
#include
<corsika/modules/BetheBlochPDG.hpp>
#include
<corsika/modules/LongitudinalProfile.hpp>
#include
<corsika/modules/ObservationPlane.hpp>
#include
<corsika/modules/StackInspector.hpp>
#include
<corsika/modules/TrackWriter.hpp>
#include
<corsika/modules/PROPOSAL.hpp>
#include
<corsika/modules/ParticleCut.hpp>
#include
<corsika/modules/Pythia8.hpp>
#include
<corsika/modules/QGSJetII.hpp>
#include
<corsika/modules/Sibyll.hpp>
#include
<corsika/modules/Sophia.hpp>
#include
<corsika/modules/StackInspector.hpp>
#include
<corsika/modules/TrackWriter.hpp>
#include
<corsika/modules/UrQMD.hpp>
#include
<corsika/modules/PROPOSAL.hpp>
#include
<corsika/modules/QGSJetII.hpp>
#include
<corsika/setup/SetupStack.hpp>
#include
<corsika/setup/SetupTrajectory.hpp>
#include
<CLI/App.hpp>
#include
<CLI/Formatter.hpp>
#include
<CLI/Config.hpp>
#include
<CLI/Formatter.hpp>
#include
<iomanip>
#include
<iostream>
...
...
@@ -118,6 +119,7 @@ void registerRandomStreams(int seed) {
RNGManager
<>::
getInstance
().
registerRandomStream
(
"cascade"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"qgsjet"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"sibyll"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"sophia"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"pythia"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"urqmd"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"proposal"
);
...
...
@@ -350,12 +352,12 @@ int main(int argc, char** argv) {
// decaySibyll.printDecayConfig();
// energy threshold for high energy hadronic model. Affects LE/HE switch for
hadron
// interactions and the hadronic photon model in proposal
// energy threshold for high energy hadronic model. Affects LE/HE switch for
//
hadron
interactions and the hadronic photon model in proposal
HEPEnergyType
heHadronModelThreshold
=
63.1
_GeV
;
corsika
::
proposal
::
Interaction
emCascade
(
env
,
sibyll
.
getHadronInteractionModel
(),
heHadronModelThreshold
);
corsika
::
sophia
::
InteractionModel
sophia
;
corsika
::
proposal
::
Interaction
emCascade
(
env
,
sophia
,
sibyll
.
getHadronInteractionModel
(),
heHadronModelThreshold
);
// NOT possible right now, due to interface difference for PROPOSAL:
// InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
...
...
@@ -398,7 +400,8 @@ int main(int argc, char** argv) {
cut
,
trackWriter
,
observationLevel
,
profile
);
/* === END: SETUP PROCESS LIST === */
// create the cascade object using the default stack and tracking implementation
// create the cascade object using the default stack and tracking
// implementation
setup
::
Tracking
tracking
;
setup
::
Stack
<
EnvType
>
stack
;
Cascade
EAS
(
env
,
tracking
,
sequence
,
output
,
stack
);
...
...
examples/vertical_EAS.cpp
View file @
322427c4
...
...
@@ -14,47 +14,48 @@
// to include it first...
#include
<corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include
<corsika/framework/process/ProcessSequence.hpp>
#include
<corsika/framework/process/SwitchProcessSequence.hpp>
#include
<corsika/framework/process/InteractionCounter.hpp>
#include
<corsika/framework/geometry/Plane.hpp>
#include
<corsika/framework/geometry/Sphere.hpp>
#include
<corsika/framework/geometry/PhysicalGeometry.hpp>
#include
<corsika/framework/core/Logging.hpp>
#include
<corsika/framework/core/Cascade.hpp>
#include
<corsika/framework/core/EnergyMomentumOperations.hpp>
#include
<corsika/framework/core/Logging.hpp>
#include
<corsika/framework/core/PhysicalUnits.hpp>
#include
<corsika/framework/core/Cascade.hpp>
#include
<corsika/framework/utility/SaveBoostHistogram.hpp>
#include
<corsika/framework/utility/CorsikaFenv.hpp>
#include
<corsika/framework/geometry/PhysicalGeometry.hpp>
#include
<corsika/framework/geometry/Plane.hpp>
#include
<corsika/framework/geometry/Sphere.hpp>
#include
<corsika/framework/process/InteractionCounter.hpp>
#include
<corsika/framework/process/ProcessSequence.hpp>
#include
<corsika/framework/process/SwitchProcessSequence.hpp>
#include
<corsika/framework/random/RNGManager.hpp>
#include
<corsika/framework/utility/CorsikaFenv.hpp>
#include
<corsika/framework/utility/SaveBoostHistogram.hpp>
#include
<corsika/output/OutputManager.hpp>
#include
<corsika/modules/writers/SubWriter.hpp>
#include
<corsika/modules/writers/EnergyLossWriter.hpp>
#include
<corsika/modules/writers/LongitudinalWriter.hpp>
#include
<corsika/modules/writers/SubWriter.hpp>
#include
<corsika/output/OutputManager.hpp>
#include
<corsika/media/CORSIKA7Atmospheres.hpp>
#include
<corsika/media/Environment.hpp>
#include
<corsika/media/FlatExponential.hpp>
#include
<corsika/media/GeomagneticModel.hpp>
#include
<corsika/media/HomogeneousMedium.hpp>
#include
<corsika/media/IMagneticFieldModel.hpp>
#include
<corsika/media/NuclearComposition.hpp>
#include
<corsika/media/MediumPropertyModel.hpp>
#include
<corsika/media/
UniformMagneticField
.hpp>
#include
<corsika/media/
NuclearComposition
.hpp>
#include
<corsika/media/ShowerAxis.hpp>
#include
<corsika/media/
CORSIKA7Atmospheres
.hpp>
#include
<corsika/media/
UniformMagneticField
.hpp>
#include
<corsika/modules/BetheBlochPDG.hpp>
#include
<corsika/modules/Epos.hpp>
#include
<corsika/modules/LongitudinalProfile.hpp>
#include
<corsika/modules/ObservationPlane.hpp>
#include
<corsika/modules/StackInspector.hpp>
#include
<corsika/modules/TrackWriter.hpp>
#include
<corsika/modules/PROPOSAL.hpp>
#include
<corsika/modules/ParticleCut.hpp>
#include
<corsika/modules/Pythia8.hpp>
#include
<corsika/modules/Sibyll.hpp>
#include
<corsika/modules/Sophia.hpp>
#include
<corsika/modules/StackInspector.hpp>
#include
<corsika/modules/TrackWriter.hpp>
#include
<corsika/modules/UrQMD.hpp>
#include
<corsika/modules/Epos.hpp>
#include
<corsika/modules/PROPOSAL.hpp>
#include
<corsika/setup/SetupStack.hpp>
#include
<corsika/setup/SetupTrajectory.hpp>
...
...
@@ -85,6 +86,7 @@ using Particle = setup::Stack<EnvType>::particle_type;
void
registerRandomStreams
(
int
seed
)
{
RNGManager
<>::
getInstance
().
registerRandomStream
(
"cascade"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"sibyll"
);
// RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager
<>::
getInstance
().
registerRandomStream
(
"pythia"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"urqmd"
);
RNGManager
<>::
getInstance
().
registerRandomStream
(
"proposal"
);
...
...
@@ -217,8 +219,8 @@ int main(int argc, char** argv) {
// construct the continuous energy loss model
BetheBlochPDG
<
SubWriter
<
decltype
(
dEdX
)
>>
emContinuous
{
dEdX
};
// construct a particle cut - cuts are set to values close to reality, put
higher
// values for faster runs
// construct a particle cut - cuts are set to values close to reality, put
//
higher
values for faster runs
ParticleCut
<
SubWriter
<
decltype
(
dEdX
)
>>
cut
{
2
_MeV
,
2
_MeV
,
2
_GeV
,
300
_MeV
,
true
,
dEdX
};
// setup longitudinal profile
...
...
@@ -238,8 +240,10 @@ int main(int argc, char** argv) {
HEPEnergyType
heThresholdNN
=
60
_GeV
;
// PROPOSAL is disabled for this example
// corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
// heThresholdNN); corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
// corsika::sophia::InteractionModel sophia;
// corsika::proposal::Interaction emCascade(env, sophia,
// sibyll.getHadronInteractionModel(), heThresholdNN);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
// emContinuous(env, dEdX);
corsika
::
pythia8
::
Decay
decayPythia
;
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment