update CONEXSourceCut to latest CONEX interface

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());
 }

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