diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f7a9687233e33e6aef241e6420a772398573fac8..4132fae2f6f0d4885e976cdb65733595e9597ca9 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -42,88 +42,81 @@ 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;
+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
- {
+ 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 )
+ 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 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;
+ 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
+ void defineEmParticles() const {
+ // create bool array identifying em particles
}
- bool isInvisible( Code pCode ) const
- {
+ 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;
+ 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
- {
+ double MinStepLength(Particle& p, setup::Trajectory&) const {
const Code pid = p.GetPID();
- if( isEmParticle( pid ) || isInvisible( pid ) ){
+ if (isEmParticle(pid) || isInvisible(pid)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
return 0.;
- } else {
+ } else {
double next_step = std::numeric_limits<double>::infinity();
cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
return next_step;
@@ -131,7 +124,7 @@ public:
}
template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&p, setup::Trajectory&t, Stack&s) const {
+ 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;
@@ -150,68 +143,55 @@ public:
// p.Delete();
// return EProcessReturn::eParticleAbsorbed;
// }
- return EProcessReturn::eOk;
+ return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const
- {
+ void DoDiscrete(Particle& p, Stack& s) const {
cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
const Code pid = p.GetPID();
- if( isEmParticle( pid ) ){
+ if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += p.GetEnergy();
- fEmCount += 1;
+ fEmCount += 1;
p.Delete();
- }
- else if ( isInvisible( pid ) ){
+ } else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += p.GetEnergy();
- fInvCount += 1;
+ fInvCount += 1;
p.Delete();
- }
- else if( isBelowEnergyCut( p ) ){
+ } else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy();
- fCount += 1;
+ fCount += 1;
p.Delete();
}
}
- void Init()
- {
- fEmEnergy = 0. * 1_GeV;
- fEmCount = 0;
+ void Init() {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
fInvEnergy = 0. * 1_GeV;
- fInvCount = 0;
- fEnergy = 0. * 1_GeV;
- //defineEmParticles();
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ // defineEmParticles();
}
- void ShowResults(){
+ 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;
+ << " 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 GetInvEnergy() { return fInvEnergy; }
- EnergyType GetCutEnergy()
- {
- return fEnergy;
- }
+ EnergyType GetCutEnergy() { return fEnergy; }
- EnergyType GetEmEnergy()
- {
- return fEmEnergy;
- }
+ EnergyType GetEmEnergy() { return fEmEnergy; }
private:
};
@@ -225,7 +205,7 @@ public:
Sibyll is hadronic generator
only hadrons decay
*/
- // set particles unstable
+ // set particles unstable
corsika::process::sibyll::setHadronsUnstable();
// make tracked particles stable
std::cout << "ProcessSplit: setting tracked hadrons stable.." << std::endl;
@@ -236,34 +216,30 @@ public:
ds.NewParticle().SetPID(Code::KMinus);
ds.NewParticle().SetPID(Code::K0Long);
ds.NewParticle().SetPID(Code::K0Short);
-
+
for (auto& p : ds) {
int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
// set particle stable by setting table value negative
- // cout << "ProcessSplit: doDiscrete: setting " << p.GetPID() << "(" << s_id << ")"
- // << " stable in Sibyll .." << endl;
- s_csydec_.idb[s_id - 1] = (-1) * abs( s_csydec_.idb[s_id - 1] );
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
p.Delete();
}
-
}
-
-
+
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
// coordinate system, get global frame of reference
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
+
const Code corsikaBeamId = p.GetPID();
-
+
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
- int kBeam = process::sibyll::GetSibyllXSCode( corsikaBeamId );
+ int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
- bool kInteraction = process::sibyll::CanInteract( 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..
@@ -272,53 +248,56 @@ public:
// FOR NOW: assume target is oxygen
int kTarget = 16;
- // proton mass in units of energy
- const EnergyType proton_mass_en = 0.93827_GeV ;
-
- EnergyType Etot = p.GetEnergy() + kTarget * proton_mass_en;
- super_stupid::MomentumVector Ptot(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+ EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+ super_stupid::MomentumVector Ptot(
+ rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
// FOR NOW: assume target is at rest
- super_stupid::MomentumVector pTarget(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+ super_stupid::MomentumVector pTarget(
+ rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
Ptot += p.GetMomentum();
Ptot += pTarget;
// calculate cm. energy
- EnergyType sqs = sqrt( Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared );
+ EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared);
double Ecm = sqs / 1_GeV;
-
- std::cout << "ProcessSplit: " << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
- << " beam can interact:" << kBeam<< endl
- << " beam XS code:" << kBeam << endl
- << " beam pid:" << p.GetPID() << endl
- << " target mass number:" << kTarget << std::endl;
+
+ std::cout << "ProcessSplit: "
+ << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
+ << " beam can interact:" << kBeam << endl
+ << " beam XS code:" << kBeam << endl
+ << " beam pid:" << p.GetPID() << endl
+ << " target mass number:" << kTarget << std::endl;
double next_step;
- if(kInteraction){
-
- double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
+ if (kInteraction) {
+
+ double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
-
- if(kTarget==1)
- sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+
+ if (kTarget == 1)
+ sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
else
- sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy );
-
- std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ sib_sigma_hnuc_(kBeam, kTarget, Ecm, 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;
+ 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;
+ 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 );
+ 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;
+ 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()
+ } else
next_step = std::numeric_limits<double>::infinity();
-
+
/*
what are the units of the output? slant depth or 3space length?
@@ -336,172 +315,152 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- cout << "ProcessSplit: " << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
- if( process::sibyll::CanInteract( p.GetPID() ) ){
+ 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
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
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)
+
+ /*
+ 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 ;
+ int kTarget = 1; // env.GetTargetParticle().GetPID();
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;
-
+ const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
+ 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;
+ 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)
*/
- // 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;
+ EnergyType E = p.GetEnergy();
+ EnergyType Etot = E + Etarget;
// 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: low en. particle, skipping.." << std::endl;
- // std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
- //fEnergy += p.GetEnergy();
- //p.Delete();
- //fCount++;
+ EnergyType Ecm = sqrt(Etot * Etot -
+ Ptot.squaredNorm() *
+ si::constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
+ /*
+ 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 = 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: low en. particle, skipping.." << std::endl;
} 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();
-
- 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();
- }
- //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
+ // 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();
+
+ for (int j = 0; j < 3; ++j) {
+ p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
+ gammaBetaComponents[j] /
+ si::constants::c;
+ en_lab -=
+ (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+ }
+
+ // add to corsika stack
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy(en_lab);
+ pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+
+ 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));
+ Ptot_final += pnew.GetMomentum();
+ }
+ // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV *
+ // si::constants::c << endl;
}
- }
+ }
}
void Init() {
fCount = 0;
- // define reference frame? --> defines boosts between corsika stack and model stack.
-
- // initialize random numbers for sibyll
- // FOR NOW USE SIBYLL INTERNAL !!!
- // rnd_ini_();
-
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
;
const std::string str_name = "s_rndm";
rmng.RegisterRandomStream(str_name);
- // // corsika::random::RNG srng;
- // auto srng = rmng.GetRandomStream("s_rndm");
-
// test random number generator
std::cout << "ProcessSplit: "
<< " test sequence of random numbers." << std::endl;
@@ -514,9 +473,9 @@ public:
setTrackedParticlesStable();
}
-
int GetCount() { return fCount; }
EnergyType GetEnergy() { return fEnergy; }
+
private:
};
@@ -527,11 +486,10 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
-
int main() {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
+
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
@@ -545,20 +503,18 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_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);
+ 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 );
+ particle.SetPID(Code::Proton);
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV << "GeV, particles below energy threshold =" << 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;
+ << (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;
}
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 2275cdf491ae34433f8afba51041d867d301b8e2..a25316553da2182f6da9b2668f7f67336ede9a84 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -6,8 +6,8 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/units/PhysicalUnits.h>
#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
using namespace std;
using namespace corsika::stack;
@@ -20,33 +20,33 @@ public:
void Init();
void Clear() { s_plist_.np = 0; }
-
- int GetSize() const { return s_plist_.np; }
-#warning check actual capacity of sibyll stack
+
+ 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 EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
- void SetMomentum(const int i, const super_stupid::MomentumVector& 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;
+ 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]; }
-
+
EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
-
- super_stupid::MomentumVector GetMomentum(const int i) const
- {
+
+ super_stupid::MomentumVector GetMomentum(const int i) const {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
- 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);
+ 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];
@@ -65,13 +65,19 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v ); }
+ void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
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()); }
- void SetMomentum(const super_stupid::MomentumVector& v) { GetStackData().SetMomentum(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());
+ }
+ void SetMomentum(const super_stupid::MomentumVector& v) {
+ GetStackData().SetMomentum(GetIndex(), v);
+ }
};
typedef Stack<SibStackData, ParticleInterface> SibStack;
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 8a9e9d6f63db1e8bd55b96e461273a9d5acd9955..e6198cffd0bbdf3c443aa79b0874a82f11d318b8 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -54,7 +54,7 @@ namespace corsika::units::si {
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
using MomentumType = phys::units::quantity<momentum_d, double>;
-
+
} // end namespace corsika::units::si
/**
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index e3597436f31f5a6e8270ee764bb72a2debde2b33..31c05bff4934049055d3dba5394f6c3183a43190 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -3,35 +3,36 @@
#include <corsika/process/ContinuousProcess.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/cascade/SibStack.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/setup/SetupTrajectory.h>
-//using namespace corsika::particles;
+// using namespace corsika::particles;
namespace corsika::process {
namespace sibyll {
-
- void setHadronsUnstable(){
+ void setHadronsUnstable() {
// name? also makes EM particles stable
-
+
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. corsika::particles::ListOfParticles()
std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
// make ALL particles unstable, then set EM stable
- for(auto &p: corsika2sibyll){
- //std::cout << (int)p << std::endl;
- const int sibCode = (int)p;
- // skip unknown and antiparticles
- if( sibCode< 1) continue;
- //std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll( static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
- s_csydec_.idb[ sibCode - 1 ] = abs( s_csydec_.idb[ sibCode - 1 ] );
- //std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] << std::endl;
+ for (auto& p : corsika2sibyll) {
+ // std::cout << (int)p << std::endl;
+ const int sibCode = (int)p;
+ // skip unknown and antiparticles
+ if (sibCode < 1) continue;
+ // std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
+ // static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
+ s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
+ // std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
+ // std::endl;
}
- // set Leptons and Proton and Neutron stable
+ // set Leptons and Proton and Neutron stable
// use stack to loop over particles
setup::Stack ds;
ds.NewParticle().SetPID(corsika::particles::Code::Proton);
@@ -46,109 +47,109 @@ namespace corsika::process {
ds.NewParticle().SetPID(corsika::particles::Code::NuMuBar);
for (auto& p : ds) {
- int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- // set particle stable by setting table value negative
- // cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
- // << " stable in Sibyll .." << endl;
- s_csydec_.idb[ s_id - 1 ] = (-1) * abs(s_csydec_.idb[ s_id - 1 ]);
- p.Delete();
- }
-
+ int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // set particle stable by setting table value negative
+ // cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
+ // << " stable in Sibyll .." << endl;
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ p.Delete();
+ }
}
-
class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
public:
ProcessDecay() {}
void Init() {
- //setHadronsUnstable();
+ // setHadronsUnstable();
}
- void setAllStable(){
- // name? also makes EM particles stable
-
- // loop over all particles in sibyll
- // should be changed to loop over human readable list
- // i.e. corsika::particles::ListOfParticles()
- for(auto &p: corsika2sibyll){
- //std::cout << (int)p << std::endl;
- const int sibCode = (int)p;
- // skip unknown and antiparticles
- if( sibCode< 1) continue;
- std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll( static_cast<SibyllCode> ( sibCode ) ) << " stable" << std::endl;
- s_csydec_.idb[ sibCode - 1 ] = -1 * abs( s_csydec_.idb[ sibCode - 1 ] );
- std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] << std::endl;
- }
+ void setAllStable() {
+ // name? also makes EM particles stable
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. corsika::particles::ListOfParticles()
+ for (auto& p : corsika2sibyll) {
+ // std::cout << (int)p << std::endl;
+ const int sibCode = (int)p;
+ // skip unknown and antiparticles
+ if (sibCode < 1) continue;
+ std::cout << "Sibyll: Decay: setting "
+ << ConvertFromSibyll(static_cast<SibyllCode>(sibCode)) << " stable"
+ << std::endl;
+ s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
+ std::cout << "decay table value: " << s_csydec_.idb[sibCode - 1] << std::endl;
+ }
}
friend void setHadronsUnstable();
-
+
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
- EnergyType E = p.GetEnergy();
- MassType m = corsika::particles::GetMass(p.GetPID());
- // env.GetDensity();
-
- const MassDensityType density = 1.25e-3 * kilogram / ( 1_cm * 1_cm * 1_cm );
-
- const double gamma = E / m / constants::cSquared;
-
- TimeType t0 = GetLifetime( p.GetPID() );
- cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
- cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
- cout << "ProcessDecay: MinStep: density: " << density << endl;
- // return as column density
- const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
- cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
- int a = 1;
- const double x = -x0 * log(s_rndm_(a));
- cout << "ProcessDecay: next decay: " << x << endl;
- return x;
+ corsika::units::hep::EnergyType E = p.GetEnergy();
+ corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
+ // env.GetDensity();
+
+ const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
+
+ const double gamma = E / m;
+
+ TimeType t0 = GetLifetime(p.GetPID());
+ cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
+ cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
+ cout << "ProcessDecay: MinStep: density: " << density << endl;
+ // return as column density
+ const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
+ cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
+ int a = 1;
+ const double x = -x0 * log(s_rndm_(a));
+ cout << "ProcessDecay: next decay: " << x << endl;
+ return x;
}
-
+
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
- SibStack ss;
- ss.Clear();
- // copy particle to sibyll stack
- auto pin = ss.NewParticle();
- 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 << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
- decsib_();
- // print output
- int print_unit = 6;
- sib_list_( print_unit );
- // copy particles from sibyll stack to corsika
- int i = -1;
- for (auto &psib: ss){
- ++i;
- // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
- if( abs(s_plist_.llist[ i ]) > 100 ) continue;
- // add to corsika stack
- //cout << "decay product: " << process::sibyll::ConvertFromSibyll( psib.GetPID() ) << endl;
- auto pnew = s.NewParticle();
- pnew.SetEnergy( psib.GetEnergy() );
- pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
- pnew.SetMomentum( psib.GetMomentum() );
- }
- // empty sibyll stack
- ss.Clear();
+ SibStack ss;
+ ss.Clear();
+ // copy particle to sibyll stack
+ auto pin = ss.NewParticle();
+ 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 << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
+ decsib_();
+ // print output
+ int print_unit = 6;
+ sib_list_(print_unit);
+ // copy particles from sibyll stack to corsika
+ int i = -1;
+ for (auto& psib : ss) {
+ ++i;
+ // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
+ if (abs(s_plist_.llist[i]) > 100) continue;
+ // add to corsika stack
+ // cout << "decay product: " << process::sibyll::ConvertFromSibyll(
+ // psib.GetPID() ) << endl;
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy(psib.GetEnergy());
+ pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+ pnew.SetMomentum(psib.GetMomentum());
+ }
+ // empty sibyll stack
+ ss.Clear();
}
-
+
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
+ return EProcessReturn::eOk;
}
-
};
- }
-}
+ } // namespace sibyll
+} // namespace corsika::process
#endif