diff --git a/Processes/QGSJetII/Interaction.cc b/Processes/QGSJetII/Interaction.cc
index 17329a2fea426ee61381013e31f98410e22c88be..98dd4b1948eb3555d0ba315e198fdfe44f5fe4be 100644
--- a/Processes/QGSJetII/Interaction.cc
+++ b/Processes/QGSJetII/Interaction.cc
@@ -21,9 +21,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/random/UniformRealDistribution.h>
-
+#include <random>
#include <sstream>
#include <string>
#include <tuple>
@@ -55,9 +53,6 @@ namespace corsika::process::qgsjetII {
Interaction::~Interaction() { cout << "QgsjetII::Interaction n=" << count_ << endl; }
void Interaction::Init() {
-
- using random::RNGManager;
-
// initialize QgsjetII
if (!initialized_) {
qgset_();
@@ -76,7 +71,7 @@ namespace corsika::process::qgsjetII {
if (process::qgsjetII::CanInteract(beamId)) {
- const int iBeam = process::qgsjetII::GetQgsjetIIXSCode(beamId);
+ const int iBeam = process::qgsjetII::GetQgsjetIIXSCodeRaw(beamId);
int iTarget = 1;
if (particles::IsNucleus(targetId)) {
iTarget = targetA;
@@ -253,56 +248,47 @@ namespace corsika::process::qgsjetII {
mediumComposition.SampleTarget(cross_section_of_components, rng_);
cout << "Interaction: target selected: " << targetCode << endl;
- int targetQgsCode = -1;
- 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.");
-
- int projQgsCode = 1;
- if (particles::IsNucleus(corsikaBeamId)) projQgsCode = vP.GetNuclearA();
- cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " "
- << corsikaBeamId << endl;
- if (projQgsCode > maxMassNumber_ || projQgsCode < 1)
- throw std::runtime_error("QgsjetII target outside range.");
-
- // beam id for qgsjetII
- QgsjetIICode qgsjet_beam_code;
- if (corsikaBeamId == particles::Code::Nucleus) {
- std::array<QgsjetIICode, 2> constexpr nucleons = {QgsjetIICode::Proton,
- QgsjetIICode::Neutron};
+ int targetMassNumber = 1; // proton
+ if (particles::IsNucleus(targetCode)) { // nucleus
+ targetMassNumber = particles::GetNucleusA(targetCode);
+ if (targetMassNumber > maxMassNumber_)
+ throw std::runtime_error("QgsjetII target mass outside range.");
+ } else {
+ if (targetCode != particles::Proton::GetCode())
+ throw std::runtime_error("QgsjetII Taget not possible.");
+ }
+ cout << "Interaction: target qgsjetII code/A: " << targetMassNumber << endl;
+
+ int projectileMassNumber = 1; // "1" means "hadron"
+ QgsjetIIHadronType qgsjet_hadron_type =
+ process::qgsjetII::GetQgsjetIIHadronType(corsikaBeamId);
+ if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
+ projectileMassNumber = vP.GetNuclearA();
+ if (projectileMassNumber > maxMassNumber_)
+ throw std::runtime_error("QgsjetII projectile mass outside range.");
+ std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
+ QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1);
- qgsjet_beam_code = nucleons[select(rng_)];
- } else { // it is an "elementary" particle
- qgsjet_beam_code = process::qgsjetII::ConvertToQgsjetII(corsikaBeamId);
- // from conex
- if (qgsjet_beam_code == QgsjetIICode::Pi0 or
- qgsjet_beam_code == QgsjetIICode::Rho0) { // replace pi0 or rho0 with pi+/pi-
- // in alternating sequence
- qgsjet_beam_code = alternate_;
- alternate_ = (alternate_ == QgsjetIICode::PiPlus ? QgsjetIICode::PiMinus
- : QgsjetIICode::PiPlus);
+ qgsjet_hadron_type = nucleons[select(rng_)];
+ } else {
+ // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
+ if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
+ qgsjet_hadron_type = alternate_;
+ alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
+ ? QgsjetIIHadronType::PiMinusType
+ : QgsjetIIHadronType::PiPlusType);
}
- // replace lambda by neutron
- if (qgsjet_beam_code == QgsjetIICode::Lambda0)
- qgsjet_beam_code = QgsjetIICode::Neutron;
- else if (qgsjet_beam_code == QgsjetIICode::Lambda0Bar)
- qgsjet_beam_code = QgsjetIICode::AntiNeutron;
- // else if (abs(qgsjet_beam_code)>6) -> throw
}
+ cout << "Interaction: projectile qgsjetII code/A: " << projectileMassNumber << " "
+ << corsikaBeamId << endl;
- int qgsjet_beam_code_int = static_cast<QgsjetIICodeIntType>(qgsjet_beam_code);
+ int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << endl;
count_++;
- qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode,
- targetQgsCode);
- // this is from CRMC, is this REALLY needed ???
- qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode,
- targetQgsCode);
+ qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
+ targetMassNumber);
qgconf_();
// bookkeeping
diff --git a/Processes/QGSJetII/Interaction.h b/Processes/QGSJetII/Interaction.h
index d06a078430590390070de5969b2425bb44a87d68..7df95c1517a75722c2f2e7c19cd21baa50694678 100644
--- a/Processes/QGSJetII/Interaction.h
+++ b/Processes/QGSJetII/Interaction.h
@@ -26,7 +26,8 @@ namespace corsika::process::qgsjetII {
std::string data_path_;
int count_ = 0;
bool initialized_ = false;
- QgsjetIICode alternate_ = QgsjetIICode::PiPlus; // for pi0, rho0 projectiles
+ QgsjetIIHadronType alternate_ =
+ QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
public:
Interaction(const std::string& dataPath = "");
diff --git a/Processes/QGSJetII/ParticleConversion.h b/Processes/QGSJetII/ParticleConversion.h
index 1a89658666deda150d5131c6bc0c01e06951174a..5828e82f583e24d2680fd9d0da9ade929b2aa455 100644
--- a/Processes/QGSJetII/ParticleConversion.h
+++ b/Processes/QGSJetII/ParticleConversion.h
@@ -17,9 +17,44 @@
namespace corsika::process::qgsjetII {
+ /**
+ These are the possible secondaries produced by QGSJetII
+ */
enum class QgsjetIICode : int8_t;
using QgsjetIICodeIntType = std::underlying_type<QgsjetIICode>::type;
+ /**
+ These are the possible projectile for which QGSJetII knwos cross section
+ */
+ enum class QgsjetIIXSClass : int8_t {
+ CannotInteract = 0,
+ LightMesons = 1,
+ Baryons = 2,
+ Kaons = 3,
+ };
+ using QgsjetIIXSClassIntType = std::underlying_type<QgsjetIIXSClass>::type;
+
+ /**
+ These are the only possible projectile types in QGSJetII
+ */
+ enum class QgsjetIIHadronType : int8_t {
+ UndefinedType = 0,
+ PiPlusType = +1,
+ PiMinusType = -1,
+ ProtonType = +2,
+ AntiProtonType = -2,
+ NeutronType = +3,
+ AntiNeutronType = -3,
+ KaonPlusType = +4,
+ KaonMinusType = -4,
+ Kaon0LType = +5,
+ Kaon0SType = -5,
+ // special codes, not in QGSJetII
+ NucleusType,
+ NeutralLightMesonType
+ };
+ using QgsjetIIHadronTypeIntType = std::underlying_type<QgsjetIIHadronType>::type;
+
#include <corsika/process/qgsjetII/Generated.inc>
QgsjetIICode constexpr ConvertToQgsjetII(corsika::particles::Code pCode) {
@@ -38,17 +73,27 @@ namespace corsika::process::qgsjetII {
return corsikaCode;
}
- int constexpr ConvertToQgsjetIIRaw(corsika::particles::Code pCode) {
- return static_cast<int>(ConvertToQgsjetII(pCode));
+ QgsjetIICodeIntType constexpr ConvertToQgsjetIIRaw(corsika::particles::Code pCode) {
+ return static_cast<QgsjetIICodeIntType>(ConvertToQgsjetII(pCode));
}
- int constexpr GetQgsjetIIXSCode(corsika::particles::Code pCode) {
- if (pCode == corsika::particles::Code::Nucleus) return 2;
+ QgsjetIIXSClass constexpr GetQgsjetIIXSCode(corsika::particles::Code pCode) {
+ // if (pCode == corsika::particles::Code::Nucleus)
+ // static_cast(QgsjetIIXSClassIntType>();
return corsika2qgsjetIIXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
}
+ QgsjetIIXSClassIntType constexpr GetQgsjetIIXSCodeRaw(corsika::particles::Code pCode) {
+ return static_cast<QgsjetIIXSClassIntType>(GetQgsjetIIXSCode(pCode));
+ }
+
bool constexpr CanInteract(corsika::particles::Code pCode) {
- return (GetQgsjetIIXSCode(pCode) > 0) && (ConvertToQgsjetIIRaw(pCode) <= 5);
+ return GetQgsjetIIXSCode(pCode) != QgsjetIIXSClass::CannotInteract;
+ }
+
+ QgsjetIIHadronType constexpr GetQgsjetIIHadronType(corsika::particles::Code pCode) {
+ return corsika2qgsjetIIHadronType[static_cast<corsika::particles::CodeIntType>(
+ pCode)];
}
} // namespace corsika::process::qgsjetII
diff --git a/Processes/QGSJetII/code_generator.py b/Processes/QGSJetII/code_generator.py
index 3b7113ec513e28fcb823fb1693cf83f6a6b6d545..d8f9a46937c01cd2576731bd7fad002da9eb77bb 100755
--- a/Processes/QGSJetII/code_generator.py
+++ b/Processes/QGSJetII/code_generator.py
@@ -14,14 +14,67 @@ import pickle, sys, itertools
# loads the pickled particle_db (which is an OrderedDict)
+# definition of particle_db dict is: "name", "antiName", "pdg", "mass", "electric_charge", "lifetime", "ngc_code", "isNucleus", "isHadron"
def load_particledb(filename):
with open(filename, "rb") as f:
particle_db = pickle.load(f)
return particle_db
+def set_default_qgsjetII_definition(particle_db):
+ for identifier, pData in particle_db.items():
+ # the cross-section types
+ xsType = "CannotInteract"
+ hadronType = "UndefinedType"
+ if (pData['isNucleus']):
+ xsType = "Baryons"
+ hadronType = "NucleusType"
+ elif (pData['isHadron']):
+ pdg = abs(pData['pdg'])
+ anti = pData['pdg'] < 0
+ isBaryon = (1000 <= pdg < 4000)
+ charge = pData['electric_charge']
+ if (pdg>=100 and pdg<300 and pdg!=130): # light mesons
+ xsType = "LightMesons"
+ if (charge==0):
+ hadronType = "NeutralLightMesonType"
+ else:
+ if (charge>0):
+ hadronType = "PiPlusType"
+ else:
+ hadronType = "PiMinusType"
+ elif ((pdg>=300 and pdg<400) or pdg in [130, 10313, 10323]): # kanos
+ xsType = "Kaons"
+ if (charge>0):
+ hadronType = "KaonPlusType"
+ else:
+ hadronType = "KaonMinusType"
+ if (charge==0):
+ hadronType = "Kaon0SType"
+ if (pdg == 130):
+ hadronType = "Kaon0LType"
+ elif (pdg == 310):
+ hadronType = "Kaon0SType"
+ elif (isBaryon or pData['isNucleus']): # baryons
+ xsType = "Baryons"
+ if (charge==0):
+ if (anti):
+ hadronType = "AntiNeutronType"
+ else:
+ hadronType = "NeutronType"
+ else:
+ if (charge>0):
+ hadronType = "ProtonType"
+ else:
+ hadronType = "AntiProtonType"
+ # all othe not-captured cased are hopefully irrelevant
+
-#
+ pData['qgsjetII_xsType'] = xsType
+ pData['qgsjetII_hadronType'] = hadronType
+
+
+# read the qgsjet-codes data file
def read_qgsjetII_codes(filename, particle_db):
with open(filename) as f:
for line in f:
@@ -33,13 +86,11 @@ def read_qgsjetII_codes(filename, particle_db):
identifier, model_code, xsType = line.split()
try:
particle_db[identifier]["qgsjetII_code"] = int(model_code)
- particle_db[identifier]["qgsjetII_xsType"] = int(xsType)
+ particle_db[identifier]["qgsjetII_xsType"] = xsType
except KeyError as e:
raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
-
-
# generates the enum to access qgsjetII particles by readable names
def generate_qgsjetII_enum(particle_db):
output = "enum class QgsjetIICode : int8_t {\n"
@@ -50,7 +101,6 @@ def generate_qgsjetII_enum(particle_db):
return output
-
# generates the look-up table to convert corsika codes to qgsjetII codes
def generate_corsika2qgsjetII(particle_db):
string = "std::array<QgsjetIICodeIntType, {:d}> constexpr corsika2qgsjetII = {{\n".format(len(particle_db))
@@ -61,13 +111,22 @@ def generate_corsika2qgsjetII(particle_db):
return string
-
# generates the look-up table to convert corsika codes to qgsjetII codes
def generate_corsika2qgsjetII_xsType(particle_db):
- string = "std::array<int, {:d}> constexpr corsika2qgsjetIIXStype = {{\n".format(len(particle_db))
+ string = "std::array<QgsjetIIXSClass, {:d}> constexpr corsika2qgsjetIIXStype = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
- modelCodeXS = pData.get("qgsjetII_xsType", -1)
- string += " {:d}, // {:s}\n".format(modelCodeXS, identifier if modelCodeXS else identifier + " (not implemented in QGSJETII)")
+ modelCodeXS = pData.get("qgsjetII_xsType", "CannotInteract")
+ string += " QgsjetIIXSClass::{:s}, // {:s}\n".format(modelCodeXS, identifier if modelCodeXS else identifier + " (not implemented in QGSJETII)")
+ string += "};\n"
+ return string
+
+
+# generates the look-up table to convert corsika codes to qgsjetII codes
+def generate_corsika2qgsjetII_hadronType(particle_db):
+ string = "std::array<QgsjetIIHadronType, {:d}> constexpr corsika2qgsjetIIHadronType = {{\n".format(len(particle_db))
+ for identifier, pData in particle_db.items():
+ modelCode = pData.get("qgsjetII_hadronType", "UndefinedType")
+ string += " QgsjetIIHadronType::{:s}, // {:s}\n".format(modelCode, identifier if modelCode else identifier + " (not implemented in QGSJETII)")
string += "};\n"
return string
@@ -98,6 +157,7 @@ def generate_qgsjetII2corsika(particle_db) :
string += "};\n"
return string
+
if __name__ == "__main__":
if len(sys.argv) != 3:
print("usage: {:s} <particle_db.pkl> <qgsjetII_codes.dat>".format(sys.argv[0]), file=sys.stderr)
@@ -107,6 +167,7 @@ if __name__ == "__main__":
particle_db = load_particledb(sys.argv[1])
read_qgsjetII_codes(sys.argv[2], particle_db)
+ set_default_qgsjetII_definition(particle_db)
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
@@ -114,3 +175,4 @@ if __name__ == "__main__":
print(generate_corsika2qgsjetII(particle_db), file=f)
print(generate_qgsjetII2corsika(particle_db), file=f)
print(generate_corsika2qgsjetII_xsType(particle_db), file=f)
+ print(generate_corsika2qgsjetII_hadronType(particle_db), file=f)
diff --git a/Processes/QGSJetII/qgsjet-II-04-codes.dat b/Processes/QGSJetII/qgsjet-II-04-codes.dat
index 3550eed83ce3f847153b8e79f85646af0816a7e4..ff0595bf1173d91cb1b732c2f99c0b9c8264d233 100644
--- a/Processes/QGSJetII/qgsjet-II-04-codes.dat
+++ b/Processes/QGSJetII/qgsjet-II-04-codes.dat
@@ -1,52 +1,37 @@
# input file for particle conversion to/from QGSJet
-# the format of this file is: "corsika-identifier" "qgsjet-id"
-
-# class 0 (cannot interact)
-Electron 11 0
-Positron -11 0
-
-# class 1
-Pi0 0 1
-PiPlus 1 1
-PiMinus -1 1
-# in crmc: all particles from 100 to 999 ???
-Eta 10 1
-Rho0 -10 1
-
-# class 2
-#Nucleus 40 2
-Neutron 3 2
-AntiNeutron -3 2
-Proton 2 2
-AntiProton -2 2
-# in crmc: 1000 to 9999 ???
-Lambda0 6 2
-Lambda0Bar -6 2
-LambdaCPlus 9 2
-LambdaCMinusBar -9 2
-
-# class 3
-K0Long -5 3
-K0 5 3
-K0Bar -5 3
-K0Short 5 3
-# ambiguity between the K0/b and K0s/l
-KPlus 4 3
-KMinus -4 3
-
-# class 4
-D0 8 4
-D0Bar 8 4
-DPlus 7 4
-DMinus -7 4
-#DS+/- (340)
-#etac (440)
-#j/psi (441)
-#h_1c0
-#psi'
-#Xi_0c0
-#Xi_1c0
-#Xi_2c0
-
-
+# the format of this file is: "corsika-identifier" "qgsjet-id" "xs-class"
+
+# Note, we list here only the particles, which are produced by
+# QGSJetII as secondaries. There is additional mapping of corsika
+# particles on QGSJetII "xs-class"es and "hadron-type"s, which are
+# programmed in code_generator.py
+
+# class 1 (-> as pions)
+Pi0 0 LightMesons
+PiPlus 1 LightMesons
+PiMinus -1 LightMesons
+Eta 10 LightMesons
+Rho0 -10 LightMesons
+
+# class 2 (-> as proton)
+Neutron 3 Baryons
+AntiNeutron -3 Baryons
+Proton 2 Baryons
+AntiProton -2 Baryons
+Lambda0 6 Baryons
+Lambda0Bar -6 Baryons
+LambdaCPlus 9 Baryons
+LambdaCMinusBar -9 Baryons
+
+# class 3 (-> as kaon)
+K0Short 5 Kaons
+K0Long -5 Kaons
+KPlus 4 Kaons
+KMinus -4 Kaons
+
+# class 4 (-> charmed mesons, not in qgsjetII)
+D0 8 Charmed
+D0Bar -8 Charmed
+DPlus 7 Charmed
+DMinus -7 Charmed
diff --git a/Processes/QGSJetII/qgsjet-II-04.f b/Processes/QGSJetII/qgsjet-II-04.f
index 40b2a4ccbc8da8506f5a6bafec36f190bec9b7bb..a3efe5977f9037e6e80eb3ec44864c71e29a6e2a 100644
--- a/Processes/QGSJetII/qgsjet-II-04.f
+++ b/Processes/QGSJetII/qgsjet-II-04.f
@@ -2051,7 +2051,7 @@ c216 format(2x,'qgaini: integrated cut Pomeron eikonals')
c=============================================================================
subroutine qgini(e0n,icp0,iap,iat)
-c-----------------------------------------------------------------------------
+c---------------------s--------------------------------------------------------
c additional initialization procedure
c e0n - interaction energy (per hadron/nucleon),
c icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~),
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index 4f284e056992a9987ae4cd70b4a5f471f5caca6a..8258b3c5eae4c52ab94633adae74f48a673a3d4b 100644
--- a/Processes/QGSJetII/testQGSJetII.cc
+++ b/Processes/QGSJetII/testQGSJetII.cc
@@ -49,10 +49,14 @@ TEST_CASE("QgsjetII", "[processes]") {
SECTION("cross-section type") {
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) == 2);
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) == 3);
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) == 2);
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) == 1);
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) ==
+ process::qgsjetII::QgsjetIIXSClass::Baryons);
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) ==
+ process::qgsjetII::QgsjetIIXSClass::Kaons);
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) ==
+ process::qgsjetII::QgsjetIIXSClass::Baryons);
+ REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) ==
+ process::qgsjetII::QgsjetIIXSClass::LightMesons);
}
}