IAP GITLAB

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

added particle cut process

parent 87948022
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 @ 853e9aad
......@@ -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;
}
Framework/Cascade/Cascade.h
+ 2
− 0
View file @ 853e9aad
......@@ -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);
......
Processes/Sibyll/ProcessDecay.h
+ 9
− 6
View file @ 853e9aad
......@@ -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>
......
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