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] 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 69b7e92f..7b4ab0ad 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 d38bf0da..c10ea7d8 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