diff --git a/Processes/QGSJetII/Interaction.cc b/Processes/QGSJetII/Interaction.cc
index ddd96c4bec45ee38a26792020e72e796dbb59c4d..b3a3ad690042e9b4609923d33496334860a35bf2 100644
--- a/Processes/QGSJetII/Interaction.cc
+++ b/Processes/QGSJetII/Interaction.cc
@@ -14,22 +14,22 @@
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/process/qgsjetII/ParticleConversion.h>
-#include <corsika/process/qgsjetII/QGSJetIIStack.h>
#include <corsika/process/qgsjetII/QGSJetIIFragmentsStack.h>
+#include <corsika/process/qgsjetII/QGSJetIIStack.h>
#include <corsika/process/qgsjetII/qgsjet-II-04.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h>
-#include <tuple>
-#include <string>
#include <sstream>
+#include <string>
+#include <tuple>
using std::cout;
using std::endl;
-using std::tuple;
-using std::string;
using std::ostringstream;
+using std::string;
+using std::tuple;
using namespace corsika;
using namespace corsika::setup;
@@ -39,18 +39,17 @@ using Track = Trajectory;
namespace corsika::process::qgsjetII {
- Interaction::Interaction(const string& dataPath) : data_path_(dataPath) {
- if (dataPath=="") {
+ Interaction::Interaction(const string& dataPath)
+ : data_path_(dataPath) {
+ if (dataPath == "") {
if (std::getenv("CORSIKA_DATA")) {
- data_path_ = string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
- cout << "Searching for QGSJetII data tables in " << data_path_ << endl;
+ data_path_ = string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
+ cout << "Searching for QGSJetII data tables in " << data_path_ << endl;
}
}
}
- Interaction::~Interaction() {
- cout << "QgsjetII::Interaction n=" << count_ << endl;
- }
+ Interaction::~Interaction() { cout << "QgsjetII::Interaction n=" << count_ << endl; }
void Interaction::Init() {
@@ -69,42 +68,44 @@ namespace corsika::process::qgsjetII {
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
const units::si::HEPEnergyType CoMenergy,
- const unsigned int Abeam,
- const unsigned int Atarget) const {
+ const unsigned int Abeam,
+ const unsigned int Atarget) const {
using namespace units::si;
double sigProd = std::numeric_limits<double>::infinity();
-
+
if (process::qgsjetII::CanInteract(BeamId)) {
const int iBeam = process::qgsjetII::GetQgsjetIIXSCode(BeamId);
if (!IsValidCoMEnergy(CoMenergy)) {
- ostringstream txt;
- txt << "Interaction: GetCrossSection: CoM energy outside range for QgsjetII! CoMenergy=" << CoMenergy;
- throw std::runtime_error(txt.str().c_str());
+ ostringstream txt;
+ txt << "Interaction: GetCrossSection: CoM energy outside range for QgsjetII! "
+ "CoMenergy="
+ << CoMenergy;
+ throw std::runtime_error(txt.str().c_str());
}
int iTarget = 1;
if (particles::IsNucleus(TargetId)) {
- iTarget = Atarget;
- if (iTarget > maxMassNumber_ || iTarget <= 0) {
- std::ostringstream txt;
- txt << "QgsjetII target outside range. iTarget=" << iTarget;
- throw std::runtime_error(txt.str().c_str());
- }
+ iTarget = Atarget;
+ if (iTarget > maxMassNumber_ || iTarget <= 0) {
+ std::ostringstream txt;
+ txt << "QgsjetII target outside range. iTarget=" << iTarget;
+ throw std::runtime_error(txt.str().c_str());
+ }
}
int iProjectile = 1;
if (particles::IsNucleus(BeamId)) {
- iProjectile = Abeam;
- if (iProjectile > maxMassNumber_ || iProjectile <= 0)
- throw std::runtime_error(
- "QgsjetII target outside range. ");
+ iProjectile = Abeam;
+ if (iProjectile > maxMassNumber_ || iProjectile <= 0)
+ throw std::runtime_error("QgsjetII target outside range. ");
}
-
+
const double dEcm = CoMenergy / 1_GeV / iProjectile; // CoM energy (per nucleon)
-
- cout << "QgsjetII::GetCrossSection " << dEcm << " " << iBeam << " " << iProjectile << " " << iTarget << endl;
+
+ cout << "QgsjetII::GetCrossSection " << dEcm << " " << iBeam << " " << iProjectile
+ << " " << iTarget << endl;
sigProd = qgsect_(dEcm, iBeam, iProjectile, iTarget);
}
-
+
return std::make_tuple(sigProd * 1_mb, 0_mb);
}
@@ -125,7 +126,7 @@ namespace corsika::process::qgsjetII {
// beam particles for qgsjetII : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = process::qgsjetII::CanInteract(corsikaBeamId);
-
+
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
@@ -136,8 +137,7 @@ namespace corsika::process::qgsjetII {
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
- const HEPEnergyType ECoM = sqrt(
- (Elab + pTotLabNorm) * (Elab - pTotLabNorm));
+ const HEPEnergyType ECoM = sqrt((Elab + pTotLabNorm) * (Elab - pTotLabNorm));
cout << "Interaction: LambdaInt: \n"
<< " input energy: " << vP.GetEnergy() / 1_GeV << endl
@@ -147,9 +147,8 @@ namespace corsika::process::qgsjetII {
if (kInteraction && IsValidCoMEnergy(ECoM)) {
int Abeam = 0;
- if (particles::IsNucleus(vP.GetPID()))
- Abeam = vP.GetNuclearA();
-
+ if (particles::IsNucleus(vP.GetPID())) Abeam = vP.GetNuclearA();
+
// get target from environment
/*
the target should be defined by the Environment,
@@ -163,10 +162,11 @@ namespace corsika::process::qgsjetII {
si::CrossSectionType weightedProdCrossSection = mediumComposition.WeightedSum(
[=](particles::Code targetID) -> si::CrossSectionType {
- int Atarget = 0;
- if (corsika::particles::IsNucleus(targetID))
- Atarget = particles::GetNucleusA(targetID);
- return std::get<0>(GetCrossSection(corsikaBeamId, targetID, ECoM, Abeam, Atarget));
+ int Atarget = 0;
+ if (corsika::particles::IsNucleus(targetID))
+ Atarget = particles::GetNucleusA(targetID);
+ return std::get<0>(
+ GetCrossSection(corsikaBeamId, targetID, ECoM, Abeam, Atarget));
});
cout << "Interaction: "
@@ -205,7 +205,7 @@ namespace corsika::process::qgsjetII {
<< process::qgsjetII::CanInteract(corsikaBeamId) << endl;
if (process::qgsjetII::CanInteract(corsikaBeamId)) {
-
+
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -219,15 +219,14 @@ namespace corsika::process::qgsjetII {
const auto eTargetLab = 0_GeV + constants::nucleonMass;
const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
-
+
// define projectile
HEPEnergyType const eProjectileLab = vP.GetEnergy();
auto const pProjectileLab = vP.GetMomentum();
-
+
int Abeam = 0;
- if (particles::IsNucleus(corsikaBeamId))
- Abeam = vP.GetNuclearA();
-
+ if (particles::IsNucleus(corsikaBeamId)) Abeam = vP.GetNuclearA();
+
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
<< endl;
@@ -237,7 +236,7 @@ namespace corsika::process::qgsjetII {
// artificial momentum, very small momentum, high mass, but direction of primary
HEPMomentumType pShort = 0.000001_eV;
const FourVector PprojLab(100_TeV, pProjectileLab * pShort / pProjectileLab.norm());
-
+
// define target kinematics in lab frame
// define boost to and from CoM frame
// CoM frame definition in QgsjetII projectile: +z
@@ -266,10 +265,11 @@ namespace corsika::process::qgsjetII {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- int Atarget = 0;
- if (corsika::particles::IsNucleus(targetId))
- Atarget = particles::GetNucleusA(targetId);
- const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm, Abeam, Atarget);
+ int Atarget = 0;
+ if (corsika::particles::IsNucleus(targetId))
+ Atarget = particles::GetNucleusA(targetId);
+ const auto [sigProd, sigEla] =
+ GetCrossSection(corsikaBeamId, targetId, Ecm, Abeam, Atarget);
[[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
}
@@ -283,38 +283,38 @@ namespace corsika::process::qgsjetII {
allowed air in atmosphere also contains some Argon.
*/
int targetQgsCode = -1;
- if (particles::IsNucleus(targetCode)) targetQgsCode = particles::GetNucleusA(targetCode);
+ if (particles::IsNucleus(targetCode))
+ targetQgsCode = particles::GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) targetQgsCode = 1;
cout << "Interaction: target qgsjetII code/A: " << targetQgsCode << endl;
if (targetQgsCode > maxMassNumber_ || targetQgsCode < 1)
- throw std::runtime_error(
- "QgsjetII target outside range.");
-
+ throw std::runtime_error("QgsjetII target outside range.");
+
int projQgsCode = 1;
if (particles::IsNucleus(corsikaBeamId)) projQgsCode = vP.GetNuclearA();
- cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " " << corsikaBeamId << endl;
+ cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " "
+ << corsikaBeamId << endl;
if (projQgsCode > maxMassNumber_ || projQgsCode < 1)
- throw std::runtime_error(
- "QgsjetII target outside range.");
+ throw std::runtime_error("QgsjetII target outside range.");
// beam id for qgsjetII
int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
if (corsikaBeamId != particles::Code::Nucleus) {
- kBeam = process::qgsjetII::ConvertToQgsjetIIRaw(corsikaBeamId);
- // from conex
- if (kBeam==0) { // replace pi0 or rho0 with pi+/pi-
- static int select = 1;
- kBeam = select;
- select *= -1;
- }
- // replace lambda by neutron
- if (kBeam==6)
- kBeam = 3;
- else if (kBeam==-6)
- kBeam = -3;
- // else if (abs(kBeam)>6) -> throw
+ kBeam = process::qgsjetII::ConvertToQgsjetIIRaw(corsikaBeamId);
+ // from conex
+ if (kBeam == 0) { // replace pi0 or rho0 with pi+/pi-
+ static int select = 1;
+ kBeam = select;
+ select *= -1;
+ }
+ // replace lambda by neutron
+ if (kBeam == 6)
+ kBeam = 3;
+ else if (kBeam == -6)
+ kBeam = -3;
+ // else if (abs(kBeam)>6) -> throw
}
-
+
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl;
@@ -328,114 +328,114 @@ namespace corsika::process::qgsjetII {
throw std::runtime_error("energy too low for QGSJETII");
} else {
count_++;
- auto Elab = eProjectileLab/projQgsCode;
- qgini_(Elab/1_GeV, kBeam, projQgsCode, targetQgsCode);
- qgini_(Elab/1_GeV, kBeam, projQgsCode, targetQgsCode);
- qgconf_();
+ auto Elab = eProjectileLab / projQgsCode;
+ qgini_(Elab / 1_GeV, kBeam, projQgsCode, targetQgsCode);
+ qgini_(Elab / 1_GeV, kBeam, projQgsCode, targetQgsCode);
+ qgconf_();
- // bookkeeping
- MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV;
+ // bookkeeping
+ MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV;
- // fragments
+ // fragments
QGSJetIIFragmentsStack qfs;
- for (auto& fragm : qfs) {
- particles::Code idFragm = particles::Code::Nucleus;
- int A = fragm.GetFragmentSize();
- int Z = 0;
- switch (A) {
- case 1:
- { // proton/neutron
- idFragm = particles::Code::Proton;
- MomentumVector Plab_nucleon(rootCS, {0.0_GeV, 0.0_GeV, sqrt((Elab+particles::Proton::GetMass())*(Elab-particles::Proton::GetMass()))});
- auto const PlabRot = boost.fromCoM(FourVector(Elab, Plab_nucleon));
- cout << "secondary fragment>" //<< static_cast<int>(psec.GetPID())
- << " " << idFragm
- << " Elab=" << Elab << " " << endl;
- auto pnew = vP.AddSecondary(
- tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- geometry::Point, units::si::TimeType>{
- idFragm,
- PlabRot.GetTimeLikeComponent(), PlabRot.GetSpaceLikeComponents(),
- pOrig, tOrig});
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
- }
- break;
- case 2: // deuterium
- Z = 1;
- break;
- case 3: // tritium
- Z = 1;
- break;
- case 4: // helium
- Z = 2;
- break;
- default: // nucleus
- {
- Z = int(A / 2.15 + 0.7);
- }
- }
-
- if (idFragm == particles::Code::Nucleus) {
- MomentumVector Plab_nucleus(rootCS, {0.0_GeV, 0.0_GeV, sqrt((Elab+constants::nucleonMass*A)*(Elab-constants::nucleonMass*A))});
- auto const PlabRot = boost.fromCoM(FourVector(Elab*A, Plab_nucleus));
- cout << "secondary fragment>" //<< static_cast<int>(psec.GetPID())
- << " " << idFragm
- << " Elab=" << Elab
- << " A=" << A
- << " Z=" << Z
- << endl;
- auto pnew = vP.AddSecondary(
- tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- geometry::Point, units::si::TimeType, unsigned short, unsigned short>{
- idFragm,
- PlabRot.GetTimeLikeComponent(), PlabRot.GetSpaceLikeComponents(),
- pOrig, tOrig, A, Z});
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
- }
- }
-
- // secondaries
- QGSJetIIStack qs;
+ for (auto& fragm : qfs) {
+ particles::Code idFragm = particles::Code::Nucleus;
+ int A = fragm.GetFragmentSize();
+ int Z = 0;
+ switch (A) {
+ case 1: { // proton/neutron
+ idFragm = particles::Code::Proton;
+ MomentumVector Plab_nucleon(rootCS,
+ {0.0_GeV, 0.0_GeV,
+ sqrt((Elab + particles::Proton::GetMass()) *
+ (Elab - particles::Proton::GetMass()))});
+ auto const PlabRot = boost.fromCoM(FourVector(Elab, Plab_nucleon));
+ cout << "secondary fragment>" //<< static_cast<int>(psec.GetPID())
+ << " " << idFragm << " Elab=" << Elab << " " << endl;
+ auto pnew = vP.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{
+ idFragm, PlabRot.GetTimeLikeComponent(),
+ PlabRot.GetSpaceLikeComponents(), pOrig, tOrig});
+ Plab_final += pnew.GetMomentum();
+ Elab_final += pnew.GetEnergy();
+ } break;
+ case 2: // deuterium
+ Z = 1;
+ break;
+ case 3: // tritium
+ Z = 1;
+ break;
+ case 4: // helium
+ Z = 2;
+ break;
+ default: // nucleus
+ {
+ Z = int(A / 2.15 + 0.7);
+ }
+ }
+
+ if (idFragm == particles::Code::Nucleus) {
+ MomentumVector Plab_nucleus(rootCS,
+ {0.0_GeV, 0.0_GeV,
+ sqrt((Elab + constants::nucleonMass * A) *
+ (Elab - constants::nucleonMass * A))});
+ auto const PlabRot = boost.fromCoM(FourVector(Elab * A, Plab_nucleus));
+ cout << "secondary fragment>" //<< static_cast<int>(psec.GetPID())
+ << " " << idFragm << " Elab=" << Elab << " A=" << A << " Z=" << Z
+ << endl;
+ auto pnew = vP.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+ geometry::Point, units::si::TimeType, unsigned short,
+ unsigned short>{idFragm, PlabRot.GetTimeLikeComponent(),
+ PlabRot.GetSpaceLikeComponents(), pOrig, tOrig, A,
+ Z});
+ Plab_final += pnew.GetMomentum();
+ Elab_final += pnew.GetEnergy();
+ }
+ }
+
+ // secondaries
+ QGSJetIIStack qs;
for (auto& psec : qs) {
// // skip particles that have decayed in QgsjetII
- //if (psec.HasDecayed()) continue;
+ // if (psec.HasDecayed()) continue;
- auto const Plab = psec.GetMomentum();
- auto const Elab = psec.GetEnergy();
+ auto const Plab = psec.GetMomentum();
+ auto const Elab = psec.GetEnergy();
// transform energy to lab. frame
- //auto const pCoM = psec.GetMomentum();
- //HEPEnergyType const eCoM = psec.GetEnergy();
+ // auto const pCoM = psec.GetMomentum();
+ // HEPEnergyType const eCoM = psec.GetEnergy();
auto const PlabRot = boost.fromCoM(FourVector(Elab, Plab));
- cout << "secondary hadron>" //<< static_cast<int>(psec.GetPID())
- << " " << process::qgsjetII::ConvertFromQgsjetII(psec.GetPID())
- << " Elab=" << Elab << " " << endl;
+ cout << "secondary hadron>" //<< static_cast<int>(psec.GetPID())
+ << " " << process::qgsjetII::ConvertFromQgsjetII(psec.GetPID())
+ << " Elab=" << Elab << " " << endl;
// add to corsika stack
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()),
- PlabRot.GetTimeLikeComponent(), PlabRot.GetSpaceLikeComponents(),
- //Elab, Plab,
- pOrig, tOrig});
-
+ PlabRot.GetTimeLikeComponent(), PlabRot.GetSpaceLikeComponents(),
+ // Elab, Plab,
+ pOrig, tOrig});
+
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
- //Ecm_final += psec.GetEnergy();
+ // Ecm_final += psec.GetEnergy();
}
cout << "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/ << endl
<< "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
- << ", N_wounded,targ=" << QGSJetIIFragmentsStackData::GetWoundedNucleonsTarget()
- << ", N_wounded,proj=" << QGSJetIIFragmentsStackData::GetWoundedNucleonsProjectile()
- << ", N_fragm,proj=" << qfs.GetSize()
- << endl;
+ << ", N_wounded,targ="
+ << QGSJetIIFragmentsStackData::GetWoundedNucleonsTarget()
+ << ", N_wounded,proj="
+ << QGSJetIIFragmentsStackData::GetWoundedNucleonsProjectile()
+ << ", N_fragm,proj=" << qfs.GetSize() << endl;
}
}
return process::EProcessReturn::eOk;
diff --git a/Processes/QGSJetII/Interaction.h b/Processes/QGSJetII/Interaction.h
index ee386086b5c22922ddce3291d78a3fed952ebde9..8404a406135b65e819999e32c67779bd4d546e89 100644
--- a/Processes/QGSJetII/Interaction.h
+++ b/Processes/QGSJetII/Interaction.h
@@ -16,8 +16,8 @@
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
-#include <tuple>
#include <string>
+#include <tuple>
namespace corsika::process::qgsjetII {
@@ -28,7 +28,7 @@ namespace corsika::process::qgsjetII {
bool initialized_ = false;
public:
- Interaction(const std::string& dataPath="");
+ Interaction(const std::string& dataPath = "");
~Interaction();
void Init();
@@ -47,9 +47,8 @@ namespace corsika::process::qgsjetII {
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code, const corsika::particles::Code,
- const corsika::units::si::HEPEnergyType,
- const unsigned int Abeam=0,
- const unsigned int Atarget=0) const;
+ const corsika::units::si::HEPEnergyType, const unsigned int Abeam = 0,
+ const unsigned int Atarget = 0) const;
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const;
diff --git a/Processes/QGSJetII/ParticleConversion.cc b/Processes/QGSJetII/ParticleConversion.cc
index f177c541256876e339ca439c4181bedf3e118bcb..a18d5673723e573034048c929abbb75e9a08a1ac 100644
--- a/Processes/QGSJetII/ParticleConversion.cc
+++ b/Processes/QGSJetII/ParticleConversion.cc
@@ -12,4 +12,3 @@
#include <corsika/process/qgsjetII/ParticleConversion.h>
using namespace corsika::process::qgsjetII;
-
diff --git a/Processes/QGSJetII/ParticleConversion.h b/Processes/QGSJetII/ParticleConversion.h
index 08b8bdb0b2fabc480b8617b54670e9965cc8fa81..382626cf7146111073c2cbefb8b3f9673ed36c73 100644
--- a/Processes/QGSJetII/ParticleConversion.h
+++ b/Processes/QGSJetII/ParticleConversion.h
@@ -43,13 +43,12 @@ namespace corsika::process::qgsjetII {
}
int constexpr GetQgsjetIIXSCode(corsika::particles::Code pCode) {
- if (pCode==corsika::particles::Code::Nucleus)
- return 2;
+ if (pCode == corsika::particles::Code::Nucleus) return 2;
return corsika2qgsjetIIXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
}
bool constexpr CanInteract(corsika::particles::Code pCode) {
- return (GetQgsjetIIXSCode(pCode) > 0) && (ConvertToQgsjetIIRaw(pCode)<=5);
+ return (GetQgsjetIIXSCode(pCode) > 0) && (ConvertToQgsjetIIRaw(pCode) <= 5);
}
} // namespace corsika::process::qgsjetII
diff --git a/Processes/QGSJetII/QGSJetIIFragmentsStack.h b/Processes/QGSJetII/QGSJetIIFragmentsStack.h
index 96e338042c809eefb90735f09a0740cea2bc9ead..958327703958324645c209b938a8a1c67b46c001 100644
--- a/Processes/QGSJetII/QGSJetIIFragmentsStack.h
+++ b/Processes/QGSJetII/QGSJetIIFragmentsStack.h
@@ -25,15 +25,18 @@ namespace corsika::process::qgsjetII {
public:
void Init();
void Dump() const {}
-
- void Clear() { qgarr13_.nsf = 0; qgarr55_.nwt = 0; }
+
+ void Clear() {
+ qgarr13_.nsf = 0;
+ qgarr55_.nwt = 0;
+ }
unsigned int GetSize() const { return qgarr13_.nsf; }
unsigned int GetCapacity() const { return iapmax; }
static unsigned int GetWoundedNucleonsTarget() { return qgarr55_.nwt; }
static unsigned int GetWoundedNucleonsProjectile() { return qgarr55_.nwp; }
-
- int GetFragmentSize(const unsigned int i) const { return qgarr13_.iaf[i]; }
+
+ int GetFragmentSize(const unsigned int i) const { return qgarr13_.iaf[i]; }
void SetFragmentSize(const unsigned int i, const int v) { qgarr13_.iaf[i] = v; }
void Copy(const unsigned int i1, const unsigned int i2) {
@@ -57,26 +60,20 @@ namespace corsika::process::qgsjetII {
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetParticleData(const int vSize) {
- SetFragmentSize(vSize);
- }
+ void SetParticleData(const int vSize) { SetFragmentSize(vSize); }
void SetParticleData(FragmentsInterface<StackIteratorInterface>& /*parent*/,
const int vSize) {
SetFragmentSize(vSize);
}
- void SetFragmentSize(const int v) {
- GetStackData().SetFragmentSize(GetIndex(), v);
- }
-
- double GetFragmentSize() const {
- return GetStackData().GetFragmentSize(GetIndex());
- }
+ void SetFragmentSize(const int v) { GetStackData().SetFragmentSize(GetIndex(), v); }
+ double GetFragmentSize() const { return GetStackData().GetFragmentSize(GetIndex()); }
};
- typedef corsika::stack::Stack<QGSJetIIFragmentsStackData, FragmentsInterface> QGSJetIIFragmentsStack;
+ typedef corsika::stack::Stack<QGSJetIIFragmentsStackData, FragmentsInterface>
+ QGSJetIIFragmentsStack;
} // end namespace corsika::process::qgsjetII
diff --git a/Processes/QGSJetII/QGSJetIIStack.h b/Processes/QGSJetII/QGSJetIIStack.h
index d907c27f3163ff97282dfbb13d4f940386c9d1a1..5f7ea3550e9f1e4dae03f78d2133b446b1c19327 100644
--- a/Processes/QGSJetII/QGSJetIIStack.h
+++ b/Processes/QGSJetII/QGSJetIIStack.h
@@ -28,7 +28,11 @@ namespace corsika::process::qgsjetII {
void Init();
void Dump() const {}
- void Clear() { qgarr12_.nsp = 0; qgarr13_.nsf = 0; qgarr55_.nwt = 0; }
+ void Clear() {
+ qgarr12_.nsp = 0;
+ qgarr13_.nsf = 0;
+ qgarr55_.nwt = 0;
+ }
unsigned int GetSize() const { return qgarr12_.nsp; }
unsigned int GetCapacity() const { return nptmax; }
@@ -58,15 +62,15 @@ namespace corsika::process::qgsjetII {
using namespace corsika::units::si;
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- QuantityVector<hepmomentum_d> components = {
- qgarr14_.esp[i][2] * 1_GeV, qgarr14_.esp[i][3] * 1_GeV, qgarr14_.esp[i][1] * 1_GeV};
+ QuantityVector<hepmomentum_d> components = {qgarr14_.esp[i][2] * 1_GeV,
+ qgarr14_.esp[i][3] * 1_GeV,
+ qgarr14_.esp[i][1] * 1_GeV};
return MomentumVector(rootCS, components);
}
void Copy(const unsigned int i1, const unsigned int i2) {
qgarr14_.ich[i2] = qgarr14_.ich[i1];
- for (unsigned int i = 0; i < 4; ++i)
- qgarr14_.esp[i2][i] = qgarr14_.esp[i1][i];
+ for (unsigned int i = 0; i < 4; ++i) qgarr14_.esp[i2][i] = qgarr14_.esp[i1][i];
}
void Swap(const unsigned int i1, const unsigned int i2) {
diff --git a/Processes/QGSJetII/qgsjet-II-04.cc b/Processes/QGSJetII/qgsjet-II-04.cc
index 17cb65b07827157460f4d28366586d7ddeefd422..2e716cf6231a5dfd710006fdf39dece200069be4 100644
--- a/Processes/QGSJetII/qgsjet-II-04.cc
+++ b/Processes/QGSJetII/qgsjet-II-04.cc
@@ -12,20 +12,19 @@
#include <corsika/random/RNGManager.h>
-#include <random>
#include <iostream>
+#include <random>
datadir::datadir(const std::string& dir) {
- if (dir.length()>130) {
- std::cerr << "QGSJetII error, will cut datadir \"" << dir << "\" to 130 characters: " << std::endl;
+ if (dir.length() > 130) {
+ std::cerr << "QGSJetII error, will cut datadir \"" << dir
+ << "\" to 130 characters: " << std::endl;
}
- int i=0;
- for (i=0;i<std::min(130,int(dir.length())); ++i)
- data[i] = dir[i];
- data[i+0] = ' ';
- data[i+1] = '\0';
+ int i = 0;
+ for (i = 0; i < std::min(130, int(dir.length())); ++i) data[i] = dir[i];
+ data[i + 0] = ' ';
+ data[i + 1] = '\0';
}
-
double qgran_(int&) {
static corsika::random::RNG& rng =
diff --git a/Processes/QGSJetII/qgsjet-II-04.h b/Processes/QGSJetII/qgsjet-II-04.h
index 8376d3b234a1acb0b5eb8751552b0de94b34c09f..94dd9ef971ac53d199dbee49f2011ee04d38ba08 100644
--- a/Processes/QGSJetII/qgsjet-II-04.h
+++ b/Processes/QGSJetII/qgsjet-II-04.h
@@ -20,98 +20,93 @@
extern "C" {
- // data memory layout
-
- extern struct {
- int nsp;
- } qgarr12_;
-
- const int nptmax = 95000;
- const int iapmax = 208;
-
- extern struct {
- double esp[nptmax][4];
- int ich[nptmax];
- } qgarr14_;
-
- extern struct {
- // c nsf - number of secondary fragments;
- // c iaf(i) - mass of the i-th fragment
- int nsf;
- int iaf[iapmax];
- } qgarr13_;
-
- extern struct {
- int nwt;
- int nwp;
- } qgarr55_;
-
-
- /**
- Small helper class to provide a data-directory name in the format qgsjetII expects
- */
- class datadir {
- private:
- datadir operator=(const std::string& dir);
- datadir operator=(const datadir&);
- public:
- datadir(const std::string& dir);
- char data[132];
- };
-
- // functions
- void qgset_();
- void qgaini_(const char* datdir); // Note: there is a length limiation 132 from fortran-qgsjet here
-
- /**
- @function qgini_
-
- additional initialization procedure per event
-
- @parameter e0n - interaction energy (per hadron/nucleon),
- @parameter icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), +-4 - K+-, +-5 - K_l/s),
- @parameter iap - projectile mass number (1 - for a hadron),
- @parameter iat - target mass number
- */
- void qgini_(const double& e0n, const int& icp0, const int& iap, const int& iat);
-
-
- /**
- @function qgconf_
-
- generate one event configuration
- */
- void qgconf_();
-
-
- /**
- @function qgsect_
-
- hadron-nucleus (hadron-nucleus) particle production cross section
-
- @parameter e0n lab. energy per projectile nucleon (hadron)
- @parameter icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
- @parameter iap projectile mass number (1=<iap<=iapmax),
- @parameter iat target mass number (1=<iat<=iapmax)
- */
- double qgsect_(const double& e0n, const int& icz, const int& iap0, const int& iat0);
-
-
- /**
- @function qgran
-
- link to random number generation
- */
- double qgran_(int&);
-
-
- /**
- dummy function from CRMC
- */
- void lzmaopenfile_(const char* name, int length);
- void lzmaclosefile_();
- void lzmafillarray_(const double& dum, const int& idum);
-
+// data memory layout
+
+extern struct { int nsp; } qgarr12_;
+
+const int nptmax = 95000;
+const int iapmax = 208;
+
+extern struct {
+ double esp[nptmax][4];
+ int ich[nptmax];
+} qgarr14_;
+
+extern struct {
+ // c nsf - number of secondary fragments;
+ // c iaf(i) - mass of the i-th fragment
+ int nsf;
+ int iaf[iapmax];
+} qgarr13_;
+
+extern struct {
+ int nwt;
+ int nwp;
+} qgarr55_;
+
+/**
+ Small helper class to provide a data-directory name in the format qgsjetII expects
+ */
+class datadir {
+private:
+ datadir operator=(const std::string& dir);
+ datadir operator=(const datadir&);
+
+public:
+ datadir(const std::string& dir);
+ char data[132];
+};
+
+// functions
+void qgset_();
+void qgaini_(
+ const char* datdir); // Note: there is a length limiation 132 from fortran-qgsjet here
+
+/**
+ @function qgini_
+
+ additional initialization procedure per event
+
+ @parameter e0n - interaction energy (per hadron/nucleon),
+ @parameter icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), +-4 - K+-, +-5 -
+ K_l/s),
+ @parameter iap - projectile mass number (1 - for a hadron),
+ @parameter iat - target mass number
+*/
+void qgini_(const double& e0n, const int& icp0, const int& iap, const int& iat);
+
+/**
+ @function qgconf_
+
+ generate one event configuration
+*/
+void qgconf_();
+
+/**
+ @function qgsect_
+
+ hadron-nucleus (hadron-nucleus) particle production cross section
+
+ @parameter e0n lab. energy per projectile nucleon (hadron)
+ @parameter icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
+ @parameter iap projectile mass number (1=<iap<=iapmax),
+ @parameter iat target mass number (1=<iat<=iapmax)
+ */
+double qgsect_(const double& e0n, const int& icz, const int& iap0, const int& iat0);
+
+/**
+ @function qgran
+
+ link to random number generation
+ */
+double qgran_(int&);
+
+/**
+ dummy function from CRMC
+ */
+void lzmaopenfile_(const char* name, int length);
+void lzmaclosefile_();
+void lzmafillarray_(const double& dum, const int& idum);
}
#endif
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index c91c6c5fa5e6f7d6b198f1f520b802f236c7d067..e084fc88852428e8d6a9baaf6ee214bf4ee89e15 100644
--- a/Processes/QGSJetII/testQGSJetII.cc
+++ b/Processes/QGSJetII/testQGSJetII.cc
@@ -26,8 +26,8 @@ using namespace corsika::process::qgsjetII;
TEST_CASE("QgsjetII", "[processes]") {
SECTION("QgsjetII -> Corsika") {
- REQUIRE(particles::PiPlus::GetCode() ==
- process::qgsjetII::ConvertFromQgsjetII(process::qgsjetII::QgsjetIICode::PiPlus));
+ REQUIRE(particles::PiPlus::GetCode() == process::qgsjetII::ConvertFromQgsjetII(
+ process::qgsjetII::QgsjetIICode::PiPlus));
}
SECTION("Corsika -> QgsjetII") {
@@ -41,7 +41,7 @@ TEST_CASE("QgsjetII", "[processes]") {
REQUIRE(process::qgsjetII::CanInteract(particles::Proton::GetCode()));
REQUIRE(process::qgsjetII::CanInteract(particles::Code::KPlus));
REQUIRE(process::qgsjetII::CanInteract(particles::Nucleus::GetCode()));
- //REQUIRE(process::qgsjetII::CanInteract(particles::Helium::GetCode()));
+ // REQUIRE(process::qgsjetII::CanInteract(particles::Helium::GetCode()));
REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::EtaC::GetCode()));
REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::SigmaC0::GetCode()));
@@ -106,10 +106,11 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType, unsigned int, unsigned int>{
- particles::Code::Nucleus, E0, plab, pos, 0_ns, 16, 8});
+ auto particle =
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned int, unsigned int>{
+ particles::Code::Nucleus, E0, plab, pos, 0_ns, 16, 8});
// corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
// particles::Code::PiPlus, E0, plab, pos, 0_ns});