diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 92c52b3a2eb71a92be1d2a79d54f8968d34768f0..03ebe5e69cf0737e78445956c7d2b6350969c48a 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -162,5 +162,5 @@ int main() { << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl; - cut.Reset(); + cut.Reset(); } diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc index 32180ffc527a3bdcc04d02cb842b767e53909e1e..5840bc597219d7bfcebb957884dc78b26f512b7c 100644 --- a/Documentation/Examples/em_shower.cc +++ b/Documentation/Examples/em_shower.cc @@ -159,15 +159,16 @@ int main(int argc, char** argv) { cut.ShowResults(); em_continuous.ShowResults(); observationLevel.ShowResults(); - cout << "Cascade energy cut: " << EAS.GetEnergyCut()/1_GeV << " GeV" << endl; - const HEPEnergyType Efinal = - cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround() + EAS.GetEnergyCut(); + cout << "Cascade energy cut: " << EAS.GetEnergyCut() / 1_GeV << " GeV" << endl; + const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + + observationLevel.GetEnergyGround() + EAS.GetEnergyCut(); cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; observationLevel.Reset(); cut.Reset(); em_continuous.Reset(); - + auto const hists = proposalCounted.GetHistogram(); hists.saveLab("inthist_lab.txt"); hists.saveCMS("inthist_cms.txt"); diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index 8b3b8fbc4d5434919818206a3b32e63a122eba5a..9e833addea3f599f6835bd511062540fe7de58cf 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -22,12 +22,12 @@ #include <corsika/process/observation_plane/ObservationPlane.h> #include <corsika/process/on_shell_check/OnShellCheck.h> #include <corsika/process/particle_cut/ParticleCut.h> +#include <corsika/process/proposal/ContinuousProcess.h> +#include <corsika/process/proposal/Interaction.h> #include <corsika/process/pythia/Decay.h> #include <corsika/process/sibyll/Decay.h> #include <corsika/process/sibyll/Interaction.h> #include <corsika/process/sibyll/NuclearInteraction.h> -#include <corsika/process/proposal/ContinuousProcess.h> -#include <corsika/process/proposal/Interaction.h> #include <corsika/process/switch_process/SwitchProcess.h> #include <corsika/process/tracking_line/TrackingLine.h> #include <corsika/process/urqmd/UrQMD.h> @@ -207,8 +207,8 @@ int main(int argc, char** argv) { 55_GeV); auto decaySequence = decayPythia << decaySibyll; - auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << em_continuous - << cut << longprof << observationLevel; + auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted + << em_continuous << cut << longprof << observationLevel; // define air shower object, run simulation tracking_line::TrackingLine tracking; @@ -220,13 +220,13 @@ int main(int argc, char** argv) { EAS.Run(); - cut.ShowResults(); em_continuous.ShowResults(); observationLevel.ShowResults(); - cout << "Cascade energy cut: " << EAS.GetEnergyCut()/1_GeV << " GeV" << endl; - const HEPEnergyType Efinal = - cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround() + EAS.GetEnergyCut(); + cout << "Cascade energy cut: " << EAS.GetEnergyCut() / 1_GeV << " GeV" << endl; + const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + + observationLevel.GetEnergyGround() + EAS.GetEnergyCut(); cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; observationLevel.Reset(); @@ -234,7 +234,7 @@ int main(int argc, char** argv) { em_continuous.Reset(); auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() + - urqmdCounted.GetHistogram() + proposalCounted.GetHistogram(); + urqmdCounted.GetHistogram() + proposalCounted.GetHistogram(); hists.saveLab("inthist_lab.txt"); hists.saveCMS("inthist_cms.txt"); diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h index 227d9fe50ecbf8756e9728ca1ec53f90307ab1d0..42b1beca0af1791e057183de49efe4d2710bd192 100644 --- a/Environment/NuclearComposition.h +++ b/Environment/NuclearComposition.h @@ -18,8 +18,6 @@ #include <stdexcept> #include <vector> - - namespace corsika::environment { class NuclearComposition { std::vector<float> const fNumberFractions; //!< relative fractions of number density @@ -29,7 +27,7 @@ namespace corsika::environment { double const fAvgMassNumber; std::size_t hash_; - + template <class AConstIterator, class BConstIterator> class WeightProviderIterator { AConstIterator fAIter; @@ -134,21 +132,17 @@ namespace corsika::environment { // Note: when this class ever modifies its internal data, the hash // must be updated, too! size_t hash() const { return hash_; } - + private: void updateHash() { std::vector<std::size_t> hashes; - for (float ifrac : GetFractions()) - hashes.push_back(std::hash<float>{}(ifrac)); + for (float ifrac : GetFractions()) hashes.push_back(std::hash<float>{}(ifrac)); for (corsika::particles::Code icode : GetComponents()) - hashes.push_back(std::hash<int>{}(static_cast<int>(icode))); + hashes.push_back(std::hash<int>{}(static_cast<int>(icode))); std::size_t h = std::hash<double>{}(GetAverageMassNumber()); - for (std::size_t ih : hashes) - h = h ^ (ih<<1); + for (std::size_t ih : hashes) h = h ^ (ih << 1); hash_ = h; } - }; } // namespace corsika::environment - diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 21fad265637112c38920d0a2e91b05090a0afd36..105cd94465e6fd1d4933d74a39a612e9407455fc 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -86,8 +86,8 @@ namespace corsika::cascade { , fStack(stack) , energy_cut_(0 * corsika::units::si::electronvolt) {} - corsika::units::si::HEPEnergyType GetEnergyCut() const { return energy_cut_; } - + corsika::units::si::HEPEnergyType GetEnergyCut() const { return energy_cut_; } + /** * set the nodes for all particles on the stack according to their numerical * position @@ -212,11 +212,11 @@ namespace corsika::cascade { // apply all continuous processes on particle + track process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step); - + if (status == process::EProcessReturn::eParticleAbsorbed) { std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " " << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl; - energy_cut_ += vParticle.GetEnergy(); + energy_cut_ += vParticle.GetEnergy(); vParticle.Delete(); return; } @@ -320,10 +320,10 @@ namespace corsika::cascade { TProcessList& fProcessSequence; TStack& fStack; corsika::random::RNG& fRNG = - corsika::random::RNGManager::GetInstance().GetRandomStream("cascade"); - + corsika::random::RNGManager::GetInstance().GetRandomStream("cascade"); + corsika::units::si::HEPEnergyType energy_cut_; - + }; // namespace corsika::cascade } // namespace corsika::cascade diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc index e10d82e59cfe759f4ffc8a83541d08109317d46d..0b2855f4ebc6f343ca10ff01b9ec066f5d4865da 100644 --- a/Framework/ProcessSequence/testProcessSequence.cc +++ b/Framework/ProcessSequence/testProcessSequence.cc @@ -145,11 +145,9 @@ public: }; class Decay1 : public DecayProcess<Decay1> { - int fV = 0; public: - Decay1(const int v) - : fV(v) { + Decay1(const int v) { cout << "Decay1()" << endl; globalCount++; } diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc index 62790e9f50d38951555b3adf6bf16ecddaea342b..69dcd7877e49b760980b15056c657a0468bf546f 100644 --- a/Processes/EnergyLoss/EnergyLoss.cc +++ b/Processes/EnergyLoss/EnergyLoss.cc @@ -31,7 +31,7 @@ using SetupTrack = corsika::setup::Trajectory; using namespace corsika::process::energy_loss; EnergyLoss::EnergyLoss(environment::ShowerAxis const& shower_axis, - corsika::units::si::HEPEnergyType emCut) + corsika::units::si::HEPEnergyType emCut) : shower_axis_(shower_axis) , emCut_(emCut) , profile_(int(shower_axis.maximumX() / dX_) + 1) {} @@ -178,7 +178,7 @@ process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack co << " E=" << E / 1_GeV << "GeV, Ekin=" << Ekin / 1_GeV << ", Enew=" << Enew / 1_GeV << "GeV" << endl; auto status = process::EProcessReturn::eOk; - if (E<emCut_) { + if (E < emCut_) { Enew = emCut_; status = process::EProcessReturn::eParticleAbsorbed; } @@ -195,7 +195,7 @@ LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle, } auto constexpr dX = 1_g / square(1_cm); - auto const dEdX = - TotalEnergyLoss(vParticle, dX) / dX; // dE > 0 + auto const dEdX = -TotalEnergyLoss(vParticle, dX) / dX; // dE > 0 //~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass(); // in any case: never go below 0.99*emCut_ This needs to be @@ -207,8 +207,8 @@ LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle, // afterwards. // const auto energy = vParticle.GetEnergy(); - auto energy_lim = std::max(0.9*energy, 0.99*emCut_); - + auto energy_lim = std::max(0.9 * energy, 0.99 * emCut_); + auto const maxGrammage = (energy - energy_lim) / dEdX; return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack, diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h index 2e34650353491c6379f53ad8e288db3b41179607..89d20bf05a09346eaf825ce647d990f316ee7006 100644 --- a/Processes/EnergyLoss/EnergyLoss.h +++ b/Processes/EnergyLoss/EnergyLoss.h @@ -27,9 +27,10 @@ namespace corsika::process::energy_loss { units::si::square(1e-2 * units::si::meter)); void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew); - + public: - EnergyLoss(environment::ShowerAxis const& showerAxis, corsika::units::si::HEPEnergyType emCut); + EnergyLoss(environment::ShowerAxis const& showerAxis, + corsika::units::si::HEPEnergyType emCut); process::EProcessReturn DoContinuous(setup::Stack::ParticleType&, setup::Trajectory const&); @@ -52,7 +53,7 @@ namespace corsika::process::energy_loss { return 10_g / square(1_cm); }); // profile binning - private: + private: environment::ShowerAxis const& shower_axis_; corsika::units::si::HEPEnergyType emCut_; std::vector<units::si::HEPEnergyType> profile_; // longitudinal profile diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc index be331f2232bf58b202fa66df75d8f9490f637416..cb7ee95bd6241c942189417da34107a6e19efbd0 100644 --- a/Processes/ObservationPlane/ObservationPlane.cc +++ b/Processes/ObservationPlane/ObservationPlane.cc @@ -35,7 +35,7 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous( if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) { return process::EProcessReturn::eOk; } - + const auto energy = particle.GetEnergy(); outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' ' << energy / 1_eV << ' ' @@ -67,10 +67,10 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&, void ObservationPlane::ShowResults() const { std::cout << " ******************************" << std::endl - << " ObservationPlane: " << std::endl; + << " ObservationPlane: " << std::endl; std::cout << " energy in ground (GeV) : " << energy_ground_ / 1_GeV << std::endl - << " no. of particles in ground : " << count_ground_ << std::endl - << " ******************************" << std::endl; + << " no. of particles in ground : " << count_ground_ << std::endl + << " ******************************" << std::endl; } void ObservationPlane::Reset() { diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h index 96d2e77d463673f7364170013818b0359fa599eb..7223defb1748b5c67a04969316725fd09c65b8d9 100644 --- a/Processes/ObservationPlane/ObservationPlane.h +++ b/Processes/ObservationPlane/ObservationPlane.h @@ -39,7 +39,7 @@ namespace corsika::process::observation_plane { void ShowResults() const; void Reset(); corsika::units::si::HEPEnergyType GetEnergyGround() const { return energy_ground_; } - + private: geometry::Plane const plane_; std::ofstream outputStream_; diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc index 05e543674631439333b18325da12ccf4b8b0e539..a503623634fb9f132797ce55049b767a0f5a0162 100644 --- a/Processes/ParticleCut/ParticleCut.cc +++ b/Processes/ParticleCut/ParticleCut.cc @@ -82,7 +82,8 @@ namespace corsika::process { return EProcessReturn::eOk; } - ParticleCut::ParticleCut(const units::si::HEPEnergyType eCut, bool discardEm, bool discardInv) + ParticleCut::ParticleCut(const units::si::HEPEnergyType eCut, bool discardEm, + bool discardInv) : eCut_(eCut) , discardEm_(discardEm) , discardInv_(discardInv) { @@ -95,14 +96,13 @@ namespace corsika::process { } void ParticleCut::ShowResults() const { - cout << " ******************************" << endl - << " ParticleCut: " << endl; - if (discardEm_) - cout << " energy in em. component (GeV) : " << emEnergy_ / 1_GeV << endl - << " no. of em. particles removed : " << emCount_ << endl; + cout << " ******************************" << endl << " ParticleCut: " << endl; + if (discardEm_) + cout << " energy in em. component (GeV) : " << emEnergy_ / 1_GeV << endl + << " no. of em. particles removed : " << emCount_ << endl; if (discardInv_) - cout << " energy in inv. component (GeV) : " << invEnergy_ / 1_GeV << endl - << " no. of inv. particles removed : " << invCount_ << endl; + cout << " energy in inv. component (GeV) : " << invEnergy_ / 1_GeV << endl + << " no. of inv. particles removed : " << invCount_ << endl; cout << " energy below energy threshold (GeV): " << energy_ / 1_GeV << endl << " ******************************" << endl; } @@ -114,6 +114,6 @@ namespace corsika::process { invCount_ = 0; energy_ = 0_GeV; } - + } // namespace particle_cut } // namespace corsika::process diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc index 81f2efa7b6698c374829b543f6301d6675ee3fb0..d856aac96052a177aaf089e44de31043e3f3bef0 100644 --- a/Processes/Proposal/ContinuousProcess.cc +++ b/Processes/Proposal/ContinuousProcess.cc @@ -36,7 +36,8 @@ namespace corsika::process::proposal { // saved in the calc map by a key build out of a hash of composed of the component and // particle code. auto disp = PROPOSAL::make_displacement(c, true); - auto scatter = PROPOSAL::make_scattering("highland", particle[code], media.at(comp.hash())); + auto scatter = + PROPOSAL::make_scattering("highland", particle[code], media.at(comp.hash())); calc[std::make_pair(comp.hash(), code)] = std::make_tuple(std::move(disp), std::move(scatter)); } @@ -80,8 +81,7 @@ namespace corsika::process::proposal { // scattering auto vec = corsika::geometry::QuantityVector( final_dir.GetX() * E_f, final_dir.GetY() * E_f, final_dir.GetZ() * E_f); - vP.SetMomentum( - corsika::stack::MomentumVector(vP_dir.GetCoordinateSystem(), vec)); + vP.SetMomentum(corsika::stack::MomentumVector(vP_dir.GetCoordinateSystem(), vec)); } template <> diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc index 623556b909e401381b7ba09941380c30b6e6f4ff..064e146ea3cd245f03219f5ca51f880e20f41184 100644 --- a/Processes/Proposal/Interaction.cc +++ b/Processes/Proposal/Interaction.cc @@ -36,15 +36,15 @@ namespace corsika::process::proposal { // interpolate the crosssection for given media and energy cut. These may // take some minutes if you have to build the tables and cannot read the // from disk - auto c = p_cross->second(media.at(&comp), emCut_); + auto c = p_cross->second(media.at(comp.hash()), emCut_); // Look which interactions take place and build the corresponding // interaction and secondarie builder. The interaction integral will // interpolated too and saved in the calc map by a key build out of a hash // of composed of the component and particle code. auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c); - calc[std::make_pair(&comp, code)] = std::make_tuple( - PROPOSAL::make_secondaries(inter_types, particle[code], media.at(&comp)), + calc[std::make_pair(comp.hash(), code)] = std::make_tuple( + PROPOSAL::make_secondaries(inter_types, particle[code], media.at(comp.hash())), PROPOSAL::make_interaction(c, true)); } @@ -67,7 +67,8 @@ namespace corsika::process::proposal { // Read how much random numbers are required to calculate the secondaries. // Calculate the secondaries and deploy them on the corsika stack. - auto rnd = vector<double>(get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type)); + auto rnd = + vector<double>(get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type)); for (auto& it : rnd) it = distr(fRNG); auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm, vP.GetPosition().GetY() / 1_cm, @@ -78,8 +79,8 @@ namespace corsika::process::proposal { d.GetZ().magnitude()); auto loss = make_tuple(static_cast<int>(type), point, direction, v * vP.GetEnergy() / 1_MeV, 0.); - auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries(vP.GetEnergy() / 1_MeV, - loss, *comp_ptr, rnd); + auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries( + vP.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd); for (auto& s : sec) { auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV; auto vec = corsika::geometry::QuantityVector( diff --git a/Processes/Proposal/ProposalProcessBase.cc b/Processes/Proposal/ProposalProcessBase.cc index 3d1bd1e980987dcefc49813661172fee4433bb7e..240c59d77f9977233e67e4017002407dfd3c2b52 100644 --- a/Processes/Proposal/ProposalProcessBase.cc +++ b/Processes/Proposal/ProposalProcessBase.cc @@ -64,7 +64,7 @@ namespace corsika::process::proposal { "table directory. "); } } - + size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept { return p.first ^ std::hash<particles::Code>{}(p.second); } diff --git a/Processes/Proposal/ProposalProcessBase.h b/Processes/Proposal/ProposalProcessBase.h index 82d8330226f666fd1ed962749b2c7e29e7803f60..17c6b22d6c2142ca6ddd3b136ca47b71429dc6aa 100644 --- a/Processes/Proposal/ProposalProcessBase.h +++ b/Processes/Proposal/ProposalProcessBase.h @@ -1,3 +1,11 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + #pragma once #include <PROPOSAL/PROPOSAL.h> diff --git a/Processes/QGSJetII/QGSJetIIStack.h b/Processes/QGSJetII/QGSJetIIStack.h index 8e3c87ddc9581464428809e86f8cce1dde87623a..04685ab7cf0a3efc89f6c953b6c3f0354758003c 100644 --- a/Processes/QGSJetII/QGSJetIIStack.h +++ b/Processes/QGSJetII/QGSJetIIStack.h @@ -86,7 +86,7 @@ namespace corsika::process::qgsjetII { public: void SetParticleData(const int vID, const corsika::units::si::HEPEnergyType vE, const MomentumVector& vP, - const corsika::units::si::HEPMassType vM) { + const corsika::units::si::HEPMassType) { SetPID(vID); SetEnergy(vE); SetMomentum(vP); @@ -95,7 +95,7 @@ namespace corsika::process::qgsjetII { void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/, const int vID, const corsika::units::si::HEPEnergyType vE, const MomentumVector& vP, - const corsika::units::si::HEPMassType vM) { + const corsika::units::si::HEPMassType) { SetPID(vID); SetEnergy(vE); SetMomentum(vP);