IAP GITLAB

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

added particle cut process

parent 820f0123
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 260
− 76
View file @ a90b32d1
...@@ -39,6 +39,183 @@ using namespace corsika::setup; ...@@ -39,6 +39,183 @@ using namespace corsika::setup;
using namespace std; using namespace std;
static int fCount = 0; 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> { class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public: public:
...@@ -141,14 +318,14 @@ public: ...@@ -141,14 +318,14 @@ public:
next_step = -int_length * log(s_rndm_(a)); next_step = -int_length * log(s_rndm_(a));
}else }else
#warning define infinite interaction length? then we can skip the test in DoDiscrete() #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? what are the units of the output? slant depth or 3space length?
*/ */
std::cout << "ProcessSplit: " std::cout << "ProcessSplit: "
<< "next step (g/cm2): " << next_step << std::endl; << "next interaction (g/cm2): " << next_step << std::endl;
return next_step; return next_step;
} }
...@@ -160,7 +337,7 @@ public: ...@@ -160,7 +337,7 @@ public:
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const { 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() ) ){ if( process::sibyll::CanInteract( p.GetPID() ) ){
cout << "defining coordinates" << endl; cout << "defining coordinates" << endl;
// coordinate system, get global frame of reference // coordinate system, get global frame of reference
...@@ -231,81 +408,82 @@ public: ...@@ -231,81 +408,82 @@ public:
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 ) { if (E < 8.5_GeV || Ecm < 10_GeV ) {
std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl; std::cout << "ProcessSplit: " << " DoDiscrete: low en. particle, skipping.." << std::endl;
p.Delete(); // std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
fCount++; //fEnergy += p.GetEnergy();
} else { //p.Delete();
// Sibyll does not know about units.. //fCount++;
double sqs = Ecm / 1_GeV; } else {
// running sibyll, filling stack // Sibyll does not know about units..
sibyll_( kBeam, kTarget, sqs); double sqs = Ecm / 1_GeV;
// running decays // running sibyll, filling stack
setTrackedParticlesStable(); sibyll_( kBeam, kTarget, sqs);
decsib_(); // running decays
// print final state setTrackedParticlesStable();
int print_unit = 6; decsib_();
sib_list_( print_unit ); // print final state
int print_unit = 6;
// delete current particle sib_list_( print_unit );
p.Delete();
// delete current particle
// add particles from sibyll to stack p.Delete();
// link to sibyll stack
SibStack ss; // add particles from sibyll to stack
// link to sibyll stack
// SibStack does not know about momentum yet so we need counter to access momentum array in Sibyll SibStack ss;
int i = -1;
super_stupid::MomentumVector Ptot_final(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second}); // SibStack does not know about momentum yet so we need counter to access momentum array in Sibyll
for (auto &psib: ss){ int i = -1;
++i; super_stupid::MomentumVector Ptot_final(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
// skip particles that have decayed in Sibyll for (auto &psib: ss){
if( abs(s_plist_.llist[ i ]) > 100 ) continue; ++i;
// skip particles that have decayed in Sibyll
if( abs(s_plist_.llist[ i ]) > 100 ) continue;
//transform energy to lab. frame, primitve //transform energy to lab. frame, primitve
// compute beta_vec * p_vec // compute beta_vec * p_vec
// arbitrary Lorentz transformation based on sibyll routines // arbitrary Lorentz transformation based on sibyll routines
const auto gammaBetaComponents = gambet.GetComponents(); const auto gammaBetaComponents = gambet.GetComponents();
const auto pSibyllComponents = psib.GetMomentum().GetComponents(); const auto pSibyllComponents = psib.GetMomentum().GetComponents();
EnergyType en_lab = 0. * 1_GeV; EnergyType en_lab = 0. * 1_GeV;
MomentumType p_lab_components[3]; MomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma; en_lab = psib.GetEnergy() * gamma;
EnergyType pnorm = 0. * 1_GeV; EnergyType pnorm = 0. * 1_GeV;
for(int j=0; j<3; ++j) for(int j=0; j<3; ++j)
pnorm += ( pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c ) / ( gamma + 1.); pnorm += ( pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c ) / ( gamma + 1.);
pnorm += psib.GetEnergy(); pnorm += psib.GetEnergy();
for(int j=0; j<3; ++j){ for(int j=0; j<3; ++j){
p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] / si::constants::c; 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 // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
// << " gb: " << gammaBetaComponents[j] << endl; // << " gb: " << gammaBetaComponents[j] << endl;
en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c; en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
} }
// const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * 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 // cout << " en cm (GeV): " << psib.GetEnergy() / 1_GeV << endl
// << " en lab (GeV): " << en_lab / 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 // 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; // << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
// add to corsika stack // add to corsika stack
auto pnew = s.NewParticle(); auto pnew = s.NewParticle();
pnew.SetEnergy( en_lab ); pnew.SetEnergy( en_lab );
pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) ); pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
//cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl; //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], corsika::geometry::QuantityVector<momentum_d> p_lab_c{ p_lab_components[0],
p_lab_components[1], p_lab_components[1],
p_lab_components[2]}; p_lab_components[2]};
pnew.SetMomentum( super_stupid::MomentumVector( rootCS, p_lab_c) ); pnew.SetMomentum( super_stupid::MomentumVector( rootCS, p_lab_c) );
//cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl; //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_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; //cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
Ptot_final += pnew.GetMomentum(); 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() { void Init() {
...@@ -339,7 +517,7 @@ public: ...@@ -339,7 +517,7 @@ public:
int GetCount() { return fCount; } int GetCount() { return fCount; }
EnergyType GetEnergy() { return fEnergy; }
private: private:
}; };
...@@ -364,7 +542,8 @@ int main() { ...@@ -364,7 +542,8 @@ int main() {
ProcessSplit p1; ProcessSplit p1;
corsika::process::sibyll::ProcessDecay p2; corsika::process::sibyll::ProcessDecay p2;
const auto sequence = p0 + p1 + p2; ProcessEMCut p3;
const auto sequence = p0 + p1 + p2 + p3;
setup::Stack stack; setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack); corsika::cascade::Cascade EAS(tracking, sequence, stack);
...@@ -383,5 +562,10 @@ int main() { ...@@ -383,5 +562,10 @@ int main() {
EAS.Init(); EAS.Init();
EAS.Run(); 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;
} }
Framework/Cascade/Cascade.h
+ 2
− 0
View file @ a90b32d1
...@@ -60,6 +60,8 @@ namespace corsika::cascade { ...@@ -60,6 +60,8 @@ namespace corsika::cascade {
} }
void Step(Particle& particle) { void Step(Particle& particle) {
std::cout << "+++++++++++++++++++++++++++++++" << std::endl;
std::cout << "Cascade: starting step.." << std::endl;
corsika::setup::Trajectory step = fTracking.GetTrack(particle); corsika::setup::Trajectory step = fTracking.GetTrack(particle);
fProcesseList.MinStepLength(particle, step); fProcesseList.MinStepLength(particle, step);
......
Processes/Sibyll/ProcessDecay.h
+ 9
− 6
View file @ a90b32d1
...@@ -100,7 +100,10 @@ namespace corsika::process { ...@@ -100,7 +100,10 @@ namespace corsika::process {
// return as column density // return as column density
const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm; const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
cout << "ProcessDecay: MinStep: x0: " << x0 << endl; 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> template <typename Particle, typename Stack>
...@@ -112,14 +115,16 @@ namespace corsika::process { ...@@ -112,14 +115,16 @@ namespace corsika::process {
pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) ); pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) );
pin.SetEnergy( p.GetEnergy() ); pin.SetEnergy( p.GetEnergy() );
pin.SetMomentum( p.GetMomentum() ); pin.SetMomentum( p.GetMomentum() );
// remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable // set all particles/hadrons unstable
setHadronsUnstable(); setHadronsUnstable();
// call sibyll decay // call sibyll decay
std::cout << "calling Sibyll decay routine.." << std::endl; std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
decsib_(); decsib_();
// print output // print output
//int print_unit = 6; int print_unit = 6;
//sib_list_( print_unit ); sib_list_( print_unit );
// copy particles from sibyll stack to corsika // copy particles from sibyll stack to corsika
int i = -1; int i = -1;
for (auto &psib: ss){ for (auto &psib: ss){
...@@ -135,8 +140,6 @@ namespace corsika::process { ...@@ -135,8 +140,6 @@ namespace corsika::process {
} }
// empty sibyll stack // empty sibyll stack
ss.Clear(); ss.Clear();
// remove original particle from stack
p.Delete();
} }
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
......
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