From 853e9aad4a22fb6afa4ab41a63fea4ea935f49d3 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sat, 8 Dec 2018 13:56:34 +0000
Subject: [PATCH] added particle cut process
---
Documentation/Examples/cascade_example.cc | 336 +++++++++++++++++-----
Framework/Cascade/Cascade.h | 2 +
Processes/Sibyll/ProcessDecay.h | 15 +-
3 files changed, 271 insertions(+), 82 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 5c631f5a0..f7a968723 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -38,6 +38,183 @@ using namespace corsika::random;
using namespace std;
static int fCount = 0;
+static EnergyType fEnergy = 0. * 1_GeV;
+
+// FOR NOW: global static variables for ParticleCut process
+// this is just wrong...
+static EnergyType fEmEnergy;
+static int fEmCount;
+
+static EnergyType fInvEnergy;
+static int fInvCount;
+
+
+class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
+public:
+ ProcessEMCut() {}
+ template <typename Particle>
+ bool isBelowEnergyCut( Particle &p ) const
+ {
+ // FOR NOW: center-of-mass energy hard coded
+ const EnergyType Ecm = sqrt( 2. * p.GetEnergy() * 0.93827_GeV );
+ if( p.GetEnergy() < 50_GeV || Ecm < 10_GeV )
+ return true;
+ else
+ return false;
+ }
+
+ bool isEmParticle( Code pCode ) const
+ {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch( pCode ){
+ case Code::Electron :
+ is_em = true;
+ break;
+ case Code::Gamma :
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ void defineEmParticles() const
+ {
+ // create bool array identifying em particles
+ }
+
+ bool isInvisible( Code pCode ) const
+ {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch( pCode ){
+ case Code::NuE :
+ is_inv = true;
+ break;
+ case Code::NuEBar :
+ is_inv = true;
+ break;
+ case Code::NuMu :
+ is_inv = true;
+ break;
+ case Code::NuMuBar :
+ is_inv = true;
+ break;
+ case Code::MuPlus :
+ is_inv = true;
+ break;
+ case Code::MuMinus :
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+
+ template <typename Particle>
+ double MinStepLength(Particle& p, setup::Trajectory&) const
+ {
+ const Code pid = p.GetPID();
+ if( isEmParticle( pid ) || isInvisible( pid ) ){
+ cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+ return 0.;
+ } else {
+ double next_step = std::numeric_limits<double>::infinity();
+ cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+ return next_step;
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&p, setup::Trajectory&t, Stack&s) const {
+ // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
+ // cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
+ // cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
+ // const Code pid = p.GetPID();
+ // if( isEmParticle( pid ) ){
+ // cout << "removing em. particle..." << endl;
+ // fEmEnergy += p.GetEnergy();
+ // fEmCount += 1;
+ // p.Delete();
+ // return EProcessReturn::eParticleAbsorbed;
+ // }
+ // if ( isInvisible( pid ) ){
+ // cout << "removing inv. particle..." << endl;
+ // fInvEnergy += p.GetEnergy();
+ // fInvCount += 1;
+ // p.Delete();
+ // return EProcessReturn::eParticleAbsorbed;
+ // }
+ return EProcessReturn::eOk;
+ }
+
+ template <typename Particle, typename Stack>
+ void DoDiscrete(Particle& p, Stack& s) const
+ {
+ cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
+ const Code pid = p.GetPID();
+ if( isEmParticle( pid ) ){
+ cout << "removing em. particle..." << endl;
+ fEmEnergy += p.GetEnergy();
+ fEmCount += 1;
+ p.Delete();
+ }
+ else if ( isInvisible( pid ) ){
+ cout << "removing inv. particle..." << endl;
+ fInvEnergy += p.GetEnergy();
+ fInvCount += 1;
+ p.Delete();
+ }
+ else if( isBelowEnergyCut( p ) ){
+ cout << "removing low en. particle..." << endl;
+ fEnergy += p.GetEnergy();
+ fCount += 1;
+ p.Delete();
+ }
+ }
+
+ void Init()
+ {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ //defineEmParticles();
+ }
+
+ void ShowResults(){
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " ******************************" << endl;
+ }
+
+ EnergyType GetInvEnergy()
+ {
+ return fInvEnergy;
+ }
+
+ EnergyType GetCutEnergy()
+ {
+ return fEnergy;
+ }
+
+ EnergyType GetEmEnergy()
+ {
+ return fEmEnergy;
+ }
+
+private:
+};
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
@@ -140,14 +317,14 @@ public:
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;
+ next_step = std::numeric_limits<double>::infinity();
/*
what are the units of the output? slant depth or 3space length?
*/
std::cout << "ProcessSplit: "
- << "next step (g/cm2): " << next_step << std::endl;
+ << "next interaction (g/cm2): " << next_step << std::endl;
return next_step;
}
@@ -159,7 +336,7 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
+ cout << "ProcessSplit: " << "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
@@ -230,81 +407,82 @@ public:
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
- setTrackedParticlesStable();
- 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;
-
- // 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;
- // skip particles that have decayed in Sibyll
- if( abs(s_plist_.llist[ i ]) > 100 ) continue;
+ if (E < 8.5_GeV || Ecm < 10_GeV ) {
+ std::cout << "ProcessSplit: " << " DoDiscrete: low en. particle, skipping.." << std::endl;
+ // std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
+ //fEnergy += p.GetEnergy();
+ //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
+ setTrackedParticlesStable();
+ 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;
+
+ // 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;
+ // skip particles that have decayed in Sibyll
+ if( abs(s_plist_.llist[ i ]) > 100 ) continue;
- //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();
+ //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;
+ 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;
+ // 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();
+ 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;
}
- //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
- }
- }else
- p.Delete();
+ }
}
void Init() {
@@ -338,7 +516,7 @@ public:
int GetCount() { return fCount; }
-
+ EnergyType GetEnergy() { return fEnergy; }
private:
};
@@ -358,7 +536,8 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
corsika::process::sibyll::ProcessDecay p2;
- const auto sequence = p0 + p1 + p2;
+ ProcessEMCut p3;
+ const auto sequence = p0 + p1 + p2 + p3;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
@@ -376,5 +555,10 @@ int main() {
particle.SetPID( Code::Proton );
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
+ cout << "Result: E0=" << E0 / 1_GeV << "GeV, particles below energy threshold =" << p1.GetCount() << endl;
+ cout << "total energy below threshold (GeV): " << p1.GetEnergy() / 1_GeV << std::endl;
+ p3.ShowResults();
+ cout << "total energy (GeV): "
+ << ( p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy() ) / 1_GeV
+ << endl;
}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index e46960dc9..c166517fa 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -60,6 +60,8 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
+ std::cout << "+++++++++++++++++++++++++++++++" << std::endl;
+ std::cout << "Cascade: starting step.." << std::endl;
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
fProcesseList.MinStepLength(particle, step);
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 449a281b4..e3597436f 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -100,7 +100,10 @@ namespace corsika::process {
// return as column density
const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
- return x0;
+ int a = 1;
+ const double x = -x0 * log(s_rndm_(a));
+ cout << "ProcessDecay: next decay: " << x << endl;
+ return x;
}
template <typename Particle, typename Stack>
@@ -112,14 +115,16 @@ namespace corsika::process {
pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) );
pin.SetEnergy( p.GetEnergy() );
pin.SetMomentum( p.GetMomentum() );
+ // remove original particle from corsika stack
+ p.Delete();
// set all particles/hadrons unstable
setHadronsUnstable();
// call sibyll decay
- std::cout << "calling Sibyll decay routine.." << std::endl;
+ std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
decsib_();
// print output
- //int print_unit = 6;
- //sib_list_( print_unit );
+ int print_unit = 6;
+ sib_list_( print_unit );
// copy particles from sibyll stack to corsika
int i = -1;
for (auto &psib: ss){
@@ -135,8 +140,6 @@ namespace corsika::process {
}
// empty sibyll stack
ss.Clear();
- // remove original particle from stack
- p.Delete();
}
template <typename Particle, typename Stack>
--
GitLab