IAP GITLAB

Skip to content
Snippets Groups Projects
Commit c4b82ab0 authored by Felix Riehn's avatar Felix Riehn
Browse files

clean and clang

parent 3c4f8281
No related branches found
No related tags found
1 merge request!65Resolve "add sibyll process NuclearInteraction"
......@@ -250,14 +250,15 @@ int main() {
ProcessCut cut(20_GeV);
// corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// corsika::process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env);
// corsika::process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
//auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
......@@ -265,15 +266,16 @@ int main() {
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Proton;
const HEPEnergyType E0 = 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
const HEPEnergyType E0 =
100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.;
double phi = 0.;
{
auto elab2plab = []( HEPEnergyType Elab, HEPMassType m){
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab( E0, corsika::particles::GetMass( beamCode ) );
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab(E0, corsika::particles::GetMass(beamCode));
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
......
......@@ -67,7 +67,7 @@ namespace corsika::stack {
protected:
/** @name Access to underlying stack data
@{
@{
*/
auto& GetStackData() { return GetIterator().GetStackData(); }
const auto& GetStackData() const { return GetIterator().GetStackData(); }
......@@ -75,7 +75,7 @@ namespace corsika::stack {
const auto& GetStack() const { return GetIterator().GetStack(); }
///@}
/**
/**
* return the index number of the underlying iterator object
*/
int GetIndex() const { return GetIterator().GetIndex(); }
......
......@@ -12,7 +12,7 @@
#ifndef _include_Stack_h__
#define _include_Stack_h__
#include <corsika/stack/StackIteratorInterface.h>
#include <corsika/stack/StackIteratorInterface.h>
#include <stdexcept>
......@@ -29,7 +29,7 @@ namespace corsika::stack {
Important: ParticleInterface must inherit from ParticleBase !
*/
template <typename>
class ParticleInterface; // forward decl
......
......@@ -112,8 +112,8 @@ namespace corsika::stack {
}
public:
/** @name Iterator interface
@{
/** @name Iterator interface
@{
*/
StackIteratorInterface& operator++() {
++fIndex;
......
......@@ -32,7 +32,7 @@ namespace corsika::process::sibyll {
int fCount = 0;
int fNucCount = 0;
bool fInitialized = false;
public:
Interaction(corsika::environment::Environment const& env)
: fEnvironment(env) {}
......@@ -48,17 +48,19 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize Sibyll
if(!fInitialized){
sibyll_ini_();
fInitialized = true;
if (!fInitialized) {
sibyll_ini_();
fInitialized = true;
}
}
bool WasInitialized(){ return fInitialized;}
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType, int> GetCrossSection(
const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy) {
bool WasInitialized() { return fInitialized; }
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
int>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy) {
using namespace corsika::units::si;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
......@@ -70,14 +72,14 @@ namespace corsika::process::sibyll {
if (iTarget > 18 || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == corsika::particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
iTarget = 1;
iTarget = 1;
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity();
}
return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
}
......@@ -192,12 +194,12 @@ namespace corsika::process::sibyll {
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl;
if(corsika::particles::IsNucleus(corsikaBeamId)){
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
if (corsika::particles::IsNucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
}
if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
......@@ -273,7 +275,8 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, sigEla, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
const auto [sigProd, sigEla, nNuc] =
GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
int ideleteme = nNuc; // to avoid not used warning in array binding
ideleteme = ideleteme; // to avoid not used warning in array binding
......@@ -306,7 +309,7 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << std::endl;
throw std::runtime_error("energy too low for SIBYLL");
throw std::runtime_error("energy too low for SIBYLL");
// p.Delete(); delete later... different process
} else {
fCount++;
......
This diff is collapsed.
......@@ -14,42 +14,37 @@
extern "C" {
// nuclib common, NUClear Multiple Scattering
/*
COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
+ ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
+ ,JJAEL(IAMAX), JJBEL(IAMAX)
*/
extern struct {
double b, bmax;
int ntry, na, nb, ni, nael, nbel;
int jja[56], jjb[56], jjint[56][56], jjael[56], jjbel[56];
} cnucms_;
/*
nuclib common, nuclear FRAGMENTS
COMMON /FRAGMENTS/ PPP(3,60)
*/
extern struct {
double ppp[60][3];
} fragments_;
// NUCLIB
// subroutine to initiate nuclib
void nuc_nuc_ini_();
// subroutine to sample nuclear interaction structure
void int_nuc_( const int&, const int&, const double&, const double&);
// subroutine to sample nuclear fragments
void fragm_(const int&, const int&, const int&, const double&, int&, int*);
void signuc_(const int&, const double&, double&);
// nuclib common, NUClear Multiple Scattering
/*
COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
+ ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
+ ,JJAEL(IAMAX), JJBEL(IAMAX)
*/
extern struct {
double b, bmax;
int ntry, na, nb, ni, nael, nbel;
int jja[56], jjb[56], jjint[56][56], jjael[56], jjbel[56];
} cnucms_;
/*
nuclib common, nuclear FRAGMENTS
COMMON /FRAGMENTS/ PPP(3,60)
*/
extern struct { double ppp[60][3]; } fragments_;
// NUCLIB
// subroutine to initiate nuclib
void nuc_nuc_ini_();
// subroutine to sample nuclear interaction structure
void int_nuc_(const int&, const int&, const double&, const double&);
// subroutine to sample nuclear fragments
void fragm_(const int&, const int&, const int&, const double&, int&, int*);
void signuc_(const int&, const double&, double&);
}
#endif
......@@ -70,7 +70,6 @@ extern struct {
int lun;
} s_debug_;
// lund random generator setup
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
......@@ -80,7 +79,6 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
// subroutine to SET DECAYS
void dec_ini_();
......@@ -96,7 +94,7 @@ void decsib_();
// interaction length
// double fpni_(double&, int&);
void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&, double&);
void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&, double&);
void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
double&);
......@@ -107,7 +105,5 @@ double get_sibyll_mass2(int&);
// phojet random generator setup
void pho_rndin_(int&, int&, int&, int&);
}
#endif
......@@ -136,12 +136,12 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
Interaction hmodel(env);
......
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