IAP GITLAB

Commit 234500f5 authored by Matthieu Carrere's avatar Matthieu Carrere Committed by Ralf Ulrich

Replace vertical_EAS.cpp

parent 284271b7
......@@ -91,6 +91,10 @@ void registerRandomStreams(const int seed) {
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
// argv : 1.number of nucleons, 2.number of protons,
// 3.total energy in GeV, 4.number of showers,
// 5.seed (0 by default to generate random values for all)
int main(int argc, char** argv) {
corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
......@@ -108,7 +112,11 @@ int main(int argc, char** argv) {
feenableexcept(FE_INVALID);
int seed = 0;
if (argc > 4) seed = std::stoi(std::string(argv[4]));
int number_showers= std::stoi(std::string(argv[4]));
if (argc > 5){
seed = std::stoi(std::string(argv[5]));
}
// initialize random number sequence(s)
registerRandomStreams(seed);
......@@ -150,173 +158,189 @@ int main(int argc, char** argv) {
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
unsigned short const A = std::stoi(std::string(argv[1]));
Code beamCode;
HEPEnergyType mass;
unsigned short Z = 0;
if (A > 0) {
beamCode = Code::Nucleus;
Z = std::stoi(std::string(argv[2]));
mass = get_nucleus_mass(A, Z);
} else {
int pdg = std::stoi(std::string(argv[2]));
beamCode = convert_from_PDG(PDGCode(pdg));
mass = get_mass(beamCode);
}
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 0_km + builder.getEarthRadius();
auto const injectionHeight = 111.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
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;
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
} else {
if (A == 1) {
if (Z == 1) {
stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
} else if (Z == 0) {
stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
string output_dir="example_outputs";
for(int i_shower=1;i_shower<number_showers+1;i_shower++){
// directory for outputs
string labHist_dir=output_dir+"/inthist_lab_verticalEAS_"+to_string(i_shower)+".npz";
string cMSHist_dir=output_dir+"/inthist_cms_verticalEAS_"+to_string(i_shower)+".npz";
string longprof_dir=output_dir+"/longprof_verticalEAS_"+to_string(i_shower)+".txt";
string tracks_dir=output_dir+"/tracks_"+to_string(i_shower)+".dat";
string particles_dir=output_dir+"/particles_"+to_string(i_shower)+".dat";
std::cout << std::endl;
std::cout << "Shower "<<i_shower<<"/"<<number_showers<<std::endl;
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
unsigned short const A = std::stoi(std::string(argv[1]));
Code beamCode;
HEPEnergyType mass;
unsigned short Z = 0;
if (A > 0) {
beamCode = Code::Nucleus;
Z = std::stoi(std::string(argv[2]));
mass = get_nucleus_mass(A, Z);
} else {
int pdg = std::stoi(std::string(argv[2]));
beamCode = convert_from_PDG(PDGCode(pdg));
mass = get_mass(beamCode);
}
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 0_km + builder.getEarthRadius();
auto const injectionHeight = 111.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
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;
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
} else {
if (A == 1) {
if (Z == 1) {
stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
} else if (Z == 0) {
stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
} else {
std::cerr << "illegal parameters" << std::endl;
return EXIT_FAILURE;
}
} else {
std::cerr << "illegal parameters" << std::endl;
return EXIT_FAILURE;
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
}
} else {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
}
// we make the axis much longer than the inj-core distance since the
// profile will go beyond the core, depending on zenith angle
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
<< std::endl;
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
// setup processes, decays and interactions
// corsika::qgsjetII::Interaction qgsjet;
corsika::sibyll::Interaction sibyll;
InteractionCounter sibyllCounted(sibyll);
corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
InteractionCounter sibyllNucCounted(sibyllNuc);
corsika::pythia8::Decay decayPythia;
// use sibyll decay routine for decays of particles unknown to pythia
corsika::sibyll::Decay decaySibyll{{
Code::N1440Plus,
Code::N1440MinusBar,
Code::N1440_0,
Code::N1440_0Bar,
Code::N1710Plus,
Code::N1710MinusBar,
Code::N1710_0,
Code::N1710_0Bar,
Code::Pi1300Plus,
Code::Pi1300Minus,
Code::Pi1300_0,
Code::KStar0_1430_0,
Code::KStar0_1430_0Bar,
Code::KStar0_1430_Plus,
Code::KStar0_1430_MinusBar,
}};
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter(tracks_dir);
LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
particles_dir);
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
StackInspector<setup::Stack> stackInspect(50000, false, E0);
// assemble all processes into an ordered process list
struct EnergySwitch {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
SwitchResult operator()(const Particle& p) {
if (p.getEnergy() < cutE_)
return SwitchResult::First;
else
return SwitchResult::Second;
}
};
auto hadronSequence = make_select(
urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted), EnergySwitch(55_GeV));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
auto sequence =
make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
emContinuous, cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.forceInteraction();
EAS.run();
cut.showResults();
emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
cut.getEmEnergy() + emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
save_hist(hists.labHist(), labHist_dir, true);
save_hist(hists.CMSHist(), cMSHist_dir, true);
longprof.save(longprof_dir);
}
// we make the axis much longer than the inj-core distance since the
// profile will go beyond the core, depending on zenith angle
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
<< std::endl;
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
// setup processes, decays and interactions
// corsika::qgsjetII::Interaction qgsjet;
corsika::sibyll::Interaction sibyll;
InteractionCounter sibyllCounted(sibyll);
corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
InteractionCounter sibyllNucCounted(sibyllNuc);
corsika::pythia8::Decay decayPythia;
// use sibyll decay routine for decays of particles unknown to pythia
corsika::sibyll::Decay decaySibyll{{
Code::N1440Plus,
Code::N1440MinusBar,
Code::N1440_0,
Code::N1440_0Bar,
Code::N1710Plus,
Code::N1710MinusBar,
Code::N1710_0,
Code::N1710_0Bar,
Code::Pi1300Plus,
Code::Pi1300Minus,
Code::Pi1300_0,
Code::KStar0_1430_0,
Code::KStar0_1430_0Bar,
Code::KStar0_1430_Plus,
Code::KStar0_1430_MinusBar,
}};
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter("tracks.dat");
LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
"particles.dat");
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
StackInspector<setup::Stack> stackInspect(50000, false, E0);
// assemble all processes into an ordered process list
struct EnergySwitch {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
SwitchResult operator()(const Particle& p) {
if (p.getEnergy() < cutE_)
return SwitchResult::First;
else
return SwitchResult::Second;
}
};
auto hadronSequence = make_select(
urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted), EnergySwitch(55_GeV));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
auto sequence =
make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
emContinuous, cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.forceInteraction();
EAS.run();
cut.showResults();
emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
cut.getEmEnergy() + emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
longprof.save("longprof_verticalEAS.txt");
}
Markdown is supported
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