Merge pull request #28814 from civanch/restore_custom_physics

Clean-up SIMP particles simulation
cms-sw · Jan 30, 2020 · a4239e5 · a4239e5
2 parents 325bccb + 3c87f42
commit a4239e5
Show file tree

Hide file tree

Showing 7 changed files with 92 additions and 1,190 deletions.
diff --git a/SimG4Core/CustomPhysics/interface/CMSQGSPSIMPBuilder.h b/SimG4Core/CustomPhysics/interface/CMSQGSPSIMPBuilder.h
@@ -3,44 +3,27 @@
 
 #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"
+#include "G4QGSModel.hh"
 
-class CMSQGSPSIMPBuilder {
-public:
-  CMSQGSPSIMPBuilder(G4bool quasiElastic = false);
-  virtual ~CMSQGSPSIMPBuilder();
+class CMSSIMPInelasticProcess;
+class G4QGSParticipants;
+class G4QGSMFragmentation;
+class G4ExcitedStringDecay;
 
+class CMSQGSPSIMPBuilder {
 public:
-  virtual void Build(G4HadronElasticProcess* aP);
-  virtual void Build(G4HadronFissionProcess* aP);
-  virtual void Build(G4HadronCaptureProcess* aP);
-  virtual void Build(CMSSIMPInelasticProcess* aP);
+  CMSQGSPSIMPBuilder();
+  ~CMSQGSPSIMPBuilder();
 
-  void SetMinEnergy(G4double aM) { theMin = aM; }
+  void Build(CMSSIMPInelasticProcess* aP);
 
 private:
-  G4TheoFSGenerator* theModel;
-  G4PreCompoundModel* thePreEquilib;
-  G4GeneratorPrecompoundInterface* theCascade;
   G4QGSModel<G4QGSParticipants>* theStringModel;
   G4ExcitedStringDecay* theStringDecay;
-  G4QuasiElasticChannel* theQuasiElastic;
-
   G4QGSMFragmentation* theQGSM;
-
-  G4double theMin;
 };
 
 #endif
diff --git a/SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h b/SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h
@@ -13,155 +13,15 @@ class CMSSIMPInelasticProcess : public G4HadronicProcess {
 
   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
diff --git a/SimG4Core/CustomPhysics/interface/CMSSIMPInelasticXS.h b/SimG4Core/CustomPhysics/interface/CMSSIMPInelasticXS.h
@@ -3,17 +3,8 @@
 
 #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 G4NeutronInelasticXS;
 
 class CMSSIMPInelasticXS : public G4VCrossSectionDataSet {
 public:
@@ -34,36 +25,17 @@ class CMSSIMPInelasticXS : public G4VCrossSectionDataSet {
                               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;
-
+  G4NeutronInelasticXS* nXsection;
+  const G4ParticleDefinition* neutron;
   G4bool isInitialized;
-
-  static const G4int amin[MAXZINEL];
-  static const G4int amax[MAXZINEL];
 };
 
 #endif
diff --git a/SimG4Core/CustomPhysics/src/CMSQGSPSIMPBuilder.cc b/SimG4Core/CustomPhysics/src/CMSQGSPSIMPBuilder.cc
@@ -1,57 +1,37 @@
 
 #include "SimG4Core/CustomPhysics/interface/CMSQGSPSIMPBuilder.h"
+#include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h"
+
 #include "G4SystemOfUnits.hh"
 #include "G4ParticleDefinition.hh"
-#include "G4ParticleTable.hh"
-#include "G4ProcessManager.hh"
-#include "G4ExcitationHandler.hh"
-
-#include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticXS.h"
-#include "SimG4Core/CustomPhysics/interface/CMSSIMP.h"
-
-CMSQGSPSIMPBuilder::CMSQGSPSIMPBuilder(G4bool quasiElastic) {
-  theMin = 12 * GeV;
-
-  theModel = new G4TheoFSGenerator("QGSP");
-
+#include "G4TheoFSGenerator.hh"
+#include "G4PreCompoundModel.hh"
+#include "G4GeneratorPrecompoundInterface.hh"
+#include "G4QGSParticipants.hh"
+#include "G4QGSMFragmentation.hh"
+#include "G4ExcitedStringDecay.hh"
+
+CMSQGSPSIMPBuilder::CMSQGSPSIMPBuilder() {
   theStringModel = new G4QGSModel<G4QGSParticipants>;
   theStringDecay = new G4ExcitedStringDecay(theQGSM = new G4QGSMFragmentation);
   theStringModel->SetFragmentationModel(theStringDecay);
-
-  theCascade = new G4GeneratorPrecompoundInterface;
-  thePreEquilib = new G4PreCompoundModel(new G4ExcitationHandler);
-  theCascade->SetDeExcitation(thePreEquilib);
-
-  theModel->SetTransport(theCascade);
-  theModel->SetHighEnergyGenerator(theStringModel);
-  if (quasiElastic) {
-    theQuasiElastic = new G4QuasiElasticChannel;
-    theModel->SetQuasiElasticChannel(theQuasiElastic);
-  } else {
-    theQuasiElastic = nullptr;
-  }
 }
 
 CMSQGSPSIMPBuilder::~CMSQGSPSIMPBuilder() {
   delete theStringDecay;
   delete theStringModel;
-  delete thePreEquilib;
-  delete theCascade;
-  if (theQuasiElastic)
-    delete theQuasiElastic;
-  delete theModel;
   delete theQGSM;
 }
 
-void CMSQGSPSIMPBuilder::Build(G4HadronElasticProcess*) {}
-
-void CMSQGSPSIMPBuilder::Build(G4HadronFissionProcess*) {}
-
-void CMSQGSPSIMPBuilder::Build(G4HadronCaptureProcess*) {}
-
 void CMSQGSPSIMPBuilder::Build(CMSSIMPInelasticProcess* aP) {
-  theModel->SetMinEnergy(theMin);
-  theModel->SetMaxEnergy(100 * TeV);
+  G4GeneratorPrecompoundInterface* theCascade = new G4GeneratorPrecompoundInterface;
+  G4PreCompoundModel* thePreEquilib = new G4PreCompoundModel();
+  theCascade->SetDeExcitation(thePreEquilib);
+
+  G4TheoFSGenerator* theModel = new G4TheoFSGenerator("QGSP");
+  theModel->SetTransport(theCascade);
+  theModel->SetHighEnergyGenerator(theStringModel);
+  theModel->SetMinEnergy(0.0);
+  theModel->SetMaxEnergy(100 * CLHEP::TeV);
   aP->RegisterMe(theModel);
-  aP->AddDataSet(new CMSSIMPInelasticXS());
 }