IAP GITLAB

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

added new test section, final work

parent babba847
No related branches found
No related tags found
No related merge requests found
Showing
with 398 additions and 245 deletions
......@@ -36,7 +36,7 @@ using namespace std;
static int fCount = 0;
class ProcessSplit : public corsika::process::DiscreteProcess<ProcessSplit> {
class ProcessSplit : public corsika::process::InteractionProcess<ProcessSplit> {
public:
ProcessSplit() {}
......@@ -86,14 +86,13 @@ public:
*/
return int_length;
//
//int a = 0;
//const double next_step = -int_length * log(s_rndm_(a));
//std::cout << "ProcessSplit: "
// int a = 0;
// const double next_step = -int_length * log(s_rndm_(a));
// std::cout << "ProcessSplit: "
// << "next step (g/cm2): " << next_step << std::endl;
//return next_step;
// return next_step;
}
template <typename Particle, typename Track, typename Stack>
EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
// corsika::utls::ignore(p);
......@@ -101,8 +100,8 @@ public:
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
cout << "DoDiscrete: " << p.GetPID() << " interaction? "
void DoInteraction(Particle& p, Stack& s) const {
cout << "DoInteraction: " << p.GetPID() << " interaction? "
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
......@@ -127,11 +126,11 @@ public:
int kTarget = 1; // p.GetPID();
std::cout << "ProcessSplit: "
<< " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
<< " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
<< std::endl;
if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "ProcessSplit: "
<< " DoDiscrete: dropping particle.." << std::endl;
<< " DoInteraction: dropping particle.." << std::endl;
p.Delete();
fCount++;
} else {
......
......@@ -13,9 +13,9 @@
#define _include_corsika_cascade_Cascade_h_
#include <corsika/process/ProcessReturn.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <type_traits>
......@@ -55,49 +55,88 @@ namespace corsika::cascade {
}
}
void Step(Particle& particle) {
void Step(Particle& particle) {
// get access to random number generator
static corsika::random::RNG& rmng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
const double total_inv_lambda = fProcesseList.GetTotalInverseInteractionLength(particle, step);
// determine combined total interaction length (inverse)
const double total_inv_lambda =
fProcesseList.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length
// ....
static corsika::random::RNG& rmng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
const double sample_step = rmng() / (double)rmng.max();
const double next_step = -log(sample_step) / total_inv_lambda;
const double next_interact = -log(sample_step) / total_inv_lambda;
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
// determine the maximum geometric step length
const double distance_max = fProcesseList.MaxStepLength(particle, step);
std::cout << "distance_max=" << distance_max << std::endl;
// determine combined total inverse decay time
const double total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
const double min_step = fProcesseList.MinStepLength(particle, step);
// convert next_step from grammage to length
// sample random exponential decay time
const double sample_decay = rmng() / (double)rmng.max();
const double next_decay = -log(sample_decay) / total_inv_lifetime;
std::cout << "total_inv_lifetime=" << total_inv_lifetime
<< ", next_decay=" << next_decay << std::endl;
// convert next_step from grammage to length [m]
// Environment::GetDistance(step, next_step);
const double distance_interact = next_interact;
// ....
// update step with actual min(min_step, next_step)
// convert next_decay from time to length [m]
// Environment::GetDistance(step, next_decay);
const double distance_decay = next_decay;
// ....
// take minimum of geometry, interaction, decay for next step
const double distance_decay_interact = std::min(next_decay, next_interact);
const double distance_next = std::min(distance_decay_interact, distance_max);
/// here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.GetPosition(1));
// .... also update time, momentum, direction, ...
// apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
fProcesseList.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
// fStack.Delete(particle); // TODO: check if this is really needed
} else {
// check if min_step < random_step -> skip discrete if rndm>min_step/random_step
if (min_step<next_step) {
const double p_next = min_step/next_step;
const double sample_skip = rmng() / (double)rmng.max();
if (sample_skip>p_next) {
return;// corsika::process::EProcessReturn::eOk;
}
}
const double sample_process = rmng() / (double)rmng.max();
double inv_lambda_count = 0;
fProcesseList.SelectDiscrete(particle, fStack, total_inv_lambda,
sample_process, inv_lambda_count);
std::cout << "select process " << (distance_decay_interact < distance_max)
<< std::endl;
// check if geometry limits step, then quit this step
if (distance_decay_interact < distance_max) {
// check weather decay or interaction limits this step
if (distance_decay > distance_interact) {
std::cout << "collide" << std::endl;
const double actual_inv_length =
fProcesseList.GetTotalInverseInteractionLength(particle, step);
const double sample_process = rmng() / (double)rmng.max();
double inv_lambda_count = 0;
fProcesseList.SelectInteraction(particle, fStack, actual_inv_length,
sample_process, inv_lambda_count);
} else {
std::cout << "decay" << std::endl;
const double actual_decay_time =
fProcesseList.GetTotalInverseLifetime(particle);
const double sample_process = rmng() / (double)rmng.max();
double inv_decay_count = 0;
fProcesseList.SelectDecay(particle, fStack, actual_decay_time, sample_process,
inv_decay_count);
}
}
}
}
......
......@@ -41,23 +41,17 @@ using namespace std;
static int fCount = 0;
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
public:
ProcessSplit() {}
template <typename Particle, typename T>
double MinStepLength(Particle&, T&) const { return 1; }
template <typename Particle, typename T>
double GetInteractionLength(Particle&, T&) const { return 1; }
template <typename Particle, typename T, typename Stack>
EProcessReturn DoContinuous(Particle&, T&, Stack&) const {
return EProcessReturn::eOk;
double MaxStepLength(Particle&, T&) const {
return 1;
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
template <typename Particle, typename T, typename Stack>
void DoContinuous(Particle& p, T&, Stack& s) const {
EnergyType E = p.GetEnergy();
if (E < 85_MeV) {
p.Delete();
......@@ -89,15 +83,15 @@ TEST_CASE("Cascade", "[Cascade]") {
rmng.RegisterRandomStream(str_name);
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
......@@ -109,7 +103,7 @@ TEST_CASE("Cascade", "[Cascade]") {
particle.SetTime(0_ns);
EAS.Init();
EAS.Run();
/*
SECTION("sectionTwo") {
for (int i = 0; i < 0; ++i) {
......@@ -119,7 +113,7 @@ TEST_CASE("Cascade", "[Cascade]") {
particle.SetEnergy(E0);
EAS.Init();
EAS.Run();
// cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
}
}
......
......@@ -29,27 +29,29 @@ namespace corsika::process {
struct BaseProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/*
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
template <typename Particle, typename Track, typename Stack>
inline EProcessReturn DoContinuous(Particle&, Track&, Stack&) const; // {}
*/
/*
template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const {
return 1./GetRef().GetInteractionLength(p, t);
}
};
*/
};
/*
template<template<typename, typename> class T, typename A, typename B>
typename BaseProcess< T<A, B> >::is_process_sequence
typename BaseProcess< T<A, B> >::is_process_sequence
{
static const bool value = true;
};
*/
/*
template <typename T>
struct is_base {
......
......@@ -12,7 +12,8 @@ set (
CORSIKAprocesssequence_HEADERS
BaseProcess.h
ContinuousProcess.h
DiscreteProcess.h
InteractionProcess.h
DecayProcess.h
ProcessSequence.h
ProcessReturn.h
)
......
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_decayprocess_h_
#define _include_corsika_decayprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
namespace corsika::process {
/**
\class DecayProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type DecayProcess<T>
*/
template <typename derived>
struct DecayProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/// here starts the interface-definition part
// -> enforce derived to implement DoDecay...
template <typename Particle, typename Stack>
inline EProcessReturn DoDecay(Particle&, Stack&) const;
template <typename Particle>
inline double GetLifetime(Particle& p) const;
template <typename Particle>
inline double GetInverseLifetime(Particle& p) const {
return 1. / GetRef().GetLifetime(p);
}
};
} // namespace corsika::process
#endif
......@@ -36,13 +36,15 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce derived to implement DoDiscrete...
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
inline EProcessReturn DoDiscrete(Particle&, Stack&) const;
template <typename Particle, typename Track>
inline double GetInteractionLength(Particle& p, Track& t) const;
template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const {
return 1./GetRef().GetInteractionLength(p, t);
return 1. / GetRef().GetInteractionLength(p, t);
}
};
} // namespace corsika::process
......
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_interactionprocess_h_
#define _include_corsika_interactionprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
namespace corsika::process {
/**
\class InteractionProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type InteractionProcess<T>
*/
template <typename derived>
struct InteractionProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/// here starts the interface-definition part
// -> enforce derived to implement DoInteraction...
template <typename Particle, typename Stack>
inline EProcessReturn DoInteraction(Particle&, Stack&) const;
template <typename Particle, typename Track>
inline double GetInteractionLength(Particle& p, Track& t) const;
template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const {
return 1. / GetRef().GetInteractionLength(p, t);
}
};
} // namespace corsika::process
#endif
......@@ -24,7 +24,8 @@ namespace corsika::process {
eOk = 1,
eParticleAbsorbed = 2,
eInteracted = 3,
};
eDecayed = 4,
};
} // namespace corsika::process
#endif
......@@ -14,11 +14,13 @@
#include <corsika/process/BaseProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/DecayProcess.h>
//#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <limits>
#include <cmath>
#include <limits>
//#include <corsika/setup/SetupTrajectory.h>
// using corsika::setup::Trajectory;
......@@ -137,24 +139,15 @@ namespace corsika::process {
https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
//template <typename T1, typename T2>
//class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> >;
// this is a marker to track which BaseProcess is also a ProcessSequence
template <typename T>
struct is_process_sequence {
static const bool value = false;
};
template <typename T>
struct is_process_sequence {
static const bool value = false;
};
template <typename T1, typename T2>
class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
/* template <>
struct BaseProcess<ProcessSequence<T1, T2> >::is_process_sequence<true> {
static const bool value = true;
};*/
public:
const T1& A;
const T2& B;
......@@ -169,27 +162,29 @@ namespace corsika::process {
template <typename Particle, typename Track, typename Stack>
inline EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
// A.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
A.DoContinuous(p, t, s);
}
if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
// B.DoContinuous(std::forward<Particle>(p), t, std::forward<Stack>(s));
if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
is_process_sequence<T2>::value) {
B.DoContinuous(p, t, s);
}
return ret;
}
template <typename Particle, typename Track>
inline double MinStepLength(Particle& p, Track& track) const {
double min_length = std::numeric_limits<double>::infinity();
if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
min_length = std::min(min_length, A.MinStepLength(p,track));
}
if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
min_length = std::min(min_length, B.MinStepLength(p,track));
}
return min_length;
inline double MaxStepLength(Particle& p, Track& track) const {
double max_length = std::numeric_limits<double>::infinity();
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
max_length = std::min(max_length, A.MaxStepLength(p, track));
}
if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
is_process_sequence<T2>::value) {
max_length = std::min(max_length, B.MaxStepLength(p, track));
}
return max_length;
}
template <typename Particle, typename Track>
......@@ -205,64 +200,116 @@ namespace corsika::process {
template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const {
double tot = 0;
if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
tot += A.GetInverseInteractionLength(p, t);
}
if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value ||
is_process_sequence<T2>::value) {
tot += B.GetInverseInteractionLength(p, t);
}
return tot;
}
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle& p, Stack& s) const {
if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
A.DoDiscrete(p, s);
inline EProcessReturn SelectInteraction(Particle& p, Stack& s,
const double lambda_inv_tot,
const double rndm_select,
double& lambda_inv_count) const {
if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectInteraction(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += A.GetInverseInteractionLength(p, s);
// check if we should execute THIS process and then EXIT
if (rndm_select < lambda_inv_count / lambda_inv_tot) {
A.DoInteraction(p, s);
return EProcessReturn::eInteracted;
}
} // end branch A
if constexpr (is_process_sequence<T2>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
B.SelectInteraction(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += B.GetInverseInteractionLength(p, s);
// check if we should execute THIS process and then EXIT
if (rndm_select < lambda_inv_count / lambda_inv_tot) {
B.DoInteraction(p, s);
return EProcessReturn::eInteracted;
}
} // end branch A
return EProcessReturn::eOk;
}
template <typename Particle>
inline double GetTotalLifetime(Particle& p) const {
return 1. / GetInverseLifetime(p);
}
template <typename Particle>
inline double GetTotalInverseLifetime(Particle& p) const {
return GetInverseLifetime(p);
}
template <typename Particle>
inline double GetInverseLifetime(Particle& p) const {
double tot = 0;
if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
tot += A.GetInverseLifetime(p);
}
if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
B.DoDiscrete(p, s);
if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value ||
is_process_sequence<T2>::value) {
tot += B.GetInverseLifetime(p);
}
return EProcessReturn::eOk;
return tot;
}
// select decay process
template <typename Particle, typename Stack>
inline EProcessReturn SelectDiscrete(Particle& p, Stack& s,
const double lambda_inv_tot,
const double rndm_select,
double& lambda_inv_count) const {
inline EProcessReturn SelectDecay(Particle& p, Stack& s, const double decay_inv_tot,
const double rndm_select,
double& decay_inv_count) const {
if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret = A.SelectDiscrete(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
// if A did suceed, stop routine
if (ret != EProcessReturn::eOk) {
return ret;
}
} else if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += A.GetInverseInteractionLength(p, s);
// check if we should execute THIS process and then EXIT
if (rndm_select<lambda_inv_count/lambda_inv_tot) {
A.DoDiscrete(p, s);
return EProcessReturn::eInteracted;
}
} // end branch A
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectDecay(p, s, decay_inv_count, rndm_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += A.GetInverseLifetime(p);
// check if we should execute THIS process and then EXIT
if (rndm_select < decay_inv_count / decay_inv_tot) {
A.DoDecay(p, s);
return EProcessReturn::eDecayed;
}
} // end branch A
if constexpr (is_process_sequence<T2>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret = B.SelectDiscrete(p, s, lambda_inv_count, rndm_select, lambda_inv_count);
// if A did suceed, stop routine
if (ret != EProcessReturn::eOk) {
return ret;
}
} else if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += B.GetInverseInteractionLength(p, s);
// check if we should execute THIS process and then EXIT
if (rndm_select<lambda_inv_count/lambda_inv_tot) {
B.DoDiscrete(p, s);
return EProcessReturn::eInteracted;
}
} // end branch A
// if A is a process sequence --> check inside
const EProcessReturn ret =
B.SelectDecay(p, s, decay_inv_count, rndm_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += B.GetInverseLifetime(p);
// check if we should execute THIS process and then EXIT
if (rndm_select < decay_inv_count / decay_inv_tot) {
B.DoDecay(p, s);
return EProcessReturn::eDecayed;
}
} // end branch B
return EProcessReturn::eOk;
}
......@@ -273,8 +320,8 @@ namespace corsika::process {
}
};
/// the +operator assembles many BaseProcess, ContinuousProcess, and
/// DiscreteProcess objects into a ProcessSequence, all combinatorics
/// the + operator assembles many BaseProcess, ContinuousProcess, and
/// InteractionProcess objects into a ProcessSequence, all combinatorics
/// must be allowed, this is why we define a macro to define all
/// combinations here:
......@@ -285,93 +332,26 @@ namespace corsika::process {
}
OPSEQ(BaseProcess, BaseProcess)
OPSEQ(BaseProcess, DiscreteProcess)
OPSEQ(BaseProcess, InteractionProcess)
OPSEQ(BaseProcess, ContinuousProcess)
OPSEQ(BaseProcess, DecayProcess)
OPSEQ(ContinuousProcess, BaseProcess)
OPSEQ(ContinuousProcess, DiscreteProcess)
OPSEQ(ContinuousProcess, InteractionProcess)
OPSEQ(ContinuousProcess, ContinuousProcess)
OPSEQ(DiscreteProcess, BaseProcess)
OPSEQ(DiscreteProcess, DiscreteProcess)
OPSEQ(DiscreteProcess, ContinuousProcess)
/*
template <typename T1>
struct depth_lhs
{
static const int num = 0;
};
// terminating condition
template <typename T1, typename T2>
struct depth_lhs< Sequence<T1,T2> >
{
// try to expand the left node (T1) which might be a Sequence type
static const int num = 1 + depth_lhs<T1>::num;
};
*/
/*
template <typename T1>
struct mat_ptrs
{
static const int num = 0;
inline static void
get_ptrs(const Process** ptrs, const T1& X)
{
ptrs[0] = reinterpret_cast<const Process*>(&X);
}
};
template <typename T1, typename T2>
struct mat_ptrs< Sequence<T1,T2> >
{
static const int num = 1 + mat_ptrs<T1>::num;
inline static void
get_ptrs(const Process** in_ptrs, const Sequence<T1,T2>& X)
{
// traverse the left node
mat_ptrs<T1>::get_ptrs(in_ptrs, X.A);
// get address of the matrix on the right node
in_ptrs[num] = reinterpret_cast<const Process*>(&X.B);
}
};
*/
/*
template<typename T1, typename T2>
const Process&
Process::operator=(const Sequence<T1,T2>& X)
{
int N = 1 + depth_lhs< Sequence<T1,T2> >::num;
const Process* ptrs[N];
mat_ptrs< Sequence<T1,T2> >::get_ptrs(ptrs, X);
int r = ptrs[0]->rows;
int c = ptrs[0]->cols;
// ... check that all matrices have the same size ...
set_size(r, c);
for(int j=0; j<r*c; ++j)
{
double sum = ptrs[0]->data[j];
for(int i=1; i<N; ++i)
{
sum += ptrs[i]->data[j];
}
data[j] = sum;
}
return *this;
}
*/
template<template<typename, typename> class T, typename A, typename B>
struct is_process_sequence< T<A,B> >
{
static const bool value = true;
};
OPSEQ(ContinuousProcess, DecayProcess)
OPSEQ(InteractionProcess, BaseProcess)
OPSEQ(InteractionProcess, InteractionProcess)
OPSEQ(InteractionProcess, ContinuousProcess)
OPSEQ(InteractionProcess, DecayProcess)
OPSEQ(DecayProcess, BaseProcess)
OPSEQ(DecayProcess, InteractionProcess)
OPSEQ(DecayProcess, ContinuousProcess)
OPSEQ(DecayProcess, DecayProcess)
template <template <typename, typename> class T, typename A, typename B>
struct is_process_sequence<T<A, B> > {
static const bool value = true;
};
} // namespace corsika::process
......
......@@ -66,7 +66,7 @@ public:
}
};
class Process1 : public DiscreteProcess<Process1> {
class Process1 : public InteractionProcess<Process1> {
public:
Process1(const int v)
: fV(v) {}
......@@ -76,7 +76,7 @@ public:
globalCount++;
}
template <typename D, typename S>
inline EProcessReturn DoDiscrete(D& d, S&) const {
inline EProcessReturn DoInteraction(D& d, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
return EProcessReturn::eOk;
}
......@@ -84,7 +84,7 @@ public:
int fV;
};
class Process2 : public DiscreteProcess<Process2> {
class Process2 : public InteractionProcess<Process2> {
int fV = 0;
public:
......@@ -96,8 +96,8 @@ public:
globalCount++;
}
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
cout << "Process2::DoDiscrete" << endl;
inline EProcessReturn DoInteraction(Particle&, Stack&) const {
cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
......@@ -107,7 +107,7 @@ public:
}
};
class Process3 : public DiscreteProcess<Process3> {
class Process3 : public InteractionProcess<Process3> {
int fV = 0;
public:
......@@ -119,8 +119,8 @@ public:
globalCount++;
}
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const {
cout << "Process3::DoDiscrete" << endl;
inline EProcessReturn DoInteraction(Particle&, Stack&) const {
cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
......@@ -148,7 +148,28 @@ public:
}
// inline double MinStepLength(D& d) {
template <typename Particle, typename Stack>
EProcessReturn DoDiscrete(Particle&, Stack&) const {
EProcessReturn DoInteraction(Particle&, Stack&) const {
return EProcessReturn::eOk;
}
};
class Decay1 : public DecayProcess<Decay1> {
int fV = 0;
public:
Decay1(const int v)
: fV(v) {}
void Init() {
cout << "Decay1::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename Particle>
double GetLifetime(Particle&) const {
return 1;
}
template <typename Particle, typename Stack>
EProcessReturn DoDecay(Particle&, Stack&) const {
return EProcessReturn::eOk;
}
};
......@@ -192,7 +213,21 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
double tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
}
SECTION("lifetime") {
ContinuousProcess1 cp1(0);
Process2 m2(1);
Process3 m3(2);
Decay1 d3(2);
DummyStack s;
const auto sequence2 = cp1 + m2 + m3 + d3;
double tot = sequence2.GetTotalLifetime(s);
double tot_inv = sequence2.GetTotalInverseLifetime(s);
cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl;
}
SECTION("sectionTwo") {
ContinuousProcess1 cp1(0);
......@@ -213,14 +248,14 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
sequence2.DoContinuous(p, t, s);
cout << "-->dodisc" << endl;
sequence2.DoDiscrete(p, s);
// sequence2.DoInteraction(p, s);
cout << "-->done" << endl;
const int nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
for (int i = 0; i < nLoop; ++i) {
sequence2.DoContinuous(p, t, s);
sequence2.DoDiscrete(p, s);
// sequence2.DoInteraction(p, s);
}
for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
cout << "done" << endl;
......
......@@ -67,9 +67,7 @@ namespace corsika::stack {
IncrementSize();
return StackIterator(*this, GetSize() - 1);
}
void Copy(StackIterator& a, StackIterator& b) {
Copy(a.GetIndex(), b.GetIndex());
}
void Copy(StackIterator& a, StackIterator& b) { Copy(a.GetIndex(), b.GetIndex()); }
/// delete this particle
void Delete(StackIterator& p) {
if (GetSize() == 0) { /*error*/
......
......@@ -17,8 +17,8 @@
#include <corsika/setup/SetupTrajectory.h>
#include <limits>
#include <iostream>
#include <limits>
using namespace std;
using namespace corsika;
......@@ -57,7 +57,7 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
}
template <typename Stack>
double StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
double StackInspector<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
return std::numeric_limits<double>::infinity();
}
......
......@@ -36,7 +36,7 @@ namespace corsika::process {
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
// template <typename Particle>
double MinStepLength(Particle&, corsika::setup::Trajectory&) const;
double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
private:
bool fReport;
......
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