diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 34b55476d9a4f869bd9a29c51688d0f49e680c37..0df3f0d6499228b06865019a13904d26da59ade8 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 7e9bd3bac3bc743bea8fcb0c2b0d1d401e16efae..57a3f9c16b099c6d00e71b405223f6c3f7eb41c5 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;