IAP GITLAB

Skip to content
Snippets Groups Projects
Commit da7bee3a authored by ralfulrich's avatar ralfulrich
Browse files

use new SecondariesProcess, added ObservationLevel (passive version) for tests

parent ce37da29
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/CMakeLists.txt
+ 24
− 1
View file @ da7bee3a
...@@ -31,7 +31,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits ...@@ -31,7 +31,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
ProcessPythia ProcessPythia
CORSIKAcascade CORSIKAcascade
ProcessEnergyLoss ProcessEnergyLoss
ProcessStackInspector # ProcessStackInspector
ProcessTrackWriter ProcessTrackWriter
ProcessTrackingLine ProcessTrackingLine
ProcessHadronicElasticModel ProcessHadronicElasticModel
...@@ -45,6 +45,29 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits ...@@ -45,6 +45,29 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
install (TARGETS cascade_example DESTINATION share/examples) install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_example) CORSIKA_ADD_TEST (cascade_example)
add_executable (vertical_EAS vertical_EAS.cc)
target_compile_options(vertical_EAS PRIVATE -g) # do not skip asserts
target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia
CORSIKAcascade
ProcessEnergyLoss
# ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS vertical_EAS DESTINATION share/examples)
CORSIKA_ADD_TEST (vertical_EAS)
add_executable (staticsequence_example staticsequence_example.cc) add_executable (staticsequence_example staticsequence_example.cc)
target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
target_link_libraries (staticsequence_example target_link_libraries (staticsequence_example
......
Documentation/Examples/cascade_example.cc
+ 63
− 45
View file @ da7bee3a
...@@ -58,7 +58,7 @@ using namespace corsika::environment; ...@@ -58,7 +58,7 @@ using namespace corsika::environment;
using namespace std; using namespace std;
using namespace corsika::units::si; using namespace corsika::units::si;
class ProcessCut : public process::ContinuousProcess<ProcessCut> { class ProcessCut : public process::SecondariesProcess<ProcessCut> {
HEPEnergyType fECut; HEPEnergyType fECut;
...@@ -152,47 +152,38 @@ public: ...@@ -152,47 +152,38 @@ public:
return is_inv; return is_inv;
} }
template <typename Particle> template <typename TSecondaries>
LengthType MaxStepLength(Particle& p, setup::Trajectory&) const { EProcessReturn DoSecondaries(TSecondaries& vS) {
cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl; auto p = vS.begin();
cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl; while (p != vS.end()) {
const Code pid = p.GetPID(); const Code pid = p.GetPID();
if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) { HEPEnergyType energy = p.GetEnergy();
cout << "ProcessCut: MinStep: next cut: " << 0. << endl; cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
return 0_m; << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
} else { << endl;
LengthType next_step = 1_m * std::numeric_limits<double>::infinity(); if (isEmParticle(pid)) {
cout << "ProcessCut: MinStep: next cut: " << next_step << endl; cout << "removing em. particle..." << endl;
return next_step; fEmEnergy += energy;
} fEmCount += 1;
} p.Delete();
} else if (isInvisible(pid)) {
template <typename Particle, typename Stack> cout << "removing inv. particle..." << endl;
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) { fInvEnergy += energy;
const Code pid = p.GetPID(); fInvCount += 1;
HEPEnergyType energy = p.GetEnergy(); p.Delete();
cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy } else if (isBelowEnergyCut(p)) {
<< ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl; cout << "removing low en. particle..." << endl;
EProcessReturn ret = EProcessReturn::eOk; fEnergy += energy;
if (isEmParticle(pid)) { p.Delete();
cout << "removing em. particle..." << endl; } else if (p.GetTime() > 10_ms) {
fEmEnergy += energy; cout << "removing OLD particle..." << endl;
fEmCount += 1; fEnergy += energy;
// p.Delete(); p.Delete();
ret = EProcessReturn::eParticleAbsorbed; } else {
} else if (isInvisible(pid)) { ++p; // next entry in SecondaryView
cout << "removing inv. particle..." << endl; }
fInvEnergy += energy;
fInvCount += 1;
// p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += energy;
// p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} }
return ret; return EProcessReturn::eOk;
} }
void Init() { void Init() {
...@@ -220,6 +211,31 @@ public: ...@@ -220,6 +211,31 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; } HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
}; };
class ObservationLevel : public process::ContinuousProcess<ObservationLevel> {
LengthType fHeight;
public:
ObservationLevel(const LengthType vHeight)
: fHeight(vHeight) {}
template <typename Particle>
LengthType MaxStepLength(Particle&, setup::Trajectory&) const {
return 1_m * std::numeric_limits<double>::infinity();
}
template <typename TParticle, typename TTrack>
EProcessReturn DoContinuous(TParticle&, TTrack& vT) {
if ((vT.GetPosition(0).GetZ() <= fHeight && vT.GetPosition(1).GetZ() > fHeight) ||
(vT.GetPosition(0).GetZ() > fHeight && vT.GetPosition(1).GetZ() <= fHeight)) {
cout << "OBSERVED " << endl;
return EProcessReturn::eParticleAbsorbed;
}
return EProcessReturn::eOk;
}
void Init() {}
};
// //
// The example main program for a particle cascade // The example main program for a particle cascade
// //
...@@ -252,7 +268,7 @@ int main() { ...@@ -252,7 +268,7 @@ int main() {
// setup processes, decays and interactions // setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env); tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
stack_inspector::StackInspector<setup::Stack> stackInspect(true); // stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = { const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus, particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
...@@ -265,6 +281,7 @@ int main() { ...@@ -265,6 +281,7 @@ int main() {
// random::RNGManager::GetInstance().RegisterRandomStream("pythia"); // random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::pythia::Decay decay(trackedHadrons); // process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV); ProcessCut cut(20_GeV);
ObservationLevel obsLevel(1400_m);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction // process::HadronicElasticModel::HadronicElasticInteraction
...@@ -277,8 +294,8 @@ int main() { ...@@ -277,8 +294,8 @@ int main() {
// auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut << // auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut <<
// trackWriter; auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << // trackWriter; auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss <<
// cut << trackWriter; // cut << trackWriter;
// auto sequence = sibyll << sibyllNuc << decay << eLoss << cut; auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << obsLevel;
auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut; // auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n"; // << "\n";
...@@ -309,7 +326,8 @@ int main() { ...@@ -309,7 +326,8 @@ int main() {
cout << "input particle: " << beamCode << endl; 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,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{ units::si::TimeType, unsigned short, unsigned short>{
......
Documentation/Examples/staticsequence_example.cc
+ 3
− 4
View file @ da7bee3a
...@@ -88,12 +88,11 @@ void modular() { ...@@ -88,12 +88,11 @@ void modular() {
auto sequence = m1 << m2 << m3 << m4; auto sequence = m1 << m2 << m3 << m4;
DummyData p; DummyData particle;
DummyStack s; DummyTrajectory track;
DummyTrajectory t;
const int n = 1000; const int n = 1000;
for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); } for (int i = 0; i < n; ++i) { sequence.DoContinuous(particle, track); }
for (int i = 0; i < nData; ++i) { for (int i = 0; i < nData; ++i) {
// cout << p.p[i] << endl; // cout << p.p[i] << endl;
......
Documentation/Examples/vertical_EAS.cc
+ 28
− 30
View file @ da7bee3a
...@@ -159,27 +159,28 @@ public: ...@@ -159,27 +159,28 @@ public:
const Code pid = p.GetPID(); const Code pid = p.GetPID();
HEPEnergyType energy = p.GetEnergy(); HEPEnergyType energy = p.GetEnergy();
cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
<< ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl; << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
<< endl;
if (isEmParticle(pid)) { if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl; cout << "removing em. particle..." << endl;
fEmEnergy += energy; fEmEnergy += energy;
fEmCount += 1; fEmCount += 1;
p.Delete(); p.Delete();
} else if (isInvisible(pid)) { } else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl; cout << "removing inv. particle..." << endl;
fInvEnergy += energy; fInvEnergy += energy;
fInvCount += 1; fInvCount += 1;
p.Delete(); p.Delete();
} else if (isBelowEnergyCut(p)) { } else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl; cout << "removing low en. particle..." << endl;
fEnergy += energy; fEnergy += energy;
p.Delete(); p.Delete();
} else if (p.GetTime()>10_ms) { } else if (p.GetTime() > 10_ms) {
cout << "removing OLD particle..." << endl; cout << "removing OLD particle..." << endl;
fEnergy += energy; fEnergy += energy;
p.Delete(); p.Delete();
} else { } else {
++p; // next entry in SecondaryView ++p; // next entry in SecondaryView
} }
} }
return EProcessReturn::eOk; return EProcessReturn::eOk;
...@@ -210,33 +211,29 @@ public: ...@@ -210,33 +211,29 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; } HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
}; };
class ObservationLevel : public process::ContinuousProcess<ObservationLevel> { class ObservationLevel : public process::ContinuousProcess<ObservationLevel> {
LengthType fHeight; LengthType fHeight;
public: public:
ObservationLevel(const LengthType vHeight) ObservationLevel(const LengthType vHeight)
: fHeight(vHeight) {} : fHeight(vHeight) {}
template <typename Particle> template <typename Particle>
LengthType MaxStepLength(Particle&, setup::Trajectory&) const { LengthType MaxStepLength(Particle&, setup::Trajectory&) const {
return 1_m * std::numeric_limits<double>::infinity(); return 1_m * std::numeric_limits<double>::infinity();
} }
template <typename TParticle, typename TTrack> template <typename TParticle, typename TTrack>
EProcessReturn DoContinuous(TParticle&, TTrack& vT) { EProcessReturn DoContinuous(TParticle&, TTrack& vT) {
if ((vT.GetPosition(0).GetZ()<=fHeight && if ((vT.GetPosition(0).GetZ() <= fHeight && vT.GetPosition(1).GetZ() > fHeight) ||
vT.GetPosition(1).GetZ()>fHeight) || (vT.GetPosition(0).GetZ() > fHeight && vT.GetPosition(1).GetZ() <= fHeight)) {
(vT.GetPosition(0).GetZ()>fHeight && cout << "OBSERVED " << endl;
vT.GetPosition(1).GetZ()<=fHeight)) { return EProcessReturn::eParticleAbsorbed;
cout << "OBSERVED " << endl;
return EProcessReturn::eParticleAbsorbed;
} }
return EProcessReturn::eOk; return EProcessReturn::eOk;
} }
void Init() {} void Init() {}
}; };
// //
...@@ -285,7 +282,7 @@ int main() { ...@@ -285,7 +282,7 @@ int main() {
// process::pythia::Decay decay(trackedHadrons); // process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV); ProcessCut cut(20_GeV);
ObservationLevel obsLevel(1400_m); ObservationLevel obsLevel(1400_m);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction // process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env); // hadronicElastic(env);
...@@ -329,7 +326,8 @@ int main() { ...@@ -329,7 +326,8 @@ int main() {
cout << "input particle: " << beamCode << endl; 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, 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe Point pos(rootCS, 0_m, 0_m,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{ units::si::TimeType, unsigned short, unsigned short>{
......
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