From 361efe9659d2b0b4e545d2af11c31889ce97244e Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 10:01:59 +0100
Subject: [PATCH 01/39] added comment
---
CMakeLists.txt | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 524f0ff9c..08fba109b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,6 +7,9 @@ project (
LANGUAGES CXX
)
+# as long as there still are modules using it:
+enable_language (Fortran)
+
# ignore many irrelevant Up-to-date messages during install
set (CMAKE_INSTALL_MESSAGE LAZY)
@@ -21,7 +24,6 @@ set (CMAKE_CXX_STANDARD 17)
enable_testing ()
set (CTEST_OUTPUT_ON_FAILURE 1)
-ENABLE_LANGUAGE (Fortran)
# unit testing coverage, does not work yet
#include (CodeCoverage)
--
GitLab
From 97c66f82981890f279858eb6132eda1c5a456f43 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 10:02:16 +0100
Subject: [PATCH 02/39] replaced compiler warning with TODO
---
Stack/SuperStupidStack/SuperStupidStack.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 2b03c5443..a3b98250c 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -121,7 +121,7 @@ namespace corsika::stack {
void IncrementSize() {
fDataPID.push_back(Code::Unknown);
fDataE.push_back(0 * joule);
-#warning this here makes no sense: see issue #48
+ //#TODO this here makes no sense: see issue #48
auto const dummyCS = corsika::geometry::CoordinateSystem::CreateRootCS();
fMomentum.push_back(MomentumVector(dummyCS,
{0 * newton_second, 0 * newton_second, 0 * newton_second}));
--
GitLab
From f9544013af6fa1d136cdf2639c73fd00dc626b7b Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 10:02:30 +0100
Subject: [PATCH 03/39] removed debug printout
---
Processes/Sibyll/code_generator.py | 2 --
1 file changed, 2 deletions(-)
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index 77c1f5dd8..a593a1f2c 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -159,8 +159,6 @@ if __name__ == "__main__":
pythia_db = load_pythiadb(sys.argv[1])
read_sibyll_codes(sys.argv[2], pythia_db)
- print (str(pythia_db))
-
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
print(generate_sibyll_enum(pythia_db), file=f)
--
GitLab
From 0ef0a7259ef1444b72021bc6e228f0a2ea8eb6e5 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 11:28:16 +0000
Subject: [PATCH 04/39] added momentum to sibyll stack, implemented boost
between sibyll stack and corsika stack
---
Documentation/Examples/cascade_example.cc | 220 ++++++++++++++++------
Framework/Cascade/SibStack.h | 27 ++-
2 files changed, 185 insertions(+), 62 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 69b7e92f0..7b4ab0ad4 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -23,6 +23,7 @@
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
+
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
@@ -42,9 +43,13 @@ public:
template <typename Particle>
double MinStepLength(Particle& p) const {
+ const Code corsikaBeamId = p.GetPID();
+
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
- int kBeam = 1;
+ int kBeam = process::sibyll::GetSibyllXSCode( corsikaBeamId );
+
+ bool kInteraction = process::sibyll::CanInteract( corsikaBeamId );
/*
the target should be defined by the Environment,
@@ -53,30 +58,44 @@ public:
*/
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
- int kTarget = 1;
- double beamEnergy = p.GetEnergy() / 1_GeV;
- std::cout << "ProcessSplit: " << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
- double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
- double dumdif[3];
-
- if(kTarget==1)
- sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
- else
- sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+ int kTarget = 16;
+ double beamEnergy = p.GetEnergy() / 1_GeV;
+#warning boost to cm. still missing, sibyll cross section input is cm. energy!
+
+ std::cout << "ProcessSplit: " << "MinStep: input en: " << beamEnergy
+ << " beam can interact:" << kBeam
+ << " beam XS code:" << kBeam
+ << " beam pid:" << p.GetPID()
+ << " target mass number:" << kTarget << std::endl;
+
+ double next_step;
+ if(kInteraction){
+
+ double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
+ double dumdif[3];
+
+ if(kTarget==1)
+ sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+ else
+ sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+
+ std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ CrossSectionType sig = prodCrossSection * 1_mbarn;
+ std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+ const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
+ std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl;
+ // calculate interaction length in medium
+ double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
+ // pick random step lenth
+ std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
+ // add exponential sampling
+ int a = 0;
+ next_step = -int_length * log(s_rndm_(a));
+ }else
+#warning define infinite interaction length? then we can skip the test in DoDiscrete()
+ next_step = 1.e8;
- std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
-
- const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
- std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl;
- // calculate interaction length in medium
- double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
- // pick random step lenth
- std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
- // add exponential sampling
- int a = 0;
- const double next_step = -int_length * log(s_rndm_(a));
/*
what are the units of the output? slant depth or 3space length?
@@ -95,28 +114,75 @@ public:
void DoDiscrete(Particle& p, Stack& s) const {
cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
if( process::sibyll::CanInteract( p.GetPID() ) ){
-
- // 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)
- */
- EnergyType E = p.GetEnergy();
- EnergyType Ecm = sqrt( 2. * E * 0.93827_GeV );
+ cout << "defining coordinates" << endl;
+ // coordinate system, get global frame of reference
+ CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+
+ 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; //p.GetPID();
- int kBeam = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+ // proton mass in units of energy
+ const EnergyType proton_mass_en = 0.93827_GeV ; //0.93827_GeV / si::constants::cSquared ;
+
+ cout << "defining target momentum.." << endl;
+ // FOR NOW: target is always at rest
+ const EnergyType Etarget = 0. * 1_GeV + proton_mass_en;
+ const auto pTarget = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+ cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
+ // const auto pBeam = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+ // cout << "beam momentum: " << pBeam.GetComponents() << endl;
+ cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV * si::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)
+ */
+ // cout << "defining total energy" << endl;
+ // total energy: E_beam + E_target
+ // in lab. frame: E_beam + m_target*c**2
+ EnergyType E = p.GetEnergy();
+ EnergyType Etot = E + Etarget;
+ // cout << "tot. energy: " << Etot / 1_GeV << endl;
+ // cout << "defining total momentum" << endl;
+ // total momentum
+ super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+ // cout << "tot. momentum: " << Ptot.GetComponents() / 1_GeV * si::constants::c << endl;
+ // cout << "inv. mass.." << endl;
+ // invariant mass, i.e. cm. energy
+ EnergyType Ecm = sqrt( Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared ); //sqrt( 2. * E * 0.93827_GeV );
+ // cout << "inv. mass: " << Ecm / 1_GeV << endl;
+ // cout << "boost parameters.." << endl;
+ /*
+ 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 = ( E + proton_mass * si::constants::cSquared ) / Ecm ;
+ // const double gambet = sqrt( E * E - proton_mass * proton_mass ) / Ecm;
+ const double gamma = Etot / Ecm ;
+ const auto gambet = Ptot / (Ecm / si::constants::c );
+
+ std::cout << "ProcessSplit: " << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "ProcessSplit: " << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- // FOR NOW: set target to proton
- int kTarget = 1; //p.GetPID();
- std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
+ std::cout << "ProcessSplit: " << " DoDiscrete: 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;
p.Delete();
@@ -138,27 +204,53 @@ public:
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
- /*
- 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
- in general transformation is rotation + boost
- */
- const EnergyType proton_mass = 0.93827_GeV;
- const double gamma = ( E + proton_mass ) / ( Ecm );
- const double gambet = sqrt( E * E - proton_mass * proton_mass ) / Ecm;
// SibStack does not know about momentum yet so we need counter to access momentum array in Sibyll
int i = -1;
- for (auto &p: ss){
+ super_stupid::MomentumVector Ptot_final(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+ for (auto &psib: ss){
++i;
- //transform to lab. frame, primitve
- const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
+ //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] * si::constants::c ) / ( gamma + 1.);
+ pnorm += psib.GetEnergy();
+
+ for(int j=0; j<3; ++j){
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
+ // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
+ // << " gb: " << gammaBetaComponents[j] << endl;
+ en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+ }
+ // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
+ // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
+ // << " en lab (GeV): " << en_lab / 1_GeV << endl;
+ // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
+ // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
+
// add to corsika stack
auto pnew = s.NewParticle();
- pnew.SetEnergy( en_lab * 1_GeV );
- pnew.SetPID( process::sibyll::ConvertFromSibyll( p.GetPID() ) );
+ pnew.SetEnergy( en_lab );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+ //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+
+ 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) );
+ //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+ //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ Ptot_final += pnew.GetMomentum();
}
+ //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
}
}else
p.Delete();
@@ -226,6 +318,12 @@ double s_rndm_(int &)
int main(){
+ // coordinate system, get global frame of reference
+ CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+
+ QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+ Point pOrig(rootCS, coordinates);
+
stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
@@ -236,8 +334,14 @@ int main(){
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV;
+ EnergyType E0 = 500_GeV;
+ MomentumType P0 = sqrt( E0*E0 - 0.93827_GeV * 0.93827_GeV ) / si::constants::c;
+ auto plab = super_stupid::MomentumVector(rootCS,
+ 0. * 1_GeV / si::constants::c,
+ 0. * 1_GeV / si::constants::c,
+ P0);
particle.SetEnergy(E0);
+ particle.SetMomentum(plab);
particle.SetPID( Code::Proton );
EAS.Init();
EAS.Run();
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index d38bf0da9..c10ea7d8d 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -10,6 +10,8 @@
using namespace std;
using namespace corsika::stack;
+using namespace corsika::units;
+using namespace corsika::geometry;
class SibStackData {
@@ -19,14 +21,30 @@ class SibStackData {
void Clear() { s_plist_.np = 0; }
int GetSize() const { return s_plist_.np; }
+#warning check actual capacity of sibyll stack
int GetCapacity() const { return 8000; }
void SetId(const int i, const int v) { s_plist_.llist[i]=v; }
- void SetEnergy(const int i, const double v) { s_plist_.p[3][i]=v; }
-
+ 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)
+ {
+ auto tmp = v.GetComponents();
+ for(int idx=0; idx<3; ++idx)
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;
+ }
+
int GetId(const int i) const { return s_plist_.llist[i]; }
- double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
+
+ EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+
+ super_stupid::MomentumVector GetMomentum(const int i) const
+ {
+ CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
+ super_stupid::MomentumVector v1(rootCS,components);
+ return v1;
+ }
void Copy(const int i1, const int i2) {
s_plist_.llist[i1] = s_plist_.llist[i2];
@@ -45,9 +63,10 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
- double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ 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()); }
};
--
GitLab
From c4ec6492e13281f105bbe8708bed2b4dedee2daf Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 10:01:59 +0100
Subject: [PATCH 05/39] added comment
---
CMakeLists.txt | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2a5cae812..1ac8449f0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,6 +7,7 @@ project (
LANGUAGES CXX
)
+# as long as there still are modules using it:
enable_language (Fortran)
# ignore many irrelevant Up-to-date messages during install
@@ -38,7 +39,6 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
-
# unit testing coverage, does not work yet
#include (CodeCoverage)
##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')
--
GitLab
From 1ea5333d573ed6f627921aef48519f6313cd6dfb Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 10:02:30 +0100
Subject: [PATCH 06/39] removed debug printout
---
Processes/Sibyll/code_generator.py | 2 --
1 file changed, 2 deletions(-)
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index 77c1f5dd8..a593a1f2c 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -159,8 +159,6 @@ if __name__ == "__main__":
pythia_db = load_pythiadb(sys.argv[1])
read_sibyll_codes(sys.argv[2], pythia_db)
- print (str(pythia_db))
-
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
print(generate_sibyll_enum(pythia_db), file=f)
--
GitLab
From b69eef7f6a33d2d10ef6445f9d4591e5d6075298 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 11:28:16 +0000
Subject: [PATCH 07/39] added momentum to sibyll stack, implemented boost
between sibyll stack and corsika stack
---
Documentation/Examples/cascade_example.cc | 301 ++++++++++++++--------
Framework/Cascade/SibStack.h | 42 ++-
2 files changed, 228 insertions(+), 115 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 34b55476d..18c8c65bc 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -24,6 +24,7 @@
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
+
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
@@ -43,46 +44,59 @@ public:
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
+ const Code corsikaBeamId = p.GetPID();
+
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
- int kBeam = 1;
+ 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 = 1;
- double beamEnergy = p.GetEnergy() / 1_GeV;
- std::cout << "ProcessSplit: "
- << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
- double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
-
- if (kTarget == 1)
- sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
- else
- sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy);
-
- std::cout << "ProcessSplit: "
- << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "ProcessSplit: "
- << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
-
- const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
- std::cout << "ProcessSplit: "
- << "nucleon mass " << nucleon_mass << std::endl;
- // calculate interaction length in medium
- double int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
- // pick random step lenth
- std::cout << "ProcessSplit: "
- << "interaction length (g/cm2): " << int_length << std::endl;
- // add exponential sampling
- int a = 0;
- const double next_step = -int_length * log(s_rndm_(a));
+ int kTarget = 16;
+ double beamEnergy = p.GetEnergy() / 1_GeV;
+#warning boost to cm. still missing, sibyll cross section input is cm. energy!
+
+ std::cout << "ProcessSplit: " << "MinStep: input en: " << beamEnergy
+ << " beam can interact:" << kBeam
+ << " beam XS code:" << kBeam
+ << " beam pid:" << p.GetPID()
+ << " target mass number:" << kTarget << std::endl;
+
+ double next_step;
+ if(kInteraction){
+
+ double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
+ double dumdif[3];
+
+ if(kTarget==1)
+ sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+ else
+ sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+
+ std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ CrossSectionType sig = prodCrossSection * 1_mbarn;
+ std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+ const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
+ std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl;
+ // calculate interaction length in medium
+ double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
+ // pick random step lenth
+ std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
+ // add exponential sampling
+ int a = 0;
+ next_step = -int_length * log(s_rndm_(a));
+ }else
+#warning define infinite interaction length? then we can skip the test in DoDiscrete()
+ next_step = 1.e8;
+
/*
what are the units of the output? slant depth or 3space length?
@@ -100,79 +114,147 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- cout << "DoDiscrete: " << p.GetPID() << " interaction? "
- << process::sibyll::CanInteract(p.GetPID()) << endl;
- if (process::sibyll::CanInteract(p.GetPID())) {
-
+ cout << "DoDiscrete: " << 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 = CoordinateSystem::CreateRootCS();
+
+ 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; //p.GetPID();
+
+ // proton mass in units of energy
+ const EnergyType proton_mass_en = 0.93827_GeV ; //0.93827_GeV / si::constants::cSquared ;
+
+ cout << "defining target momentum.." << endl;
+ // FOR NOW: target is always at rest
+ const EnergyType Etarget = 0. * 1_GeV + proton_mass_en;
+ const auto pTarget = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+ cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
+ // const auto pBeam = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+ // cout << "beam momentum: " << pBeam.GetComponents() << endl;
+ cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV * si::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)
+ 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)
*/
- EnergyType E = p.GetEnergy();
- EnergyType Ecm = sqrt(2. * E * 0.93827_GeV);
-
- int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // cout << "defining total energy" << endl;
+ // total energy: E_beam + E_target
+ // in lab. frame: E_beam + m_target*c**2
+ EnergyType E = p.GetEnergy();
+ EnergyType Etot = E + Etarget;
+ // cout << "tot. energy: " << Etot / 1_GeV << endl;
+ // cout << "defining total momentum" << endl;
+ // total momentum
+ super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+ // cout << "tot. momentum: " << Ptot.GetComponents() / 1_GeV * si::constants::c << endl;
+ // cout << "inv. mass.." << endl;
+ // invariant mass, i.e. cm. energy
+ EnergyType Ecm = sqrt( Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared ); //sqrt( 2. * E * 0.93827_GeV );
+ // cout << "inv. mass: " << Ecm / 1_GeV << endl;
+ // cout << "boost parameters.." << endl;
+ /*
+ 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 = ( E + proton_mass * si::constants::cSquared ) / Ecm ;
+ // const double gambet = sqrt( E * E - proton_mass * proton_mass ) / Ecm;
+ const double gamma = Etot / Ecm ;
+ const auto gambet = Ptot / (Ecm / si::constants::c );
+
+ std::cout << "ProcessSplit: " << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "ProcessSplit: " << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+
+
+ std::cout << "ProcessSplit: " << " DoDiscrete: 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;
+ p.Delete();
+ fCount++;
+ } else {
+ // Sibyll does not know about units..
+ double sqs = Ecm / 1_GeV;
+ // running sibyll, filling stack
+ sibyll_( kBeam, kTarget, sqs);
+ // running decays
+ //decsib_();
+ // print final state
+ int print_unit = 6;
+ sib_list_( print_unit );
+
+ // delete current particle
+ p.Delete();
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- // FOR NOW: set target to proton
- int kTarget = 1; // p.GetPID();
-
- std::cout << "ProcessSplit: "
- << " DoDiscrete: 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;
- p.Delete();
- fCount++;
- } else {
- // Sibyll does not know about units..
- double sqs = Ecm / 1_GeV;
- // running sibyll, filling stack
- sibyll_(kBeam, kTarget, sqs);
- // running decays
- // 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;
- /*
- 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
- in general transformation is rotation + boost
- */
- const EnergyType proton_mass = 0.93827_GeV;
- const double gamma = (E + proton_mass) / (Ecm);
- const double gambet = sqrt(E * E - proton_mass * proton_mass) / Ecm;
-
- // SibStack does not know about momentum yet so we need counter to access momentum
- // array in Sibyll
- int i = -1;
- for (auto& p : ss) {
- ++i;
- // transform to lab. frame, primitve
- const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetEnergy(en_lab * 1_GeV);
- pnew.SetPID(process::sibyll::ConvertFromSibyll(p.GetPID()));
- }
+ // 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;
+ //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] * si::constants::c ) / ( gamma + 1.);
+ pnorm += psib.GetEnergy();
+
+ for(int j=0; j<3; ++j){
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
+ // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
+ // << " gb: " << gammaBetaComponents[j] << endl;
+ en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+ }
+ // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
+ // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
+ // << " en lab (GeV): " << en_lab / 1_GeV << endl;
+ // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
+ // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
+
+ // add to corsika stack
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy( en_lab );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+ //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+
+ 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) );
+ //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+ //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ Ptot_final += pnew.GetMomentum();
}
- } else
+ //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
+ }
+ }else
p.Delete();
}
@@ -238,8 +320,14 @@ double s_rndm_(int&) {
int main() {
- tracking_line::TrackingLine<setup::Stack> tracking;
- stack_inspector::StackInspector<setup::Stack> p0(true);
+ // coordinate system, get global frame of reference
+ CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+
+ QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+ Point pOrig(rootCS, coordinates);
+
+ stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+
ProcessSplit p1;
const auto sequence = p0 + p1;
setup::Stack stack;
@@ -248,9 +336,16 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV;
+ EnergyType E0 = 500_GeV;
+ MomentumType P0 = sqrt( E0*E0 - 0.93827_GeV * 0.93827_GeV ) / si::constants::c;
+ auto plab = super_stupid::MomentumVector(rootCS,
+ 0. * 1_GeV / si::constants::c,
+ 0. * 1_GeV / si::constants::c,
+ P0);
particle.SetEnergy(E0);
- particle.SetPID(Code::Proton);
+ particle.SetMomentum(plab);
+ particle.SetPID( Code::Proton );
+
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 7e9bd3bac..c18cb3e33 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -10,6 +10,8 @@
using namespace std;
using namespace corsika::stack;
+using namespace corsika::units;
+using namespace corsika::geometry;
class SibStackData {
@@ -17,16 +19,33 @@ public:
void Init();
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
int GetCapacity() const { return 8000; }
- void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
- void SetEnergy(const int i, const double v) { s_plist_.p[3][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 SetMomentum(const int i, const super_stupid::MomentumVector& v)
+ {
+ auto tmp = v.GetComponents();
+ for(int idx=0; idx<3; ++idx)
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;
+ }
+
int GetId(const int i) const { return s_plist_.llist[i]; }
- double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
-
+
+ EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+
+ super_stupid::MomentumVector GetMomentum(const int i) const
+ {
+ CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
+ super_stupid::MomentumVector v1(rootCS,components);
+ return v1;
+ }
+
void Copy(const int i1, const int i2) {
s_plist_.llist[i1] = s_plist_.llist[i2];
s_plist_.p[3][i1] = s_plist_.p[3][i2];
@@ -46,12 +65,11 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
public:
void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
- double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ 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()));
- }
+ 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()); }
+
};
typedef Stack<SibStackData, ParticleInterface> SibStack;
--
GitLab
From f8a177e2de31394585735a36af99d869d2f267af Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 12:36:03 +0000
Subject: [PATCH 08/39] added cm energy calculation for minsteplength in sibyll
process
---
Documentation/Examples/cascade_example.cc | 33 ++++++++++++++++-------
1 file changed, 24 insertions(+), 9 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 18c8c65bc..5cdcf73e8 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -44,6 +44,9 @@ public:
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
const Code corsikaBeamId = p.GetPID();
// beam particles for sibyll : 1, 2, 3 for p, pi, k
@@ -60,13 +63,24 @@ public:
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
int kTarget = 16;
- double beamEnergy = p.GetEnergy() / 1_GeV;
-#warning boost to cm. still missing, sibyll cross section input is cm. energy!
+
+ // proton mass in units of energy
+ const EnergyType proton_mass_en = 0.93827_GeV ;
+
+ EnergyType Etot = p.GetEnergy() + kTarget * proton_mass_en;
+ 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() * si::constants::cSquared );
+ double Ecm = sqs / 1_GeV;
- std::cout << "ProcessSplit: " << "MinStep: input en: " << beamEnergy
- << " beam can interact:" << kBeam
- << " beam XS code:" << kBeam
- << " beam pid:" << p.GetPID()
+ std::cout << "ProcessSplit: " << "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 next_step;
@@ -76,9 +90,9 @@ public:
double dumdif[3];
if(kTarget==1)
- sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+ sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
else
- sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+ sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy );
std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
CrossSectionType sig = prodCrossSection * 1_mbarn;
@@ -291,6 +305,7 @@ public:
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);
@@ -336,7 +351,7 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 500_GeV;
+ EnergyType E0 = 100_GeV;
MomentumType P0 = sqrt( E0*E0 - 0.93827_GeV * 0.93827_GeV ) / si::constants::c;
auto plab = super_stupid::MomentumVector(rootCS,
0. * 1_GeV / si::constants::c,
--
GitLab
From 7e1ca43a3204c39c2832dc94baea798610600e4a Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 5 Dec 2018 08:39:21 +0000
Subject: [PATCH 09/39] fixed rebase
---
Documentation/Examples/cascade_example.cc | 11 ++++++-----
Framework/Cascade/SibStack.h | 2 +-
2 files changed, 7 insertions(+), 6 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 5cdcf73e8..5697e0d57 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -121,7 +121,7 @@ public:
}
template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
// corsika::utls::ignore(p);
return EProcessReturn::eOk;
}
@@ -132,7 +132,7 @@ public:
if( process::sibyll::CanInteract( p.GetPID() ) ){
cout << "defining coordinates" << endl;
// coordinate system, get global frame of reference
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point pOrig(rootCS, coordinates);
@@ -336,12 +336,13 @@ double s_rndm_(int&) {
int main() {
// coordinate system, get global frame of reference
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point pOrig(rootCS, coordinates);
-
- stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
+
+ tracking_line::TrackingLine<setup::Stack> tracking;
+ stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index c18cb3e33..57a3f9c16 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -40,7 +40,7 @@ public:
super_stupid::MomentumVector GetMomentum(const int i) const
{
- CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
super_stupid::MomentumVector v1(rootCS,components);
return v1;
--
GitLab
From d03f06af2443013efa1c4003af92bd9c60266359 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 08:20:41 +0000
Subject: [PATCH 10/39] added proto decay process
---
Documentation/Examples/cascade_example.cc | 19 +++++++++++++++++++
1 file changed, 19 insertions(+)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 5697e0d57..f1663654a 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -333,6 +333,25 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
+class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+public:
+ ProcessDecay() {}
+ void Init() {}
+ template <typename Particle>
+ double MinStepLength(Particle& p) const {
+ }
+ template <typename Particle, typename Stack>
+ void DoDiscrete(Particle& p, Stack& s) const {
+ }
+
+ template <typename Particle, typename Trajectory, typename Stack>
+ EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+};
+
+
int main() {
// coordinate system, get global frame of reference
--
GitLab
From a0286a682accba00e8d66906da5bda5d3d568c78 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 17:43:21 +0000
Subject: [PATCH 11/39] added mass density type
---
Framework/Units/PhysicalUnits.h | 9 +++++++++
1 file changed, 9 insertions(+)
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index e97003eee..b51634796 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -36,6 +36,15 @@ 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 MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
+
+ // defining momentum you suckers
+ // dimensions, i.e. composition in base SI dimensions
+ using momentum_d = phys::units::dimensions< 1, 1, -1 >;
+ // defining the unit of momentum, so far newton-meter, maybe go to HEP?
+ constexpr phys::units::quantity< momentum_d > newton_second { meter * kilogram / second };
+ // defining the type
+
using MomentumType = phys::units::quantity<momentum_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
--
GitLab
From fa0f9b7a6ee111c5e7a230c3886fc489b929ea29 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 17:44:17 +0000
Subject: [PATCH 12/39] added primitive decay step length
---
Documentation/Examples/cascade_example.cc | 37 +++++++++++++++++++++--
1 file changed, 34 insertions(+), 3 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f1663654a..3d5625232 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -337,12 +337,42 @@ class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
public:
ProcessDecay() {}
void Init() {}
+
template <typename Particle>
double MinStepLength(Particle& p) const {
+ EnergyType E = p.GetEnergy();
+ MassType m = corsika::particles::GetMass(p.GetPID());
+ // env.GetDensity();
+ const MassDensityType density = 1.e3 * kilogram / ( 1_cm * 1_cm * 1_cm );
+
+ const double gamma = E / m / constants::cSquared;
+ // lifetimes not implemented yet
+ TimeType t0;
+ switch( p.GetPID() ){
+ case Code::PiPlus :
+ t0 = 1.e-5 * 1_s;
+ break;
+
+ case Code::KPlus :
+ t0 = 1.e-5 * 1_s;
+ break;
+
+ default:
+ t0 = 1.e8 * 1_s;
+ break;
+ }
+ 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>
+
+ template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- }
+ }
template <typename Particle, typename Trajectory, typename Stack>
EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
@@ -364,7 +394,8 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
- const auto sequence = p0 + p1;
+ ProcessDecay p2;
+ const auto sequence = p0 + p1 + p2;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
--
GitLab
From 4366cb4926aa12b7fbc96b631f026c9597a09bb7 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 17:50:29 +0000
Subject: [PATCH 13/39] added more realistic numbers for decay
---
Documentation/Examples/cascade_example.cc | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 3d5625232..cec989d69 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -343,14 +343,14 @@ public:
EnergyType E = p.GetEnergy();
MassType m = corsika::particles::GetMass(p.GetPID());
// env.GetDensity();
- const MassDensityType density = 1.e3 * 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 / constants::cSquared;
// lifetimes not implemented yet
TimeType t0;
switch( p.GetPID() ){
case Code::PiPlus :
- t0 = 1.e-5 * 1_s;
+ t0 = 2.6e-8 * 1_s;
break;
case Code::KPlus :
--
GitLab
From ece530f9f0a5e35ff066b931a3e910cd57e93aad Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sun, 2 Dec 2018 11:05:34 +0000
Subject: [PATCH 14/39] fixed rebase
---
Framework/Units/PhysicalUnits.h | 9 ---------
1 file changed, 9 deletions(-)
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index b51634796..54b7d3371 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -37,15 +37,6 @@ namespace corsika::units::si {
using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
-
- // defining momentum you suckers
- // dimensions, i.e. composition in base SI dimensions
- using momentum_d = phys::units::dimensions< 1, 1, -1 >;
- // defining the unit of momentum, so far newton-meter, maybe go to HEP?
- constexpr phys::units::quantity< momentum_d > newton_second { meter * kilogram / second };
- // defining the type
-
- using MomentumType = phys::units::quantity<momentum_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
} // end namespace corsika::units::si
--
GitLab
From 94e892c8b2f23d5853514262d3bc400debbc1a20 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 12:56:08 +0000
Subject: [PATCH 15/39] fixed rebase
---
Documentation/Examples/cascade_example.cc | 6 +++---
Framework/Units/PhysicalUnits.h | 3 ++-
2 files changed, 5 insertions(+), 4 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index cec989d69..da916db78 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -339,7 +339,7 @@ public:
void Init() {}
template <typename Particle>
- double MinStepLength(Particle& p) const {
+ double MinStepLength(Particle& p, setup::Trajectory&) const {
EnergyType E = p.GetEnergy();
MassType m = corsika::particles::GetMass(p.GetPID());
// env.GetDensity();
@@ -374,8 +374,8 @@ public:
void DoDiscrete(Particle& p, Stack& s) const {
}
- template <typename Particle, typename Trajectory, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
return EProcessReturn::eOk;
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 54b7d3371..dfddec627 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -38,7 +38,8 @@ namespace corsika::units::si {
using MassType = phys::units::quantity<phys::units::mass_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
/**
--
GitLab
From 4043ac5b635fb56620d5b72801055b9ef48c6ee2 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 13:13:24 +0000
Subject: [PATCH 16/39] moved decay into sibyll process
---
Documentation/Examples/cascade_example.cc | 52 ++----------------
Processes/Sibyll/CMakeLists.txt | 1 +
Processes/Sibyll/ProcessDecay.h | 65 +++++++++++++++++++++++
3 files changed, 69 insertions(+), 49 deletions(-)
create mode 100644 Processes/Sibyll/ProcessDecay.h
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index da916db78..7bab3ecda 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -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;
@@ -333,54 +335,6 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
-class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
-public:
- ProcessDecay() {}
- void Init() {}
-
- 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;
- // lifetimes not implemented yet
- TimeType t0;
- switch( p.GetPID() ){
- case Code::PiPlus :
- t0 = 2.6e-8 * 1_s;
- break;
-
- case Code::KPlus :
- t0 = 1.e-5 * 1_s;
- break;
-
- default:
- t0 = 1.e8 * 1_s;
- break;
- }
- 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 {
- }
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
- }
-
-};
-
int main() {
@@ -394,7 +348,7 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
- ProcessDecay p2;
+ corsika::process::sibyll::ProcessDecay p2;
const auto sequence = p0 + p1 + p2;
setup::Stack stack;
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index f2e4a9fee..a0985a1ad 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -20,6 +20,7 @@ set (
set (
MODEL_HEADERS
ParticleConversion.h
+ ProcessDecay.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
new file mode 100644
index 000000000..20f493607
--- /dev/null
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -0,0 +1,65 @@
+#ifndef _include_ProcessDecay_h_
+#define _include_ProcessDecay_h_
+
+#include <corsika/process/ContinuousProcess.h>
+
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+
+//using namespace corsika::particles;
+
+namespace corsika::process {
+
+ namespace sibyll {
+
+class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+public:
+ ProcessDecay() {}
+ void Init() {}
+
+ 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;
+ // lifetimes not implemented yet
+ TimeType t0;
+ switch( p.GetPID() ){
+ case corsika::particles::Code::PiPlus :
+ t0 = 2.6e-8 * 1_s;
+ break;
+
+ case corsika::particles::Code::KPlus :
+ t0 = 1.e-5 * 1_s;
+ break;
+
+ default:
+ t0 = 1.e8 * 1_s;
+ break;
+ }
+ 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 {
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+};
+ }
+}
+
+#endif
--
GitLab
From 309963036df7a287119a9b0163697f787909f0a5 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 5 Dec 2018 08:17:39 +0000
Subject: [PATCH 17/39] added negative charge pions and kaons
---
Processes/Sibyll/ProcessDecay.h | 9 +++++++++
1 file changed, 9 insertions(+)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 20f493607..4989ab90a 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -32,9 +32,17 @@ public:
t0 = 2.6e-8 * 1_s;
break;
+ case corsika::particles::Code::PiMinus :
+ t0 = 2.6e-8 * 1_s;
+ break;
+
case corsika::particles::Code::KPlus :
t0 = 1.e-5 * 1_s;
break;
+
+ case corsika::particles::Code::KMinus :
+ t0 = 1.e-5 * 1_s;
+ break;
default:
t0 = 1.e8 * 1_s;
@@ -51,6 +59,7 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
+
}
template <typename Particle, typename Stack>
--
GitLab
From dd6c3424768ed67e50fd595c11c9113ba0185591 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 5 Dec 2018 10:36:39 +0000
Subject: [PATCH 18/39] added decay configuration routines, turned on decay for
sibyll interaction
---
Documentation/Examples/cascade_example.cc | 58 +++++----
Processes/Sibyll/ProcessDecay.h | 145 ++++++++++++++--------
2 files changed, 131 insertions(+), 72 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 7bab3ecda..e9a4a5ed3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -43,6 +43,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 {
@@ -211,7 +240,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 );
@@ -228,6 +258,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
@@ -300,29 +333,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:
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 4989ab90a..ecae303e7 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -12,62 +12,107 @@ namespace corsika::process {
namespace sibyll {
-class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
-public:
- ProcessDecay() {}
- void Init() {}
-
- 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 );
+
+ 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();
+ }
+
+ }
+
- const double gamma = E / m / constants::cSquared;
- // lifetimes not implemented yet
- TimeType t0;
- switch( p.GetPID() ){
- case corsika::particles::Code::PiPlus :
- t0 = 2.6e-8 * 1_s;
- break;
-
- case corsika::particles::Code::PiMinus :
- t0 = 2.6e-8 * 1_s;
- break;
+ class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+ public:
+ ProcessDecay() {}
+ void Init() {
+ //setHadronsUnstable();
+ }
- case corsika::particles::Code::KPlus :
- t0 = 1.e-5 * 1_s;
- break;
+ 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;
+ }
+ }
- case corsika::particles::Code::KMinus :
- t0 = 1.e-5 * 1_s;
- break;
+ friend void setHadronsUnstable();
- default:
- t0 = 1.e8 * 1_s;
- break;
- }
- 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 {
+ 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 );
- }
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
- }
+ 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 {
+
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+ };
}
}
--
GitLab
From ba25b159ee3e3e45a2decf8a643531a7ac9040d4 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 6 Dec 2018 10:59:14 +0000
Subject: [PATCH 19/39] added actual decay to ProcessDecay
---
Framework/Cascade/SibStack.h | 4 +++-
Processes/Sibyll/ProcessDecay.h | 34 ++++++++++++++++++++++++++++++++-
2 files changed, 36 insertions(+), 2 deletions(-)
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 57a3f9c16..2275cdf49 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -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); }
};
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index ecae303e7..ab11fc7e4 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -5,6 +5,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/cascade/SibStack.h>
//using namespace corsika::particles;
@@ -104,7 +105,38 @@ namespace corsika::process {
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>
--
GitLab
From e1bfd98d100b085d016165b2418665144c28db1d Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 6 Dec 2018 11:01:35 +0000
Subject: [PATCH 20/39] suppressed output
---
Processes/Sibyll/ProcessDecay.h | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index ab11fc7e4..449a281b4 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -118,8 +118,8 @@ namespace corsika::process {
std::cout << "calling Sibyll decay routine.." << std::endl;
decsib_();
// print output
- int print_unit = 6;
- sib_list_( print_unit );
+ //int print_unit = 6;
+ //sib_list_( print_unit );
// copy particles from sibyll stack to corsika
int i = -1;
for (auto &psib: ss){
@@ -127,7 +127,7 @@ namespace corsika::process {
// 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;
+ //cout << "decay product: " << process::sibyll::ConvertFromSibyll( psib.GetPID() ) << endl;
auto pnew = s.NewParticle();
pnew.SetEnergy( psib.GetEnergy() );
pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
--
GitLab
From a90b32d15805cf5392e9c0746eb75a1fc9b223c3 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sat, 8 Dec 2018 13:56:34 +0000
Subject: [PATCH 21/39] added particle cut process
---
Documentation/Examples/cascade_example.cc | 336 +++++++++++++++++-----
Framework/Cascade/Cascade.h | 2 +
Processes/Sibyll/ProcessDecay.h | 15 +-
3 files changed, 271 insertions(+), 82 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 367b7700c..7bd286b11 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -39,6 +39,183 @@ using namespace corsika::setup;
using namespace std;
static int fCount = 0;
+static EnergyType fEnergy = 0. * 1_GeV;
+
+// FOR NOW: global static variables for ParticleCut process
+// this is just wrong...
+static EnergyType fEmEnergy;
+static int fEmCount;
+
+static EnergyType fInvEnergy;
+static int fInvCount;
+
+
+class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
+public:
+ ProcessEMCut() {}
+ template <typename Particle>
+ bool isBelowEnergyCut( Particle &p ) const
+ {
+ // FOR NOW: center-of-mass energy hard coded
+ const EnergyType Ecm = sqrt( 2. * p.GetEnergy() * 0.93827_GeV );
+ if( p.GetEnergy() < 50_GeV || Ecm < 10_GeV )
+ return true;
+ else
+ return false;
+ }
+
+ bool isEmParticle( Code pCode ) const
+ {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch( pCode ){
+ case Code::Electron :
+ is_em = true;
+ break;
+ case Code::Gamma :
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ void defineEmParticles() const
+ {
+ // create bool array identifying em particles
+ }
+
+ bool isInvisible( Code pCode ) const
+ {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch( pCode ){
+ case Code::NuE :
+ is_inv = true;
+ break;
+ case Code::NuEBar :
+ is_inv = true;
+ break;
+ case Code::NuMu :
+ is_inv = true;
+ break;
+ case Code::NuMuBar :
+ is_inv = true;
+ break;
+ case Code::MuPlus :
+ is_inv = true;
+ break;
+ case Code::MuMinus :
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+
+ template <typename Particle>
+ double MinStepLength(Particle& p, setup::Trajectory&) const
+ {
+ const Code pid = p.GetPID();
+ if( isEmParticle( pid ) || isInvisible( pid ) ){
+ cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+ return 0.;
+ } else {
+ double next_step = std::numeric_limits<double>::infinity();
+ cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+ return next_step;
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&p, setup::Trajectory&t, Stack&s) const {
+ // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
+ // 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;
+ }
+
+ template <typename Particle, typename Stack>
+ void DoDiscrete(Particle& p, Stack& s) const
+ {
+ cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
+ const Code pid = p.GetPID();
+ if( isEmParticle( pid ) ){
+ cout << "removing em. particle..." << endl;
+ fEmEnergy += p.GetEnergy();
+ fEmCount += 1;
+ p.Delete();
+ }
+ else if ( isInvisible( pid ) ){
+ cout << "removing inv. particle..." << endl;
+ fInvEnergy += p.GetEnergy();
+ fInvCount += 1;
+ p.Delete();
+ }
+ else if( isBelowEnergyCut( p ) ){
+ cout << "removing low en. particle..." << endl;
+ fEnergy += p.GetEnergy();
+ fCount += 1;
+ p.Delete();
+ }
+ }
+
+ void Init()
+ {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ //defineEmParticles();
+ }
+
+ void ShowResults(){
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " ******************************" << endl;
+ }
+
+ EnergyType GetInvEnergy()
+ {
+ return fInvEnergy;
+ }
+
+ EnergyType GetCutEnergy()
+ {
+ return fEnergy;
+ }
+
+ EnergyType GetEmEnergy()
+ {
+ return fEmEnergy;
+ }
+
+private:
+};
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
@@ -141,14 +318,14 @@ public:
next_step = -int_length * log(s_rndm_(a));
}else
#warning define infinite interaction length? then we can skip the test in DoDiscrete()
- next_step = 1.e8;
+ next_step = std::numeric_limits<double>::infinity();
/*
what are the units of the output? slant depth or 3space length?
*/
std::cout << "ProcessSplit: "
- << "next step (g/cm2): " << next_step << std::endl;
+ << "next interaction (g/cm2): " << next_step << std::endl;
return next_step;
}
@@ -160,7 +337,7 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
+ cout << "ProcessSplit: " << "DoDiscrete: " << 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
@@ -231,81 +408,82 @@ public:
std::cout << "ProcessSplit: " << " DoDiscrete: 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;
- p.Delete();
- fCount++;
- } 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_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;
+ if (E < 8.5_GeV || Ecm < 10_GeV ) {
+ std::cout << "ProcessSplit: " << " DoDiscrete: low en. particle, skipping.." << std::endl;
+ // std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
+ //fEnergy += p.GetEnergy();
+ //p.Delete();
+ //fCount++;
+ } 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_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] * si::constants::c ) / ( gamma + 1.);
- pnorm += psib.GetEnergy();
+ //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] * si::constants::c ) / ( gamma + 1.);
+ pnorm += psib.GetEnergy();
- for(int j=0; j<3; ++j){
- p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
- // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
- // << " gb: " << gammaBetaComponents[j] << endl;
- en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
- }
- // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
- // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
- // << " en lab (GeV): " << en_lab / 1_GeV << endl;
- // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
- // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
+ for(int j=0; j<3; ++j){
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
+ // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
+ // << " gb: " << gammaBetaComponents[j] << endl;
+ en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+ }
+ // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
+ // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
+ // << " en lab (GeV): " << en_lab / 1_GeV << endl;
+ // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
+ // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetEnergy( en_lab );
- pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
- //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+ // add to corsika stack
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy( en_lab );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+ //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
- 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) );
- //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
- //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
- //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
- Ptot_final += pnew.GetMomentum();
+ 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) );
+ //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+ //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ Ptot_final += pnew.GetMomentum();
+ }
+ //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
}
- //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
- }
- }else
- p.Delete();
+ }
}
void Init() {
@@ -339,7 +517,7 @@ public:
int GetCount() { return fCount; }
-
+ EnergyType GetEnergy() { return fEnergy; }
private:
};
@@ -364,7 +542,8 @@ int main() {
ProcessSplit p1;
corsika::process::sibyll::ProcessDecay p2;
- const auto sequence = p0 + p1 + p2;
+ ProcessEMCut p3;
+ const auto sequence = p0 + p1 + p2 + p3;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
@@ -383,5 +562,10 @@ int main() {
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
+ cout << "Result: E0=" << E0 / 1_GeV << "GeV, particles below energy threshold =" << p1.GetCount() << endl;
+ cout << "total energy below threshold (GeV): " << p1.GetEnergy() / 1_GeV << std::endl;
+ p3.ShowResults();
+ cout << "total energy (GeV): "
+ << ( p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy() ) / 1_GeV
+ << endl;
}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index e46960dc9..c166517fa 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -60,6 +60,8 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
+ std::cout << "+++++++++++++++++++++++++++++++" << std::endl;
+ std::cout << "Cascade: starting step.." << std::endl;
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
fProcesseList.MinStepLength(particle, step);
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 449a281b4..e3597436f 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -100,7 +100,10 @@ namespace corsika::process {
// 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;
+ int a = 1;
+ const double x = -x0 * log(s_rndm_(a));
+ cout << "ProcessDecay: next decay: " << x << endl;
+ return x;
}
template <typename Particle, typename Stack>
@@ -112,14 +115,16 @@ namespace corsika::process {
pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) );
pin.SetEnergy( p.GetEnergy() );
pin.SetMomentum( p.GetMomentum() );
+ // remove original particle from corsika stack
+ p.Delete();
// set all particles/hadrons unstable
setHadronsUnstable();
// call sibyll decay
- std::cout << "calling Sibyll decay routine.." << std::endl;
+ std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
decsib_();
// print output
- //int print_unit = 6;
- //sib_list_( print_unit );
+ int print_unit = 6;
+ sib_list_( print_unit );
// copy particles from sibyll stack to corsika
int i = -1;
for (auto &psib: ss){
@@ -135,8 +140,6 @@ namespace corsika::process {
}
// empty sibyll stack
ss.Clear();
- // remove original particle from stack
- p.Delete();
}
template <typename Particle, typename Stack>
--
GitLab
From 892e8365a05e00d8a097918678899c57e13961f6 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 10:02:30 +0100
Subject: [PATCH 22/39] removed debug printout
---
Processes/Sibyll/code_generator.py | 2 --
1 file changed, 2 deletions(-)
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index 77c1f5dd8..a593a1f2c 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -159,8 +159,6 @@ if __name__ == "__main__":
pythia_db = load_pythiadb(sys.argv[1])
read_sibyll_codes(sys.argv[2], pythia_db)
- print (str(pythia_db))
-
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
print(generate_sibyll_enum(pythia_db), file=f)
--
GitLab
From 08c06a19602f0fe9bc2a27a8a8fd69b5e52985be Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 11:28:16 +0000
Subject: [PATCH 23/39] added momentum to sibyll stack, implemented boost
between sibyll stack and corsika stack
---
Documentation/Examples/cascade_example.cc | 294 ++++++++++++++--------
Framework/Cascade/SibStack.h | 42 +++-
2 files changed, 222 insertions(+), 114 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 34b55476d..0df3f0d64 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -24,6 +24,7 @@
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
+
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
@@ -43,46 +44,59 @@ public:
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
+ const Code corsikaBeamId = p.GetPID();
+
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
- int kBeam = 1;
+ 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 = 1;
- double beamEnergy = p.GetEnergy() / 1_GeV;
- std::cout << "ProcessSplit: "
- << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
- double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
-
- if (kTarget == 1)
- sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
- else
- sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy);
-
- std::cout << "ProcessSplit: "
- << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "ProcessSplit: "
- << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
-
- const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
- std::cout << "ProcessSplit: "
- << "nucleon mass " << nucleon_mass << std::endl;
- // calculate interaction length in medium
- double int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
- // pick random step lenth
- std::cout << "ProcessSplit: "
- << "interaction length (g/cm2): " << int_length << std::endl;
- // add exponential sampling
- int a = 0;
- const double next_step = -int_length * log(s_rndm_(a));
+ int kTarget = 16;
+ double beamEnergy = p.GetEnergy() / 1_GeV;
+#warning boost to cm. still missing, sibyll cross section input is cm. energy!
+
+ std::cout << "ProcessSplit: " << "MinStep: input en: " << beamEnergy
+ << " beam can interact:" << kBeam
+ << " beam XS code:" << kBeam
+ << " beam pid:" << p.GetPID()
+ << " target mass number:" << kTarget << std::endl;
+
+ double next_step;
+ if(kInteraction){
+
+ double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
+ double dumdif[3];
+
+ if(kTarget==1)
+ sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+ else
+ sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+
+ std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ CrossSectionType sig = prodCrossSection * 1_mbarn;
+ std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+ const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
+ std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl;
+ // calculate interaction length in medium
+ double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
+ // pick random step lenth
+ std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
+ // add exponential sampling
+ int a = 0;
+ next_step = -int_length * log(s_rndm_(a));
+ }else
+#warning define infinite interaction length? then we can skip the test in DoDiscrete()
+ next_step = 1.e8;
+
/*
what are the units of the output? slant depth or 3space length?
@@ -93,86 +107,154 @@ public:
}
template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
// corsika::utls::ignore(p);
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- cout << "DoDiscrete: " << p.GetPID() << " interaction? "
- << process::sibyll::CanInteract(p.GetPID()) << endl;
- if (process::sibyll::CanInteract(p.GetPID())) {
-
+ cout << "DoDiscrete: " << 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; //p.GetPID();
+
+ // proton mass in units of energy
+ const EnergyType proton_mass_en = 0.93827_GeV ; //0.93827_GeV / si::constants::cSquared ;
+
+ cout << "defining target momentum.." << endl;
+ // FOR NOW: target is always at rest
+ const EnergyType Etarget = 0. * 1_GeV + proton_mass_en;
+ const auto pTarget = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+ cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
+ // const auto pBeam = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+ // cout << "beam momentum: " << pBeam.GetComponents() << endl;
+ cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV * si::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)
+ 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)
*/
- EnergyType E = p.GetEnergy();
- EnergyType Ecm = sqrt(2. * E * 0.93827_GeV);
-
- int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // cout << "defining total energy" << endl;
+ // total energy: E_beam + E_target
+ // in lab. frame: E_beam + m_target*c**2
+ EnergyType E = p.GetEnergy();
+ EnergyType Etot = E + Etarget;
+ // cout << "tot. energy: " << Etot / 1_GeV << endl;
+ // cout << "defining total momentum" << endl;
+ // total momentum
+ super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+ // cout << "tot. momentum: " << Ptot.GetComponents() / 1_GeV * si::constants::c << endl;
+ // cout << "inv. mass.." << endl;
+ // invariant mass, i.e. cm. energy
+ EnergyType Ecm = sqrt( Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared ); //sqrt( 2. * E * 0.93827_GeV );
+ // cout << "inv. mass: " << Ecm / 1_GeV << endl;
+ // cout << "boost parameters.." << endl;
+ /*
+ 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 = ( E + proton_mass * si::constants::cSquared ) / Ecm ;
+ // const double gambet = sqrt( E * E - proton_mass * proton_mass ) / Ecm;
+ const double gamma = Etot / Ecm ;
+ const auto gambet = Ptot / (Ecm / si::constants::c );
+
+ std::cout << "ProcessSplit: " << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "ProcessSplit: " << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+
+
+ std::cout << "ProcessSplit: " << " DoDiscrete: 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;
+ p.Delete();
+ fCount++;
+ } else {
+ // Sibyll does not know about units..
+ double sqs = Ecm / 1_GeV;
+ // running sibyll, filling stack
+ sibyll_( kBeam, kTarget, sqs);
+ // running decays
+ //decsib_();
+ // print final state
+ int print_unit = 6;
+ sib_list_( print_unit );
+
+ // delete current particle
+ p.Delete();
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- // FOR NOW: set target to proton
- int kTarget = 1; // p.GetPID();
-
- std::cout << "ProcessSplit: "
- << " DoDiscrete: 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;
- p.Delete();
- fCount++;
- } else {
- // Sibyll does not know about units..
- double sqs = Ecm / 1_GeV;
- // running sibyll, filling stack
- sibyll_(kBeam, kTarget, sqs);
- // running decays
- // 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;
- /*
- 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
- in general transformation is rotation + boost
- */
- const EnergyType proton_mass = 0.93827_GeV;
- const double gamma = (E + proton_mass) / (Ecm);
- const double gambet = sqrt(E * E - proton_mass * proton_mass) / Ecm;
-
- // SibStack does not know about momentum yet so we need counter to access momentum
- // array in Sibyll
- int i = -1;
- for (auto& p : ss) {
- ++i;
- // transform to lab. frame, primitve
- const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetEnergy(en_lab * 1_GeV);
- pnew.SetPID(process::sibyll::ConvertFromSibyll(p.GetPID()));
- }
+ // 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;
+ //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] * si::constants::c ) / ( gamma + 1.);
+ pnorm += psib.GetEnergy();
+
+ for(int j=0; j<3; ++j){
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
+ // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
+ // << " gb: " << gammaBetaComponents[j] << endl;
+ en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+ }
+ // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
+ // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
+ // << " en lab (GeV): " << en_lab / 1_GeV << endl;
+ // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
+ // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
+
+ // add to corsika stack
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy( en_lab );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+ //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+
+ 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) );
+ //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+ //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ Ptot_final += pnew.GetMomentum();
}
- } else
+ //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
+ }
+ }else
p.Delete();
}
@@ -238,6 +320,8 @@ double s_rndm_(int&) {
int main() {
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
@@ -248,9 +332,15 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV;
+ EnergyType E0 = 500_GeV;
+ MomentumType P0 = sqrt( E0*E0 - 0.93827_GeV * 0.93827_GeV ) / si::constants::c;
+ auto plab = super_stupid::MomentumVector(rootCS,
+ 0. * 1_GeV / si::constants::c,
+ 0. * 1_GeV / si::constants::c,
+ P0);
particle.SetEnergy(E0);
- particle.SetPID(Code::Proton);
+ particle.SetMomentum(plab);
+ particle.SetPID( Code::Proton );
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 7e9bd3bac..57a3f9c16 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -10,6 +10,8 @@
using namespace std;
using namespace corsika::stack;
+using namespace corsika::units;
+using namespace corsika::geometry;
class SibStackData {
@@ -17,16 +19,33 @@ public:
void Init();
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
int GetCapacity() const { return 8000; }
- void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
- void SetEnergy(const int i, const double v) { s_plist_.p[3][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 SetMomentum(const int i, const super_stupid::MomentumVector& v)
+ {
+ auto tmp = v.GetComponents();
+ for(int idx=0; idx<3; ++idx)
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;
+ }
+
int GetId(const int i) const { return s_plist_.llist[i]; }
- double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
-
+
+ EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+
+ super_stupid::MomentumVector GetMomentum(const int i) const
+ {
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+ corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
+ super_stupid::MomentumVector v1(rootCS,components);
+ return v1;
+ }
+
void Copy(const int i1, const int i2) {
s_plist_.llist[i1] = s_plist_.llist[i2];
s_plist_.p[3][i1] = s_plist_.p[3][i2];
@@ -46,12 +65,11 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
public:
void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
- double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ 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()));
- }
+ 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()); }
+
};
typedef Stack<SibStackData, ParticleInterface> SibStack;
--
GitLab
From 8fd298e29e990b2c6678fadde88e7e23856b5f67 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 2 Dec 2018 10:01:59 +0100
Subject: [PATCH 24/39] added comment
---
CMakeLists.txt | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2a5cae812..1ac8449f0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,6 +7,7 @@ project (
LANGUAGES CXX
)
+# as long as there still are modules using it:
enable_language (Fortran)
# ignore many irrelevant Up-to-date messages during install
@@ -38,7 +39,6 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
-
# unit testing coverage, does not work yet
#include (CodeCoverage)
##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')
--
GitLab
From f21be6305a8c0b2727cd7d8d506c66c294e98782 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 12:36:03 +0000
Subject: [PATCH 25/39] added cm energy calculation for minsteplength in sibyll
process
---
Documentation/Examples/cascade_example.cc | 33 ++++++++++++++++-------
1 file changed, 24 insertions(+), 9 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 0df3f0d64..abb635ac6 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -44,6 +44,9 @@ public:
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
const Code corsikaBeamId = p.GetPID();
// beam particles for sibyll : 1, 2, 3 for p, pi, k
@@ -60,13 +63,24 @@ public:
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
int kTarget = 16;
- double beamEnergy = p.GetEnergy() / 1_GeV;
-#warning boost to cm. still missing, sibyll cross section input is cm. energy!
+
+ // proton mass in units of energy
+ const EnergyType proton_mass_en = 0.93827_GeV ;
+
+ EnergyType Etot = p.GetEnergy() + kTarget * proton_mass_en;
+ 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() * si::constants::cSquared );
+ double Ecm = sqs / 1_GeV;
- std::cout << "ProcessSplit: " << "MinStep: input en: " << beamEnergy
- << " beam can interact:" << kBeam
- << " beam XS code:" << kBeam
- << " beam pid:" << p.GetPID()
+ std::cout << "ProcessSplit: " << "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 next_step;
@@ -76,9 +90,9 @@ public:
double dumdif[3];
if(kTarget==1)
- sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+ sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
else
- sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+ sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy );
std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
CrossSectionType sig = prodCrossSection * 1_mbarn;
@@ -291,6 +305,7 @@ public:
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);
@@ -332,7 +347,7 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 500_GeV;
+ EnergyType E0 = 100_GeV;
MomentumType P0 = sqrt( E0*E0 - 0.93827_GeV * 0.93827_GeV ) / si::constants::c;
auto plab = super_stupid::MomentumVector(rootCS,
0. * 1_GeV / si::constants::c,
--
GitLab
From c4572c99ce8333908de5045d8354c421c155d2d1 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 08:20:41 +0000
Subject: [PATCH 26/39] added proto decay process
---
Documentation/Examples/cascade_example.cc | 19 +++++++++++++++++++
1 file changed, 19 insertions(+)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index abb635ac6..f208cf296 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -333,6 +333,25 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
+class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+public:
+ ProcessDecay() {}
+ void Init() {}
+ template <typename Particle>
+ double MinStepLength(Particle& p) const {
+ }
+ template <typename Particle, typename Stack>
+ void DoDiscrete(Particle& p, Stack& s) const {
+ }
+
+ template <typename Particle, typename Trajectory, typename Stack>
+ EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+};
+
+
int main() {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
--
GitLab
From 2a862a95c3ed95a1910471245bd92f9611dc5c02 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 17:43:21 +0000
Subject: [PATCH 27/39] added mass density type
---
Framework/Units/PhysicalUnits.h | 9 +++++++++
1 file changed, 9 insertions(+)
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 8975a6e99..0bcaa7a45 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -51,6 +51,15 @@ 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 MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
+
+ // defining momentum you suckers
+ // dimensions, i.e. composition in base SI dimensions
+ using momentum_d = phys::units::dimensions< 1, 1, -1 >;
+ // defining the unit of momentum, so far newton-meter, maybe go to HEP?
+ constexpr phys::units::quantity< momentum_d > newton_second { meter * kilogram / second };
+ // defining the type
+
using MomentumType = phys::units::quantity<momentum_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
--
GitLab
From b93be13c3a5efbcc6ad6f68461294c94a67efd55 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 17:44:17 +0000
Subject: [PATCH 28/39] added primitive decay step length
---
Documentation/Examples/cascade_example.cc | 37 +++++++++++++++++++++--
1 file changed, 34 insertions(+), 3 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f208cf296..b44863358 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -337,12 +337,42 @@ class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
public:
ProcessDecay() {}
void Init() {}
+
template <typename Particle>
double MinStepLength(Particle& p) const {
+ EnergyType E = p.GetEnergy();
+ MassType m = corsika::particles::GetMass(p.GetPID());
+ // env.GetDensity();
+ const MassDensityType density = 1.e3 * kilogram / ( 1_cm * 1_cm * 1_cm );
+
+ const double gamma = E / m / constants::cSquared;
+ // lifetimes not implemented yet
+ TimeType t0;
+ switch( p.GetPID() ){
+ case Code::PiPlus :
+ t0 = 1.e-5 * 1_s;
+ break;
+
+ case Code::KPlus :
+ t0 = 1.e-5 * 1_s;
+ break;
+
+ default:
+ t0 = 1.e8 * 1_s;
+ break;
+ }
+ 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>
+
+ template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- }
+ }
template <typename Particle, typename Trajectory, typename Stack>
EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
@@ -359,7 +389,8 @@ int main() {
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
- const auto sequence = p0 + p1;
+ ProcessDecay p2;
+ const auto sequence = p0 + p1 + p2;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
--
GitLab
From 14fa58d6eee6d43de04655cce97e065fb5afe1ee Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 30 Nov 2018 17:50:29 +0000
Subject: [PATCH 29/39] added more realistic numbers for decay
---
Documentation/Examples/cascade_example.cc | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index b44863358..2152ec810 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -343,14 +343,14 @@ public:
EnergyType E = p.GetEnergy();
MassType m = corsika::particles::GetMass(p.GetPID());
// env.GetDensity();
- const MassDensityType density = 1.e3 * 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 / constants::cSquared;
// lifetimes not implemented yet
TimeType t0;
switch( p.GetPID() ){
case Code::PiPlus :
- t0 = 1.e-5 * 1_s;
+ t0 = 2.6e-8 * 1_s;
break;
case Code::KPlus :
--
GitLab
From b4a890b18507e8edb7d7dd39556c15ef9f8866f8 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sun, 2 Dec 2018 11:05:34 +0000
Subject: [PATCH 30/39] fixed rebase
---
Framework/Units/PhysicalUnits.h | 9 ---------
1 file changed, 9 deletions(-)
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 0bcaa7a45..0f6c82c2e 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -52,15 +52,6 @@ namespace corsika::units::si {
using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
-
- // defining momentum you suckers
- // dimensions, i.e. composition in base SI dimensions
- using momentum_d = phys::units::dimensions< 1, 1, -1 >;
- // defining the unit of momentum, so far newton-meter, maybe go to HEP?
- constexpr phys::units::quantity< momentum_d > newton_second { meter * kilogram / second };
- // defining the type
-
- using MomentumType = phys::units::quantity<momentum_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
} // end namespace corsika::units::si
--
GitLab
From efbdd4ca978fb717d9842850d155279536e97873 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 12:56:08 +0000
Subject: [PATCH 31/39] fixed rebase
---
Documentation/Examples/cascade_example.cc | 6 +++---
Framework/Units/PhysicalUnits.h | 3 ++-
2 files changed, 5 insertions(+), 4 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 2152ec810..a37ced7d1 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -339,7 +339,7 @@ public:
void Init() {}
template <typename Particle>
- double MinStepLength(Particle& p) const {
+ double MinStepLength(Particle& p, setup::Trajectory&) const {
EnergyType E = p.GetEnergy();
MassType m = corsika::particles::GetMass(p.GetPID());
// env.GetDensity();
@@ -374,8 +374,8 @@ public:
void DoDiscrete(Particle& p, Stack& s) const {
}
- template <typename Particle, typename Trajectory, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
return EProcessReturn::eOk;
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 0f6c82c2e..8a9e9d6f6 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -53,7 +53,8 @@ namespace corsika::units::si {
using MassType = phys::units::quantity<phys::units::mass_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
/**
--
GitLab
From bbb8f0ecbba6007ffa2afed1bb401f3ba9142449 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 4 Dec 2018 13:13:24 +0000
Subject: [PATCH 32/39] moved decay into sibyll process
---
Documentation/Examples/cascade_example.cc | 52 ++----------------
Processes/Sibyll/CMakeLists.txt | 1 +
Processes/Sibyll/ProcessDecay.h | 65 +++++++++++++++++++++++
3 files changed, 69 insertions(+), 49 deletions(-)
create mode 100644 Processes/Sibyll/ProcessDecay.h
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index a37ced7d1..fdc706949 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -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;
@@ -333,54 +335,6 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
-class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
-public:
- ProcessDecay() {}
- void Init() {}
-
- 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;
- // lifetimes not implemented yet
- TimeType t0;
- switch( p.GetPID() ){
- case Code::PiPlus :
- t0 = 2.6e-8 * 1_s;
- break;
-
- case Code::KPlus :
- t0 = 1.e-5 * 1_s;
- break;
-
- default:
- t0 = 1.e8 * 1_s;
- break;
- }
- 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 {
- }
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
- }
-
-};
-
int main() {
@@ -389,7 +343,7 @@ int main() {
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
- ProcessDecay p2;
+ corsika::process::sibyll::ProcessDecay p2;
const auto sequence = p0 + p1 + p2;
setup::Stack stack;
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index f2e4a9fee..a0985a1ad 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -20,6 +20,7 @@ set (
set (
MODEL_HEADERS
ParticleConversion.h
+ ProcessDecay.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
new file mode 100644
index 000000000..20f493607
--- /dev/null
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -0,0 +1,65 @@
+#ifndef _include_ProcessDecay_h_
+#define _include_ProcessDecay_h_
+
+#include <corsika/process/ContinuousProcess.h>
+
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+
+//using namespace corsika::particles;
+
+namespace corsika::process {
+
+ namespace sibyll {
+
+class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+public:
+ ProcessDecay() {}
+ void Init() {}
+
+ 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;
+ // lifetimes not implemented yet
+ TimeType t0;
+ switch( p.GetPID() ){
+ case corsika::particles::Code::PiPlus :
+ t0 = 2.6e-8 * 1_s;
+ break;
+
+ case corsika::particles::Code::KPlus :
+ t0 = 1.e-5 * 1_s;
+ break;
+
+ default:
+ t0 = 1.e8 * 1_s;
+ break;
+ }
+ 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 {
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+};
+ }
+}
+
+#endif
--
GitLab
From b92396ee84f45be12cd2aad19079ff893d6ecd94 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 5 Dec 2018 08:17:39 +0000
Subject: [PATCH 33/39] added negative charge pions and kaons
---
Processes/Sibyll/ProcessDecay.h | 9 +++++++++
1 file changed, 9 insertions(+)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 20f493607..4989ab90a 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -32,9 +32,17 @@ public:
t0 = 2.6e-8 * 1_s;
break;
+ case corsika::particles::Code::PiMinus :
+ t0 = 2.6e-8 * 1_s;
+ break;
+
case corsika::particles::Code::KPlus :
t0 = 1.e-5 * 1_s;
break;
+
+ case corsika::particles::Code::KMinus :
+ t0 = 1.e-5 * 1_s;
+ break;
default:
t0 = 1.e8 * 1_s;
@@ -51,6 +59,7 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
+
}
template <typename Particle, typename Stack>
--
GitLab
From 3c9b91067dbd1f8ed9196a7c6be9c23e1eb34fb6 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 5 Dec 2018 10:36:39 +0000
Subject: [PATCH 34/39] added decay configuration routines, turned on decay for
sibyll interaction
---
Documentation/Examples/cascade_example.cc | 58 +++++----
Processes/Sibyll/ProcessDecay.h | 145 ++++++++++++++--------
2 files changed, 131 insertions(+), 72 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index fdc706949..5c631f5a0 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -43,6 +43,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 {
@@ -211,7 +240,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 );
@@ -228,6 +258,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
@@ -300,29 +333,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:
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 4989ab90a..ecae303e7 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -12,62 +12,107 @@ namespace corsika::process {
namespace sibyll {
-class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
-public:
- ProcessDecay() {}
- void Init() {}
-
- 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 );
+
+ 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();
+ }
+
+ }
+
- const double gamma = E / m / constants::cSquared;
- // lifetimes not implemented yet
- TimeType t0;
- switch( p.GetPID() ){
- case corsika::particles::Code::PiPlus :
- t0 = 2.6e-8 * 1_s;
- break;
-
- case corsika::particles::Code::PiMinus :
- t0 = 2.6e-8 * 1_s;
- break;
+ class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+ public:
+ ProcessDecay() {}
+ void Init() {
+ //setHadronsUnstable();
+ }
- case corsika::particles::Code::KPlus :
- t0 = 1.e-5 * 1_s;
- break;
+ 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;
+ }
+ }
- case corsika::particles::Code::KMinus :
- t0 = 1.e-5 * 1_s;
- break;
+ friend void setHadronsUnstable();
- default:
- t0 = 1.e8 * 1_s;
- break;
- }
- 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 {
+ 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 );
- }
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
- }
+ 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 {
+
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+ };
}
}
--
GitLab
From d10704f0e3716ddc6c569f12635345ed198ad88b Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 6 Dec 2018 10:59:14 +0000
Subject: [PATCH 35/39] added actual decay to ProcessDecay
---
Framework/Cascade/SibStack.h | 4 +++-
Processes/Sibyll/ProcessDecay.h | 34 ++++++++++++++++++++++++++++++++-
2 files changed, 36 insertions(+), 2 deletions(-)
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 57a3f9c16..2275cdf49 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -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); }
};
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index ecae303e7..ab11fc7e4 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -5,6 +5,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/cascade/SibStack.h>
//using namespace corsika::particles;
@@ -104,7 +105,38 @@ namespace corsika::process {
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>
--
GitLab
From 879480221679ae2051f1b2384c9f49317fb742fb Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 6 Dec 2018 11:01:35 +0000
Subject: [PATCH 36/39] suppressed output
---
Processes/Sibyll/ProcessDecay.h | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index ab11fc7e4..449a281b4 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -118,8 +118,8 @@ namespace corsika::process {
std::cout << "calling Sibyll decay routine.." << std::endl;
decsib_();
// print output
- int print_unit = 6;
- sib_list_( print_unit );
+ //int print_unit = 6;
+ //sib_list_( print_unit );
// copy particles from sibyll stack to corsika
int i = -1;
for (auto &psib: ss){
@@ -127,7 +127,7 @@ namespace corsika::process {
// 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;
+ //cout << "decay product: " << process::sibyll::ConvertFromSibyll( psib.GetPID() ) << endl;
auto pnew = s.NewParticle();
pnew.SetEnergy( psib.GetEnergy() );
pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
--
GitLab
From 853e9aad4a22fb6afa4ab41a63fea4ea935f49d3 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sat, 8 Dec 2018 13:56:34 +0000
Subject: [PATCH 37/39] added particle cut process
---
Documentation/Examples/cascade_example.cc | 336 +++++++++++++++++-----
Framework/Cascade/Cascade.h | 2 +
Processes/Sibyll/ProcessDecay.h | 15 +-
3 files changed, 271 insertions(+), 82 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 5c631f5a0..f7a968723 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -38,6 +38,183 @@ using namespace corsika::random;
using namespace std;
static int fCount = 0;
+static EnergyType fEnergy = 0. * 1_GeV;
+
+// FOR NOW: global static variables for ParticleCut process
+// this is just wrong...
+static EnergyType fEmEnergy;
+static int fEmCount;
+
+static EnergyType fInvEnergy;
+static int fInvCount;
+
+
+class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
+public:
+ ProcessEMCut() {}
+ template <typename Particle>
+ bool isBelowEnergyCut( Particle &p ) const
+ {
+ // FOR NOW: center-of-mass energy hard coded
+ const EnergyType Ecm = sqrt( 2. * p.GetEnergy() * 0.93827_GeV );
+ if( p.GetEnergy() < 50_GeV || Ecm < 10_GeV )
+ return true;
+ else
+ return false;
+ }
+
+ bool isEmParticle( Code pCode ) const
+ {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch( pCode ){
+ case Code::Electron :
+ is_em = true;
+ break;
+ case Code::Gamma :
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ void defineEmParticles() const
+ {
+ // create bool array identifying em particles
+ }
+
+ bool isInvisible( Code pCode ) const
+ {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch( pCode ){
+ case Code::NuE :
+ is_inv = true;
+ break;
+ case Code::NuEBar :
+ is_inv = true;
+ break;
+ case Code::NuMu :
+ is_inv = true;
+ break;
+ case Code::NuMuBar :
+ is_inv = true;
+ break;
+ case Code::MuPlus :
+ is_inv = true;
+ break;
+ case Code::MuMinus :
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+
+ template <typename Particle>
+ double MinStepLength(Particle& p, setup::Trajectory&) const
+ {
+ const Code pid = p.GetPID();
+ if( isEmParticle( pid ) || isInvisible( pid ) ){
+ cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+ return 0.;
+ } else {
+ double next_step = std::numeric_limits<double>::infinity();
+ cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+ return next_step;
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&p, setup::Trajectory&t, Stack&s) const {
+ // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
+ // 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;
+ }
+
+ template <typename Particle, typename Stack>
+ void DoDiscrete(Particle& p, Stack& s) const
+ {
+ cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
+ const Code pid = p.GetPID();
+ if( isEmParticle( pid ) ){
+ cout << "removing em. particle..." << endl;
+ fEmEnergy += p.GetEnergy();
+ fEmCount += 1;
+ p.Delete();
+ }
+ else if ( isInvisible( pid ) ){
+ cout << "removing inv. particle..." << endl;
+ fInvEnergy += p.GetEnergy();
+ fInvCount += 1;
+ p.Delete();
+ }
+ else if( isBelowEnergyCut( p ) ){
+ cout << "removing low en. particle..." << endl;
+ fEnergy += p.GetEnergy();
+ fCount += 1;
+ p.Delete();
+ }
+ }
+
+ void Init()
+ {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ //defineEmParticles();
+ }
+
+ void ShowResults(){
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " ******************************" << endl;
+ }
+
+ EnergyType GetInvEnergy()
+ {
+ return fInvEnergy;
+ }
+
+ EnergyType GetCutEnergy()
+ {
+ return fEnergy;
+ }
+
+ EnergyType GetEmEnergy()
+ {
+ return fEmEnergy;
+ }
+
+private:
+};
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
@@ -140,14 +317,14 @@ public:
next_step = -int_length * log(s_rndm_(a));
}else
#warning define infinite interaction length? then we can skip the test in DoDiscrete()
- next_step = 1.e8;
+ next_step = std::numeric_limits<double>::infinity();
/*
what are the units of the output? slant depth or 3space length?
*/
std::cout << "ProcessSplit: "
- << "next step (g/cm2): " << next_step << std::endl;
+ << "next interaction (g/cm2): " << next_step << std::endl;
return next_step;
}
@@ -159,7 +336,7 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
+ cout << "ProcessSplit: " << "DoDiscrete: " << 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
@@ -230,81 +407,82 @@ public:
std::cout << "ProcessSplit: " << " DoDiscrete: 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;
- p.Delete();
- fCount++;
- } 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_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;
+ if (E < 8.5_GeV || Ecm < 10_GeV ) {
+ std::cout << "ProcessSplit: " << " DoDiscrete: low en. particle, skipping.." << std::endl;
+ // std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
+ //fEnergy += p.GetEnergy();
+ //p.Delete();
+ //fCount++;
+ } 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_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] * si::constants::c ) / ( gamma + 1.);
- pnorm += psib.GetEnergy();
+ //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] * si::constants::c ) / ( gamma + 1.);
+ pnorm += psib.GetEnergy();
- for(int j=0; j<3; ++j){
- p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
- // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
- // << " gb: " << gammaBetaComponents[j] << endl;
- en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
- }
- // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
- // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
- // << " en lab (GeV): " << en_lab / 1_GeV << endl;
- // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
- // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
+ for(int j=0; j<3; ++j){
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
+ // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
+ // << " gb: " << gammaBetaComponents[j] << endl;
+ en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+ }
+ // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
+ // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
+ // << " en lab (GeV): " << en_lab / 1_GeV << endl;
+ // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
+ // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetEnergy( en_lab );
- pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
- //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+ // add to corsika stack
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy( en_lab );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+ //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
- 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) );
- //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
- //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
- //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
- Ptot_final += pnew.GetMomentum();
+ 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) );
+ //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+ //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+ Ptot_final += pnew.GetMomentum();
+ }
+ //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
}
- //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
- }
- }else
- p.Delete();
+ }
}
void Init() {
@@ -338,7 +516,7 @@ public:
int GetCount() { return fCount; }
-
+ EnergyType GetEnergy() { return fEnergy; }
private:
};
@@ -358,7 +536,8 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
corsika::process::sibyll::ProcessDecay p2;
- const auto sequence = p0 + p1 + p2;
+ ProcessEMCut p3;
+ const auto sequence = p0 + p1 + p2 + p3;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
@@ -376,5 +555,10 @@ int main() {
particle.SetPID( Code::Proton );
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
+ cout << "Result: E0=" << E0 / 1_GeV << "GeV, particles below energy threshold =" << p1.GetCount() << endl;
+ cout << "total energy below threshold (GeV): " << p1.GetEnergy() / 1_GeV << std::endl;
+ p3.ShowResults();
+ cout << "total energy (GeV): "
+ << ( p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy() ) / 1_GeV
+ << endl;
}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index e46960dc9..c166517fa 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -60,6 +60,8 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
+ std::cout << "+++++++++++++++++++++++++++++++" << std::endl;
+ std::cout << "Cascade: starting step.." << std::endl;
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
fProcesseList.MinStepLength(particle, step);
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 449a281b4..e3597436f 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -100,7 +100,10 @@ namespace corsika::process {
// 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;
+ int a = 1;
+ const double x = -x0 * log(s_rndm_(a));
+ cout << "ProcessDecay: next decay: " << x << endl;
+ return x;
}
template <typename Particle, typename Stack>
@@ -112,14 +115,16 @@ namespace corsika::process {
pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) );
pin.SetEnergy( p.GetEnergy() );
pin.SetMomentum( p.GetMomentum() );
+ // remove original particle from corsika stack
+ p.Delete();
// set all particles/hadrons unstable
setHadronsUnstable();
// call sibyll decay
- std::cout << "calling Sibyll decay routine.." << std::endl;
+ std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
decsib_();
// print output
- //int print_unit = 6;
- //sib_list_( print_unit );
+ int print_unit = 6;
+ sib_list_( print_unit );
// copy particles from sibyll stack to corsika
int i = -1;
for (auto &psib: ss){
@@ -135,8 +140,6 @@ namespace corsika::process {
}
// empty sibyll stack
ss.Clear();
- // remove original particle from stack
- p.Delete();
}
template <typename Particle, typename Stack>
--
GitLab
From 72bd4ba21aea367f3f6077b79620e2c2bec28a58 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sat, 8 Dec 2018 19:35:45 +0000
Subject: [PATCH 38/39] clean and clang
---
Documentation/Examples/cascade_example.cc | 516 ++++++++++------------
Framework/Cascade/SibStack.h | 54 ++-
Framework/Units/PhysicalUnits.h | 2 +-
Processes/Sibyll/ProcessDecay.h | 197 +++++----
4 files changed, 366 insertions(+), 403 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f7a968723..4132fae2f 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -42,88 +42,81 @@ static EnergyType fEnergy = 0. * 1_GeV;
// FOR NOW: global static variables for ParticleCut process
// this is just wrong...
-static EnergyType fEmEnergy;
-static int fEmCount;
-
-static EnergyType fInvEnergy;
-static int fInvCount;
+static EnergyType fEmEnergy;
+static int fEmCount;
+static EnergyType fInvEnergy;
+static int fInvCount;
class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
public:
ProcessEMCut() {}
template <typename Particle>
- bool isBelowEnergyCut( Particle &p ) const
- {
+ bool isBelowEnergyCut(Particle& p) const {
// FOR NOW: center-of-mass energy hard coded
- const EnergyType Ecm = sqrt( 2. * p.GetEnergy() * 0.93827_GeV );
- if( p.GetEnergy() < 50_GeV || Ecm < 10_GeV )
+ const EnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+ if (p.GetEnergy() < 50_GeV || Ecm < 10_GeV)
return true;
else
return false;
}
- bool isEmParticle( Code pCode ) const
- {
+ bool isEmParticle(Code pCode) const {
bool is_em = false;
// FOR NOW: switch
- switch( pCode ){
- case Code::Electron :
- is_em = true;
- break;
- case Code::Gamma :
- is_em = true;
- break;
- default:
- break;
+ switch (pCode) {
+ case Code::Electron:
+ is_em = true;
+ break;
+ case Code::Gamma:
+ is_em = true;
+ break;
+ default:
+ break;
}
return is_em;
}
- void defineEmParticles() const
- {
- // create bool array identifying em particles
+ void defineEmParticles() const {
+ // create bool array identifying em particles
}
- bool isInvisible( Code pCode ) const
- {
+ bool isInvisible(Code pCode) const {
bool is_inv = false;
// FOR NOW: switch
- switch( pCode ){
- case Code::NuE :
- is_inv = true;
- break;
- case Code::NuEBar :
- is_inv = true;
- break;
- case Code::NuMu :
- is_inv = true;
- break;
- case Code::NuMuBar :
- is_inv = true;
- break;
- case Code::MuPlus :
- is_inv = true;
- break;
- case Code::MuMinus :
- is_inv = true;
- break;
-
- default:
- break;
+ switch (pCode) {
+ case Code::NuE:
+ is_inv = true;
+ break;
+ case Code::NuEBar:
+ is_inv = true;
+ break;
+ case Code::NuMu:
+ is_inv = true;
+ break;
+ case Code::NuMuBar:
+ is_inv = true;
+ break;
+ case Code::MuPlus:
+ is_inv = true;
+ break;
+ case Code::MuMinus:
+ is_inv = true;
+ break;
+
+ default:
+ break;
}
return is_inv;
}
-
template <typename Particle>
- double MinStepLength(Particle& p, setup::Trajectory&) const
- {
+ double MinStepLength(Particle& p, setup::Trajectory&) const {
const Code pid = p.GetPID();
- if( isEmParticle( pid ) || isInvisible( pid ) ){
+ if (isEmParticle(pid) || isInvisible(pid)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
return 0.;
- } else {
+ } else {
double next_step = std::numeric_limits<double>::infinity();
cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
return next_step;
@@ -131,7 +124,7 @@ public:
}
template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&p, setup::Trajectory&t, Stack&s) const {
+ EProcessReturn DoContinuous(Particle& p, setup::Trajectory& t, Stack& s) const {
// cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
// cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
// cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
@@ -150,68 +143,55 @@ public:
// p.Delete();
// return EProcessReturn::eParticleAbsorbed;
// }
- return EProcessReturn::eOk;
+ return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const
- {
+ void DoDiscrete(Particle& p, Stack& s) const {
cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
const Code pid = p.GetPID();
- if( isEmParticle( pid ) ){
+ if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += p.GetEnergy();
- fEmCount += 1;
+ fEmCount += 1;
p.Delete();
- }
- else if ( isInvisible( pid ) ){
+ } else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += p.GetEnergy();
- fInvCount += 1;
+ fInvCount += 1;
p.Delete();
- }
- else if( isBelowEnergyCut( p ) ){
+ } else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy();
- fCount += 1;
+ fCount += 1;
p.Delete();
}
}
- void Init()
- {
- fEmEnergy = 0. * 1_GeV;
- fEmCount = 0;
+ void Init() {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
fInvEnergy = 0. * 1_GeV;
- fInvCount = 0;
- fEnergy = 0. * 1_GeV;
- //defineEmParticles();
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ // defineEmParticles();
}
- void ShowResults(){
+ void ShowResults() {
cout << " ******************************" << endl
- << " ParticleCut: " << endl
- << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
- << " no. of em. particles injected: " << fEmCount << endl
- << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
- << " no. of inv. particles injected: " << fInvCount << endl
- << " ******************************" << endl;
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " ******************************" << endl;
}
- EnergyType GetInvEnergy()
- {
- return fInvEnergy;
- }
+ EnergyType GetInvEnergy() { return fInvEnergy; }
- EnergyType GetCutEnergy()
- {
- return fEnergy;
- }
+ EnergyType GetCutEnergy() { return fEnergy; }
- EnergyType GetEmEnergy()
- {
- return fEmEnergy;
- }
+ EnergyType GetEmEnergy() { return fEmEnergy; }
private:
};
@@ -225,7 +205,7 @@ public:
Sibyll is hadronic generator
only hadrons decay
*/
- // set particles unstable
+ // set particles unstable
corsika::process::sibyll::setHadronsUnstable();
// make tracked particles stable
std::cout << "ProcessSplit: setting tracked hadrons stable.." << std::endl;
@@ -236,34 +216,30 @@ public:
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] );
+ 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 {
// coordinate system, get global frame of reference
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
+
const 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 );
+ int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
- bool kInteraction = process::sibyll::CanInteract( 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..
@@ -272,53 +248,56 @@ public:
// FOR NOW: assume target is oxygen
int kTarget = 16;
- // proton mass in units of energy
- const EnergyType proton_mass_en = 0.93827_GeV ;
-
- EnergyType Etot = p.GetEnergy() + kTarget * proton_mass_en;
- super_stupid::MomentumVector Ptot(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+ 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});
+ 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() * si::constants::cSquared );
+ EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared);
double Ecm = sqs / 1_GeV;
-
- std::cout << "ProcessSplit: " << "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;
+
+ std::cout << "ProcessSplit: "
+ << "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 next_step;
- if(kInteraction){
-
- double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
+ 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 );
+
+ 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 << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
+
+ std::cout << "ProcessSplit: "
+ << "MinStep: sibyll return: " << prodCrossSection << std::endl;
CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+ std::cout << "ProcessSplit: "
+ << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
- std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl;
+ std::cout << "ProcessSplit: "
+ << "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
- double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
+ double int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
// pick random step lenth
- std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
+ std::cout << "ProcessSplit: "
+ << "interaction length (g/cm2): " << int_length << std::endl;
// add exponential sampling
int a = 0;
next_step = -int_length * log(s_rndm_(a));
- }else
-#warning define infinite interaction length? then we can skip the test in DoDiscrete()
+ } else
next_step = std::numeric_limits<double>::infinity();
-
+
/*
what are the units of the output? slant depth or 3space length?
@@ -336,172 +315,152 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- cout << "ProcessSplit: " << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
- if( process::sibyll::CanInteract( p.GetPID() ) ){
+ cout << "ProcessSplit: "
+ << "DoDiscrete: " << 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)
+
+ /*
+ 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; //p.GetPID();
-
- // proton mass in units of energy
- const EnergyType proton_mass_en = 0.93827_GeV ; //0.93827_GeV / si::constants::cSquared ;
+ int kTarget = 1; // env.GetTargetParticle().GetPID();
cout << "defining target momentum.." << endl;
// FOR NOW: target is always at rest
- const EnergyType Etarget = 0. * 1_GeV + proton_mass_en;
- const auto pTarget = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
- cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
- // const auto pBeam = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
- // cout << "beam momentum: " << pBeam.GetComponents() << endl;
- cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
-
+ const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
+ const auto pTarget = super_stupid::MomentumVector(
+ rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c,
+ 0. * 1_GeV / si::constants::c);
+ cout << "target momentum (GeV/c): "
+ << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
+ cout << "beam momentum (GeV/c): "
+ << p.GetMomentum().GetComponents() / 1_GeV * si::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)
+ 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)
*/
- // cout << "defining total energy" << endl;
// total energy: E_beam + E_target
// in lab. frame: E_beam + m_target*c**2
- EnergyType E = p.GetEnergy();
- EnergyType Etot = E + Etarget;
- // cout << "tot. energy: " << Etot / 1_GeV << endl;
- // cout << "defining total momentum" << endl;
+ EnergyType E = p.GetEnergy();
+ EnergyType Etot = E + Etarget;
// total momentum
super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
- // cout << "tot. momentum: " << Ptot.GetComponents() / 1_GeV * si::constants::c << endl;
- // cout << "inv. mass.." << endl;
// invariant mass, i.e. cm. energy
- EnergyType Ecm = sqrt( Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared ); //sqrt( 2. * E * 0.93827_GeV );
- // cout << "inv. mass: " << Ecm / 1_GeV << endl;
- // cout << "boost parameters.." << endl;
- /*
- 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 = ( E + proton_mass * si::constants::cSquared ) / Ecm ;
- // const double gambet = sqrt( E * E - proton_mass * proton_mass ) / Ecm;
- const double gamma = Etot / Ecm ;
- const auto gambet = Ptot / (Ecm / si::constants::c );
-
- std::cout << "ProcessSplit: " << " DoDiscrete: gamma:" << gamma << endl;
- std::cout << "ProcessSplit: " << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
-
- int kBeam = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
-
-
- std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
- if (E < 8.5_GeV || Ecm < 10_GeV ) {
- std::cout << "ProcessSplit: " << " DoDiscrete: low en. particle, skipping.." << std::endl;
- // std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
- //fEnergy += p.GetEnergy();
- //p.Delete();
- //fCount++;
+ EnergyType Ecm = sqrt(Etot * Etot -
+ Ptot.squaredNorm() *
+ si::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 / si::constants::c);
+
+ std::cout << "ProcessSplit: "
+ << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "ProcessSplit: "
+ << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+
+ std::cout << "ProcessSplit: "
+ << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
+ << std::endl;
+ if (E < 8.5_GeV || Ecm < 10_GeV) {
+ std::cout << "ProcessSplit: "
+ << " DoDiscrete: low en. particle, skipping.." << std::endl;
} 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_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] * si::constants::c ) / ( gamma + 1.);
- pnorm += psib.GetEnergy();
-
- for(int j=0; j<3; ++j){
- p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c;
- // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
- // << " gb: " << gammaBetaComponents[j] << endl;
- en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
- }
- // const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
- // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
- // << " en lab (GeV): " << en_lab / 1_GeV << endl;
- // cout << " pz cm (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
- // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
-
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetEnergy( en_lab );
- pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
- //cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
-
- 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) );
- //cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
- //cout << "s_cm (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
- //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
- Ptot_final += pnew.GetMomentum();
- }
- //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
+ // 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_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] * si::constants::c) /
+ (gamma + 1.);
+ pnorm += psib.GetEnergy();
+
+ for (int j = 0; j < 3; ++j) {
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
+ gammaBetaComponents[j] /
+ si::constants::c;
+ en_lab -=
+ (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::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 *
+ // si::constants::c << endl;
}
- }
+ }
}
void Init() {
fCount = 0;
- // define reference frame? --> defines boosts between corsika stack and model stack.
-
- // initialize random numbers for sibyll
- // FOR NOW USE SIBYLL INTERNAL !!!
- // rnd_ini_();
-
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
;
const std::string str_name = "s_rndm";
rmng.RegisterRandomStream(str_name);
- // // corsika::random::RNG srng;
- // auto srng = rmng.GetRandomStream("s_rndm");
-
// test random number generator
std::cout << "ProcessSplit: "
<< " test sequence of random numbers." << std::endl;
@@ -514,9 +473,9 @@ public:
setTrackedParticlesStable();
}
-
int GetCount() { return fCount; }
EnergyType GetEnergy() { return fEnergy; }
+
private:
};
@@ -527,11 +486,10 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
-
int main() {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
+
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
@@ -545,20 +503,18 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
- MomentumType P0 = sqrt( E0*E0 - 0.93827_GeV * 0.93827_GeV ) / si::constants::c;
- auto plab = super_stupid::MomentumVector(rootCS,
- 0. * 1_GeV / si::constants::c,
- 0. * 1_GeV / si::constants::c,
- P0);
+ MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV) / si::constants::c;
+ auto plab = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c,
+ 0. * 1_GeV / si::constants::c, P0);
particle.SetEnergy(E0);
particle.SetMomentum(plab);
- particle.SetPID( Code::Proton );
+ particle.SetPID(Code::Proton);
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV << "GeV, particles below energy threshold =" << p1.GetCount() << endl;
+ cout << "Result: E0=" << E0 / 1_GeV
+ << "GeV, particles below energy threshold =" << p1.GetCount() << endl;
cout << "total energy below threshold (GeV): " << p1.GetEnergy() / 1_GeV << std::endl;
p3.ShowResults();
cout << "total energy (GeV): "
- << ( p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy() ) / 1_GeV
- << endl;
+ << (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;
}
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 2275cdf49..a25316553 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -6,8 +6,8 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/units/PhysicalUnits.h>
#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
using namespace std;
using namespace corsika::stack;
@@ -20,33 +20,33 @@ public:
void Init();
void Clear() { s_plist_.np = 0; }
-
- int GetSize() const { return s_plist_.np; }
-#warning check actual capacity of sibyll stack
+
+ int GetSize() const { return s_plist_.np; }
+#warning check actual capacity of sibyll stack
int GetCapacity() const { return 8000; }
-
- 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 SetMomentum(const int i, const super_stupid::MomentumVector& 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 SetMomentum(const int i, const super_stupid::MomentumVector& v) {
auto tmp = v.GetComponents();
- for(int idx=0; idx<3; ++idx)
- s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;
+ for (int idx = 0; idx < 3; ++idx)
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;
}
-
+
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; }
-
- super_stupid::MomentumVector GetMomentum(const int i) const
- {
+
+ super_stupid::MomentumVector GetMomentum(const int i) const {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
- corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
- super_stupid::MomentumVector v1(rootCS,components);
+ corsika::geometry::QuantityVector<momentum_d> components{
+ s_plist_.p[0][i] * 1_GeV / si::constants::c,
+ s_plist_.p[1][i] * 1_GeV / si::constants::c,
+ s_plist_.p[2][i] * 1_GeV / si::constants::c};
+ super_stupid::MomentumVector v1(rootCS, components);
return v1;
}
-
+
void Copy(const int i1, const int i2) {
s_plist_.llist[i1] = s_plist_.llist[i2];
s_plist_.p[3][i1] = s_plist_.p[3][i2];
@@ -65,13 +65,19 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetIndex;
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()); }
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); }
-
+ 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);
+ }
};
typedef Stack<SibStackData, ParticleInterface> SibStack;
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 8a9e9d6f6..e6198cffd 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -54,7 +54,7 @@ namespace corsika::units::si {
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
/**
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index e3597436f..31c05bff4 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -3,35 +3,36 @@
#include <corsika/process/ContinuousProcess.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/cascade/SibStack.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/setup/SetupTrajectory.h>
-//using namespace corsika::particles;
+// using namespace corsika::particles;
namespace corsika::process {
namespace sibyll {
-
- void setHadronsUnstable(){
+ 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;
+ 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
+ // set Leptons and Proton and Neutron stable
// use stack to loop over particles
setup::Stack ds;
ds.NewParticle().SetPID(corsika::particles::Code::Proton);
@@ -46,109 +47,109 @@ namespace corsika::process {
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();
- }
-
+ 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();
+ // 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;
- }
+ 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;
- int a = 1;
- const double x = -x0 * log(s_rndm_(a));
- cout << "ProcessDecay: next decay: " << x << endl;
- return x;
+ corsika::units::hep::EnergyType E = p.GetEnergy();
+ corsika::units::hep::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;
+
+ 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;
+ int a = 1;
+ const double x = -x0 * log(s_rndm_(a));
+ cout << "ProcessDecay: next decay: " << x << endl;
+ return x;
}
-
+
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() );
- // remove original particle from corsika stack
- p.Delete();
- // set all particles/hadrons unstable
- setHadronsUnstable();
- // call sibyll decay
- std::cout << "ProcessDecay: 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();
+ 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());
+ // remove original particle from corsika stack
+ p.Delete();
+ // set all particles/hadrons unstable
+ setHadronsUnstable();
+ // call sibyll decay
+ std::cout << "ProcessDecay: 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();
}
-
+
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
+ return EProcessReturn::eOk;
}
-
};
- }
-}
+ } // namespace sibyll
+} // namespace corsika::process
#endif
--
GitLab
From 1ab99a5f878f5c8c5e72304ea7b343e08a0c7ccc Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 11 Dec 2018 16:42:45 +0100
Subject: [PATCH 39/39] a few fixes
---
Documentation/Examples/cascade_example.cc | 9 ++++++---
Processes/Sibyll/ProcessDecay.h | 1 +
Processes/StackInspector/StackInspector.cc | 11 ++++++-----
Processes/StackInspector/StackInspector.h | 1 +
4 files changed, 14 insertions(+), 8 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f34f09449..e6e2fe483 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -125,7 +125,7 @@ public:
}
template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle& p, setup::Trajectory& t, Stack& s) const {
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
// cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
// cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
// cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
@@ -148,7 +148,7 @@ public:
}
template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const {
+ void DoDiscrete(Particle& p, Stack&) const {
cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
const Code pid = p.GetPID();
if (isEmParticle(pid)) {
@@ -499,7 +499,7 @@ int main() {
ProcessSplit p1;
corsika::process::sibyll::ProcessDecay p2;
ProcessEMCut p3;
- const auto sequence = p0 + p1 + p2 + p3;
+ const auto sequence = /*p0 +*/ p1 + p2 + p3;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
@@ -513,6 +513,9 @@ int main() {
particle.SetEnergy(E0);
particle.SetMomentum(plab);
particle.SetPID(Code::Proton);
+ particle.SetTime(0_ns);
+ Point p(rootCS, 0_m, 0_m, 0_m);
+ particle.SetPosition(p);
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 31c05bff4..e828ca9ae 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -95,6 +95,7 @@ namespace corsika::process {
const double gamma = E / m;
TimeType t0 = GetLifetime(p.GetPID());
+ cout << "ProcessDecay: code: " << (p.GetPID()) << endl;
cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
cout << "ProcessDecay: MinStep: density: " << density << endl;
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index b0acefc4a..949cfd025 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -26,7 +26,7 @@ using namespace corsika::process::stack_inspector;
template <typename Stack>
StackInspector<Stack>::StackInspector(const bool aReport)
- : fReport(aReport) {}
+ : fReport(aReport), fCountStep(0) {}
template <typename Stack>
StackInspector<Stack>::~StackInspector() {}
@@ -34,7 +34,6 @@ StackInspector<Stack>::~StackInspector() {}
template <typename Stack>
process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
Stack& s) const {
- static int countStep = 0;
if (!fReport) return EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
EnergyType Etot = 0_GeV;
@@ -49,8 +48,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
<< " pos=" << pos << endl;
}
- countStep++;
- cout << "StackInspector: nStep=" << countStep << " stackSize=" << s.GetSize()
+ fCountStep++;
+ cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()
<< " Estack=" << Etot / 1_GeV << " GeV" << endl;
return EProcessReturn::eOk;
}
@@ -61,7 +60,9 @@ void StackInspector<Stack>::MinStepLength(Particle&, setup::Trajectory&) const {
}
template <typename Stack>
-void StackInspector<Stack>::Init() {}
+void StackInspector<Stack>::Init() {
+ fCountStep = 0;
+}
#include <corsika/setup/SetupStack.h>
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 7b68d92d3..dddb7e4f6 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -40,6 +40,7 @@ namespace corsika::process {
private:
bool fReport;
+ mutable int fCountStep = 0;
};
} // namespace stack_inspector
--
GitLab