diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 1b428066812c63c2a1cad080e0010446fb69ac08..6dda883951276c029cd0413351ff0a25c0075097 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -218,8 +218,9 @@ int main() {
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
corsika::environment::NuclearComposition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Nitrogen,corsika::particles::Code::Oxygen},
- std::vector<float>{(float)1.-fox, fox}));
+ std::vector<corsika::particles::Code>{corsika::particles::Code::Nitrogen,
+ corsika::particles::Code::Oxygen},
+ std::vector<float>{(float)1. - fox, fox}));
universe.AddChild(std::move(theMedium));
@@ -272,8 +273,8 @@ int main() {
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
- const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
- cout << "total energy (GeV): "
- << Efinal / 1_GeV << endl
- << "relative difference (%): " << (Efinal / E0 - 1. ) * 100 << endl;
+ const HEPEnergyType Efinal =
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cout << "total energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
}
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index 4c9a03f366a52c68b96e3ffe90bca8421b8f7810..9344d200922ad54a8584110c2b411b58fed14990 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -16,8 +16,8 @@
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
/**
* a homogeneous medium
@@ -42,29 +42,30 @@ namespace corsika::environment {
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
- corsika::particles::Code const& GetTarget( std::vector<corsika::units::si::CrossSectionType> &sigma) const override {
+ corsika::particles::Code const& GetTarget(
+ std::vector<corsika::units::si::CrossSectionType>& sigma) const override {
using namespace corsika::units::si;
- int i=-1;
+ int i = -1;
corsika::units::si::CrossSectionType total_weighted_sigma = 0._mbarn;
- std::vector<float> fractions;
- for(auto w: fNuclComp.GetFractions() ){
- i++;
- std::cout << "HomogeneousMedium: fraction: " << w << std::endl;
- total_weighted_sigma += w * sigma[i];
- fractions.push_back( w * sigma[i] / 1_mbarn );
+ std::vector<float> fractions;
+ for (auto w : fNuclComp.GetFractions()) {
+ i++;
+ std::cout << "HomogeneousMedium: fraction: " << w << std::endl;
+ total_weighted_sigma += w * sigma[i];
+ fractions.push_back(w * sigma[i] / 1_mbarn);
}
-
- for(auto f: fractions){
- f = f / ( total_weighted_sigma / 1_mbarn );
- std::cout << "HomogeneousMedium: reweighted fraction: " << f << std::endl;
+
+ for (auto f : fractions) {
+ f = f / (total_weighted_sigma / 1_mbarn);
+ std::cout << "HomogeneousMedium: reweighted fraction: " << f << std::endl;
}
- std::discrete_distribution channelDist( fractions.begin(), fractions.end() );
+ std::discrete_distribution channelDist(fractions.begin(), fractions.end());
static corsika::random::RNG& kRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
const int ichannel = channelDist(kRNG);
return fNuclComp.GetComponents()[ichannel];
}
-
+
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::LengthType pTo) const override {
@@ -76,7 +77,7 @@ namespace corsika::environment {
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
- }
+ }
};
} // namespace corsika::environment
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index 24393e805eb46a90febe4ac568c2949334f88bfe..3dabc3f983e187c0d88371262e38abc14720e02c 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -39,7 +39,8 @@ namespace corsika::environment {
virtual NuclearComposition const& GetNuclearComposition() const = 0;
- virtual corsika::particles::Code const& GetTarget( std::vector<corsika::units::si::CrossSectionType> &) const = 0;
+ virtual corsika::particles::Code const& GetTarget(
+ std::vector<corsika::units::si::CrossSectionType>&) const = 0;
};
} // namespace corsika::environment
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 632c57848fc734861cd2dab29325a0c0961206ff..6af26b0abdbc9cac073d1c833e085c4c8295ce78 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -138,7 +138,8 @@ namespace corsika::process {
template <typename Particle, typename Stack>
EProcessReturn SelectInteraction(
- Particle& p, Stack& s, [[maybe_unused]]corsika::units::si::InverseGrammageType lambda_select,
+ Particle& p, Stack& s,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if constexpr (is_process_sequence<T1type>::value) {
@@ -204,9 +205,10 @@ namespace corsika::process {
// select decay process
template <typename Particle, typename Stack>
- EProcessReturn SelectDecay(Particle& p, Stack& s,
- [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
- corsika::units::si::InverseTimeType& decay_inv_count) {
+ EProcessReturn SelectDecay(
+ Particle& p, Stack& s,
+ [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
+ corsika::units::si::InverseTimeType& decay_inv_count) {
if constexpr (is_process_sequence<T1>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret = A.SelectDecay(p, s, decay_select, decay_inv_count);
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index a798f318db53e994ba7a1edff828e9010dec979d..e7ed3165dc46bd50e193298f8a63956becf028f4 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -169,7 +169,7 @@ namespace corsika::process {
using corsika::geometry::Point;
using namespace corsika::units::si;
- // TODO: this should be done in a central, common place. Not here..
+ // TODO: this should be done in a central, common place. Not here..
#ifndef CORSIKA_OSX
feenableexcept(FE_INVALID);
#endif
@@ -185,7 +185,8 @@ namespace corsika::process {
pin.SetMomentum(p.GetMomentum());
// setting particle mass with Corsika values, may be inconsistent with sibyll
// internal values
- // TODO: #warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
+ // TODO: #warning setting particle mass with Corsika values, may be inconsistent
+ // with sibyll internal values
pin.SetMass(corsika::particles::GetMass(pCode));
// remember position
Point decayPoint = p.GetPosition();
@@ -219,7 +220,7 @@ namespace corsika::process {
// empty sibyll stack
ss.Clear();
- // TODO: this should be done in a central, common place. Not here..
+ // TODO: this should be done in a central, common place. Not here..
#ifndef CORSIKA_OSX
fedisableexcept(FE_INVALID);
#endif
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index b8709216a75bd7a827fea1b063520fa615805553..003d22a5f1bd5d97373be79ab0cbdf3980d48560 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -14,15 +14,15 @@
#include <corsika/process/InteractionProcess.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
-#include <corsika/utl/COMBoost.h>
-#include <corsika/particles/ParticleProperties.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/NuclearComposition.h>
+#include <corsika/utl/COMBoost.h>
namespace corsika::process::sibyll {
@@ -30,8 +30,10 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
+
public:
- Interaction(corsika::environment::Environment const& env) : fEnvironment(env) { }
+ Interaction(corsika::environment::Environment const& env)
+ : fEnvironment(env) {}
~Interaction() {
std::cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount
<< std::endl;
@@ -45,38 +47,37 @@ namespace corsika::process::sibyll {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm");
-
+
// initialize Sibyll
sibyll_ini_();
}
- std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(const corsika::particles::Code &BeamId, const corsika::particles::Code & TargetId, const corsika::units::si::HEPEnergyType& CoMenergy)
- {
+ std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
+ const corsika::particles::Code& BeamId, const corsika::particles::Code& TargetId,
+ const corsika::units::si::HEPEnergyType& CoMenergy) {
using namespace corsika::units::si;
double sigProd, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
const double dEcm = CoMenergy / 1_GeV;
- if(corsika::particles::IsNucleus( TargetId )){
- const int iTarget = corsika::particles::GetNucleusA( TargetId );
- if(iTarget>18||iTarget==0)
- throw std::runtime_error("Sibyll target outside range. Only nuclei with A<18 are allowed.");
- sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
- return std::make_tuple( sigProd * 1_mbarn, iTarget );
+ if (corsika::particles::IsNucleus(TargetId)) {
+ const int iTarget = corsika::particles::GetNucleusA(TargetId);
+ if (iTarget > 18 || iTarget == 0)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 are allowed.");
+ sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
+ return std::make_tuple(sigProd * 1_mbarn, iTarget);
+ } else {
+ if (TargetId == corsika::particles::Proton::GetCode()) {
+ sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
+ return std::make_tuple(sigProd * 1_mbarn, 1);
+ } else
+ // no interaction in sibyll possible, return infinite cross section? or throw?
+ sigProd = std::numeric_limits<double>::infinity();
+ return std::make_tuple(sigProd * 1_mbarn, 0);
}
- else
- {
- if (TargetId == corsika::particles::Proton::GetCode()){
- sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
- return std::make_tuple( sigProd * 1_mbarn, 1 );
- }
- else
- // no interaction in sibyll possible, return infinite cross section? or throw?
- sigProd = std::numeric_limits<double>::infinity();
- return std::make_tuple( sigProd * 1_mbarn, 0 );
- }
}
-
+
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) {
@@ -95,9 +96,9 @@ namespace corsika::process::sibyll {
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
+
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
+ corsika::particles::Neutron::GetMass());
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
@@ -117,42 +118,48 @@ namespace corsika::process::sibyll {
if (kInteraction && Elab >= 8.5_GeV && ECoM >= 10_GeV) {
- // get target from environment
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
- // determine average interaction length
- // weighted sum
- int i =-1;
- double avgTargetMassNumber = 0.;
- si::CrossSectionType weightedProdCrossSection = 0_mbarn;
- // get weights of components from environment/medium
- const auto w = mediumComposition.GetFractions();
- // loop over components in medium
- for (auto targetId: mediumComposition.GetComponents() ){
- i++;
- cout << "Interaction: get interaction length for target: " << targetId << endl;
-
- auto const [ productionCrossSection, numberOfNucleons ] = GetCrossSection( corsikaBeamId, targetId, ECoM);
-
- std::cout << "Interaction: "
- << " IntLength: sibyll return (mb): " << productionCrossSection / 1_mbarn << std::endl;
- weightedProdCrossSection += w[i] * productionCrossSection;
- avgTargetMassNumber += w[i] * numberOfNucleons;
- }
- cout << "Interaction: "
- << "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mbarn
- << endl;
-
- // calculate interaction length in medium
-#warning check interaction length units
- GrammageType int_length = avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
- std::cout << "Interaction: "
- << "interaction length (g/cm2): " << int_length / ( 0.001_kg ) * 1_cm * 1_cm << std::endl;
+ // get target from environment
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+ */
+ const auto currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ // determine average interaction length
+ // weighted sum
+ int i = -1;
+ double avgTargetMassNumber = 0.;
+ si::CrossSectionType weightedProdCrossSection = 0_mbarn;
+ // get weights of components from environment/medium
+ const auto w = mediumComposition.GetFractions();
+ // loop over components in medium
+ for (auto targetId : mediumComposition.GetComponents()) {
+ i++;
+ cout << "Interaction: get interaction length for target: " << targetId << endl;
+
+ auto const [productionCrossSection, numberOfNucleons] =
+ GetCrossSection(corsikaBeamId, targetId, ECoM);
+
+ std::cout << "Interaction: "
+ << " IntLength: sibyll return (mb): "
+ << productionCrossSection / 1_mbarn << std::endl;
+ weightedProdCrossSection += w[i] * productionCrossSection;
+ avgTargetMassNumber += w[i] * numberOfNucleons;
+ }
+ cout << "Interaction: "
+ << "IntLength: weighted CrossSection (mb): "
+ << weightedProdCrossSection / 1_mbarn << endl;
+
+ // calculate interaction length in medium
+#warning check interaction length units
+ GrammageType int_length =
+ avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
+ std::cout << "Interaction: "
+ << "interaction length (g/cm2): "
+ << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
return int_length;
}
@@ -173,142 +180,102 @@ namespace corsika::process::sibyll {
const auto corsikaBeamId = p.GetPID();
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
- << process::sibyll::CanInteract( corsikaBeamId ) << endl;
- if (process::sibyll::CanInteract(corsikaBeamId) ) {
+ << process::sibyll::CanInteract(corsikaBeamId) << endl;
+ if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ // position and time of interaction, not used in Sibyll
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
- // sibyll CS has z along particle momentum
- // FOR NOW: hard coded z-axis for corsika frame
- QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m};
- QuantityVector<length_d> const yAxis{0_m, 1_m, 0_m};
- auto rotation_angles = [](MomentumVector const& pin) {
- const auto p = pin.GetComponents();
- const auto th = acos(p[2] / p.norm());
- const auto ph = atan2(
- p[1] / 1_GeV, p[0] / 1_GeV); // acos( p[0] / sqrt(p[0]*p[0]+p[1]*p[1] ) );
- return std::make_tuple(th, ph);
- };
- // auto pt = []( MomentumVector &p ){
- // return sqrt(p.GetComponents()[0] * p.GetComponents()[0] +
- // p.GetComponents()[1] * p.GetComponents()[1]);
- // };
- // double theta = acos( p.GetMomentum().GetComponents()[2] /
- // p.GetMomentum().norm());
- auto const [theta, phi] = rotation_angles(p.GetMomentum());
- cout << "ProcessSibyll: zenith angle between sibyllCS and rootCS: "
- << theta / M_PI * 180. << endl;
- cout << "ProcessSibyll: azimuth angle between sibyllCS and rootCS: "
- << phi / M_PI * 180. << endl;
- // double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
- const CoordinateSystem tempCS = rootCS.rotate(zAxis, phi);
- const CoordinateSystem sibyllCS = tempCS.rotate(yAxis, theta);
-
-
- // FOR NOW: target is always at rest
+ // define target
+ // for Sibyll is always a single nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
- const auto Etarget = 0_GeV + nucleon_mass;
- const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
- cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
- cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
- << endl;
- cout << "beam momentum in sibyll frame (GeV/c): "
- << p.GetMomentum().GetComponents(sibyllCS) / 1_GeV << endl;
+ corsika::particles::Neutron::GetMass());
+ // FOR NOW: target is always at rest
+ const auto eTargetLab = 0_GeV + nucleon_mass;
+ const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
- cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
- cout << "time: " << tOrig << endl;
+ // define projectile
+ auto const pProjectileLab = p.GetMomentum();
+ HEPEnergyType const eProjectileLab = p.GetEnergy();
- // 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)
- */
- // total energy: E_beam + E_target
- // in lab. frame: E_beam + m_target*c**2
- HEPEnergyType E = p.GetEnergy();
- HEPEnergyType Etot = E + Etarget;
- // total momentum
- MomentumVector Ptot = p.GetMomentum();
- // invariant mass, i.e. cm. energy
- HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+ // define target kinematics in lab frame
+ HEPMassType const targetMass = nucleon_mass;
+ // define boost to and from CoM frame
+ // CoM frame definition in Sibyll projectile: +z
+ COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
-
- // sample target mass number
- const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
- const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition();
- // get cross sections for target materials
- /*
- Here we read the cross section from the interaction model again,
- should be passed from GetInteractionLength if possible
- */
-#warning reading interaction cross section again, should not be necessary
- std::vector<si::CrossSectionType> cross_section_of_components;
- for(auto targetId: mediumComposition.GetComponents() ){
- const auto [sigProd, nNuc] = GetCrossSection( corsikaBeamId, targetId, Ecm);
- cross_section_of_components.push_back( sigProd );
- }
-
- const auto targetCode = currentNode->GetModelProperties().GetTarget( cross_section_of_components);
- cout << "Interaction: target selected: " << targetCode << endl;
- /*
- FOR NOW: allow nuclei with A<18 or protons only.
- when medium composition becomes more complex, approximations will have to be allowed
- air in atmosphere also contains some Argon.
- */
- int targetSibCode = -1;
- if( IsNucleus(targetCode))
- targetSibCode = GetNucleusA( targetCode );
- if( targetCode == corsika::particles::Proton::GetCode())
- targetSibCode = 1;
- cout << "Interaction: sibyll code: " << targetSibCode << endl;
- if(targetSibCode>18||targetSibCode<1)
- throw std::runtime_error("Sibyll target outside range. Only nuclei with A<18 or protons are allowed.");
-
- /*
- 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;
+ cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
+ << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
+ << endl;
+ cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
+ << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
+ << endl;
- std::cout << "Interaction: "
- << " DoDiscrete: gamma:" << gamma << endl;
- std::cout << "Interaction: "
- << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+ // just for show:
+ // boost projecticle
+ auto const [eProjectileCoM, pProjectileCoM] =
+ boost.toCoM(eProjectileLab, pProjectileLab);
- auto const pProjectileLab = p.GetMomentum();
- //{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
- HEPEnergyType const eProjectileLab = p.GetEnergy();
- //energy(projectileMass, pProjectileLab);
+ // boost target
+ auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
- // define target kinematics in lab frame
- HEPMassType const targetMass = nucleon_mass;
- // define boost to com frame
- COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
+ cout << "Interaction: ebeam CoM: " << eProjectileCoM / 1_GeV << endl
+ << "Interaction: pbeam CoM: " << pProjectileCoM / 1_GeV << endl;
+ cout << "Interaction: etarget CoM: " << eTargetCoM / 1_GeV << endl
+ << "Interaction: ptarget CoM: " << pTargetCoM / 1_GeV << endl;
+
+ cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
+ << endl;
+ cout << "Interaction: time: " << tOrig << endl;
- cout << "Interaction: new boost: ebeam lab: " << eProjectileLab / 1_GeV << endl
- << "Interaction: new boost: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV << endl;
+ HEPEnergyType Etot = eProjectileLab + eTargetLab;
+ MomentumVector Ptot = p.GetMomentum();
+ // invariant mass, i.e. cm. energy
+ HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
- // boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
+ // sample target mass number
+ const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
+ const auto mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ // get cross sections for target materials
+ /*
+ Here we read the cross section from the interaction model again,
+ should be passed from GetInteractionLength if possible
+ */
+#warning reading interaction cross section again, should not be necessary
+ std::vector<si::CrossSectionType> cross_section_of_components;
+ for (auto targetId : mediumComposition.GetComponents()) {
+ const auto [sigProd, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ cross_section_of_components.push_back(sigProd);
+ }
- cout << "Interaction: new boost: ebeam com: " << eProjectileCoM / 1_GeV << endl
- << "Interaction: new boost: pbeam com: " << pProjectileCoM / 1_GeV << endl;
-
+ const auto targetCode =
+ currentNode->GetModelProperties().GetTarget(cross_section_of_components);
+ cout << "Interaction: target selected: " << targetCode << endl;
+ /*
+ FOR NOW: allow nuclei with A<18 or protons only.
+ when medium composition becomes more complex, approximations will have to be
+ allowed air in atmosphere also contains some Argon.
+ */
+ int targetSibCode = -1;
+ if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
+ if (targetCode == corsika::particles::Proton::GetCode()) targetSibCode = 1;
+ cout << "Interaction: sibyll code: " << targetSibCode << endl;
+ if (targetSibCode > 18 || targetSibCode < 1)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 or protons are "
+ "allowed.");
+
+ // beam id for sibyll
const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
std::cout << "Interaction: "
- << " DoInteraction: E(GeV):" << E / 1_GeV
+ << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
- if (E < 8.5_GeV || Ecm < 10_GeV) {
+ if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
<< " DoInteraction: should have dropped particle.." << std::endl;
// p.Delete(); delete later... different process
@@ -331,70 +298,44 @@ namespace corsika::process::sibyll {
// add particles from sibyll to stack
// link to sibyll stack
- // here we need to pass projectile momentum and energy to define the local
- // sibyll frame and the boosts to the lab. frame
SibStack ss;
- // momentum array in Sibyll
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
+ HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// skip particles that have decayed in Sibyll
if (psib.HasDecayed()) continue;
- // transform energy to lab. frame, primitve
- // compute beta_vec * p_vec
- // arbitrary Lorentz transformation based on sibyll routines
- const auto gammaBetaComponents = gambet.GetComponents();
- // FOR NOW: fill vector in sibCS and then rotate into rootCS
- // can be done in SibStack by passing sibCS
- // get momentum vector in sibyllCS
- const auto pSibyllComponentsSibCS = psib.GetMomentum().GetComponents();
- // temporary vector in sibyllCS
- auto SibVector = MomentumVector(sibyllCS, pSibyllComponentsSibCS);
- // rotatate to rootCS
- const auto pSibyllComponents = SibVector.GetComponents(rootCS);
- // boost to lab. frame
- HEPEnergyType en_lab = 0. * 1_GeV;
- HEPMomentumType p_lab_components[3];
- en_lab = psib.GetEnergy() * gamma;
- HEPEnergyType pnorm = 0. * 1_GeV;
- for (int j = 0; j < 3; ++j)
- pnorm += (pSibyllComponents[j] * gammaBetaComponents[j]) / (gamma + 1.);
- pnorm += psib.GetEnergy();
-
- for (int j = 0; j < 3; ++j) {
- p_lab_components[j] =
- pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j];
- en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j];
- }
+ // transform energy to lab. frame
+ auto const pCoM = psib.GetMomentum().GetComponents();
+ HEPEnergyType const eCoM = psib.GetEnergy();
+ auto const [eLab, pLab] = boost.fromCoM(eCoM, pCoM);
// add to corsika stack
auto pnew = s.NewParticle();
- pnew.SetEnergy(en_lab);
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- corsika::geometry::QuantityVector<hepmomentum_d> p_lab_c{
- p_lab_components[0], p_lab_components[1], p_lab_components[2]};
- pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
+ pnew.SetEnergy(eLab);
+ pnew.SetMomentum(pLab);
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
+
Plab_final += pnew.GetMomentum();
- E_final += en_lab;
+ Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
}
- std::cout << "conservation (all GeV): E_final=" << E_final / 1_GeV
- << ", Ecm_final=" << Ecm_final / 1_GeV
+ std::cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV
+ << std::endl
+ << "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
<< std::endl;
}
}
return process::EProcessReturn::eOk;
}
-
+
private:
corsika::environment::Environment const& fEnvironment;
-
};
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index eb3259e792aa86899f6569b7a1247dce8acd4cf6..8c8ce3b5e69e563238fadc2416bd1f91351603ac 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -70,9 +70,9 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/particles/ParticleProperties.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
@@ -83,12 +83,12 @@ using namespace corsika::units;
TEST_CASE("SibyllInterface", "[processes]") {
- // setup environment, geometry
+ // setup environment, geometry
corsika::environment::Environment env;
auto& universe = *(env.GetUniverse());
auto theMedium = corsika::environment::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel =
@@ -109,14 +109,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
-
SECTION("InteractionInterface") {
setup::Stack stack;
auto particle = stack.NewParticle();
-
+
Interaction model(env);
-
+
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
@@ -131,7 +130,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
{
const HEPEnergyType E0 = 10_GeV;
particle.SetPID(particles::Code::Proton);
- HEPMomentumType P0 = sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
particle.SetEnergy(E0);
particle.SetMomentum(plab);
@@ -139,9 +139,9 @@ TEST_CASE("SibyllInterface", "[processes]") {
geometry::Point p(cs, 0_m, 0_m, 0_m);
particle.SetPosition(p);
}
-
+
Decay model;
-
+
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
stack);