IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 450c7e0c authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

correct units in sibyll process interface

parent 2dafac86
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 297 additions and 322 deletions
Documentation/Examples/cascade_example.cc
+ 5
− 5
View file @ 450c7e0c
...@@ -200,7 +200,6 @@ public: ...@@ -200,7 +200,6 @@ public:
private: private:
}; };
int main() { int main() {
corsika::environment::Environment env; // dummy environment corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse()); auto& universe = *(env.GetUniverse());
...@@ -224,13 +223,13 @@ int main() { ...@@ -224,13 +223,13 @@ int main() {
tracking_line::TrackingLine<setup::Stack> tracking(env); tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true); stack_inspector::StackInspector<setup::Stack> p0(true);
corsika::process::sibyll::Interaction/*<setup::Stack>,setup::Trajectory>*/ p1; corsika::process::sibyll::Interaction /*<setup::Stack>,setup::Trajectory>*/ p1;
corsika::process::sibyll::Decay p2; corsika::process::sibyll::Decay p2;
ProcessEMCut p3; ProcessEMCut p3;
const auto sequence = /*p0 +*/ p1 + p2 + p3; const auto sequence = /*p0 +*/ p1 + p2 + p3;
setup::Stack stack; setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack); corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
stack.Clear(); stack.Clear();
auto particle = stack.NewParticle(); auto particle = stack.NewParticle();
...@@ -246,8 +245,9 @@ int main() { ...@@ -246,8 +245,9 @@ int main() {
particle.SetPosition(p); particle.SetPosition(p);
EAS.Init(); EAS.Init();
EAS.Run(); EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV cout << "Result: E0="
//<< "GeV, particles below energy threshold =" << p1.GetCount() << E0 / 1_GeV
//<< "GeV, particles below energy threshold =" << p1.GetCount()
<< endl; << endl;
cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
<< std::endl; << std::endl;
......
Framework/Cascade/Cascade.h
+ 1
− 7
View file @ 450c7e0c
...@@ -18,7 +18,6 @@ ...@@ -18,7 +18,6 @@
#include <corsika/random/UniformRealDistribution.h> #include <corsika/random/UniformRealDistribution.h>
#include <corsika/setup/SetupTrajectory.h> #include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h> #include <corsika/units/PhysicalUnits.h>
#include <corsika/environment/Environment.h>
#include <type_traits> #include <type_traits>
...@@ -58,7 +57,7 @@ namespace corsika::cascade { ...@@ -58,7 +57,7 @@ namespace corsika::cascade {
void Step(Particle& particle) { void Step(Particle& particle) {
using namespace corsika::units::si; using namespace corsika::units::si;
// determine geometric tracking // determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle); corsika::setup::Trajectory step = fTracking.GetTrack(particle);
...@@ -158,12 +157,7 @@ namespace corsika::cascade { ...@@ -158,12 +157,7 @@ namespace corsika::cascade {
Stack& fStack; Stack& fStack;
corsika::environment::Environment const& fEnvironment; corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG = corsika::random::RNG& fRNG =
<<<<<<< HEAD
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade"); corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
=======
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
>>>>>>> e9023467d9ae486f2436418f2004da612ebd02a7
}; };
} // namespace corsika::cascade } // namespace corsika::cascade
......
Framework/Geometry/Line.h
+ 0
− 2
View file @ 450c7e0c
...@@ -31,8 +31,6 @@ namespace corsika::geometry { ...@@ -31,8 +31,6 @@ namespace corsika::geometry {
, v0(pV0) {} , v0(pV0) {}
Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; } Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
Point PositionFromArclength(corsika::units::si::LengthType l) const { return r0 + v0.normalized() * l; }
Point PositionFromArclength(corsika::units::si::LengthType l) const { Point PositionFromArclength(corsika::units::si::LengthType l) const {
return r0 + v0.normalized() * l; return r0 + v0.normalized() * l;
......
Framework/ProcessSequence/ProcessSequence.h
+ 1
− 1
View file @ 450c7e0c
...@@ -361,7 +361,7 @@ namespace corsika::process { ...@@ -361,7 +361,7 @@ namespace corsika::process {
OPSEQ(DecayProcess, DecayProcess) OPSEQ(DecayProcess, DecayProcess)
template <typename A, typename B> template <typename A, typename B>
struct is_process_sequence<corsika::process::ProcessSequence<A, B> > { struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
static const bool value = true; static const bool value = true;
}; };
......
Processes/Sibyll/Decay.h
+ 17
− 17
View file @ 450c7e0c
#ifndef _include_corsika_process_sibyll_decay_h_ #ifndef _include_corsika_process_sibyll_decay_h_
#define _include_corsika_process_sibyll_decay_h_ #define _include_corsika_process_sibyll_decay_h_
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/DecayProcess.h> #include <corsika/process/DecayProcess.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/setup/SetupTrajectory.h> #include <corsika/setup/SetupTrajectory.h>
#include <corsika/particles/ParticleProperties.h> #include <corsika/particles/ParticleProperties.h>
...@@ -57,8 +57,8 @@ namespace corsika::process { ...@@ -57,8 +57,8 @@ namespace corsika::process {
void setTrackedParticlesStable() { void setTrackedParticlesStable() {
/* /*
Sibyll is hadronic generator Sibyll is hadronic generator
only hadrons decay only hadrons decay
*/ */
// set particles unstable // set particles unstable
setHadronsUnstable(); setHadronsUnstable();
...@@ -71,12 +71,12 @@ namespace corsika::process { ...@@ -71,12 +71,12 @@ namespace corsika::process {
ds.NewParticle().SetPID(particles::Code::KMinus); ds.NewParticle().SetPID(particles::Code::KMinus);
ds.NewParticle().SetPID(particles::Code::K0Long); ds.NewParticle().SetPID(particles::Code::K0Long);
ds.NewParticle().SetPID(particles::Code::K0Short); ds.NewParticle().SetPID(particles::Code::K0Short);
for (auto& p : ds) { for (auto& p : ds) {
int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID()); int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
// set particle stable by setting table value negative // set particle stable by setting table value negative
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]); s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
p.Delete(); p.Delete();
} }
} }
...@@ -84,7 +84,7 @@ namespace corsika::process { ...@@ -84,7 +84,7 @@ namespace corsika::process {
public: public:
Decay() {} Decay() {}
void Init() { void Init() {
setHadronsUnstable(); setHadronsUnstable();
setTrackedParticlesStable(); setTrackedParticlesStable();
} }
...@@ -110,11 +110,11 @@ namespace corsika::process { ...@@ -110,11 +110,11 @@ namespace corsika::process {
friend void setHadronsUnstable(); friend void setHadronsUnstable();
template <typename Particle> template <typename Particle>
double GetLifetime(Particle& p) const { corsika::units::si::TimeType GetLifetime(Particle& p) const {
corsika::units::hep::EnergyType E = p.GetEnergy(); corsika::units::hep::EnergyType E = p.GetEnergy();
corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID()); corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
//const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm); // const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
const double gamma = E / m; const double gamma = E / m;
...@@ -126,11 +126,11 @@ namespace corsika::process { ...@@ -126,11 +126,11 @@ namespace corsika::process {
// return as column density // return as column density
// const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm; // const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
// cout << "Decay: MinStep: x0: " << x0 << endl; // cout << "Decay: MinStep: x0: " << x0 << endl;
const double lifetime = gamma*t0 / 1_s; corsika::units::si::TimeType const lifetime = gamma * t0;
cout << "Decay: MinStep: tau: " << lifetime << endl; cout << "Decay: MinStep: tau: " << lifetime << endl;
//int a = 1; // int a = 1;
//const double x = -x0 * log(s_rndm_(a)); // const double x = -x0 * log(s_rndm_(a));
//cout << "Decay: next decay: " << x << endl; // cout << "Decay: next decay: " << x << endl;
return lifetime; return lifetime;
} }
......
Processes/Sibyll/Interaction.h
+ 205
− 221
View file @ 450c7e0c
...@@ -6,278 +6,262 @@ ...@@ -6,278 +6,262 @@
//#include <corsika/setup/SetupStack.h> //#include <corsika/setup/SetupStack.h>
//#include <corsika/setup/SetupTrajectory.h> //#include <corsika/setup/SetupTrajectory.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h> #include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h> #include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/particles/ParticleProperties.h> #include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h> #include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
using namespace corsika; using namespace corsika;
using namespace corsika::process::sibyll; using namespace corsika::process::sibyll;
using namespace corsika::units::si; using namespace corsika::units::si;
namespace corsika::process::sibyll { namespace corsika::process::sibyll {
// template <typename Stack, typename Track> // template <typename Stack, typename Track>
//template <typename Stack> // template <typename Stack>
class Interaction : public corsika::process::InteractionProcess<Interaction> { // <Stack,Track>> { class Interaction
: public corsika::process::InteractionProcess<Interaction> { // <Stack,Track>> {
//typedef typename Stack::ParticleType Particle; // typedef typename Stack::ParticleType Particle;
//typedef typename corsika::setup::Stack::ParticleType Particle; // typedef typename corsika::setup::Stack::ParticleType Particle;
//typedef corsika::setup::Trajectory Track; // typedef corsika::setup::Trajectory Track;
public: public:
Interaction() {} Interaction() {}
~Interaction() {} ~Interaction() {}
void Init() { void Init() {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance(); corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm"); rmng.RegisterRandomStream("s_rndm");
// test random number generator // test random number generator
std::cout << "Interaction: " std::cout << "Interaction: "
<< " test sequence of random numbers." << std::endl; << " test sequence of random numbers." << std::endl;
int a = 0; int a = 0;
for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl; for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
// initialize Sibyll // initialize Sibyll
sibyll_ini_(); sibyll_ini_();
} }
//void setTrackedParticlesStable(); // void setTrackedParticlesStable();
template <typename Particle, typename Track> template <typename Particle, typename Track>
double GetInteractionLength(Particle& p, Track&) const { corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) const {
// coordinate system, get global frame of reference
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
const particles::Code corsikaBeamId = p.GetPID();
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
int kTarget = 16;
EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
super_stupid::MomentumVector Ptot(
rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
// FOR NOW: assume target is at rest
super_stupid::MomentumVector pTarget(
rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
Ptot += p.GetMomentum();
Ptot += pTarget;
// calculate cm. energy
EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * constants::cSquared);
double Ecm = sqs / 1_GeV;
std::cout << "Interaction: "
<< "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kBeam << endl
<< " beam XS code:" << kBeam << endl
<< " beam pid:" << p.GetPID() << endl
<< " target mass number:" << kTarget << std::endl;
double int_length = 0;
if (kInteraction) {
double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
if (kTarget == 1)
sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
else
sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
std::cout << "Interaction: "
<< "MinStep: sibyll return: " << prodCrossSection << std::endl;
CrossSectionType sig = prodCrossSection * 1_mbarn;
std::cout << "Interaction: "
<< "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
std::cout << "Interaction: "
<< "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
// pick random step lenth
std::cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length << std::endl;
} else
int_length = std::numeric_limits<double>::infinity();
/*
what are the units of the output? slant depth or 3space length?
*/
return int_length;
//
// int a = 0;
// const double next_step = -int_length * log(s_rndm_(a));
// std::cout << "Interaction: "
// << "next step (g/cm2): " << next_step << std::endl;
// return next_step;
}
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
cout << "ProcessSibyll: "
<< "DoInteraction: " << p.GetPID() << " interaction? "
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
cout << "defining coordinates" << endl;
// coordinate system, get global frame of reference // coordinate system, get global frame of reference
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS(); CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m}; const particles::Code corsikaBeamId = p.GetPID();
Point pOrig(rootCS, coordinates);
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/* /*
the target should be defined by the Environment, the target should be defined by the Environment,
ideally as full particle object so that the four momenta ideally as full particle object so that the four momenta
and the boosts can be defined.. and the boosts can be defined..
*/
here we need: GetTargetMassNumber() or GetTargetPID()?? // target nuclei: A < 18
GetTargetMomentum() (zero in EAS) // FOR NOW: assume target is oxygen
*/ int kTarget = 16;
// FOR NOW: set target to proton
int kTarget = 1; // env.GetTargetParticle().GetPID(); EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
super_stupid::MomentumVector Ptot(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
cout << "defining target momentum.." << endl; // FOR NOW: assume target is at rest
// FOR NOW: target is always at rest super_stupid::MomentumVector pTarget(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass(); Ptot += p.GetMomentum();
const auto pTarget = super_stupid::MomentumVector( Ptot += pTarget;
rootCS, 0. * 1_GeV / constants::c, 0. * 1_GeV / constants::c, // calculate cm. energy
0. * 1_GeV / constants::c); EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * constants::cSquared);
cout << "target momentum (GeV/c): " double Ecm = sqs / 1_GeV;
<< pTarget.GetComponents() / 1_GeV * constants::c << endl;
cout << "beam momentum (GeV/c): "
<< p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
// get energy of particle from stack
/*
stack is in GeV in lab. frame
convert to GeV in cm. frame
(assuming proton at rest as target AND
assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
*/
// total energy: E_beam + E_target
// in lab. frame: E_beam + m_target*c**2
EnergyType E = p.GetEnergy();
EnergyType Etot = E + Etarget;
// total momentum
super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
// invariant mass, i.e. cm. energy
EnergyType Ecm = sqrt(Etot * Etot -
Ptot.squaredNorm() *
constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
/*
get transformation between Stack-frame and SibStack-frame
for EAS Stack-frame is lab. frame, could be different for CRMC-mode
the transformation should be derived from the input momenta
*/
const double gamma = Etot / Ecm;
const auto gambet = Ptot / (Ecm / constants::c);
std::cout << "Interaction: " std::cout << "Interaction: "
<< " DoDiscrete: gamma:" << gamma << endl; << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
std::cout << "Interaction: " << " beam can interact:" << kBeam << endl
<< " DoDiscrete: gambet:" << gambet.GetComponents() << endl; << " beam XS code:" << kBeam << endl
<< " beam pid:" << p.GetPID() << endl
<< " target mass number:" << kTarget << std::endl;
int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID()); if (kInteraction) {
std::cout << "Interaction: " double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
<< " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV double dumdif[3];
<< std::endl;
if (E < 8.5_GeV || Ecm < 10_GeV) { if (kTarget == 1)
sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
else
sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
std::cout << "Interaction: "
<< "MinStep: sibyll return: " << prodCrossSection << std::endl;
CrossSectionType sig = prodCrossSection * 1_mbarn;
std::cout << "Interaction: "
<< "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
const MassType nucleon_mass =
0.93827_GeV / corsika::units::si::constants::cSquared;
std::cout << "Interaction: "
<< "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
GrammageType int_length = kTarget * nucleon_mass / sig;
// pick random step lenth
std::cout << "Interaction: " std::cout << "Interaction: "
<< " DoInteraction: dropping particle.." << std::endl; << "interaction length (g/cm2): " << int_length << std::endl;
p.Delete();
return int_length;
} else { } else {
// Sibyll does not know about units.. return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
double sqs = Ecm / 1_GeV; }
// running sibyll, filling stack }
sibyll_(kBeam, kTarget, sqs);
// running decays
//setTrackedParticlesStable();
decsib_();
// print final state
int print_unit = 6;
sib_list_(print_unit);
// delete current particle
p.Delete();
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
// SibStack does not know about momentum yet so we need counter to access momentum
// array in Sibyll
int i = -1;
super_stupid::MomentumVector Ptot_final(
rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
for (auto& psib : ss) {
++i;
// skip particles that have decayed in Sibyll
if (abs(s_plist_.llist[i]) > 100) continue;
// transform energy to lab. frame, primitve
// compute beta_vec * p_vec
// arbitrary Lorentz transformation based on sibyll routines
const auto gammaBetaComponents = gambet.GetComponents();
const auto pSibyllComponents = psib.GetMomentum().GetComponents();
EnergyType en_lab = 0. * 1_GeV;
MomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma;
EnergyType pnorm = 0. * 1_GeV;
for (int j = 0; j < 3; ++j)
pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * constants::c) /
(gamma + 1.);
pnorm += psib.GetEnergy();
for (int j = 0; j < 3; ++j) {
p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
gammaBetaComponents[j] /
constants::c;
en_lab -=
(-1) * pSibyllComponents[j] * gammaBetaComponents[j] * constants::c;
}
// add to corsika stack template <typename Particle, typename Stack>
auto pnew = s.NewParticle(); corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
pnew.SetEnergy(en_lab); cout << "ProcessSibyll: "
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID())); << "DoInteraction: " << p.GetPID() << " interaction? "
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
cout << "defining coordinates" << endl;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point pOrig(rootCS, coordinates);
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
here we need: GetTargetMassNumber() or GetTargetPID()??
GetTargetMomentum() (zero in EAS)
*/
// FOR NOW: set target to proton
int kTarget = 1; // env.GetTargetParticle().GetPID();
cout << "defining target momentum.." << endl;
// FOR NOW: target is always at rest
const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
const auto pTarget = super_stupid::MomentumVector(
rootCS, 0. * 1_GeV / constants::c, 0. * 1_GeV / constants::c,
0. * 1_GeV / constants::c);
cout << "target momentum (GeV/c): "
<< pTarget.GetComponents() / 1_GeV * constants::c << endl;
cout << "beam momentum (GeV/c): "
<< p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
// get energy of particle from stack
/*
stack is in GeV in lab. frame
convert to GeV in cm. frame
(assuming proton at rest as target AND
assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
*/
// total energy: E_beam + E_target
// in lab. frame: E_beam + m_target*c**2
EnergyType E = p.GetEnergy();
EnergyType Etot = E + Etarget;
// total momentum
super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
// invariant mass, i.e. cm. energy
EnergyType Ecm =
sqrt(Etot * Etot - Ptot.squaredNorm() *
constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
/*
get transformation between Stack-frame and SibStack-frame
for EAS Stack-frame is lab. frame, could be different for CRMC-mode
the transformation should be derived from the input momenta
*/
const double gamma = Etot / Ecm;
const auto gambet = Ptot / (Ecm / constants::c);
std::cout << "Interaction: "
<< " DoDiscrete: gamma:" << gamma << endl;
std::cout << "Interaction: "
<< " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
corsika::geometry::QuantityVector<momentum_d> p_lab_c{ std::cout << "Interaction: "
p_lab_components[0], p_lab_components[1], p_lab_components[2]}; << " DoInteraction: E(GeV):" << E / 1_GeV
pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c)); << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
Ptot_final += pnew.GetMomentum(); if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
<< " DoInteraction: dropping particle.." << std::endl;
p.Delete();
} else {
// Sibyll does not know about units..
double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, kTarget, sqs);
// running decays
// setTrackedParticlesStable();
decsib_();
// print final state
int print_unit = 6;
sib_list_(print_unit);
// delete current particle
p.Delete();
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
// SibStack does not know about momentum yet so we need counter to access
// momentum array in Sibyll
int i = -1;
super_stupid::MomentumVector Ptot_final(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
for (auto& psib : ss) {
++i;
// skip particles that have decayed in Sibyll
if (abs(s_plist_.llist[i]) > 100) continue;
// transform energy to lab. frame, primitve
// compute beta_vec * p_vec
// arbitrary Lorentz transformation based on sibyll routines
const auto gammaBetaComponents = gambet.GetComponents();
const auto pSibyllComponents = psib.GetMomentum().GetComponents();
EnergyType en_lab = 0. * 1_GeV;
MomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma;
EnergyType pnorm = 0. * 1_GeV;
for (int j = 0; j < 3; ++j)
pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * constants::c) /
(gamma + 1.);
pnorm += psib.GetEnergy();
for (int j = 0; j < 3; ++j) {
p_lab_components[j] = pSibyllComponents[j] -
(-1) * pnorm * gammaBetaComponents[j] / constants::c;
en_lab -=
(-1) * pSibyllComponents[j] * gammaBetaComponents[j] * constants::c;
}
// add to corsika stack
auto pnew = s.NewParticle();
pnew.SetEnergy(en_lab);
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
corsika::geometry::QuantityVector<momentum_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]};
pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
Ptot_final += pnew.GetMomentum();
}
// cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
// * constants::c << endl;
} }
// cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV *
// constants::c << endl;
} }
return process::EProcessReturn::eOk;
} }
return process::EProcessReturn::eOk;
}
}; };
} } // namespace corsika::process::sibyll
#endif #endif
Processes/Sibyll/ParticleConversion.h
+ 1
− 1
View file @ 450c7e0c
...@@ -25,7 +25,7 @@ namespace corsika::process::sibyll { ...@@ -25,7 +25,7 @@ namespace corsika::process::sibyll {
#include <corsika/process/sibyll/Generated.inc> #include <corsika/process/sibyll/Generated.inc>
bool constexpr KnownBySibyll(corsika::particles::Code pCode) { bool constexpr KnownBySibyll(corsika::particles::Code pCode) {
return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)]; return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)];
} }
......
Processes/Sibyll/SibStack.h
+ 67
− 67
View file @ 450c7e0c
#ifndef _include_sibstack_h_ #ifndef _include_sibstack_h_
#define _include_sibstack_h_ #define _include_sibstack_h_
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h> #include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/stack/Stack.h> #include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h> #include <corsika/units/PhysicalUnits.h>
...@@ -14,75 +14,75 @@ using namespace corsika::units::si; ...@@ -14,75 +14,75 @@ using namespace corsika::units::si;
using namespace corsika::geometry; using namespace corsika::geometry;
namespace corsika::process::sibyll { namespace corsika::process::sibyll {
class SibStackData {
public: class SibStackData {
void Init();
public:
void Init();
void Clear() { s_plist_.np = 0; } void Clear() { s_plist_.np = 0; }
int GetSize() const { return s_plist_.np; } int GetSize() const { return s_plist_.np; }
#warning check actual capacity of sibyll stack #warning check actual capacity of sibyll stack
int GetCapacity() const { return 8000; } int GetCapacity() const { return 8000; }
void SetId(const int i, const int v) { s_plist_.llist[i] = v; } void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; } void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
void SetMomentum(const int i, const super_stupid::MomentumVector& v) { void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
auto tmp = v.GetComponents(); auto tmp = v.GetComponents();
for (int idx = 0; idx < 3; ++idx) for (int idx = 0; idx < 3; ++idx)
s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c; s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
} }
int GetId(const int i) const { return s_plist_.llist[i]; } int GetId(const int i) const { return s_plist_.llist[i]; }
EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; } EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
super_stupid::MomentumVector GetMomentum(const int i) const { super_stupid::MomentumVector GetMomentum(const int i) const {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS(); CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
corsika::geometry::QuantityVector<momentum_d> components{ corsika::geometry::QuantityVector<momentum_d> components{
s_plist_.p[0][i] * 1_GeV / constants::c, s_plist_.p[0][i] * 1_GeV / constants::c,
s_plist_.p[1][i] * 1_GeV / constants::c, s_plist_.p[1][i] * 1_GeV / constants::c,
s_plist_.p[2][i] * 1_GeV / constants::c}; s_plist_.p[2][i] * 1_GeV / constants::c};
super_stupid::MomentumVector v1(rootCS, components); super_stupid::MomentumVector v1(rootCS, components);
return v1; return v1;
} }
void Copy(const int i1, const int i2) { void Copy(const int i1, const int i2) {
s_plist_.llist[i1] = s_plist_.llist[i2]; s_plist_.llist[i1] = s_plist_.llist[i2];
s_plist_.p[3][i1] = s_plist_.p[3][i2]; s_plist_.p[3][i1] = s_plist_.p[3][i2];
} }
protected: protected:
void IncrementSize() { s_plist_.np++; } void IncrementSize() { s_plist_.np++; }
void DecrementSize() { void DecrementSize() {
if (s_plist_.np > 0) { s_plist_.np--; } if (s_plist_.np > 0) { s_plist_.np--; }
} }
}; };
template <typename StackIteratorInterface> template <typename StackIteratorInterface>
class ParticleInterface : public ParticleBase<StackIteratorInterface> { class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetStackData; using ParticleBase<StackIteratorInterface>::GetStackData;
using ParticleBase<StackIteratorInterface>::GetIndex; using ParticleBase<StackIteratorInterface>::GetIndex;
public: public:
void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); } void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); } EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); } void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const { corsika::process::sibyll::SibyllCode GetPID() const {
return static_cast<corsika::process::sibyll::SibyllCode>( return static_cast<corsika::process::sibyll::SibyllCode>(
GetStackData().GetId(GetIndex())); GetStackData().GetId(GetIndex()));
} }
super_stupid::MomentumVector GetMomentum() const { super_stupid::MomentumVector GetMomentum() const {
return GetStackData().GetMomentum(GetIndex()); return GetStackData().GetMomentum(GetIndex());
} }
void SetMomentum(const super_stupid::MomentumVector& v) { void SetMomentum(const super_stupid::MomentumVector& v) {
GetStackData().SetMomentum(GetIndex(), v); GetStackData().SetMomentum(GetIndex(), v);
} }
}; };
typedef Stack<SibStackData, ParticleInterface> SibStack; typedef Stack<SibStackData, ParticleInterface> SibStack;
} } // namespace corsika::process::sibyll
#endif #endif
Processes/Sibyll/sibyll2.3c.cc
+ 0
− 1
View file @ 450c7e0c
...@@ -2,7 +2,6 @@ ...@@ -2,7 +2,6 @@
#include <corsika/random/RNGManager.h> #include <corsika/random/RNGManager.h>
double s_rndm_(int&) { double s_rndm_(int&) {
static corsika::random::RNG& rmng = static corsika::random::RNG& rmng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm"); corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
......
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