From f18e43c226d5a7ee8a71372bfdb910eee50903b0 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 23 Jan 2019 18:26:18 +0000
Subject: [PATCH] started new process
---
Documentation/Examples/cascade_example.cc | 4 +-
Processes/Sibyll/NuclearInteraction.h | 46 ++++++++---------------
Processes/Sibyll/sibyll2.3c.h | 5 ++-
Processes/Sibyll/testSibyll.cc | 28 +++++++++++++-
4 files changed, 47 insertions(+), 36 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index fd6055f83..3f7047375 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -270,8 +270,6 @@ int main() {
double phi = 0.;
{
- auto particle = stack.NewParticle();
- particle.SetPID(beamCode);
auto elab2plab = []( HEPEnergyType Elab, HEPMassType m){
return sqrt(Elab * Elab - m * m);
};
@@ -286,7 +284,7 @@ int main() {
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
- stack.AddParticle(Code::Proton, E0, plab, pos, 0_ns);
+ stack.AddParticle(beamCode, E0, plab, pos, 0_ns);
}
// define air shower object, run simulation
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 63ee3aa13..917cda8c3 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -46,7 +46,7 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize hadronic interaction module
- //sibyll_ini_();
+ sibyll_ini_();
// initialize nuclib
nuc_nuc_ini_();
@@ -58,17 +58,16 @@ namespace corsika::process::sibyll {
using namespace corsika::units::si;
double sigProd, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
-
- corsika::particles::Code BeamIdToUse;
- if(corsika::particles::IsNucleus(BeamId)){
-
- // TODO: use nuclib to calc. nuclear cross sections
- // FOR NOW: use proton cross section for nuclei
- BeamIdToUse = corsika::particles::Proton::GetCode();
- std::cout << "WARNING: replacing beam nucleus with proton!" << std::endl;
- } else {
- BeamIdToUse = BeamId;
+
+ if(!corsika::particles::IsNucleus(BeamId)){
+ sigProd = std::numeric_limits<double>::infinity();
+ return std::make_tuple(sigProd * 1_mbarn, 1);
}
+
+ // TODO: use nuclib to calc. nuclear cross sections
+ // FOR NOW: use proton cross section for nuclei
+ auto const BeamIdToUse = corsika::particles::Proton::GetCode();
+ std::cout << "WARNING: replacing beam nucleus with proton!" << std::endl;
const int iBeam = process::sibyll::GetSibyllXSCode(BeamIdToUse);
const double dEcm = CoMenergy / 1_GeV;
@@ -183,7 +182,7 @@ namespace corsika::process::sibyll {
}
template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
+ corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
@@ -424,33 +423,26 @@ namespace corsika::process::sibyll {
// put nuclear fragments on corsika stack
for(int j=0; j<nFragments; ++j){
- auto pnew = s.NewParticle();
// here we need the nucleonNumber to corsika Id conversion
// A = AFragments[j]
// pnew.SetPID( corsika::particles::GetCode( corsika::particles::Nucleus(A,Z) ) );
auto pCode = Nucleus( AFragments[j] );
- pnew.SetPID(pCode);
// CORSIKA 7 way
// spectators inherit momentum from original projectile
const double mass_ratio = corsika::particles::GetMass( pCode ) / corsika::particles::GetMass( corsikaProjId );
auto const Plab = PprojLab * mass_ratio;
-
- pnew.SetEnergy(Plab.GetTimeLikeComponent());
- pnew.SetMomentum(Plab.GetSpaceLikeComponents());
- pnew.SetPosition(pOrig);
- pnew.SetTime(tOrig);
+ p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
}
// add elastic nucleons to corsika stack
for(int j=0; j<nElasticNucleons; ++j){
- auto pnew = s.NewParticle();
// TODO: sample proton or neutron
auto pCode = corsika::particles::Proton::GetCode();
- pnew.SetPID( pCode );
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
@@ -458,10 +450,7 @@ namespace corsika::process::sibyll {
const double mass_ratio = corsika::particles::GetMass( pCode ) / corsika::particles::GetMass( corsikaProjId );
auto const Plab = PprojLab * mass_ratio;
- pnew.SetEnergy(Plab.GetTimeLikeComponent());
- pnew.SetMomentum(Plab.GetSpaceLikeComponents());
- pnew.SetPosition(pOrig);
- pnew.SetTime(tOrig);
+ p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
@@ -498,12 +487,7 @@ namespace corsika::process::sibyll {
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- pnew.SetEnergy(Plab.GetTimeLikeComponent());
- pnew.SetMomentum(Plab.GetSpaceLikeComponents());
- pnew.SetPosition(pOrig);
- pnew.SetTime(tOrig);
+ auto pnew = p.AddSecondary(process::sibyll::ConvertFromSibyll(psib.GetPID()), Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig);
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
index 639c5a600..df5a57609 100644
--- a/Processes/Sibyll/sibyll2.3c.h
+++ b/Processes/Sibyll/sibyll2.3c.h
@@ -101,7 +101,10 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
-
+
+// subroutine to initiate nuclib
+void nuc_nuc_ini_();
+
// subroutine to SET DECAYS
void dec_ini_();
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index d421e98de..795ba53aa 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -136,7 +136,33 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
- auto particle = stack.NewParticle();
+ const HEPEnergyType E0 = 10_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+ auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ geometry::Point pos(cs, 0_m, 0_m, 0_m);
+
+ auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
+
+ Interaction model(env);
+
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret =
+ model.DoInteraction(particle, stack);
+ [[maybe_unused]] const GrammageType length =
+ model.GetInteractionLength(particle, track);
+ }
+
+ SECTION("DecayInterface") {
+
+ setup::Stack stack;
+ const HEPEnergyType E0 = 10_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+ auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ geometry::Point pos(cs, 0_m, 0_m, 0_m);
+
+ auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
NuclearInteraction model(env);
--
GitLab