IAP GITLAB

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

implemented DoInteraction for nuclei

parent 1f3a1ac2
No related branches found
No related tags found
1 merge request!65Resolve "add sibyll process NuclearInteraction"
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 24
− 10
View file @ f64c7d4f
...@@ -125,6 +125,14 @@ public: ...@@ -125,6 +125,14 @@ public:
is_inv = true; is_inv = true;
break; break;
case Code::Neutron:
is_inv = true;
break;
case Code::AntiNeutron:
is_inv = true;
break;
default: default:
break; break;
} }
...@@ -235,33 +243,39 @@ int main() { ...@@ -235,33 +243,39 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true); stack_inspector::StackInspector<setup::Stack> p0(true);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
// corsika::process::sibyll::Interaction sibyll(env); corsika::process::sibyll::Interaction sibyll(env);
corsika::process::sibyll::NuclearInteraction sibyll(env); corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
corsika::process::sibyll::Decay decay; corsika::process::sibyll::Decay decay;
ProcessCut cut(8_GeV); ProcessCut cut(8_GeV);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); // 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"); corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list // 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() // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n"; // << "\n";
// setup particle stack, and add primary particle // setup particle stack, and add primary particle
setup::Stack stack; setup::Stack stack;
stack.Clear(); stack.Clear();
const HEPEnergyType E0 = const Code beamCode = Code::Carbon;
100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later) const HEPEnergyType E0 = 100_TeV;
// 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.; double theta = 0.;
double phi = 0.; double phi = 0.;
{ {
auto particle = stack.NewParticle(); auto particle = stack.NewParticle();
particle.SetPID(Code::Helium); particle.SetPID(beamCode);
HEPMomentumType P0 = sqrt(E0 * E0 - Helium::GetMass() * Helium::GetMass()); 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) { auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi), return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta)); -ptot * cos(theta));
......
Framework/Particles/NuclearData.xml
+ 3
− 0
View file @ f64c7d4f
...@@ -45,6 +45,9 @@ ...@@ -45,6 +45,9 @@
<particle id="1000100220" name="neon" A="22" Z="10" > <particle id="1000100220" name="neon" A="22" Z="10" >
</particle> </particle>
<particle id="1000160330" name="sulphur" A="33" Z="16" >
</particle>
<particle id="1000180400" name="argon" A="40" Z="18" > <particle id="1000180400" name="argon" A="40" Z="18" >
</particle> </particle>
......
Processes/Sibyll/NuclearInteraction.h
+ 252
− 63
View file @ f64c7d4f
...@@ -46,7 +46,7 @@ namespace corsika::process::sibyll { ...@@ -46,7 +46,7 @@ namespace corsika::process::sibyll {
using std::endl; using std::endl;
// initialize hadronic interaction module // initialize hadronic interaction module
sibyll_ini_(); //sibyll_ini_();
// initialize nuclib // initialize nuclib
nuc_nuc_ini_(); nuc_nuc_ini_();
...@@ -185,6 +185,9 @@ namespace corsika::process::sibyll { ...@@ -185,6 +185,9 @@ namespace corsika::process::sibyll {
template <typename Particle, typename Stack> template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) { corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
using namespace corsika::units; using namespace corsika::units;
using namespace corsika::utl; using namespace corsika::utl;
using namespace corsika::units::si; using namespace corsika::units::si;
...@@ -192,13 +195,19 @@ namespace corsika::process::sibyll { ...@@ -192,13 +195,19 @@ namespace corsika::process::sibyll {
using std::cout; using std::cout;
using std::endl; using std::endl;
const auto corsikaBeamId = p.GetPID(); const auto corsikaProjId = p.GetPID();
cout << "NuclearInteraction: DoInteraction: called with:" << corsikaBeamId cout << "NuclearInteraction: DoInteraction: called with:" << corsikaProjId << endl
<< "NuclearInteraction: projectile mass: " << corsika::particles::GetMass(corsikaProjId) / 1_GeV
<< endl; << endl;
if(!IsNucleus(corsikaBeamId)) if(!IsNucleus(corsikaProjId)){
// this should not happen
throw std::runtime_error("Non nuclear projectile in NUCLIB!");
return process::EProcessReturn::eOk; return process::EProcessReturn::eOk;
}
fCount++;
const CoordinateSystem& rootCS = const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
...@@ -206,43 +215,75 @@ namespace corsika::process::sibyll { ...@@ -206,43 +215,75 @@ namespace corsika::process::sibyll {
Point pOrig = p.GetPosition(); Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime(); TimeType tOrig = p.GetTime();
// beam nucleon number cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
const int kABeam = GetNucleusA(corsikaBeamId); << endl;
// TODO: verify this number !!! cout << "Interaction: time: " << tOrig << endl;
if(kABeam>56)
throw std::runtime_error("Beam nucleus too large for SIBYLL!"); // projectile nucleon number
const int kAProj = GetNucleusA(corsikaProjId);
if(kAProj>56)
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics // kinematics
// define projectile NUCLEON // define projectile nucleus
HEPEnergyType const eProjectileLab = p.GetEnergy() / kABeam; HEPEnergyType const eProjectileLab = p.GetEnergy();
auto const pProjectileLab = p.GetMomentum() / kABeam; auto const pProjectileLab = p.GetMomentum();
const FourVector PprojLab(eProjectileLab, pProjectileLab); const FourVector PprojLab(eProjectileLab, pProjectileLab);
cout << "NuclearInteraction: ebeam lab: " << eProjectileLab / 1_GeV << endl cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl
<< "NuclearInteraction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV << "NuclearInteraction: pProj lab: " << pProjectileLab.GetComponents() / 1_GeV
<< endl; << endl;
// define projectile nucleon
HEPEnergyType const eProjectileNucLab = p.GetEnergy() / kAProj;
auto const pProjectileNucLab = p.GetMomentum() / kAProj;
const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab);
cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV << endl
<< "NuclearInteraction: pProjNucleon lab: " << pProjectileNucLab.GetComponents() / 1_GeV
<< endl;
// define target // define target
// always a nucleon // always a nucleon
// for Sibyll is always a single nucleon // for Sibyll is always a single nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() + auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass()); corsika::particles::Neutron::GetMass());
// target is always at rest // target is always at rest
const auto eTargetLab = 0_GeV + nucleon_mass; const auto eTargetNucLab = 0_GeV + nucleon_mass;
const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV); const auto pTargetNucLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab); const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab);
cout << "NuclearInteraction: etarget lab: " << eTargetLab / 1_GeV << endl cout << "NuclearInteraction: etarget lab: " << eTargetNucLab / 1_GeV << endl
<< "NuclearInteraction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV << "NuclearInteraction: ptarget lab: " << pTargetNucLab.GetComponents() / 1_GeV
<< endl; << endl;
// center-of-mass energy in nucleon-nucleon frame // center-of-mass energy in nucleon-nucleon frame
auto const Ptot4 = PtargLab + PprojLab; auto const PtotNN4 = PtargNucLab + PprojNucLab;
HEPEnergyType Ecm = Ptot4.GetNorm(); HEPEnergyType EcmNN = PtotNN4.GetNorm();
cout << "NuclearInteraction: nuc-nuc cm energy: " << Ecm / 1_GeV << endl; cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
// define boost to NUCLEON-NUCLEON frame
COMBoost const boost(PprojNucLab, nucleon_mass);
// boost projecticle
auto const PprojNucCoM = boost.toCoM(PprojNucLab);
// boost target
auto const PtargNucCoM = boost.toCoM(PtargNucLab);
cout << "Interaction: ebeam CoM: " << PprojNucCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: pbeam CoM: "
<< PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargNucCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: ptarget CoM: "
<< PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
const auto beamId = corsika::particles::Proton::GetCode();
// sample target nucleon number // sample target nucleon number
//
// proton stand-in for nucleon
const auto beamId = corsika::particles::Proton::GetCode();
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig); const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
const auto& mediumComposition = const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition(); currentNode->GetModelProperties().GetNuclearComposition();
...@@ -256,7 +297,7 @@ namespace corsika::process::sibyll { ...@@ -256,7 +297,7 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) { for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i]; auto const targetId = compVec[i];
const auto [sigProd, nNuc] = GetCrossSection(beamId,targetId,Ecm); const auto [sigProd, nNuc] = GetCrossSection(beamId,targetId,EcmNN);
cross_section_of_components[i] = sigProd; cross_section_of_components[i] = sigProd;
} }
...@@ -282,59 +323,207 @@ namespace corsika::process::sibyll { ...@@ -282,59 +323,207 @@ namespace corsika::process::sibyll {
cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "<< endl; cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "<< endl;
// get nucleon-nucleon cross section // get nucleon-nucleon cross section
// (needed to determine number of scatterings) // (needed to determine number of scatterings)
double sigProd, sigEla, dum1, dum2, dum3, dum4; //
// TODO: this is an explicit dependence on sibyll
// should be changed to a call to GetCrossSection() of the had. int. implementation
// this also means GetCrossSection has to return the elastic cross section as well
double sigProd, sigEla, dum1, dum2, dum3;
double dumdif[3]; double dumdif[3];
const double sqsnuc = Ecm / 1_GeV; const double sqsnuc = EcmNN / 1_GeV;
const int iBeam = 1; const int iBeam = 1;
sib_sigma_hp_(iBeam, sqsnuc, dum1, sigEla, sigProd, dumdif, dum3, dum4); // read hadron-proton cross section table (input: hadron id, CoMenergy, output: cross sections)
sib_sigma_hp_(iBeam, sqsnuc, dum1, sigEla, sigProd, dumdif, dum2, dum3);
// sample number of interactions (output in common cnucms) // sample number of interactions (only input variables, output in common cnucms)
// nuclear multiple scattering according to glauber (r.i.p.) // nuclear multiple scattering according to glauber (r.i.p.)
int_nuc_( kATarget, kABeam, sigProd, sigEla); int_nuc_( kATarget, kAProj, sigProd, sigEla);
cout << "number of wounded nucleons in target : " << cnucms_.na << endl cout << "number of nucleons in target : " << kATarget << endl
<< "number of wounded nucleons in target : " << cnucms_.na << endl
<< "number of nucleons in projectile : " << kAProj << endl
<< "number of wounded nucleons in project. : " << cnucms_.nb << endl << "number of wounded nucleons in project. : " << cnucms_.nb << endl
<< "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl << "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl
<< "number of elastic nucleons in target : " << cnucms_.nael << endl << "number of elastic nucleons in target : " << cnucms_.nael << endl
<< "number of elastic nucleons in project. : " << cnucms_.nbel << endl; << "number of elastic nucleons in project. : " << cnucms_.nbel << endl
<< "impact parameter: " << cnucms_.b << endl;
throw std::runtime_error(" end now");
// calculate fragmentation // calculate fragmentation
// call FRAGM (IAT,IAP, NW,B, NF, IAF) cout << "calculating nuclear fragments.." << endl;
// fragm_(kATarget, kABeam, NBT, B, NF, IAF) // number of interactions
// include elastic
const int nElasticNucleons = cnucms_.nbel;
const int nInelNucleons = cnucms_.nb;
const int nIntProj = nInelNucleons + nElasticNucleons;
const double impactPar = cnucms_.b; // only needed to avoid passing common var.
int nFragments;
// number of fragments is limited to 60
int AFragments[60];
// call fragmentation routine
// input: target A, projectile A, number of int. nucleons in projectile, impact parameter (fm)
// output: nFragments, AFragments
// in addition the momenta ar stored in pf in common fragments, neglected
fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :)
if(nFragments>60)
throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
cout << "number of fragments: " << nFragments << endl;
for(int j=0; j<nFragments; ++j)
cout << "fragment: " << j << " A=" << AFragments[j]
<< " px=" << fragments_.ppp[j][0]
<< " py=" << fragments_.ppp[j][1]
<< " pz=" << fragments_.ppp[j][2]
<< endl;
// bookeeping accross nucleon-nucleon interactions
MomentumVector Plab_all(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_all = 0_GeV;
/*
C. INPUT: IAP = mass of incident nucleus
C. IAT = mass of target nucleus
C. NW = number of wounded nucleons in the beam nucleus
C. B = impact parameter in the interaction
C.
C. OUTPUT : NF = number of fragments of the spectator nucleus
C. IAF(1:NF) = mass number of each fragment
C. PF(3,60) in common block /FRAGMENTS/ contains
C. the three momentum components (MeV/c) of each
C. fragment in the projectile frame
C..............................................................
*/
// put spectators on intermediate stack auto Nucleus = []( int iA ){
//Stack nucs; // int znew = iA / 2.15 + 0.7;
corsika::particles::Code pCode;
switch( iA ){
case 2:
pCode = corsika::particles::Deuterium::GetCode();
break;
case 3:
pCode = corsika::particles::Tritium::GetCode();
break;
case 4:
pCode = corsika::particles::Helium::GetCode();
break;
case 12:
pCode = corsika::particles::Carbon::GetCode();
break;
case 14:
pCode = corsika::particles::Nitrogen::GetCode();
break;
case 16:
pCode = corsika::particles::Oxygen::GetCode();
break;
case 33:
pCode = corsika::particles::Sulphur::GetCode();
break;
case 40:
pCode = corsika::particles::Argon::GetCode();
break;
default:
pCode = corsika::particles::Proton::GetCode();
}
return pCode;
};
// add elastic nucleons to stack // put nuclear fragments on corsika stack
for(int j=0; j<nFragments; ++j){
// add inelastic interactions auto pnew = s.NewParticle();
// here we need the nucleonNumber to corsika Id conversion
// A = AFragments[j]
// pnew.SetPID( corsika::particles::GetCode( corsika::particles::Nucleus(A,Z) ) );
auto pCode = Nucleus( AFragments[j] );
pnew.SetPID(pCode);
// CORSIKA 7 way
// spectators inherit momentum from original projectile
const double mass_ratio = corsika::particles::GetMass( pCode ) / corsika::particles::GetMass( corsikaProjId );
auto const Plab = PprojLab * mass_ratio;
pnew.SetEnergy(Plab.GetTimeLikeComponent());
pnew.SetMomentum(Plab.GetSpaceLikeComponents());
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
}
// add elastic nucleons to corsika stack
for(int j=0; j<nElasticNucleons; ++j){
auto pnew = s.NewParticle();
// TODO: sample proton or neutron
auto pCode = corsika::particles::Proton::GetCode();
pnew.SetPID( pCode );
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
const double mass_ratio = corsika::particles::GetMass( pCode ) / corsika::particles::GetMass( corsikaProjId );
auto const Plab = PprojLab * mass_ratio;
pnew.SetEnergy(Plab.GetTimeLikeComponent());
pnew.SetMomentum(Plab.GetSpaceLikeComponents());
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
}
// move particles to corsika stack // add inelastic interactions
// boost for(int j=0; j<nInelNucleons; ++j){
// create nucleon-nucleus inelastic interaction
// TODO: switch to DoInteraction() of had. interaction implementation
// for now use SIBYLL directly
const int kBeamCode = process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode());
cout << "creating interaction no. "<< j << endl;
sibyll_( kBeamCode, kATarget, sqsnuc );
decsib_();
// print final state
int print_unit = 6;
sib_list_(print_unit);
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// skip particles that have decayed in Sibyll
if (psib.HasDecayed()) continue;
// // transform energy to lab. frame
auto const pCoM = psib.GetMomentum();
HEPEnergyType const eCoM = psib.GetEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
auto pnew = s.NewParticle();
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
pnew.SetEnergy(Plab.GetTimeLikeComponent());
pnew.SetMomentum(Plab.GetSpaceLikeComponents());
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
Plab_all += Plab.GetSpaceLikeComponents();
Elab_all += Plab.GetTimeLikeComponent();
}
cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV
<< endl
<< "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
<< endl;
}
cout << "accross all nucleon interactions: Etot lab: " << Elab_all / 1_GeV << endl
<< "accross all nucleon interactions: Ptot lab: " << (Plab_all / 1_GeV).GetComponents()
<< endl;
// add proton instead //throw std::runtime_error(" stop here");
auto pnew = s.NewParticle();
pnew.SetPID( corsika::particles::Proton::GetCode() );
pnew.SetMomentum( p.GetMomentum() );
pnew.SetEnergy( p.GetEnergy() );
// delete current particle // delete current particle
p.Delete(); p.Delete();
......
Processes/Sibyll/sibyll2.3c.h
+ 23
− 11
View file @ f64c7d4f
...@@ -83,8 +83,15 @@ extern struct { ...@@ -83,8 +83,15 @@ extern struct {
int jja[56], jjb[56], jjint[56][56], jjael[56], jjbel[56]; int jja[56], jjb[56], jjint[56][56], jjael[56], jjbel[56];
} cnucms_; } cnucms_;
/*
nuclib common, nuclear FRAGMENTS
COMMON /FRAGMENTS/ PPP(3,60)
*/
extern struct {
double ppp[60][3];
} fragments_;
// lund random generator setup // lund random generator setup
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_; // extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
...@@ -94,16 +101,7 @@ void sibyll_(const int&, const int&, const double&); ...@@ -94,16 +101,7 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll // subroutine to initiate sibyll
void sibyll_ini_(); void sibyll_ini_();
// 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&, int&, double&, int&, int&);
// subroutine to SET DECAYS // subroutine to SET DECAYS
void dec_ini_(); void dec_ini_();
...@@ -130,5 +128,19 @@ double get_sibyll_mass2(int&); ...@@ -130,5 +128,19 @@ double get_sibyll_mass2(int&);
// phojet random generator setup // phojet random generator setup
void pho_rndin_(int&, int&, int&, int&); void pho_rndin_(int&, int&, int&, int&);
// 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*);
} }
#endif #endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum