IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 820f0123 authored by Felix Riehn's avatar Felix Riehn
Browse files

Merge branch 'sibyll_decays' into sibyll

parents 9ca69001 e1bfd98d
No related branches found
No related tags found
No related merge requests found
......@@ -23,6 +23,8 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/ProcessDecay.h>
#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
......@@ -42,6 +44,35 @@ class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
ProcessSplit() {}
void setTrackedParticlesStable() const {
/*
Sibyll is hadronic generator
only hadrons decay
*/
// set particles unstable
corsika::process::sibyll::setHadronsUnstable();
// make tracked particles stable
std::cout << "ProcessSplit: setting tracked hadrons stable.." << std::endl;
setup::Stack ds;
ds.NewParticle().SetPID(Code::PiPlus);
ds.NewParticle().SetPID(Code::PiMinus);
ds.NewParticle().SetPID(Code::KPlus);
ds.NewParticle().SetPID(Code::KMinus);
ds.NewParticle().SetPID(Code::K0Long);
ds.NewParticle().SetPID(Code::K0Short);
for (auto& p : ds) {
int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
// set particle stable by setting table value negative
// cout << "ProcessSplit: doDiscrete: setting " << p.GetPID() << "(" << s_id << ")"
// << " stable in Sibyll .." << endl;
s_csydec_.idb[s_id - 1] = (-1) * abs( s_csydec_.idb[s_id - 1] );
p.Delete();
}
}
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
......@@ -210,7 +241,8 @@ public:
// running sibyll, filling stack
sibyll_( kBeam, kTarget, sqs);
// running decays
//decsib_();
setTrackedParticlesStable();
decsib_();
// print final state
int print_unit = 6;
sib_list_( print_unit );
......@@ -227,6 +259,9 @@ public:
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
......@@ -299,29 +334,10 @@ public:
// initialize Sibyll
sibyll_ini_();
// set particles stable / unstable
// use stack to loop over particles
setup::Stack ds;
ds.NewParticle().SetPID(Code::Proton);
ds.NewParticle().SetPID(Code::Neutron);
ds.NewParticle().SetPID(Code::PiPlus);
ds.NewParticle().SetPID(Code::PiMinus);
ds.NewParticle().SetPID(Code::Pi0);
ds.NewParticle().SetPID(Code::KPlus);
ds.NewParticle().SetPID(Code::KMinus);
ds.NewParticle().SetPID(Code::K0Long);
ds.NewParticle().SetPID(Code::K0Short);
for (auto& p : ds) {
int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
// set particle stable by setting table value negative
cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")"
<< " stable in Sibyll .." << endl;
s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1];
p.Delete();
}
setTrackedParticlesStable();
}
int GetCount() { return fCount; }
private:
......@@ -334,6 +350,7 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
int main() {
// coordinate system, get global frame of reference
......@@ -346,7 +363,8 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
corsika::process::sibyll::ProcessDecay p2;
const auto sequence = p0 + p1 + p2;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
......
......@@ -6,6 +6,7 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/stack/Stack.h>
using namespace std;
......@@ -64,11 +65,12 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v ); }
EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const { return static_cast<corsika::process::sibyll::SibyllCode> (GetStackData().GetId(GetIndex())); }
super_stupid::MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
void SetMomentum(const super_stupid::MomentumVector& v) { GetStackData().SetMomentum(GetIndex(), v); }
};
......
......@@ -36,9 +36,10 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::electric_charge_d, double>;
using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
using MomentumType = phys::units::quantity<momentum_d, double>;
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
using MomentumType = phys::units::quantity<momentum_d, double>;
} // end namespace corsika::units::si
/**
......
......@@ -20,6 +20,7 @@ set (
set (
MODEL_HEADERS
ParticleConversion.h
ProcessDecay.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
......
#ifndef _include_ProcessDecay_h_
#define _include_ProcessDecay_h_
#include <corsika/process/ContinuousProcess.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/cascade/SibStack.h>
//using namespace corsika::particles;
namespace corsika::process {
namespace sibyll {
void setHadronsUnstable(){
// name? also makes EM particles stable
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. corsika::particles::ListOfParticles()
std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
// make ALL particles unstable, then set EM stable
for(auto &p: corsika2sibyll){
//std::cout << (int)p << std::endl;
const int sibCode = (int)p;
// skip unknown and antiparticles
if( sibCode< 1) continue;
//std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll( static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
s_csydec_.idb[ sibCode - 1 ] = abs( s_csydec_.idb[ sibCode - 1 ] );
//std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] << std::endl;
}
// set Leptons and Proton and Neutron stable
// use stack to loop over particles
setup::Stack ds;
ds.NewParticle().SetPID(corsika::particles::Code::Proton);
ds.NewParticle().SetPID(corsika::particles::Code::Neutron);
ds.NewParticle().SetPID(corsika::particles::Code::Electron);
ds.NewParticle().SetPID(corsika::particles::Code::Positron);
ds.NewParticle().SetPID(corsika::particles::Code::NuE);
ds.NewParticle().SetPID(corsika::particles::Code::NuEBar);
ds.NewParticle().SetPID(corsika::particles::Code::MuMinus);
ds.NewParticle().SetPID(corsika::particles::Code::MuPlus);
ds.NewParticle().SetPID(corsika::particles::Code::NuMu);
ds.NewParticle().SetPID(corsika::particles::Code::NuMuBar);
for (auto& p : ds) {
int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
// set particle stable by setting table value negative
// cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
// << " stable in Sibyll .." << endl;
s_csydec_.idb[ s_id - 1 ] = (-1) * abs(s_csydec_.idb[ s_id - 1 ]);
p.Delete();
}
}
class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
public:
ProcessDecay() {}
void Init() {
//setHadronsUnstable();
}
void setAllStable(){
// name? also makes EM particles stable
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. corsika::particles::ListOfParticles()
for(auto &p: corsika2sibyll){
//std::cout << (int)p << std::endl;
const int sibCode = (int)p;
// skip unknown and antiparticles
if( sibCode< 1) continue;
std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll( static_cast<SibyllCode> ( sibCode ) ) << " stable" << std::endl;
s_csydec_.idb[ sibCode - 1 ] = -1 * abs( s_csydec_.idb[ sibCode - 1 ] );
std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] << std::endl;
}
}
friend void setHadronsUnstable();
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
EnergyType E = p.GetEnergy();
MassType m = corsika::particles::GetMass(p.GetPID());
// env.GetDensity();
const MassDensityType density = 1.25e-3 * kilogram / ( 1_cm * 1_cm * 1_cm );
const double gamma = E / m / constants::cSquared;
TimeType t0 = GetLifetime( p.GetPID() );
cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
cout << "ProcessDecay: MinStep: density: " << density << endl;
// return as column density
const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
return x0;
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
SibStack ss;
ss.Clear();
// copy particle to sibyll stack
auto pin = ss.NewParticle();
pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) );
pin.SetEnergy( p.GetEnergy() );
pin.SetMomentum( p.GetMomentum() );
// set all particles/hadrons unstable
setHadronsUnstable();
// call sibyll decay
std::cout << "calling Sibyll decay routine.." << std::endl;
decsib_();
// print output
//int print_unit = 6;
//sib_list_( print_unit );
// copy particles from sibyll stack to corsika
int i = -1;
for (auto &psib: ss){
++i;
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
if( abs(s_plist_.llist[ i ]) > 100 ) continue;
// add to corsika stack
//cout << "decay product: " << process::sibyll::ConvertFromSibyll( psib.GetPID() ) << endl;
auto pnew = s.NewParticle();
pnew.SetEnergy( psib.GetEnergy() );
pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
pnew.SetMomentum( psib.GetMomentum() );
}
// empty sibyll stack
ss.Clear();
// remove original particle from stack
p.Delete();
}
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
return EProcessReturn::eOk;
}
};
}
}
#endif
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