diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.cc b/Processes/CONEXSourceCut/CONEXSourceCut.cc
index c3109315c9f66784a24c68ce2bddf3c952fe4649..a312349085fe00e08c2df5123c8a8e8d65b0d208 100644
--- a/Processes/CONEXSourceCut/CONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.cc
@@ -9,6 +9,7 @@
*/
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
+#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalConstants.h>
#include <algorithm>
#include <iomanip>
@@ -27,9 +28,9 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
while (p != vS.end()) {
Code const pid = p.GetPID();
- auto const it = std::find_if(em_codes_.cbegin(), em_codes_.cend(),
+ auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
[=](auto const& p) { return pid == p.first; });
- if (it == em_codes_.cend()) {
+ if (it == egs_em_codes_.cend()) {
continue; // no EM particle
}
@@ -42,8 +43,8 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
energy / 1_MeV; // total energy, TODO: check if maybe kinetic should be used
auto coords = position.GetCoordinates(conexObservationCS_) / 1_m;
- double x = coords[0];
- double y = coords[1];
+ double x = coords[0].magnitude();
+ double y = coords[1].magnitude();
double altitude = ((position - center_).norm() - conex::earthRadius) / 1_m;
@@ -61,9 +62,9 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
double time = (p.GetTime() * units::constants::c - groundDist_) / 1_m;
// fill u,v,w momentum direction in EGS frame
- double u = direction.dot(y_sf_);
- double v = direction.dot(x_sf_);
- double w = direction.dot(showerAxis_.GetDirection());
+ double u = direction.dot(y_sf_).magnitude();
+ double v = direction.dot(x_sf_).magnitude();
+ double w = direction.dot(showerAxis_.GetDirection()).magnitude();
int iri = 2; // EGS medium air
@@ -84,21 +85,17 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
void CONEXSourceCut::Init() {}
-void CONEXSourceCut::SolveCE() {
- int zero = 0;
- int iCEmode = 1;
- conex::HadronCascade_(id, nshtot_, zero, iCEmode);
- conex::SolveMomentEquations_(zero);
-}
+void CONEXSourceCut::SolveCE() { conex::conexcascade_(); }
CONEXSourceCut::CONEXSourceCut(geometry::Point center, environment::ShowerAxis showerAxis,
units::si::LengthType groundDist,
units::si::GrammageType Xcut,
- units::si::EnergyType primaryEnergy)
+ units::si::HEPEnergyType primaryEnergy,
+ particles::PDGCode primaryID)
: center_{center}
, showerAxis_{showerAxis}
, groundDist_{groundDist}
- , conexObservation_{std::invoke([]() {
+ , conexObservationCS_{std::invoke([&]() {
auto const& c8cs = center.GetCoordinateSystem();
auto const showerCore =
showerAxis.GetStart() + showerAxis.GetDirection() * groundDist;
@@ -110,31 +107,34 @@ CONEXSourceCut::CONEXSourceCut(geometry::Point center, environment::ShowerAxis s
c8cs, intermediateCS2); // either this way or vice versa... TODO: test this!
return geometry::CoordinateSystem(c8cs, transform);
})}
- , x_sf_{geometry::Vector<dimensionless_d>{conexObservationCS_, 0., 0., 1.}
- .cross(showerAxis_.GetDirection())
- .normalized()}
+ , x_sf_{std::invoke([&]() {
+ return geometry::Vector<length_d>{conexObservationCS_, 0._m, 0._m, 1._m}
+ .cross(showerAxis_.GetDirection())
+ .normalized();
+ })}
, y_sf_{showerAxis_.GetDirection().cross(x_sf_)} {
- auto id = static_cast<int>(InitialParticle_(particles::GetPDG(pid)));
- conex::eprima_ = primaryEnergy / 1_GeV;
- nshtot = 1; // not sure about this...
- conex::fehcut_ = conex::femcut_ = ; // what?
- conex::ehcut_ = max(enymin, min(1.d10, eprima / aNbrNucl * fehcut));
- conex::emcut_ = max(enymin, min(1.d10, eprima / aNbrNucl * femcut));
-
- // set phisho, thetas
- geometry::Vector<dimensionless_d> ez{conexObservationCS_, 0., 0., 1.};
- auto const c = showerAxis_.GetDirection().dot(ez);
- auto const theta = M_PI - std::acos(c);
- conex::cxbas4_::thetas_ = theta * 180;
+
+ double eprima = primaryEnergy / 1_GeV;
+
+ // set phi, theta
+ geometry::Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, 1_m}};
+ auto const c = showerAxis_.GetDirection().dot(ez) / 1_m;
+ double theta = 180 * (M_PI - std::acos(c));
auto const showerAxisConex =
showerAxis_.GetDirection().GetComponents(conexObservationCS_);
- auto const phi = std::atan2(-showerAxisConex.GetY(), showerAxisConex.GetX());
- conex::cxbas4_::phisho_ = phi * 180;
+ double phi = 180 * std::atan2(-showerAxisConex.GetY().magnitude(),
+ showerAxisConex.GetX().magnitude());
+ double XmaxP_ = Xcut / (1_g / 1_cm / 1_cm);
+
+ int ipart = static_cast<int>(primaryID);
+ auto rng = corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
- // call ranfgt(seed)?
+ double dimpact = 0.; // valid only if shower core is fixed on the observation plane; for
+ // skimming showers an offset is needed like in CONEX
- cxthin_::ethin = 0;
+ std::array<int, 3> ioseed{static_cast<int>(rng()), static_cast<int>(rng()),
+ static_cast<int>(rng())};
- conex::cxoutput1_::XminP_ = conex::cxoutput1_::XmaxP_ = Xcut / (1_g / 1_cm / 1_cm);
+ conex::conexrun_(ipart, eprima, theta, phi, dimpact, ioseed.data());
}
diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.h b/Processes/CONEXSourceCut/CONEXSourceCut.h
index 3dced0507552bd0e89b931448e14f36153bf12b5..918f6b1f3156c5fb260cde9aef1373d39a29f824 100644
--- a/Processes/CONEXSourceCut/CONEXSourceCut.h
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.h
@@ -11,6 +11,9 @@
#ifndef _corsika_process_particle_cut_CONEXSourceCut_h_
#define _corsika_process_particle_cut_CONEXSourceCut_h_
+#include <corsika/environment/ShowerAxis.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/setup/SetupStack.h>
@@ -18,20 +21,20 @@
namespace conex {
extern "C" {
- extern struct cxcut_ {
+ extern struct {
double eecut;
double epcut;
double ehcut;
double emcut;
- };
+ } cxcut_;
- extern struct cxsubcut_ {
+ extern struct {
double feecut_;
double fehcut_;
double femcut_;
- };
+ } cxsubcut_;
- extern struct cxbas4_ {
+ extern struct {
double eprima_;
double thetas_;
double costhet_;
@@ -44,94 +47,14 @@ namespace conex {
double XminSlant;
double HGrd;
double distMaxi;
- };
-
- extern struct cxoptl_ {
- struct dptl_ {
- double px, py, pz, E, m; // 5-momentum
- double x, y;
- double h;
- double t;
- double id;
- double weight;
- double generation;
- double Xslant;
- double lateralX, lateralY;
- double slantDistance;
- };
- };
-
- extern struct cxbas6_ {
- double altitude;
- double RadGrd;
- double DistALt;
- double dphmaxi0;
- double dphlim0;
- double dphmin0;
- bool goOutGrd;
- };
-
- extern struct cxetc_ {
- int mode;
- int iwrt;
- int i1DMC;
- int iphonu;
- };
-
- extern struct cxthin_ {
- double thin;
- double ethin;
- double wtmax;
- double rthmax;
- int iothin;
- };
-
- extern struct cxoutput3_ {
- std::array<std::array<double, 4>, mxPxpro + 2> XmaxShow;
- std::array<std::array<double, 5>, mxPxpro + 2> XmaxMean;
- std::array<std::array<double, maximz>, mxPxpro> XmaxProf;
- int iXmax;
- };
-
- int constexpr mxExpro = 1;
- int constexpr mxPxpro = 10,
-
- int constexpr maximE = 281;
-
-// CONEX ifdef
-#ifdef __SAVEMEMO__
- int constexpr maximZ = 201;
-#else
- int constexpr maximZ = 3701;
-#endif
-
- extern struct cxoutput1_ {
- std::array<std::array<std::array<double, maximZ>, mxExpro>, mxPxpro + 2> XProf,
- XmeanP, XmeanP2;
- std::array<double, mxExpro> EMCutP, HaCutP;
- double XminP;
- double XmaxP;
- int nminX;
- int nmaxX;
- int mZEMHa;
- int ifout;
- int ivers;
- bool lheader;
- };
-
- /*
- common/cxoutput3/XmaxShow(4,-1:mxPxpro),XmaxMean(5,-1:mxPxpro)
- &,XmaxProf(maximz,-1:mxPxpro),iXmax
- common/cxoutput1/XProf(maximz,mxExpro,-1:mxPxpro)
- &,XmeanP(maximz,mxExpro,-1:mxPxpro)
- &,XmeanP2(maximz,mxExpro,-1:mxPxpro)
- &,EMCutP(mxExpro),HaCutP(mxExpro),XminP,XmaxP
- &,nminX,nmaxX,mZEMHa,ifout,ivers,lheader
- */
-
- int InitialParticle_(int&);
- void HadronCascade_(int&, int&, int&, int&);
- void SolveMomentEquations_(int&);
+ } cxbas4_;
+
+ // ipart,energy,theta,phi,dimpact,ioseed
+ void conexrun_(int& ipart, double& energy, double& theta, double& phi, double& dimpact,
+ int ioseed[3]);
+ void conexcascade_();
+ void hadroncascade_(int&, int&, int&, int&);
+ void solvemomentequations_(int&);
void show_(int& iqi, double& ei, double& xmi, double& ymi, double& zmi, double& dmi,
double& xi, double& yi, double& zi, double& tmi, double& ui, double& vi,
double& wi, int& iri, double& wti, int& latchi);
@@ -146,8 +69,10 @@ namespace corsika::process {
class CONEXSourceCut : public process::SecondariesProcess<CONEXSourceCut> {
public:
- CONEXSourceCut(geometry::Point const& center, environment::ShowerAxis const&,
- units::si::LengthType, units::si::GrammageType);
+ CONEXSourceCut(geometry::Point center, environment::ShowerAxis showerAxis,
+ units::si::LengthType groundDist, units::si::GrammageType Xcut,
+ units::si::HEPEnergyType primaryEnergy,
+ particles::PDGCode primaryID);
corsika::process::EProcessReturn DoSecondaries(corsika::setup::StackView&);
void Init();
@@ -157,18 +82,15 @@ namespace corsika::process {
private:
//! CONEX e.m. particle codes
static std::array<std::pair<particles::Code, int>, 3> constexpr egs_em_codes_{
- {particles::Code::Gamma, 0},
- {particles::Code::Electron, -1},
- {particles::Code::Positron, -1}};
-
- int egs_calls_{1};
- int nshtot_{1};
+ {{particles::Code::Gamma, 0},
+ {particles::Code::Electron, -1},
+ {particles::Code::Positron, -1}}};
- units::si::LengthType groundDist_; //!< length from injection point to shower core
- geometry::Point const center_; //!< center of CONEX Earth
+ geometry::Point const center_; //!< center of CONEX Earth
environment::ShowerAxis const& showerAxis_;
+ units::si::LengthType groundDist_; //!< length from injection point to shower core
geometry::CoordinateSystem const conexObservationCS_; //!< CONEX observation frame
- Vector<units::si::dimensionless_d> const x_sf_,
+ geometry::Vector<units::si::dimensionless_d> const x_sf_,
y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction
};
} // namespace conex_source_cut