IAP GITLAB

Skip to content
Snippets Groups Projects
Commit b69eef7f authored by Felix Riehn's avatar Felix Riehn
Browse files

added momentum to sibyll stack, implemented boost between sibyll stack and corsika stack

parent 1ea5333d
No related branches found
No related tags found
1 merge request!28Sibyll
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 198
− 103
View file @ b69eef7f
......@@ -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;
......
Framework/Cascade/SibStack.h
+ 30
− 12
View file @ b69eef7f
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum