IAP GITLAB

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

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

parent 7bc8214f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/CMakeLists.txt
+ 30
− 5
View file @ 979a7a7c
......@@ -31,7 +31,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
ProcessPythia
CORSIKAcascade
ProcessEnergyLoss
ProcessStackInspector
# ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
CORSIKAprocesses
......@@ -47,8 +47,8 @@ CORSIKA_ADD_TEST (cascade_example)
add_executable (boundary_example boundary_example.cc)
target_compile_options(boundary_example PRIVATE -g) # do not skip asserts
target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
ProcessTrackingLine
......@@ -70,20 +70,45 @@ target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
ProcessSibyll
ProcessPythia
CORSIKAcascade
ProcessStackInspector
ProcessEnergyLoss
# ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS cascade_proton_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_proton_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)
target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
target_link_libraries (staticsequence_example
......
Documentation/Examples/cascade_example.cc
+ 62
− 44
View file @ 979a7a7c
......@@ -54,7 +54,7 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
class ProcessCut : public process::ContinuousProcess<ProcessCut> {
class ProcessCut : public process::SecondariesProcess<ProcessCut> {
HEPEnergyType fECut;
......@@ -148,47 +148,38 @@ public:
return is_inv;
}
template <typename Particle>
LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
const Code pid = p.GetPID();
if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
return 0_m;
} else {
LengthType next_step = 1_m * std::numeric_limits<double>::infinity();
cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
return next_step;
}
}
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) {
const Code pid = p.GetPID();
HEPEnergyType energy = p.GetEnergy();
cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
<< ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
EProcessReturn ret = EProcessReturn::eOk;
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += energy;
fEmCount += 1;
// p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) {
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;
template <typename TSecondaries>
EProcessReturn DoSecondaries(TSecondaries& vS) {
auto p = vS.begin();
while (p != vS.end()) {
const Code pid = p.GetPID();
HEPEnergyType energy = p.GetEnergy();
cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
<< ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
<< endl;
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += energy;
fEmCount += 1;
p.Delete();
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += energy;
fInvCount += 1;
p.Delete();
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += energy;
p.Delete();
} else if (p.GetTime() > 10_ms) {
cout << "removing OLD particle..." << endl;
fEnergy += energy;
p.Delete();
} else {
++p; // next entry in SecondaryView
}
}
return ret;
return EProcessReturn::eOk;
}
void Init() {
......@@ -216,6 +207,31 @@ public:
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
//
......@@ -255,7 +271,7 @@ int main() {
// setup processes, decays and interactions
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 = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
......@@ -267,12 +283,13 @@ int main() {
process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
process::sibyll::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
ObservationLevel obsLevel(1400_m);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
process::EnergyLoss::EnergyLoss eLoss;
// assemble all processes into an ordered process list
auto sequence = sibyll << sibyllNuc << decay << cut << trackWriter;
auto sequence = sibyll << sibyllNuc << decay << eLoss << cut << trackWriter; // << obsLevel
// setup particle stack, and add primary particle
setup::Stack stack;
......@@ -300,7 +317,8 @@ int main() {
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);
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,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
......
Documentation/Examples/staticsequence_example.cc
+ 3
− 4
View file @ 979a7a7c
......@@ -88,12 +88,11 @@ void modular() {
auto sequence = m1 << m2 << m3 << m4;
DummyData p;
DummyStack s;
DummyTrajectory t;
DummyData particle;
DummyTrajectory track;
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) {
// cout << p.p[i] << endl;
......
Documentation/Examples/vertical_EAS.cc
+ 28
− 30
View file @ 979a7a7c
......@@ -159,27 +159,28 @@ public:
const Code pid = p.GetPID();
HEPEnergyType energy = p.GetEnergy();
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)) {
cout << "removing em. particle..." << endl;
fEmEnergy += energy;
fEmCount += 1;
p.Delete();
cout << "removing em. particle..." << endl;
fEmEnergy += energy;
fEmCount += 1;
p.Delete();
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += energy;
fInvCount += 1;
p.Delete();
cout << "removing inv. particle..." << endl;
fInvEnergy += energy;
fInvCount += 1;
p.Delete();
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += energy;
p.Delete();
} else if (p.GetTime()>10_ms) {
cout << "removing OLD particle..." << endl;
fEnergy += energy;
p.Delete();
cout << "removing low en. particle..." << endl;
fEnergy += energy;
p.Delete();
} else if (p.GetTime() > 10_ms) {
cout << "removing OLD particle..." << endl;
fEnergy += energy;
p.Delete();
} else {
++p; // next entry in SecondaryView
++p; // next entry in SecondaryView
}
}
return EProcessReturn::eOk;
......@@ -210,33 +211,29 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
class ObservationLevel : public process::ContinuousProcess<ObservationLevel> {
LengthType fHeight;
public:
ObservationLevel(const LengthType vHeight)
: fHeight(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;
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() {}
};
//
......@@ -285,7 +282,7 @@ int main() {
// process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
ObservationLevel obsLevel(1400_m);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
......@@ -329,7 +326,8 @@ int main() {
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, 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,
corsika::stack::MomentumVector, geometry::Point,
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