IAP GITLAB

Skip to content
Snippets Groups Projects
Commit ee9ec049 authored by Ralf Ulrich's avatar Ralf Ulrich
Browse files

Merge branch 'qgsjet-fixes' into 'master' 

QgsjetII fixes

There are issues with the particle ID mapping to QGSJetII

See merge request AirShowerPhysics/corsika!182
parents 76ad9926 405b093f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Processes/QGSJetII/Interaction.cc
+ 35
− 49
View file @ ee9ec049
...@@ -21,9 +21,7 @@ ...@@ -21,9 +21,7 @@
#include <corsika/setup/SetupTrajectory.h> #include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h> #include <corsika/utl/COMBoost.h>
#include <corsika/random/RNGManager.h> #include <random>
#include <corsika/random/UniformRealDistribution.h>
#include <sstream> #include <sstream>
#include <string> #include <string>
#include <tuple> #include <tuple>
...@@ -55,9 +53,6 @@ namespace corsika::process::qgsjetII { ...@@ -55,9 +53,6 @@ namespace corsika::process::qgsjetII {
Interaction::~Interaction() { cout << "QgsjetII::Interaction n=" << count_ << endl; } Interaction::~Interaction() { cout << "QgsjetII::Interaction n=" << count_ << endl; }
void Interaction::Init() { void Interaction::Init() {
using random::RNGManager;
// initialize QgsjetII // initialize QgsjetII
if (!initialized_) { if (!initialized_) {
qgset_(); qgset_();
...@@ -76,7 +71,7 @@ namespace corsika::process::qgsjetII { ...@@ -76,7 +71,7 @@ namespace corsika::process::qgsjetII {
if (process::qgsjetII::CanInteract(beamId)) { if (process::qgsjetII::CanInteract(beamId)) {
const int iBeam = process::qgsjetII::GetQgsjetIIXSCode(beamId); const int iBeam = process::qgsjetII::GetQgsjetIIXSCodeRaw(beamId);
int iTarget = 1; int iTarget = 1;
if (particles::IsNucleus(targetId)) { if (particles::IsNucleus(targetId)) {
iTarget = targetA; iTarget = targetA;
...@@ -253,56 +248,47 @@ namespace corsika::process::qgsjetII { ...@@ -253,56 +248,47 @@ namespace corsika::process::qgsjetII {
mediumComposition.SampleTarget(cross_section_of_components, rng_); mediumComposition.SampleTarget(cross_section_of_components, rng_);
cout << "Interaction: target selected: " << targetCode << endl; cout << "Interaction: target selected: " << targetCode << endl;
int targetQgsCode = -1; int targetMassNumber = 1; // proton
if (particles::IsNucleus(targetCode)) if (particles::IsNucleus(targetCode)) { // nucleus
targetQgsCode = particles::GetNucleusA(targetCode); targetMassNumber = particles::GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) targetQgsCode = 1; if (targetMassNumber > maxMassNumber_)
cout << "Interaction: target qgsjetII code/A: " << targetQgsCode << endl; throw std::runtime_error("QgsjetII target mass outside range.");
if (targetQgsCode > maxMassNumber_ || targetQgsCode < 1) } else {
throw std::runtime_error("QgsjetII target outside range."); if (targetCode != particles::Proton::GetCode())
throw std::runtime_error("QgsjetII Taget not possible.");
int projQgsCode = 1; }
if (particles::IsNucleus(corsikaBeamId)) projQgsCode = vP.GetNuclearA(); cout << "Interaction: target qgsjetII code/A: " << targetMassNumber << endl;
cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " "
<< corsikaBeamId << endl; int projectileMassNumber = 1; // "1" means "hadron"
if (projQgsCode > maxMassNumber_ || projQgsCode < 1) QgsjetIIHadronType qgsjet_hadron_type =
throw std::runtime_error("QgsjetII target outside range."); process::qgsjetII::GetQgsjetIIHadronType(corsikaBeamId);
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
// beam id for qgsjetII projectileMassNumber = vP.GetNuclearA();
QgsjetIICode qgsjet_beam_code; if (projectileMassNumber > maxMassNumber_)
if (corsikaBeamId == particles::Code::Nucleus) { throw std::runtime_error("QgsjetII projectile mass outside range.");
std::array<QgsjetIICode, 2> constexpr nucleons = {QgsjetIICode::Proton, std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
QgsjetIICode::Neutron}; QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1); std::uniform_int_distribution select(0, 1);
qgsjet_beam_code = nucleons[select(rng_)]; qgsjet_hadron_type = nucleons[select(rng_)];
} else { // it is an "elementary" particle } else {
qgsjet_beam_code = process::qgsjetII::ConvertToQgsjetII(corsikaBeamId); // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
// from conex if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
if (qgsjet_beam_code == QgsjetIICode::Pi0 or qgsjet_hadron_type = alternate_;
qgsjet_beam_code == QgsjetIICode::Rho0) { // replace pi0 or rho0 with pi+/pi- alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
// in alternating sequence ? QgsjetIIHadronType::PiMinusType
qgsjet_beam_code = alternate_; : QgsjetIIHadronType::PiPlusType);
alternate_ = (alternate_ == QgsjetIICode::PiPlus ? QgsjetIICode::PiMinus
: QgsjetIICode::PiPlus);
} }
// 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: " cout << "Interaction: "
<< " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << endl; << " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << endl;
count_++; count_++;
qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode, qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
targetQgsCode); targetMassNumber);
// this is from CRMC, is this REALLY needed ???
qgini_(projectileEnergyLab / 1_GeV, qgsjet_beam_code_int, projQgsCode,
targetQgsCode);
qgconf_(); qgconf_();
// bookkeeping // bookkeeping
......
Processes/QGSJetII/Interaction.h
+ 2
− 1
View file @ ee9ec049
...@@ -26,7 +26,8 @@ namespace corsika::process::qgsjetII { ...@@ -26,7 +26,8 @@ namespace corsika::process::qgsjetII {
std::string data_path_; std::string data_path_;
int count_ = 0; int count_ = 0;
bool initialized_ = false; bool initialized_ = false;
QgsjetIICode alternate_ = QgsjetIICode::PiPlus; // for pi0, rho0 projectiles QgsjetIIHadronType alternate_ =
QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
public: public:
Interaction(const std::string& dataPath = ""); Interaction(const std::string& dataPath = "");
......
Processes/QGSJetII/ParticleConversion.h
+ 50
− 5
View file @ ee9ec049
...@@ -17,9 +17,44 @@ ...@@ -17,9 +17,44 @@
namespace corsika::process::qgsjetII { namespace corsika::process::qgsjetII {
/**
These are the possible secondaries produced by QGSJetII
*/
enum class QgsjetIICode : int8_t; enum class QgsjetIICode : int8_t;
using QgsjetIICodeIntType = std::underlying_type<QgsjetIICode>::type; 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> #include <corsika/process/qgsjetII/Generated.inc>
QgsjetIICode constexpr ConvertToQgsjetII(corsika::particles::Code pCode) { QgsjetIICode constexpr ConvertToQgsjetII(corsika::particles::Code pCode) {
...@@ -38,17 +73,27 @@ namespace corsika::process::qgsjetII { ...@@ -38,17 +73,27 @@ namespace corsika::process::qgsjetII {
return corsikaCode; return corsikaCode;
} }
int constexpr ConvertToQgsjetIIRaw(corsika::particles::Code pCode) { QgsjetIICodeIntType constexpr ConvertToQgsjetIIRaw(corsika::particles::Code pCode) {
return static_cast<int>(ConvertToQgsjetII(pCode)); return static_cast<QgsjetIICodeIntType>(ConvertToQgsjetII(pCode));
} }
int constexpr GetQgsjetIIXSCode(corsika::particles::Code pCode) { QgsjetIIXSClass constexpr GetQgsjetIIXSCode(corsika::particles::Code pCode) {
if (pCode == corsika::particles::Code::Nucleus) return 2; // if (pCode == corsika::particles::Code::Nucleus)
// static_cast(QgsjetIIXSClassIntType>();
return corsika2qgsjetIIXStype[static_cast<corsika::particles::CodeIntType>(pCode)]; 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) { 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 } // namespace corsika::process::qgsjetII
......
Processes/QGSJetII/code_generator.py
+ 71
− 9
View file @ ee9ec049
...@@ -14,14 +14,67 @@ import pickle, sys, itertools ...@@ -14,14 +14,67 @@ import pickle, sys, itertools
# loads the pickled particle_db (which is an OrderedDict) # 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): def load_particledb(filename):
with open(filename, "rb") as f: with open(filename, "rb") as f:
particle_db = pickle.load(f) particle_db = pickle.load(f)
return particle_db 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): def read_qgsjetII_codes(filename, particle_db):
with open(filename) as f: with open(filename) as f:
for line in f: for line in f:
...@@ -33,13 +86,11 @@ def read_qgsjetII_codes(filename, particle_db): ...@@ -33,13 +86,11 @@ def read_qgsjetII_codes(filename, particle_db):
identifier, model_code, xsType = line.split() identifier, model_code, xsType = line.split()
try: try:
particle_db[identifier]["qgsjetII_code"] = int(model_code) 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: except KeyError as e:
raise Exception("Identifier '{:s}' not found in particle_db".format(identifier)) raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
# generates the enum to access qgsjetII particles by readable names # generates the enum to access qgsjetII particles by readable names
def generate_qgsjetII_enum(particle_db): def generate_qgsjetII_enum(particle_db):
output = "enum class QgsjetIICode : int8_t {\n" output = "enum class QgsjetIICode : int8_t {\n"
...@@ -50,7 +101,6 @@ def generate_qgsjetII_enum(particle_db): ...@@ -50,7 +101,6 @@ def generate_qgsjetII_enum(particle_db):
return output return output
# generates the look-up table to convert corsika codes to qgsjetII codes # generates the look-up table to convert corsika codes to qgsjetII codes
def generate_corsika2qgsjetII(particle_db): def generate_corsika2qgsjetII(particle_db):
string = "std::array<QgsjetIICodeIntType, {:d}> constexpr corsika2qgsjetII = {{\n".format(len(particle_db)) string = "std::array<QgsjetIICodeIntType, {:d}> constexpr corsika2qgsjetII = {{\n".format(len(particle_db))
...@@ -61,13 +111,22 @@ def generate_corsika2qgsjetII(particle_db): ...@@ -61,13 +111,22 @@ def generate_corsika2qgsjetII(particle_db):
return string return string
# generates the look-up table to convert corsika codes to qgsjetII codes # generates the look-up table to convert corsika codes to qgsjetII codes
def generate_corsika2qgsjetII_xsType(particle_db): 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(): for identifier, pData in particle_db.items():
modelCodeXS = pData.get("qgsjetII_xsType", -1) modelCodeXS = pData.get("qgsjetII_xsType", "CannotInteract")
string += " {:d}, // {:s}\n".format(modelCodeXS, identifier if modelCodeXS else identifier + " (not implemented in QGSJETII)") 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" string += "};\n"
return string return string
...@@ -98,6 +157,7 @@ def generate_qgsjetII2corsika(particle_db) : ...@@ -98,6 +157,7 @@ def generate_qgsjetII2corsika(particle_db) :
string += "};\n" string += "};\n"
return string return string
if __name__ == "__main__": if __name__ == "__main__":
if len(sys.argv) != 3: if len(sys.argv) != 3:
print("usage: {:s} <particle_db.pkl> <qgsjetII_codes.dat>".format(sys.argv[0]), file=sys.stderr) print("usage: {:s} <particle_db.pkl> <qgsjetII_codes.dat>".format(sys.argv[0]), file=sys.stderr)
...@@ -107,6 +167,7 @@ if __name__ == "__main__": ...@@ -107,6 +167,7 @@ if __name__ == "__main__":
particle_db = load_particledb(sys.argv[1]) particle_db = load_particledb(sys.argv[1])
read_qgsjetII_codes(sys.argv[2], particle_db) read_qgsjetII_codes(sys.argv[2], particle_db)
set_default_qgsjetII_definition(particle_db)
with open("Generated.inc", "w") as f: with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f) print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
...@@ -114,3 +175,4 @@ if __name__ == "__main__": ...@@ -114,3 +175,4 @@ if __name__ == "__main__":
print(generate_corsika2qgsjetII(particle_db), file=f) print(generate_corsika2qgsjetII(particle_db), file=f)
print(generate_qgsjetII2corsika(particle_db), file=f) print(generate_qgsjetII2corsika(particle_db), file=f)
print(generate_corsika2qgsjetII_xsType(particle_db), file=f) print(generate_corsika2qgsjetII_xsType(particle_db), file=f)
print(generate_corsika2qgsjetII_hadronType(particle_db), file=f)
Processes/QGSJetII/qgsjet-II-04-codes.dat
+ 35
− 50
View file @ ee9ec049
# input file for particle conversion to/from QGSJet # input file for particle conversion to/from QGSJet
# the format of this file is: "corsika-identifier" "qgsjet-id" # the format of this file is: "corsika-identifier" "qgsjet-id" "xs-class"
# class 0 (cannot interact) # Note, we list here only the particles, which are produced by
Electron 11 0 # QGSJetII as secondaries. There is additional mapping of corsika
Positron -11 0 # particles on QGSJetII "xs-class"es and "hadron-type"s, which are
# programmed in code_generator.py
# class 1
Pi0 0 1 # class 1 (-> as pions)
PiPlus 1 1 Pi0 0 LightMesons
PiMinus -1 1 PiPlus 1 LightMesons
# in crmc: all particles from 100 to 999 ??? PiMinus -1 LightMesons
Eta 10 1 Eta 10 LightMesons
Rho0 -10 1 Rho0 -10 LightMesons
# class 2 # class 2 (-> as proton)
#Nucleus 40 2 Neutron 3 Baryons
Neutron 3 2 AntiNeutron -3 Baryons
AntiNeutron -3 2 Proton 2 Baryons
Proton 2 2 AntiProton -2 Baryons
AntiProton -2 2 Lambda0 6 Baryons
# in crmc: 1000 to 9999 ??? Lambda0Bar -6 Baryons
Lambda0 6 2 LambdaCPlus 9 Baryons
Lambda0Bar -6 2 LambdaCMinusBar -9 Baryons
LambdaCPlus 9 2
LambdaCMinusBar -9 2 # class 3 (-> as kaon)
K0Short 5 Kaons
# class 3 K0Long -5 Kaons
K0Long -5 3 KPlus 4 Kaons
K0 5 3 KMinus -4 Kaons
K0Bar -5 3
K0Short 5 3 # class 4 (-> charmed mesons, not in qgsjetII)
# ambiguity between the K0/b and K0s/l D0 8 Charmed
KPlus 4 3 D0Bar -8 Charmed
KMinus -4 3 DPlus 7 Charmed
DMinus -7 Charmed
# 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
Processes/QGSJetII/qgsjet-II-04.f
+ 1
− 1
View file @ ee9ec049
...@@ -2051,7 +2051,7 @@ c216 format(2x,'qgaini: integrated cut Pomeron eikonals') ...@@ -2051,7 +2051,7 @@ c216 format(2x,'qgaini: integrated cut Pomeron eikonals')
   
c============================================================================= c=============================================================================
subroutine qgini(e0n,icp0,iap,iat) subroutine qgini(e0n,icp0,iap,iat)
c----------------------------------------------------------------------------- c---------------------s--------------------------------------------------------
c additional initialization procedure c additional initialization procedure
c e0n - interaction energy (per hadron/nucleon), c e0n - interaction energy (per hadron/nucleon),
c icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), c icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~),
......
Processes/QGSJetII/testQGSJetII.cc
+ 8
− 4
View file @ ee9ec049
...@@ -49,10 +49,14 @@ TEST_CASE("QgsjetII", "[processes]") { ...@@ -49,10 +49,14 @@ TEST_CASE("QgsjetII", "[processes]") {
SECTION("cross-section type") { SECTION("cross-section type") {
REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) == 2); REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) ==
REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) == 3); process::qgsjetII::QgsjetIIXSClass::Baryons);
REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) == 2); REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) ==
REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) == 1); 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);
} }
} }
......
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