IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 5bb7d683 authored by ralfulrich's avatar ralfulrich
Browse files

formating and copyrights

parent 67d74461
No related branches found
No related tags found
1 merge request!245Include proposal process rebase
Pipeline #2311 failed
Showing
with 72 additions and 69 deletions
......@@ -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();
}
......@@ -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");
......
......@@ -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");
......
......@@ -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
......@@ -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
......@@ -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++;
}
......
......@@ -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,
......
......@@ -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
......
......@@ -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() {
......
......@@ -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_;
......
......@@ -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
......@@ -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 <>
......
......@@ -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(
......
......@@ -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);
}
......
/*
* (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>
......
......@@ -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);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment