Fixed PID mapping between C8 and QGSjetII .

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

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 = "");

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

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)

Processes/QGSJetII/qgsjet-II-04-codes.dat

 # 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
 

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~),

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);
   }
 }