diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.cc b/Processes/CONEXSourceCut/CONEXSourceCut.cc index 8b90ff0ab167fe6893436b84e0f1aa667b2bb29a..77fffb0e65363aada42d45b4f5ab1c6ba7aa9a47 100644 --- a/Processes/CONEXSourceCut/CONEXSourceCut.cc +++ b/Processes/CONEXSourceCut/CONEXSourceCut.cc @@ -34,10 +34,10 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries( continue; // no EM particle } - int egs_pid = it->second; + auto const egs_pid = it->second; - addParticle(egs_pid, p.GetEnergy(), p.GetPosition(), p.GetMomentum().normalized(), - p.GetTime()); + addParticle(egs_pid, p.GetEnergy(), p.GetMass(), p.GetPosition(), + p.GetMomentum().normalized(), p.GetTime()); p.Delete(); } @@ -45,47 +45,52 @@ corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries( return corsika::process::EProcessReturn::eOk; } -void CONEXSourceCut::addParticle(int egs_pid, HEPEnergyType energy, +void CONEXSourceCut::addParticle(int egs_pid, HEPEnergyType energy, HEPEnergyType mass, geometry::Point const& position, geometry::Vector<dimensionless_d> const& direction, TimeType t) { std::cout << "position conexObs: " << position.GetCoordinates(conexObservationCS_) << std::endl; - auto coords = position.GetCoordinates(conexObservationCS_) / 1_m; - double x = coords[0].magnitude(); - double y = coords[1].magnitude(); + auto const coords = position.GetCoordinates(conexObservationCS_) / 1_m; + double const x = coords[0].magnitude(); + double const y = coords[1].magnitude(); - double altitude = ((position - center_).norm() - conex::earthRadius) / 1_m; + double const altitude = ((position - center_).norm() - conex::earthRadius) / 1_m; auto const d = position - showerCore_; // distance from core to particle projected along shower axis - double slantDistance = -d.dot(showerAxis_.GetDirection()) / 1_m; + double const slantDistance = -d.dot(showerAxis_.GetDirection()) / 1_m; // lateral coordinates in CONEX shower frame 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 const lateralX = dShowerPlane.dot(x_sf_) / 1_m; + double const lateralY = dShowerPlane.dot(y_sf_) / 1_m; - double slantX = showerAxis_.projectedX(position) * (1_cm * 1_cm / 1_g); + double const slantX = showerAxis_.projectedX(position) * (1_cm * 1_cm / 1_g); - double time = (t * units::constants::c - groundDist_) / 1_m; + double const time = (t * units::constants::c - groundDist_) / 1_m; // fill u,v,w momentum direction in EGS frame - double u = direction.dot(y_sf_).magnitude(); - double v = direction.dot(x_sf_).magnitude(); - double w = direction.dot(showerAxis_.GetDirection()).magnitude(); + double const u = direction.dot(y_sf_).magnitude(); + double const v = direction.dot(x_sf_).magnitude(); + double const w = direction.dot(showerAxis_.GetDirection()).magnitude(); - double weight = 1; + double const weight = 1; // NEEDS TO BE CHANGED WHEN WE HAVE WEIGHTS! - double E = energy / 1_GeV; + // generation, TO BE CHANGED WHEN WE HAVE THAT INFORMATION AVAILABLE + int const latchin = 1; + + double const E = energy / 1_GeV; + double const m = mass / 1_GeV; std::cout << "CONEXSourceCut: removing " << egs_pid << " " << std::scientific << energy - << std::endl; + << " GeV" << std::endl; - std::cout << "#### parameters to show_() ####" << std::endl; + std::cout << "#### parameters to cegs4_() ####" << std::endl; std::cout << "egs_pid = " << egs_pid << std::endl; std::cout << "E = " << E << std::endl; + std::cout << "m = " << m << std::endl; std::cout << "x = " << x << std::endl; std::cout << "y = " << y << std::endl; std::cout << "altitude = " << altitude << std::endl; @@ -100,11 +105,13 @@ void CONEXSourceCut::addParticle(int egs_pid, HEPEnergyType energy, conex::cxoptl_.dptl[10 - 1] = egs_pid; conex::cxoptl_.dptl[4 - 1] = E; + conex::cxoptl_.dptl[5 - 1] = m; conex::cxoptl_.dptl[6 - 1] = x; conex::cxoptl_.dptl[7 - 1] = y; conex::cxoptl_.dptl[8 - 1] = altitude; conex::cxoptl_.dptl[9 - 1] = time; conex::cxoptl_.dptl[11 - 1] = weight; + conex::cxoptl_.dptl[12 - 1] = latchin; conex::cxoptl_.dptl[13 - 1] = slantX; conex::cxoptl_.dptl[14 - 1] = lateralX; conex::cxoptl_.dptl[15 - 1] = lateralY; diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.h b/Processes/CONEXSourceCut/CONEXSourceCut.h index fec7190d3e057d2f9e63859238a7f30a1e9d9333..121d8c16369a144bac054678f1c6b9af021b308a 100644 --- a/Processes/CONEXSourceCut/CONEXSourceCut.h +++ b/Processes/CONEXSourceCut/CONEXSourceCut.h @@ -38,7 +38,7 @@ namespace corsika::process { void SolveCE(); void addParticle(int egs_pid, units::si::HEPEnergyType energy, - geometry::Point const& position, + units::si::HEPEnergyType mass, geometry::Point const& position, geometry::Vector<units::si::dimensionless_d> const& direction, units::si::TimeType t); diff --git a/Processes/CONEXSourceCut/testCONEXSourceCut.cc b/Processes/CONEXSourceCut/testCONEXSourceCut.cc index dacfc8728d07f5056b5aee36153e97f362097478..18df797fe78d62d216df5c5e91de0b1dc0b4c400 100644 --- a/Processes/CONEXSourceCut/testCONEXSourceCut.cc +++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc @@ -88,7 +88,7 @@ TEST_CASE("CONEXSourceCut") { std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " " << emPosition.GetCoordinates(rootCS) << std::endl; - conex.addParticle(0, Eem, emPosition, momentum.normalized(), 0_s); + conex.addParticle(0, Eem, 0_eV, emPosition, momentum.normalized(), 0_s); conex.SolveCE(); }