IAP GITLAB

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

added all of Felix patches

parent 19ecfbbd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 8
− 20
View file @ 92137f4c
...@@ -119,8 +119,10 @@ public: ...@@ -119,8 +119,10 @@ public:
template <typename Particle> template <typename Particle>
LengthType MaxStepLength(Particle& p, setup::Trajectory&) const { 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(); const Code pid = p.GetPID();
if (isEmParticle(pid) || isInvisible(pid)) { if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl; cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
return 0_m; return 0_m;
} else { } else {
...@@ -134,40 +136,26 @@ public: ...@@ -134,40 +136,26 @@ public:
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const { EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl; cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
const Code pid = p.GetPID(); const Code pid = p.GetPID();
EProcessReturn ret = EProcessReturn::eOk;
if (isEmParticle(pid)) { if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl; cout << "removing em. particle..." << endl;
fEmEnergy += p.GetEnergy(); fEmEnergy += p.GetEnergy();
fEmCount += 1; fEmCount += 1;
p.Delete(); p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) { } else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl; cout << "removing inv. particle..." << endl;
fInvEnergy += p.GetEnergy(); fInvEnergy += p.GetEnergy();
fInvCount += 1; fInvCount += 1;
p.Delete(); p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) { } else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl; cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy(); fEnergy += p.GetEnergy();
p.Delete(); p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} }
// cout << "ProcessCut: DoContinous: " << p.GetPID() << endl; return ret;
// cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
// cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
// const Code pid = p.GetPID();
// if( isEmParticle( pid ) ){
// cout << "removing em. particle..." << endl;
// fEmEnergy += p.GetEnergy();
// fEmCount += 1;
// p.Delete();
// return EProcessReturn::eParticleAbsorbed;
// }
// if ( isInvisible( pid ) ){
// cout << "removing inv. particle..." << endl;
// fInvEnergy += p.GetEnergy();
// fInvCount += 1;
// p.Delete();
// return EProcessReturn::eParticleAbsorbed;
// }
return EProcessReturn::eOk;
} }
void Init() { void Init() {
......
Processes/Sibyll/Decay.h
+ 26
− 8
View file @ 92137f4c
...@@ -42,13 +42,23 @@ namespace corsika::process { ...@@ -42,13 +42,23 @@ namespace corsika::process {
} }
} }
void setUnstable(const corsika::particles::Code pCode) const {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
void setStable(const corsika::particles::Code pCode) const {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
void setAllStable() const { void setAllStable() const {
// name? also makes EM particles stable
// using namespace corsika::io;
using std::cout; using std::cout;
using std::endl; using std::endl;
// name? also makes EM particles stable std::cout << "Decay: setting all particles stable.." << std::endl;
// loop over all particles in sibyll // loop over all particles in sibyll
// should be changed to loop over human readable list // should be changed to loop over human readable list
...@@ -58,11 +68,7 @@ namespace corsika::process { ...@@ -58,11 +68,7 @@ namespace corsika::process {
const int sibCode = (int)p; const int sibCode = (int)p;
// skip unknown and antiparticles // skip unknown and antiparticles
if (sibCode < 1) continue; if (sibCode < 1) continue;
const corsika::particles::Code pid =
ConvertFromSibyll(static_cast<SibyllCode>(sibCode));
std::cout << "Sibyll: Decay: setting " << pid << " stable" << std::endl;
s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]); s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
std::cout << "decay table value: " << s_csydec_.idb[sibCode - 1] << std::endl;
} }
} }
...@@ -124,6 +130,7 @@ namespace corsika::process { ...@@ -124,6 +130,7 @@ namespace corsika::process {
corsika::particles::GetLifetime(p.GetPID()); corsika::particles::GetLifetime(p.GetPID());
cout << "Decay: code: " << p.GetPID() << endl; cout << "Decay: code: " << p.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl; cout << "Decay: MinStep: t0: " << t0 << endl;
cout << "Decay: MinStep: energy: " << E / 1_GeV << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl; cout << "Decay: MinStep: gamma: " << gamma << endl;
// cout << "Decay: MinStep: density: " << density << endl; // cout << "Decay: MinStep: density: " << density << endl;
// return as column density // return as column density
...@@ -139,20 +146,29 @@ namespace corsika::process { ...@@ -139,20 +146,29 @@ namespace corsika::process {
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
void DoDecay(Particle& p, Stack& s) const { void DoDecay(Particle& p, Stack& s) const {
using corsika::geometry::Point;
using namespace corsika::units::si;
SibStack ss; SibStack ss;
ss.Clear(); ss.Clear();
// copy particle to sibyll stack // copy particle to sibyll stack
auto pin = ss.NewParticle(); auto pin = ss.NewParticle();
pin.SetPID(process::sibyll::ConvertToSibyllRaw(p.GetPID())); const corsika::particles::Code pCode = p.GetPID();
pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
pin.SetEnergy(p.GetEnergy()); pin.SetEnergy(p.GetEnergy());
pin.SetMomentum(p.GetMomentum()); pin.SetMomentum(p.GetMomentum());
// remember position
Point decayPoint = p.GetPosition();
TimeType t0 = p.GetTime();
// remove original particle from corsika stack // remove original particle from corsika stack
p.Delete(); p.Delete();
// set all particles/hadrons unstable // set all particles/hadrons unstable
setHadronsUnstable(); // setHadronsUnstable();
setUnstable(pCode);
// call sibyll decay // call sibyll decay
std::cout << "Decay: calling Sibyll decay routine.." << std::endl; std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
decsib_(); decsib_();
// reset to stable
setStable(pCode);
// print output // print output
int print_unit = 6; int print_unit = 6;
sib_list_(print_unit); sib_list_(print_unit);
...@@ -169,6 +185,8 @@ namespace corsika::process { ...@@ -169,6 +185,8 @@ namespace corsika::process {
pnew.SetEnergy(psib.GetEnergy()); pnew.SetEnergy(psib.GetEnergy());
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID())); pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
pnew.SetMomentum(psib.GetMomentum()); pnew.SetMomentum(psib.GetMomentum());
pnew.SetPosition(decayPoint);
pnew.SetTime(t0);
} }
// empty sibyll stack // empty sibyll stack
ss.Clear(); ss.Clear();
......
Processes/Sibyll/Interaction.h
+ 15
− 11
View file @ 92137f4c
...@@ -56,9 +56,8 @@ namespace corsika::process::sibyll { ...@@ -56,9 +56,8 @@ namespace corsika::process::sibyll {
// beam particles for sibyll : 1, 2, 3 for p, pi, k // beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table // read from cross section code table
int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId); const int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/* /*
the target should be defined by the Environment, the target should be defined by the Environment,
...@@ -67,7 +66,7 @@ namespace corsika::process::sibyll { ...@@ -67,7 +66,7 @@ namespace corsika::process::sibyll {
*/ */
// target nuclei: A < 18 // target nuclei: A < 18
// FOR NOW: assume target is oxygen // FOR NOW: assume target is oxygen
int kTarget = 16; const int kTarget = corsika::particles::Oxygen::GetNucleusA();
hep::EnergyType Etot = hep::EnergyType Etot =
p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass(); p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
...@@ -77,8 +76,8 @@ namespace corsika::process::sibyll { ...@@ -77,8 +76,8 @@ namespace corsika::process::sibyll {
Ptot += p.GetMomentum(); Ptot += p.GetMomentum();
Ptot += pTarget; Ptot += pTarget;
// calculate cm. energy // calculate cm. energy
hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm()); const hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
double Ecm = sqs / 1_GeV; const double Ecm = sqs / 1_GeV;
std::cout << "Interaction: " std::cout << "Interaction: "
<< "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
...@@ -149,8 +148,8 @@ namespace corsika::process::sibyll { ...@@ -149,8 +148,8 @@ namespace corsika::process::sibyll {
CoordinateSystem& rootCS = CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m}; Point pOrig = p.GetPosition();
Point pOrig(rootCS, coordinates); TimeType tOrig = p.GetTime();
/* /*
the target should be defined by the Environment, the target should be defined by the Environment,
...@@ -160,8 +159,9 @@ namespace corsika::process::sibyll { ...@@ -160,8 +159,9 @@ namespace corsika::process::sibyll {
here we need: GetTargetMassNumber() or GetTargetPID()?? here we need: GetTargetMassNumber() or GetTargetPID()??
GetTargetMomentum() (zero in EAS) GetTargetMomentum() (zero in EAS)
*/ */
// FOR NOW: set target to proton // FOR NOW: set target to oxygen
int kTarget = 1; // env.GetTargetParticle().GetPID(); const int kTarget = corsika::particles::Oxygen::
GetNucleusA(); // env.GetTargetParticle().GetPID();
cout << "defining target momentum.." << endl; cout << "defining target momentum.." << endl;
// FOR NOW: target is always at rest // FOR NOW: target is always at rest
...@@ -171,6 +171,8 @@ namespace corsika::process::sibyll { ...@@ -171,6 +171,8 @@ namespace corsika::process::sibyll {
<< pTarget.GetComponents() / 1_GeV * constants::c << endl; << pTarget.GetComponents() / 1_GeV * constants::c << endl;
cout << "beam momentum (GeV/c): " cout << "beam momentum (GeV/c): "
<< p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl; << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "time: " << tOrig << endl;
// get energy of particle from stack // get energy of particle from stack
/* /*
...@@ -209,7 +211,7 @@ namespace corsika::process::sibyll { ...@@ -209,7 +211,7 @@ namespace corsika::process::sibyll {
if (E < 8.5_GeV || Ecm < 10_GeV) { if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: " std::cout << "Interaction: "
<< " DoInteraction: dropping particle.." << std::endl; << " DoInteraction: dropping particle.." << std::endl;
p.Delete(); // p.Delete(); delete later... different process
} else { } else {
// Sibyll does not know about units.. // Sibyll does not know about units..
double sqs = Ecm / 1_GeV; double sqs = Ecm / 1_GeV;
...@@ -265,6 +267,8 @@ namespace corsika::process::sibyll { ...@@ -265,6 +267,8 @@ namespace corsika::process::sibyll {
corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{ corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]}; p_lab_components[0], p_lab_components[1], p_lab_components[2]};
pnew.SetMomentum(MomentumVector(rootCS, p_lab_c)); pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
Ptot_final += pnew.GetMomentum(); Ptot_final += pnew.GetMomentum();
} }
// cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
......
Processes/Sibyll/sibyll2.3c.h
+ 4
− 3
View file @ 92137f4c
...@@ -59,7 +59,7 @@ extern struct { ...@@ -59,7 +59,7 @@ extern struct {
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_; // extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
// sibyll main subroutine // sibyll main subroutine
void sibyll_(int&, int&, double&); void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll // subroutine to initiate sibyll
void sibyll_ini_(); void sibyll_ini_();
...@@ -79,8 +79,9 @@ void decsib_(); ...@@ -79,8 +79,9 @@ void decsib_();
// interaction length // interaction length
// double fpni_(double&, int&); // double fpni_(double&, int&);
void sib_sigma_hnuc_(int&, int&, double&, double&, double&); void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&);
void sib_sigma_hp_(int&, double&, double&, double&, double&, double*, double&, double&); void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
double&);
double s_rndm_(int&); double s_rndm_(int&);
......
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