diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index c29cfdc902ad47f104edfa60750e98120c1efe24..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,12 +86,21 @@ 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 (assuming proton at rest as target) + /* + 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 ); // 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 / 1_GeV << std::endl; @@ -83,6 +113,8 @@ public: double sqs = Ecm / 1_GeV; // running sibyll sibyll_( kBeam, kTarget, sqs); + // running decays + //decsib_(); // print final state int print_unit = 6; sib_list_( print_unit ); @@ -96,7 +128,7 @@ public: // add particles from sibyll to stack for(int i=0; i<s_plist_.np; ++i){ - //transform to lab. frame + //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 ); @@ -108,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_(); @@ -127,6 +161,8 @@ public: //initialize Sibyll sibyll_ini_(); + + // set particles stable / unstable } int GetCount() { return fCount; }