IAP GITLAB

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

some small adjustments, made for the talk in Mainz

parent bd8c3a91
No related branches found
No related tags found
No related merge requests found
...@@ -27,8 +27,20 @@ namespace corsika::cascade { ...@@ -27,8 +27,20 @@ namespace corsika::cascade {
void Run() { void Run() {
while (!fStack.IsEmpty()) { while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) { while (!fStack.IsEmpty()) {
Particle& p = *fStack.GetNextParticle(); //Particle& p = *fStack.GetNextParticle();
Step(p); EnergyType Emin;
typename Stack::StackIterator pMin(fStack, 0);
bool first = true;
for (typename Stack::StackIterator ip = fStack.begin(); ip!=fStack.end(); ++ip)
{
if (first || ip.GetEnergy()<Emin) {
first = false;
pMin = ip;
Emin = pMin.GetEnergy();
}
}
Step(pMin);
} }
// do cascade equations, which can put new particles on Stack, // do cascade equations, which can put new particles on Stack,
// thus, the double loop // thus, the double loop
......
...@@ -13,6 +13,8 @@ using namespace corsika::units; ...@@ -13,6 +13,8 @@ using namespace corsika::units;
#include <iostream> #include <iostream>
using namespace std; using namespace std;
static int fCount = 0;
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> { class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public: public:
ProcessSplit() {} ProcessSplit() {}
...@@ -28,13 +30,13 @@ public: ...@@ -28,13 +30,13 @@ public:
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const { void DoDiscrete(Particle& p, Stack& s) const {
EnergyType E = p.GetEnergy(); EnergyType E = p.GetEnergy();
if (E < 1_GeV) { if (E < 85_MeV) {
p.Delete(); p.Delete();
fCount++; fCount++;
} else { } else {
p.SetEnergy(E / 2); p.SetEnergy(E / 2);
s.NewParticle().SetEnergy(E / 2); s.NewParticle().SetEnergy(E / 2);
} }
} }
void Init() { fCount = 0; } void Init() { fCount = 0; }
...@@ -42,7 +44,6 @@ public: ...@@ -42,7 +44,6 @@ public:
int GetCount() { return fCount; } int GetCount() { return fCount; }
private: private:
mutable int fCount = 0;
}; };
class ProcessReport : public corsika::process::BaseProcess<ProcessReport> { class ProcessReport : public corsika::process::BaseProcess<ProcessReport> {
...@@ -59,17 +60,22 @@ public: ...@@ -59,17 +60,22 @@ public:
template <typename Particle, typename Trajectory, typename Stack> template <typename Particle, typename Trajectory, typename Stack>
void DoContinuous(Particle& p, Trajectory& t, Stack& s) const { void DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
static int countStep = 0;
if (!fReport) return; if (!fReport) return;
static int fCount = 0; //std::cout << "generation " << countStep << std::endl;
std::cout << "generation " << fCount << std::endl;
int i = 0; int i = 0;
EnergyType Etot = 0_GeV;
for (auto& iterP : s) { for (auto& iterP : s) {
EnergyType E = iterP.GetEnergy(); EnergyType E = iterP.GetEnergy();
std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID() Etot += E;
/* std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()
<< ", E=" << double(E / 1_GeV) << " GeV " << ", E=" << double(E / 1_GeV) << " GeV "
<< " | " << std::endl; << " | " << std::endl;
*/
} }
fCount++; countStep++;
//cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl;
cout << countStep << " " << s.GetSize() << " " << Etot/1_GeV << " " << fCount << endl;
} }
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
...@@ -79,25 +85,32 @@ public: ...@@ -79,25 +85,32 @@ public:
TEST_CASE("Cascade", "[Cascade]") { TEST_CASE("Cascade", "[Cascade]") {
ProcessReport p0(false); ProcessReport p0(true);
ProcessSplit p1; ProcessSplit p1;
const auto sequence = p0 + p1; const auto sequence = p0 + p1;
corsika::stack::super_stupid::SuperStupidStack stack; corsika::stack::super_stupid::SuperStupidStack stack;
corsika::cascade::Cascade<corsika::geometry::LineTrajectory, decltype(sequence), corsika::cascade::Cascade<corsika::geometry::LineTrajectory, decltype(sequence),
decltype(stack)> decltype(stack)>
EAS(sequence, stack); EAS(sequence, stack);
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
particle.SetEnergy(E0);
EAS.Init();
EAS.Run();
SECTION("sectionTwo") { SECTION("sectionTwo") {
for (int i = 0; i < 5; ++i) { for (int i = 0; i < 0; ++i) {
stack.Clear(); stack.Clear();
auto particle = stack.NewParticle(); auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV * pow(10, i); EnergyType E0 = 100_GeV * pow(10, i);
particle.SetEnergy(E0); particle.SetEnergy(E0);
EAS.Init(); EAS.Init();
EAS.Run(); EAS.Run();
cout << "E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl; //cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
} }
} }
} }
...@@ -63,9 +63,12 @@ namespace corsika::stack { ...@@ -63,9 +63,12 @@ namespace corsika::stack {
StackIteratorInterface(const StackIteratorInterface& mit) StackIteratorInterface(const StackIteratorInterface& mit)
: fData(mit.fData) : fData(mit.fData)
, fIndex(mit.fIndex) {} , fIndex(mit.fIndex) {}
public:
StackIteratorInterface& operator=(const StackIteratorInterface& mit) { StackIteratorInterface& operator=(const StackIteratorInterface& mit) {
fData = mit.fData; fData = mit.fData;
fIndex = mit.fIndex; fIndex = mit.fIndex;
return *this;
} }
public: public:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment