Forward port of Geant4 SIMP custom physics model to master branch

SimG4Core/CustomPhysics/data/particles_simp_1000_GeV.txt

@@ -0,0 +1,4 @@
+Block MASS   #
+#  PDG code    mass (GeV)         particle
+   9000006     1000.0          # chi
+Block

SimG4Core/CustomPhysics/data/particles_simp_100_GeV.txt

@@ -0,0 +1,4 @@
+Block MASS   #
+#  PDG code    mass                 particle
+   9000006     100.0                # chi
+Block

SimG4Core/CustomPhysics/data/particles_simp_10_GeV.txt

@@ -0,0 +1,4 @@
+Block MASS   #
+#  PDG code    mass                 particle
+   9000006     10.0                  # chi
+Block

SimG4Core/CustomPhysics/data/particles_simp_1_GeV.txt

@@ -0,0 +1,4 @@
+Block MASS   #
+#  PDG code    mass                 particle
+   9000006     1.0                  # chi
+Block

SimG4Core/CustomPhysics/data/particles_simp_200_GeV.txt

@@ -0,0 +1,4 @@
+Block MASS   #
+#  PDG code    mass                 particle
+   9000006     200.0                  # chi
+Block

SimG4Core/CustomPhysics/data/particles_simp_400_GeV.txt

@@ -0,0 +1,4 @@
+Block MASS   #
+#  PDG code    mass                 particle
+   9000006     400.0                  # chi
+Block

SimG4Core/CustomPhysics/data/particles_simp_700_GeV.txt

@@ -0,0 +1,4 @@
+Block MASS   #
+#  PDG code    mass                 particle
+   9000006     700.0                  # chi
+Block

SimG4Core/CustomPhysics/interface/CMSAntiSIMP.h

@@ -0,0 +1,20 @@
+#ifndef SimG4Core_CustomPhysics_CMSAntiSIMP_H
+#define SimG4Core_CustomPhysics_CMSAntiSIMP_H
+
+#include "globals.hh"
+#include "G4ios.hh"
+#include "G4ParticleDefinition.hh"
+
+class CMSAntiSIMP : public G4ParticleDefinition {
+private:
+  static CMSAntiSIMP* theInstance;
+  CMSAntiSIMP() {}
+  ~CMSAntiSIMP() override {}
+
+public:
+  static CMSAntiSIMP* Definition(double mass);
+  static CMSAntiSIMP* AntiSIMPDefinition(double mass);
+  static CMSAntiSIMP* AntiSIMP();
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSQGSPSIMPBuilder.h

@@ -0,0 +1,46 @@
+#ifndef SimG4Core_CustomPhysics_CMSQGSPSIMPBuilder_H
+#define SimG4Core_CustomPhysics_CMSQGSPSIMPBuilder_H
+
+#include "globals.hh"
+
+#include "G4HadronElasticProcess.hh"
+#include "G4HadronFissionProcess.hh"
+#include "G4HadronCaptureProcess.hh"
+#include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h"
+
+#include "G4TheoFSGenerator.hh"
+#include "G4PreCompoundModel.hh"
+#include "G4GeneratorPrecompoundInterface.hh"
+#include "G4QGSModel.hh"
+#include "G4QGSParticipants.hh"
+#include "G4QGSMFragmentation.hh"
+#include "G4ExcitedStringDecay.hh"
+#include "G4QuasiElasticChannel.hh"
+
+class CMSQGSPSIMPBuilder {
+public:
+  CMSQGSPSIMPBuilder(G4bool quasiElastic = false);
+  virtual ~CMSQGSPSIMPBuilder();
+
+public:
+  virtual void Build(G4HadronElasticProcess* aP);
+  virtual void Build(G4HadronFissionProcess* aP);
+  virtual void Build(G4HadronCaptureProcess* aP);
+  virtual void Build(CMSSIMPInelasticProcess* aP);
+
+  void SetMinEnergy(G4double aM) { theMin = aM; }
+
+private:
+  G4TheoFSGenerator* theModel;
+  G4PreCompoundModel* thePreEquilib;
+  G4GeneratorPrecompoundInterface* theCascade;
+  G4QGSModel<G4QGSParticipants>* theStringModel;
+  G4ExcitedStringDecay* theStringDecay;
+  G4QuasiElasticChannel* theQuasiElastic;
+
+  G4QGSMFragmentation* theQGSM;
+
+  G4double theMin;
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSIMP.h

@@ -0,0 +1,20 @@
+#ifndef SimG4Core_CustomPhysics_CMSSIMP_H
+#define SimG4Core_CustomPhysics_CMSSIMP_H
+
+#include "globals.hh"
+#include "G4ios.hh"
+#include "G4ParticleDefinition.hh"
+
+class CMSSIMP : public G4ParticleDefinition {
+private:
+  static CMSSIMP* theInstance;
+  CMSSIMP() {}
+  ~CMSSIMP() override {}
+
+public:
+  static CMSSIMP* Definition(double mass);
+  static CMSSIMP* SIMPDefinition(double mass);
+  static CMSSIMP* SIMP();
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h

@@ -0,0 +1,167 @@
+#ifndef SimG4Core_CustomPhysics_CMSSIMPInelasticProcess_H
+#define SimG4Core_CustomPhysics_CMSSIMPInelasticProcess_H
+
+#include "G4HadronicProcess.hh"
+
+class G4ParticleDefinition;
+
+class CMSSIMPInelasticProcess : public G4HadronicProcess {
+public:
+  CMSSIMPInelasticProcess(const G4String& processName = "SIMPInelastic");
+
+  ~CMSSIMPInelasticProcess() override;
+
+  G4bool IsApplicable(const G4ParticleDefinition& aParticleType) override;
+
+  // register generator of secondaries
+  void RegisterMe(G4HadronicInteraction* a);
+
+  // get cross section per element
+  inline G4double GetElementCrossSection(const G4DynamicParticle* part,
+                                         const G4Element* elm,
+                                         const G4Material* mat = nullptr) {
+    G4double x = theCrossSectionDataStore->GetCrossSection(part, elm, mat);
+    if (x < 0.0) {
+      x = 0.0;
+    }
+    return x;
+  }
+
+  // obsolete method to get cross section per element
+  inline G4double GetMicroscopicCrossSection(const G4DynamicParticle* part,
+                                             const G4Element* elm,
+                                             const G4Material* mat = nullptr) {
+    return GetElementCrossSection(part, elm, mat);
+  }
+
+  // generic PostStepDoIt recommended for all derived classes
+  G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) override;
+
+  // initialisation of physics tables and CMSSIMPInelasticProcessStore
+  void PreparePhysicsTable(const G4ParticleDefinition&) override;
+
+  // build physics tables and print out the configuration of the process
+  void BuildPhysicsTable(const G4ParticleDefinition&) override;
+
+  // dump physics tables
+  inline void DumpPhysicsTable(const G4ParticleDefinition& p) { theCrossSectionDataStore->DumpPhysicsTable(p); }
+
+  // add cross section data set
+  inline void AddDataSet(G4VCrossSectionDataSet* aDataSet) { theCrossSectionDataStore->AddDataSet(aDataSet); }
+
+  // access to the manager
+  inline G4EnergyRangeManager* GetManagerPointer() { return &theEnergyRangeManager; }
+
+  // get inverse cross section per volume
+  G4double GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*) override;
+
+  // access to the target nucleus
+  inline const G4Nucleus* GetTargetNucleus() const { return &targetNucleus; }
+
+  //  G4ParticleDefinition* GetTargetDefinition();
+  inline const G4Isotope* GetTargetIsotope() { return targetNucleus.GetIsotope(); }
+
+  void ProcessDescription(std::ostream& outFile) const override;
+
+protected:
+  // generic method to choose secondary generator
+  // recommended for all derived classes
+  inline G4HadronicInteraction* ChooseHadronicInteraction(G4double kineticEnergy,
+                                                          G4Material* aMaterial,
+                                                          G4Element* anElement) {
+    return theEnergyRangeManager.GetHadronicInteraction(kineticEnergy, aMaterial, anElement);
+  }
+
+  // access to the target nucleus
+  inline G4Nucleus* GetTargetNucleusPointer() { return &targetNucleus; }
+
+public:
+  void BiasCrossSectionByFactor(G4double aScale);
+
+  // Energy-momentum non-conservation limits and reporting
+  inline void SetEpReportLevel(G4int level) { epReportLevel = level; }
+
+  inline void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double absoluteLevel) {
+    epCheckLevels.first = relativeLevel;
+    epCheckLevels.second = absoluteLevel;
+    levelsSetByProcess = true;
+  }
+
+  inline std::pair<G4double, G4double> GetEnergyMomentumCheckLevels() const { return epCheckLevels; }
+
+  // access to the cross section data store
+  inline G4CrossSectionDataStore* GetCrossSectionDataStore() { return theCrossSectionDataStore; }
+
+  inline void MultiplyCrossSectionBy(G4double factor) { aScaleFactor = factor; }
+
+protected:
+  void DumpState(const G4Track&, const G4String&, G4ExceptionDescription&);
+
+  // obsolete method will be removed
+  inline const G4EnergyRangeManager& GetEnergyRangeManager() const { return theEnergyRangeManager; }
+
+  // obsolete method will be removed
+  inline void SetEnergyRangeManager(const G4EnergyRangeManager& value) { theEnergyRangeManager = value; }
+
+  // access to the chosen generator
+  inline G4HadronicInteraction* GetHadronicInteraction() const { return theInteraction; }
+
+  // access to the cross section data set
+  inline G4double GetLastCrossSection() { return theLastCrossSection; }
+
+  // fill result
+  void FillResult(G4HadFinalState* aR, const G4Track& aT);
+
+  // Check the result for catastrophic energy non-conservation
+  G4HadFinalState* CheckResult(const G4HadProjectile& thePro,
+                               const G4Nucleus& targetNucleus,
+                               G4HadFinalState* result) const;
+
+  // Check 4-momentum balance
+  void CheckEnergyMomentumConservation(const G4Track&, const G4Nucleus&);
+
+private:
+  G4double XBiasSurvivalProbability();
+  G4double XBiasSecondaryWeight();
+
+  // hide assignment operator as private
+  CMSSIMPInelasticProcess& operator=(const CMSSIMPInelasticProcess& right) = delete;
+  CMSSIMPInelasticProcess(const CMSSIMPInelasticProcess&) = delete;
+
+  // Set E/p conservation check levels from environment variables
+  void GetEnergyMomentumCheckEnvvars();
+
+protected:
+  G4HadProjectile thePro;
+
+  G4ParticleChange* theTotalResult;
+
+  G4int epReportLevel;
+
+private:
+  G4ParticleDefinition* theParticle;
+
+  G4EnergyRangeManager theEnergyRangeManager;
+
+  G4HadronicInteraction* theInteraction;
+
+  G4CrossSectionDataStore* theCrossSectionDataStore;
+
+  G4Nucleus targetNucleus;
+
+  bool CMSSIMPInelasticProcess_debug_flag;
+
+  // Energy-momentum checking
+  std::pair<G4double, G4double> epCheckLevels;
+  G4bool levelsSetByProcess;
+
+  std::vector<G4VLeadingParticleBiasing*> theBias;
+
+  G4double theInitialNumberOfInteractionLength;
+
+  G4double aScaleFactor;
+  G4bool xBiasOn;
+  G4double theLastCrossSection;
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSIMPInelasticXS.h

@@ -0,0 +1,69 @@
+#ifndef SimG4Core_CustomPhysics_CMSSIMPInelasticXS_H
+#define SimG4Core_CustomPhysics_CMSSIMPInelasticXS_H
+
+#include "G4VCrossSectionDataSet.hh"
+#include "globals.hh"
+#include "G4ElementData.hh"
+#include <vector>
+
+const G4int MAXZINEL = 93;
+
+class G4DynamicParticle;
+class G4ParticleDefinition;
+class G4Element;
+class G4PhysicsVector;
+class G4ComponentGGHadronNucleusXsc;
+class G4HadronNucleonXsc;
+
+class CMSSIMPInelasticXS : public G4VCrossSectionDataSet {
+public:
+  CMSSIMPInelasticXS();
+
+  ~CMSSIMPInelasticXS() override;
+
+  G4bool IsElementApplicable(const G4DynamicParticle*, G4int Z, const G4Material*) override;
+
+  G4bool IsIsoApplicable(const G4DynamicParticle*, G4int Z, G4int A, const G4Element*, const G4Material*) override;
+
+  G4double GetElementCrossSection(const G4DynamicParticle*, G4int Z, const G4Material* mat = nullptr) override;
+
+  G4double GetIsoCrossSection(const G4DynamicParticle*,
+                              G4int Z,
+                              G4int A,
+                              const G4Isotope* iso,
+                              const G4Element* elm,
+                              const G4Material* mat) override;
+
+  G4Isotope* SelectIsotope(const G4Element*, G4double kinEnergy) override;
+
+  void BuildPhysicsTable(const G4ParticleDefinition&) override;
+
+  void CrossSectionDescription(std::ostream&) const override;
+
+private:
+  void Initialise(G4int Z, G4DynamicParticle* dp = nullptr, const char* = nullptr);
+
+  G4PhysicsVector* RetrieveVector(std::ostringstream& in, G4bool warn);
+
+  G4double IsoCrossSection(G4double ekin, G4int Z, G4int A);
+
+  CMSSIMPInelasticXS& operator=(const CMSSIMPInelasticXS& right) = delete;
+  CMSSIMPInelasticXS(const CMSSIMPInelasticXS&) = delete;
+
+  G4ComponentGGHadronNucleusXsc* ggXsection;
+  G4HadronNucleonXsc* fNucleon;
+
+  const G4ParticleDefinition* proton;
+
+  G4ElementData data;
+  std::vector<G4PhysicsVector*> work;
+  std::vector<G4double> temp;
+  std::vector<G4double> coeff;
+
+  G4bool isInitialized;
+
+  static const G4int amin[MAXZINEL];
+  static const G4int amax[MAXZINEL];
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CustomPDGParser.h

@@ -16,6 +16,7 @@ class CustomPDGParser {
   static bool s_isRGlueball(int pdg);
   static bool s_isDphoton(int pdg);
   static bool s_isChargino(int pdg);
+  static bool s_isSIMP(int pdg);
   static double s_charge(int pdg);
   static double s_spin(int pdg);
   static std::vector<int> s_containedQuarks(int pdg);