IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 8f5374ea authored by ralfulrich's avatar ralfulrich
Browse files

why you did not do clan-format?

parent c5380eaf
No related branches found
No related tags found
No related merge requests found
......@@ -252,10 +252,9 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true);
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::sibyll::Interaction sibyll(env);
......@@ -274,7 +273,7 @@ int main() {
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
// auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
auto sequence = p0 << pythia << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
......@@ -285,7 +284,7 @@ int main() {
stack.Clear();
const Code beamCode = Code::Proton;
const HEPMassType mass = particles::Proton::GetMass();
const HEPEnergyType E0 = 100_GeV;
const HEPEnergyType E0 = 100_GeV;
double theta = 0.;
double phi = 0.;
......@@ -305,10 +304,10 @@ int main() {
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns});
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns});
}
// define air shower object, run simulation
......
......@@ -8,5 +8,3 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
......@@ -107,7 +107,7 @@ namespace corsika::environment {
auto const& GetVolume() const { return *fGeoVolume; }
auto const& GetModelProperties() const { return *fModelProperties; }
auto const& HasModelProperties() const { return fModelProperties.get() != nullptr; }
template <typename ModelProperties, typename... Args>
......
......@@ -33,8 +33,8 @@ namespace corsika::stack {
<b>Important:</b> ParticleInterface must inherit from ParticleBase !
*/
template <typename>
class ParticleInterface;
template <typename>
class ParticleInterface;
/**
The Stack class provides (and connects) the main particle data storage machinery.
......@@ -98,8 +98,8 @@ namespace corsika::stack {
typedef StackDataType
StackImpl; ///< this is the type of the user-provided data structure
template <typename SI>
using PIType = ParticleInterface<SI>;
template <typename SI>
using PIType = ParticleInterface<SI>;
/**
* Via the StackIteratorInterface and ConstStackIteratorInterface
......
......@@ -8,5 +8,3 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
......@@ -8,5 +8,3 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
......@@ -62,8 +62,7 @@ namespace corsika::process::HadronicElasticModel {
avgCrossSection += CrossSection(s) * fractions[i];
}
std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb"
<< std::endl;
std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb" << std::endl;
return avgCrossSection;
}();
......
......@@ -32,13 +32,11 @@ using Track = Trajectory;
namespace corsika::process::pythia {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
Interaction::Interaction(environment::Environment const& env)
: fEnvironment(env) {}
: fEnvironment(env) {}
Interaction::~Interaction() {
cout << "Pythia::Interaction n=" << fCount << endl;
}
Interaction::~Interaction() { cout << "Pythia::Interaction n=" << fCount << endl; }
void Interaction::Init() {
......@@ -49,58 +47,58 @@ namespace corsika::process::pythia {
fPythia.readString("Print:quiet = on");
// TODO: proper process initialization for MinBias needed
fPythia.readString("HardQCD:all = on");
fPythia.readString("HardQCD:all = on");
fPythia.readString("ProcessLevel:resonanceDecays = off");
fPythia.init();
// any decays in pythia? if yes need to define which particles
if(fInternalDecays){
// define which particles are passed to corsika, i.e. which particles make it into history
// even very shortlived particles like charm or pi0 are of interest here
const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Pi0,
particles::Code::KMinus, particles::Code::KPlus,
particles::Code::K0Long, particles::Code::K0Short,
particles::Code::SigmaPlus, particles::Code::SigmaMinus,
particles::Code::Lambda0,
particles::Code::Xi0, particles::Code::XiMinus,
particles::Code::OmegaMinus,
particles::Code::DPlus, particles::Code::DMinus, particles::Code::D0, particles::Code::D0Bar};
Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
if (fInternalDecays) {
// define which particles are passed to corsika, i.e. which particles make it into
// history even very shortlived particles like charm or pi0 are of interest here
const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::Pi0, particles::Code::KMinus,
particles::Code::KPlus, particles::Code::K0Long,
particles::Code::K0Short, particles::Code::SigmaPlus,
particles::Code::SigmaMinus, particles::Code::Lambda0,
particles::Code::Xi0, particles::Code::XiMinus,
particles::Code::OmegaMinus, particles::Code::DPlus,
particles::Code::DMinus, particles::Code::D0,
particles::Code::D0Bar};
Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
}
// basic initialization of cross section routines
fSigma.init( &fPythia.info, fPythia.settings, &fPythia.particleData, &fPythia.rndm );
fSigma.init(&fPythia.info, fPythia.settings, &fPythia.particleData, &fPythia.rndm);
fInitialized = true;
}
}
void Interaction::SetParticleListStable(const std::vector<particles::Code> particleList) {
for (auto p : particleList)
Interaction::SetStable( p );
void Interaction::SetParticleListStable(
const std::vector<particles::Code> particleList) {
for (auto p : particleList) Interaction::SetStable(p);
}
void Interaction::SetUnstable(const particles::Code pCode) {
cout << "Pythia::Interaction: setting " << pCode << " unstable.." << endl;
fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , true);
fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
}
void Interaction::SetStable(const particles::Code pCode) {
cout << "Pythia::Interaction: setting " << pCode << " stable.." << endl;
fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , false);
fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), false);
}
void Interaction::ConfigureLabFrameCollision( const particles::Code BeamId, const particles::Code TargetId,
const units::si::HEPEnergyType BeamEnergy )
{
void Interaction::ConfigureLabFrameCollision(
const particles::Code BeamId, const particles::Code TargetId,
const units::si::HEPEnergyType BeamEnergy) {
using namespace units::si;
// Pythia configuration of the current event
// very clumsy. I am sure this can be done better..
// set beam
// beam id for pythia
auto const pdgBeam = static_cast<int>(particles::GetPDG(BeamId));
......@@ -110,8 +108,8 @@ namespace corsika::process::pythia {
// set target
auto pdgTarget = static_cast<int>(particles::GetPDG(TargetId));
// replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
if(TargetId == particles::Code::Hydrogen)
pdgTarget = static_cast<int>( particles::GetPDG(particles::Code::Proton));
if (TargetId == particles::Code::Hydrogen)
pdgTarget = static_cast<int>(particles::GetPDG(particles::Code::Proton));
std::stringstream stTarget;
stTarget << "Beams:idB = " << pdgTarget;
fPythia.readString(stTarget.str());
......@@ -128,14 +126,15 @@ namespace corsika::process::pythia {
fPythia.init();
}
bool Interaction::CanInteract(const corsika::particles::Code pCode)
{
return pCode == corsika::particles::Code::Proton || pCode == corsika::particles::Code::Neutron
|| pCode == corsika::particles::Code::AntiProton || pCode == corsika::particles::Code::AntiNeutron
|| pCode == corsika::particles::Code::PiMinus || pCode == corsika::particles::Code::PiPlus;
}
bool Interaction::CanInteract(const corsika::particles::Code pCode) {
return pCode == corsika::particles::Code::Proton ||
pCode == corsika::particles::Code::Neutron ||
pCode == corsika::particles::Code::AntiProton ||
pCode == corsika::particles::Code::AntiNeutron ||
pCode == corsika::particles::Code::PiMinus ||
pCode == corsika::particles::Code::PiPlus;
}
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
......@@ -143,27 +142,27 @@ namespace corsika::process::pythia {
using namespace units::si;
// interaction possible in pythia?
if( TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen ){
if( CanInteract(BeamId) && ValidCoMEnergy(CoMenergy) ){
// input particle PDG
auto const pdgCodeBeam = static_cast<int>( particles::GetPDG( BeamId ));
auto const pdgCodeTarget = static_cast<int>( particles::GetPDG( TargetId ));
const double ecm = CoMenergy / 1_GeV;
// calculate cross section
fSigma.calc( pdgCodeBeam, pdgCodeTarget, ecm);
if(fSigma.hasSigmaTot()){
const double sigEla = fSigma.sigmaEl();
const double sigProd = fSigma.sigmaTot() - sigEla;
return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
} else
throw std::runtime_error("pythia cross section init failed");
if (TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen) {
if (CanInteract(BeamId) && ValidCoMEnergy(CoMenergy)) {
// input particle PDG
auto const pdgCodeBeam = static_cast<int>(particles::GetPDG(BeamId));
auto const pdgCodeTarget = static_cast<int>(particles::GetPDG(TargetId));
const double ecm = CoMenergy / 1_GeV;
// calculate cross section
fSigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
if (fSigma.hasSigmaTot()) {
const double sigEla = fSigma.sigmaEl();
const double sigProd = fSigma.sigmaTot() - sigEla;
return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
} else
throw std::runtime_error("pythia cross section init failed");
} else {
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
std::numeric_limits<double>::infinity() * 1_mb);
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
std::numeric_limits<double>::infinity() * 1_mb);
}
} else {
throw std::runtime_error("invalid target for pythia");
......@@ -206,7 +205,7 @@ namespace corsika::process::pythia {
<< " beam pid:" << p.GetPID() << endl;
// TODO: move limits into variables
if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM) ) {
if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM)) {
// get target from environment
/*
......@@ -235,19 +234,19 @@ namespace corsika::process::pythia {
elaCrossSection; // ONLY TO AVOID COMPILER WARNING
cout << "Interaction: IntLength: pythia return (mb): "
<< productionCrossSection / 1_mb << endl
<< "Interaction: IntLength: weight : " << w[i] << endl;
<< productionCrossSection / 1_mb << endl
<< "Interaction: IntLength: weight : " << w[i] << endl;
weightedProdCrossSection += w[i] * productionCrossSection;
}
cout << "Interaction: IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mb << endl
<< "Interaction: IntLength: average mass number: "
<< mediumComposition.GetAverageMassNumber() << endl;
<< "Interaction: IntLength: average mass number: "
<< mediumComposition.GetAverageMassNumber() << endl;
// calculate interaction length in medium
GrammageType const int_length =
mediumComposition.GetAverageMassNumber() * units::constants::u / weightedProdCrossSection;
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
units::constants::u / weightedProdCrossSection;
cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< endl;
......@@ -257,7 +256,7 @@ namespace corsika::process::pythia {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
/**
In this function PYTHIA is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
......@@ -271,7 +270,6 @@ namespace corsika::process::pythia {
using namespace units::si;
using namespace geometry;
const auto corsikaBeamId = p.GetPID();
cout << "Pythia::Interaction: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
......@@ -283,7 +281,7 @@ namespace corsika::process::pythia {
}
if (process::pythia::Interaction::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
......@@ -353,8 +351,7 @@ namespace corsika::process::pythia {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, sigEla] =
GetCrossSection(corsikaBeamId, targetId, Ecm);
const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy =
sigEla; // to avoid not used warning in array binding
......@@ -363,56 +360,60 @@ namespace corsika::process::pythia {
const auto corsikaTargetId =
mediumComposition.SampleTarget(cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << corsikaTargetId << endl;
if (corsikaTargetId != particles::Code::Hydrogen && corsikaTargetId != particles::Code::Neutron
&& corsikaTargetId != particles::Code::Proton )
throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
if (corsikaTargetId != particles::Code::Hydrogen &&
corsikaTargetId != particles::Code::Neutron &&
corsikaTargetId != particles::Code::Proton)
throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl;
if (eProjectileLab < 8.5_GeV || !ValidCoMEnergy(Ecm)) {
cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
throw std::runtime_error("energy too low for PYTHIA");
} else {
} else {
fCount++;
ConfigureLabFrameCollision( corsikaBeamId, corsikaTargetId, eProjectileLab );
// create event in pytia
if(!fPythia.next())
throw std::runtime_error("Pythia::DoInteraction: failed!");
ConfigureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
// create event in pytia
if (!fPythia.next()) throw std::runtime_error("Pythia::DoInteraction: failed!");
// link to pythia stack
Pythia8::Event& event = fPythia.event;
// link to pythia stack
Pythia8::Event& event = fPythia.event;
// print final state
event.list();
event.list();
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
for (int i = 0; i < event.size(); ++i){
Pythia8::Particle &pp = event[i];
for (int i = 0; i < event.size(); ++i) {
Pythia8::Particle& pp = event[i];
// skip particles that have decayed in pythia
if (!pp.isFinal()) continue;
if (!pp.isFinal()) continue;
auto const pyId = particles::ConvertFromPDG(static_cast<particles::PDGCode>(pp.id()));
auto const pyId =
particles::ConvertFromPDG(static_cast<particles::PDGCode>(pp.id()));
const MomentumVector pyPlab(rootCS, {pp.px()*1_GeV, pp.py()*1_GeV, pp.pz()*1_GeV});
const MomentumVector pyPlab(
rootCS, {pp.px() * 1_GeV, pp.py() * 1_GeV, pp.pz() * 1_GeV});
HEPEnergyType const pyEn = pp.e() * 1_GeV;
// add to corsika stack
auto pnew = p.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig, tOrig});
geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig,
tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
}
cout << "conservation (all GeV): " << "Elab_final=" << Elab_final / 1_GeV
cout << "conservation (all GeV): "
<< "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
// delete current particle
......
......@@ -37,8 +37,8 @@ namespace corsika::process::pythia {
void Init();
void SetParticleListStable(const std::vector<particles::Code>);
void SetUnstable(const corsika::particles::Code );
void SetStable(const corsika::particles::Code );
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
bool WasInitialized() { return fInitialized; }
bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
......@@ -47,8 +47,9 @@ namespace corsika::process::pythia {
}
bool CanInteract(const corsika::particles::Code);
void ConfigureLabFrameCollision(const corsika::particles::Code, const corsika::particles::Code,
const corsika::units::si::HEPEnergyType);
void ConfigureLabFrameCollision(const corsika::particles::Code,
const corsika::particles::Code,
const corsika::units::si::HEPEnergyType);
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
......@@ -68,7 +69,7 @@ namespace corsika::process::pythia {
private:
corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
Pythia8::Pythia fPythia;
Pythia8::SigmaTotal fSigma;
const bool fInternalDecays = true;
......
......@@ -9,10 +9,9 @@
* the license.
*/
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/pythia/Interaction.h>
#include <Pythia8/Pythia.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/pythia/Interaction.h>
#include <corsika/random/RNGManager.h>
......@@ -96,7 +95,7 @@ TEST_CASE("Pythia", "[processes]") {
using namespace corsika;
using namespace corsika::units::si;
TEST_CASE("pythia process"){
TEST_CASE("pythia process") {
// setup environment, geometry
environment::Environment env;
......@@ -110,7 +109,8 @@ TEST_CASE("pythia process"){
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Hydrogen}, std::vector<float>{1.}));
std::vector<particles::Code>{particles::Code::Hydrogen},
std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
......@@ -152,7 +152,7 @@ TEST_CASE("pythia process"){
}
SECTION("pythia interaction") {
setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
......@@ -164,14 +164,12 @@ TEST_CASE("pythia process"){
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::PiPlus, E0, plab, pos, 0_ns});
process::pythia::Interaction model(env);
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoInteraction(particle,
stack);
stack);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
model.GetInteractionLength(particle, track);
}
}
......@@ -290,8 +290,8 @@ namespace corsika::process::sibyll {
weightedProdCrossSection += w[i] * productionCrossSection;
}
cout << "NuclearInteraction: "
<< "IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mb << endl;
<< "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
<< endl;
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
......
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