diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.cc b/Processes/CONEXSourceCut/CONEXSourceCut.cc
index d82ecb7ead67e94ba2df3f478802a28e99450b5f..c3109315c9f66784a24c68ce2bddf3c952fe4649 100644
--- a/Processes/CONEXSourceCut/CONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.cc
@@ -9,10 +9,10 @@
*/
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
+#include <corsika/units/PhysicalConstants.h>
#include <algorithm>
#include <iomanip>
-
-using namespace std;
+#include <utility>
using namespace corsika::process::conex_source_cut;
using namespace corsika::units::si;
@@ -22,20 +22,61 @@ using namespace corsika::setup;
corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
corsika::setup::StackView& vS) {
auto p = vS.begin();
+ HEPEnergyType const energy = p.GetEnergy();
+
while (p != vS.end()) {
Code const pid = p.GetPID();
- HEPEnergyType const energy = p.GetEnergy();
- if (std::find(em_codes_.cbegin(), em_codes_.cend(), pid) != em_codes_.cend()) {
- std::cout << "CONEXSourceCut: removing " << pid << " " << std::scientific << energy
- << std::endl;
-
- p.Delete();
+ auto const it = std::find_if(em_codes_.cbegin(), em_codes_.cend(),
+ [=](auto const& p) { return pid == p.first; });
+ if (it == em_codes_.cend()) {
+ continue; // no EM particle
}
- else {
- ++p;
- }
+ auto const& position = p.GetPosition();
+ auto const direction = p.GetMomentum().normalized();
+
+ int egs_pid = it->second;
+
+ double E =
+ 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 altitude = ((position - center_).norm() - conex::earthRadius) / 1_m;
+
+ double slantDistance =
+ (position - showerAxis_.GetStart()).dot(showerAxis_.GetDirection()) / 1_m;
+
+ // lateral coordinates in CONEX shower frame
+ auto const d = position - showerAxis_.GetStart();
+ auto const dShowerPlane = d - d.parallelProjectionOnto(showerAxis_.GetDirection());
+ double lateralX = dShowerPlane.dot(x_sf_) / 1_m;
+ double lateralY = dShowerPlane.dot(y_sf_) / 1_m;
+
+ double slantX = showerAxis_.projectedX(position) * (1_cm * 1_cm / 1_g);
+
+ 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());
+
+ int iri = 2; // EGS medium air
+
+ double weight = 1;
+
+ int latchin = 1; // generation, we don't have the actual value...
+
+ std::cout << "CONEXSourceCut: removing " << pid << " " << std::scientific << energy
+ << std::endl;
+ p.Delete();
+
+ conex::show_(egs_pid, E, x, y, altitude, slantDistance, lateralX, lateralY, slantX,
+ time, u, v, w, iri, weight, latchin);
}
return corsika::process::EProcessReturn::eOk;
@@ -43,4 +84,57 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
void CONEXSourceCut::Init() {}
-CONEXSourceCut::CONEXSourceCut() {}
+void CONEXSourceCut::SolveCE() {
+ int zero = 0;
+ int iCEmode = 1;
+ conex::HadronCascade_(id, nshtot_, zero, iCEmode);
+ conex::SolveMomentEquations_(zero);
+}
+
+CONEXSourceCut::CONEXSourceCut(geometry::Point center, environment::ShowerAxis showerAxis,
+ units::si::LengthType groundDist,
+ units::si::GrammageType Xcut,
+ units::si::EnergyType primaryEnergy)
+ : center_{center}
+ , showerAxis_{showerAxis}
+ , groundDist_{groundDist}
+ , conexObservation_{std::invoke([]() {
+ auto const& c8cs = center.GetCoordinateSystem();
+ auto const showerCore =
+ showerAxis.GetStart() + showerAxis.GetDirection() * groundDist;
+ auto const translation = showerCore - center;
+ auto const intermediateCS = c8cs.translate(translation.GetComponents(c8cs));
+ auto const intermediateCS2 = intermediateCS.RotateToZ(translation);
+
+ auto const transform = geometry::CoordinateSystem::GetTransformation(
+ 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()}
+ , 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;
+
+ auto const showerAxisConex =
+ showerAxis_.GetDirection().GetComponents(conexObservationCS_);
+ auto const phi = std::atan2(-showerAxisConex.GetY(), showerAxisConex.GetX());
+ conex::cxbas4_::phisho_ = phi * 180;
+
+ // call ranfgt(seed)?
+
+ cxthin_::ethin = 0;
+
+ conex::cxoutput1_::XminP_ = conex::cxoutput1_::XmaxP_ = Xcut / (1_g / 1_cm / 1_cm);
+}
diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.h b/Processes/CONEXSourceCut/CONEXSourceCut.h
index 037451bb37d8b4d8006165c9270608c2d1981774..3dced0507552bd0e89b931448e14f36153bf12b5 100644
--- a/Processes/CONEXSourceCut/CONEXSourceCut.h
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.h
@@ -16,19 +16,160 @@
#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
+namespace conex {
+ extern "C" {
+ extern struct cxcut_ {
+ double eecut;
+ double epcut;
+ double ehcut;
+ double emcut;
+ };
+
+ extern struct cxsubcut_ {
+ double feecut_;
+ double fehcut_;
+ double femcut_;
+ };
+
+ extern struct cxbas4_ {
+ double eprima_;
+ double thetas_;
+ double costhet_;
+ double phisho_;
+ int muse;
+ int musz;
+ double c2bas;
+ double sinthet;
+ double sinphi;
+ 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&);
+ 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);
+ }
+
+ corsika::units::si::LengthType constexpr earthRadius{6371315 *
+ corsika::units::si::meter};
+} // namespace conex
+
namespace corsika::process {
namespace conex_source_cut {
class CONEXSourceCut : public process::SecondariesProcess<CONEXSourceCut> {
public:
- CONEXSourceCut();
+ CONEXSourceCut(geometry::Point const& center, environment::ShowerAxis const&,
+ units::si::LengthType, units::si::GrammageType);
corsika::process::EProcessReturn DoSecondaries(corsika::setup::StackView&);
void Init();
+ void SolveCE();
+
private:
- static std::array<particles::Code, 3> constexpr em_codes_{
- particles::Code::Gamma, particles::Code::Electron, particles::Code::Positron};
+ //! 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};
+
+ units::si::LengthType groundDist_; //!< length from injection point to shower core
+ geometry::Point const center_; //!< center of CONEX Earth
+ environment::ShowerAxis const& showerAxis_;
+ geometry::CoordinateSystem const conexObservationCS_; //!< CONEX observation frame
+ 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
} // namespace corsika::process