diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 08a55e7621e8a00481439180beb806be178ddb6c..c4d7a7e8ee24daf5cec4e1ca29941ce26d755025 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -20,6 +20,8 @@
#include <corsika/random/RNGManager.h>
#include <corsika/cascade/sibyll2.3c.h>
+//#include <corsika/units/PhysicalConstants.h>
+#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
@@ -39,21 +41,40 @@ public:
double MinStepLength(Particle& p) const {
// beam particles for sibyll : 1, 2, 3 for p, pi, k
int kBeam = 1;
+
+ /*
+ 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
- int kTarget = 0.4*16 + 0.6*14;
+ // FOR NOW: assume target is oxygen
+ int kTarget = 16;
double beamEnergy = p.GetEnergy() / 1_GeV;
std::cout << "ProcessSplit: " << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
double prodCrossSection,dummy;
+
sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+
std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
CrossSectionType sig = prodCrossSection / 1000. * barn;
- std::cout << "ProcessSplit: " << "MinStep: CrossSection= " << sig << 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;
// calculate interaction length in medium
-
+ double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
// pick random step lenth
-
- return prodCrossSection;
+ 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?
+
+ */
+ std::cout << "ProcessSplit: " << "next step (g/cm2): " << next_step << std::endl;
+ return next_step;
}
template <typename Particle, typename Trajectory, typename Stack>
@@ -65,35 +86,52 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
// get energy of particle from stack
- // stack is in GeV in lab. frame
- // convert to GeV in cm. frame
+ /*
+ 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();
- double Ecm = sqrt( 2. * E * 0.93827_GeV ) / 1_GeV ;
+ EnergyType Ecm = sqrt( 2. * E * 0.93827_GeV );
// FOR NOW: set beam to proton
int kBeam = 13; //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 << std::endl;
- if (E < 8.5_GeV || Ecm < 10. ) {
+ 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 {
- //p.SetEnergy(E / 2);
- //s.NewParticle().SetEnergy(E / 2);
-
+ // Sibyll does not know about units..
+ double sqs = Ecm / 1_GeV;
// running sibyll
- sibyll_( kBeam, kTarget, Ecm);
+ sibyll_( kBeam, kTarget, sqs);
+ // running decays
+ //decsib_();
+ // print final state
int print_unit = 6;
sib_list_( print_unit );
// delete current particle
p.Delete();
+
+ 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;
// add particles from sibyll to stack
for(int i=0; i<s_plist_.np; ++i){
-
- s.NewParticle().SetEnergy( s_plist_.p[3][i] * 1_GeV );
+ //transform to lab. frame, primitve
+ const double en_lab = gambet * s_plist_.p[2][i] + gamma * s_plist_.p[3][i];
+ // add to corsika stack
+ s.NewParticle().SetEnergy( en_lab * 1_GeV );
}
}
}
@@ -102,6 +140,8 @@ public:
{
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_();
@@ -121,6 +161,8 @@ public:
//initialize Sibyll
sibyll_ini_();
+
+ // set particles stable / unstable
}
int GetCount() { return fCount; }
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 915f2d59cf7f418255574218f1530b09f1ff37e5..c3d5824e8a1c82fba366b9409add9b543cad8b60 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -54,7 +54,7 @@ namespace corsika::units::si::constants {
// unified atomic mass unit
constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
- // millibarn
+ // barn
constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
// etc.
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index adaa004d30592edd91bf0a2ae154d8906e1841ff..9060037e37486045a8bae8fd133a875ed521e108 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -24,6 +24,31 @@ namespace phys {
} // namespace units
} // namespace phys
+namespace phys {
+ namespace units {
+ namespace literals {
+ QUANTITY_DEFINE_SCALING_LITERALS(barn, area_d,
+ magnitude(corsika::units::si::constants::barn))
+
+ // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+ } // namespace literals
+ } // namespace units
+}
+
+namespace phys {
+ namespace units {
+ namespace literals {
+ QUANTITY_DEFINE_SCALING_LITERALS(meter, length_d,
+ magnitude(corsika::units::si::constants::meter))
+
+ // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+ } // namespace literals
+ } // namespace units
+}
+
+
namespace corsika::units::si {
using namespace phys::units;
using namespace phys::units::literals;
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 3528cf67cd85cd07d3b38d9f9e23cd6c973b9057..fc844c10860c006466b73d6b2619ece646bb139b 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -54,7 +54,8 @@ TEST_CASE("PhysicalUnits", "[Units]") {
REQUIRE(1_mol / 1_amol == Approx(1e18));
REQUIRE(1_K / 1_zK == Approx(1e21));
REQUIRE(1_K / 1_yK == Approx(1e24));
-
+ // REQUIRE(1_barn / 1_mbarn == Approx(1e3));
+
REQUIRE(1_A / 1_hA == Approx(1e-2));
REQUIRE(1_m / 1_km == Approx(1e-3));
REQUIRE(1_m / 1_Mm == Approx(1e-6));
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 5717953f93296478c8692e3903e0c0b91b5f53e5..b475e9c359437164eacd0782047f53bbd84b7eb2 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -38,14 +38,15 @@ process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle
// std::cout << "generation " << countStep << std::endl;
[[maybe_unused]] int i = 0;
EnergyType Etot = 0_GeV;
+
for (auto& iterP : s) {
EnergyType E = iterP.GetEnergy();
Etot += E;
- // std::cout << " particle data: " << i++ << ", id=" << iterP << " | " << std::endl;
+ std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()<< ", en=" << iterP.GetEnergy() / 1_GeV << " | " << std::endl;
}
countStep++;
// cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl;
- cout << countStep << " " << s.GetSize() << " " << Etot / 1_GeV << " " << endl;
+ cout << "call: " << countStep << " stack size: " << s.GetSize() << " energy: " << Etot / 1_GeV << endl;
return EProcessReturn::eOk;
}