Merge pull request #43403 from civanch/phys_list_phase2

[14_0_X SIM] Update FTFP_BERT_EMN physics

SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h

@@ -30,6 +30,7 @@ class CMSEmStandardPhysicsXS : public G4VPhysicsConstructor {
   G4double fSafetyFactor;
   G4double fLambdaLimit;
   G4MscStepLimitType fStepLimitType;
+  bool fG4HepEmActive;
 };
 
 #endif

SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_EMN.cc

@@ -8,6 +8,7 @@
 #include "G4IonPhysics.hh"
 #include "G4StoppingPhysics.hh"
 #include "G4HadronElasticPhysics.hh"
+#include "G4HadronicParameters.hh"
 
 FTFPCMS_BERT_EMN::FTFPCMS_BERT_EMN(const edm::ParameterSet& p) : PhysicsList(p) {
   int ver = p.getUntrackedParameter<int>("Verbosity", 0);
@@ -28,12 +29,25 @@ FTFPCMS_BERT_EMN::FTFPCMS_BERT_EMN(const edm::ParameterSet& p) : PhysicsList(p)
     // Synchroton Radiation & GN Physics
     G4EmExtraPhysics* gn = new G4EmExtraPhysics(ver);
     RegisterPhysics(gn);
+    bool mu = p.getParameter<bool>("G4MuonPairProductionByMuon");
+    gn->MuonToMuMu(mu);
+    edm::LogVerbatim("PhysicsList") << " Muon pair production by muons: " << mu;
   }
 
   // Decays
   this->RegisterPhysics(new G4DecayPhysics(ver));
 
   if (hadPhys) {
+    bool ngen = p.getParameter<bool>("G4NeutronGeneralProcess");
+    bool bc = p.getParameter<bool>("G4BCHadronicProcess");
+    bool hn = p.getParameter<bool>("G4LightHyperNucleiTracking");
+    auto param = G4HadronicParameters::Instance();
+    param->SetEnableNeutronGeneralProcess(ngen);
+    param->SetEnableBCParticles(bc);
+    param->SetEnableHyperNuclei(hn);
+    edm::LogVerbatim("PhysicsList") << " Eneble neutron general process: " << ngen
+                                    << "\n Enable b- and c- hadron physics: " << bc
+                                    << "\n Enable light hyper-nuclei physics: " << hn;
     // Hadron Elastic scattering
     RegisterPhysics(new G4HadronElasticPhysics(ver));
 

SimG4Core/PhysicsLists/src/CMSEmStandardPhysics.cc

@@ -4,13 +4,13 @@
 
 #include "G4SystemOfUnits.hh"
 #include "G4ParticleDefinition.hh"
+#include "G4LossTableManager.hh"
 #include "G4EmParameters.hh"
 #include "G4EmBuilder.hh"
 
 #include "G4ComptonScattering.hh"
 #include "G4GammaConversion.hh"
 #include "G4PhotoElectricEffect.hh"
-#include "G4LivermorePhotoElectricModel.hh"
 
 #include "G4MscStepLimitType.hh"
 
@@ -37,14 +37,12 @@
 #include "G4PhysicsListHelper.hh"
 #include "G4BuilderType.hh"
 #include "G4GammaGeneralProcess.hh"
-#include "G4LossTableManager.hh"
 
 #include "G4ProcessManager.hh"
 #include "G4TransportationWithMsc.hh"
 
 #include "G4RegionStore.hh"
 #include "G4Region.hh"
-#include <string>
 
 CMSEmStandardPhysics::CMSEmStandardPhysics(G4int ver, const edm::ParameterSet& p)
     : G4VPhysicsConstructor("CMSEmStandard_emm") {
@@ -136,8 +134,6 @@ void CMSEmStandardPhysics::ConstructProcess() {
   // e-
   particle = G4Electron::Electron();
 
-  G4eIonisation* eioni = new G4eIonisation();
-
   G4UrbanMscModel* msc1 = new G4UrbanMscModel();
   G4WentzelVIModel* msc2 = new G4WentzelVIModel();
   msc1->SetHighEnergyLimit(highEnergyLimit);
@@ -203,13 +199,12 @@ void CMSEmStandardPhysics::ConstructProcess() {
   ssm->SetLowEnergyLimit(highEnergyLimit);
   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
 
-  ph->RegisterProcess(eioni, particle);
+  ph->RegisterProcess(new G4eIonisation(), particle);
   ph->RegisterProcess(new G4eBremsstrahlung(), particle);
   ph->RegisterProcess(ss, particle);
 
   // e+
   particle = G4Positron::Positron();
-  eioni = new G4eIonisation();
 
   msc1 = new G4UrbanMscModel();
   msc2 = new G4WentzelVIModel();
@@ -274,7 +269,7 @@ void CMSEmStandardPhysics::ConstructProcess() {
   ssm->SetLowEnergyLimit(highEnergyLimit);
   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
 
-  ph->RegisterProcess(eioni, particle);
+  ph->RegisterProcess(new G4eIonisation(), particle);
   ph->RegisterProcess(new G4eBremsstrahlung(), particle);
   ph->RegisterProcess(new G4eplusAnnihilation(), particle);
   ph->RegisterProcess(ss, particle);

SimG4Core/PhysicsLists/src/CMSEmStandardPhysicsXS.cc

@@ -1,4 +1,5 @@
 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h"
+#include "SimG4Core/PhysicsLists/interface/CMSHepEmTrackingManager.h"
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
 #include "G4SystemOfUnits.hh"
@@ -10,35 +11,25 @@
 #include "G4ComptonScattering.hh"
 #include "G4GammaConversion.hh"
 #include "G4PhotoElectricEffect.hh"
-#include "G4RayleighScattering.hh"
-#include "G4PEEffectFluoModel.hh"
-#include "G4KleinNishinaModel.hh"
-#include "G4LowEPComptonModel.hh"
-#include "G4BetheHeitler5DModel.hh"
-#include "G4LivermorePhotoElectricModel.hh"
+
+#include "G4MscStepLimitType.hh"
 
 #include "G4eMultipleScattering.hh"
 #include "G4hMultipleScattering.hh"
-#include "G4MscStepLimitType.hh"
+#include "G4eCoulombScatteringModel.hh"
+#include "G4CoulombScattering.hh"
+#include "G4WentzelVIModel.hh"
 #include "G4UrbanMscModel.hh"
 #include "G4GoudsmitSaundersonMscModel.hh"
-#include "G4DummyModel.hh"
-#include "G4WentzelVIModel.hh"
-#include "G4CoulombScattering.hh"
-#include "G4eCoulombScatteringModel.hh"
 
 #include "G4eIonisation.hh"
 #include "G4eBremsstrahlung.hh"
-#include "G4Generator2BS.hh"
-#include "G4SeltzerBergerModel.hh"
-#include "G4ePairProduction.hh"
-#include "G4UniversalFluctuation.hh"
-
 #include "G4eplusAnnihilation.hh"
 
 #include "G4hIonisation.hh"
 #include "G4ionIonisation.hh"
 
+#include "G4ParticleTable.hh"
 #include "G4Gamma.hh"
 #include "G4Electron.hh"
 #include "G4Positron.hh"
@@ -53,9 +44,6 @@
 
 #include "G4RegionStore.hh"
 #include "G4Region.hh"
-#include "G4GammaGeneralProcess.hh"
-
-#include "G4SystemOfUnits.hh"
 
 CMSEmStandardPhysicsXS::CMSEmStandardPhysicsXS(G4int ver, const edm::ParameterSet& p)
     : G4VPhysicsConstructor("CMSEmStandard_emn") {
@@ -65,12 +53,10 @@ CMSEmStandardPhysicsXS::CMSEmStandardPhysicsXS(G4int ver, const edm::ParameterSe
   param->SetDefaults();
   param->SetVerbose(ver);
   param->SetApplyCuts(true);
-  param->SetMinEnergy(100 * CLHEP::eV);
-  param->SetNumberOfBinsPerDecade(20);
   param->SetStepFunction(0.8, 1 * CLHEP::mm);
   param->SetMscRangeFactor(0.2);
   param->SetMscStepLimitType(fMinimal);
-  param->SetFluo(true);
+  param->SetFluo(false);
   param->SetUseMottCorrection(true);  // use Mott-correction for e-/e+ msc gs
   SetPhysicsType(bElectromagnetic);
   fRangeFactor = p.getParameter<double>("G4MscRangeFactor");
@@ -88,6 +74,16 @@ CMSEmStandardPhysicsXS::CMSEmStandardPhysicsXS(G4int ver, const edm::ParameterSe
   double tcut = p.getParameter<double>("G4TrackingCut") * CLHEP::MeV;
   param->SetLowestElectronEnergy(tcut);
   param->SetLowestMuHadEnergy(tcut);
+  fG4HepEmActive = p.getParameter<bool>("G4HepEmActive");
+  if (fG4HepEmActive) {
+    // At the moment, G4HepEm supports only one configuration of MSC, so use
+    // the most generic parameters everywhere.
+    param->SetMscRangeFactor(fRangeFactor);
+    param->SetMscGeomFactor(fGeomFactor);
+    param->SetMscSafetyFactor(fSafetyFactor);
+    param->SetMscLambdaLimit(fLambdaLimit);
+    param->SetMscStepLimitType(fStepLimitType);
+  }
 }
 
 void CMSEmStandardPhysicsXS::ConstructParticle() {
@@ -112,40 +108,29 @@ void CMSEmStandardPhysicsXS::ConstructProcess() {
   G4NuclearStopping* pnuc(nullptr);
 
   // high energy limit for e+- scattering models
-  G4double highEnergyLimit = G4EmParameters::Instance()->MscEnergyLimit();
+  auto param = G4EmParameters::Instance();
+  G4double highEnergyLimit = param->MscEnergyLimit();
+
   const G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
   const G4Region* bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
 
   // Add gamma EM Processes
   G4ParticleDefinition* particle = G4Gamma::Gamma();
 
-  // Photoelectric
-  G4PhotoElectricEffect* pe = new G4PhotoElectricEffect();
-  G4VEmModel* theLivermorePEModel = new G4LivermorePhotoElectricModel();
-  pe->SetEmModel(theLivermorePEModel);
-
-  // Compton scattering
-  G4ComptonScattering* cs = new G4ComptonScattering;
-  cs->SetEmModel(new G4KleinNishinaModel());
+  G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
 
-  // Gamma conversion
-  G4GammaConversion* gc = new G4GammaConversion();
-  G4VEmModel* conv = new G4BetheHeitler5DModel();
-  gc->SetEmModel(conv);
-
-  if (G4EmParameters::Instance()->GeneralProcessActive()) {
+  if (param->GeneralProcessActive()) {
     G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
-    sp->AddEmProcess(pe);
-    sp->AddEmProcess(cs);
-    sp->AddEmProcess(gc);
-    sp->AddEmProcess(new G4RayleighScattering());
+    sp->AddEmProcess(pee);
+    sp->AddEmProcess(new G4ComptonScattering());
+    sp->AddEmProcess(new G4GammaConversion());
     G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
     ph->RegisterProcess(sp, particle);
+
   } else {
-    ph->RegisterProcess(pe, particle);
-    ph->RegisterProcess(cs, particle);
-    ph->RegisterProcess(gc, particle);
-    ph->RegisterProcess(new G4RayleighScattering(), particle);
+    ph->RegisterProcess(pee, particle);
+    ph->RegisterProcess(new G4ComptonScattering(), particle);
+    ph->RegisterProcess(new G4GammaConversion(), particle);
   }
 
   // e-
@@ -228,29 +213,15 @@ void CMSEmStandardPhysicsXS::ConstructProcess() {
   ssm->SetLowEnergyLimit(highEnergyLimit);
   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
 
-  // ionisation
-  G4eIonisation* eioni = new G4eIonisation();
-
-  // bremsstrahlung
-  G4eBremsstrahlung* brem = new G4eBremsstrahlung();
-  G4SeltzerBergerModel* br1 = new G4SeltzerBergerModel();
-  G4eBremsstrahlungRelModel* br2 = new G4eBremsstrahlungRelModel();
-  br1->SetAngularDistribution(new G4Generator2BS());
-  br2->SetAngularDistribution(new G4Generator2BS());
-  brem->SetEmModel(br1);
-  brem->SetEmModel(br2);
-  br1->SetHighEnergyLimit(CLHEP::GeV);
-
-  G4ePairProduction* ee = new G4ePairProduction();
-
   // register processes
-  ph->RegisterProcess(eioni, particle);
-  ph->RegisterProcess(brem, particle);
-  ph->RegisterProcess(ee, particle);
+  ph->RegisterProcess(new G4eIonisation(), particle);
+  ph->RegisterProcess(new G4eBremsstrahlung(), particle);
   ph->RegisterProcess(ss, particle);
 
   // e+
   particle = G4Positron::Positron();
+  msc3 = nullptr;
+  msc4 = nullptr;
 
   // multiple scattering
   msc1 = new G4UrbanMscModel();
@@ -326,26 +297,18 @@ void CMSEmStandardPhysicsXS::ConstructProcess() {
   ssm->SetLowEnergyLimit(highEnergyLimit);
   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
 
-  // ionisation
-  eioni = new G4eIonisation();
-
-  // bremsstrahlung
-  brem = new G4eBremsstrahlung();
-  br1 = new G4SeltzerBergerModel();
-  br2 = new G4eBremsstrahlungRelModel();
-  br1->SetAngularDistribution(new G4Generator2BS());
-  br2->SetAngularDistribution(new G4Generator2BS());
-  brem->SetEmModel(br1);
-  brem->SetEmModel(br2);
-  br1->SetHighEnergyLimit(CLHEP::GeV);
-
   // register processes
-  ph->RegisterProcess(eioni, particle);
-  ph->RegisterProcess(brem, particle);
-  ph->RegisterProcess(ee, particle);
+  ph->RegisterProcess(new G4eIonisation(), particle);
+  ph->RegisterProcess(new G4eBremsstrahlung(), particle);
   ph->RegisterProcess(new G4eplusAnnihilation(), particle);
   ph->RegisterProcess(ss, particle);
 
+  if (fG4HepEmActive) {
+    auto* hepEmTM = new CMSHepEmTrackingManager(highEnergyLimit);
+    G4Electron::Electron()->SetTrackingManager(hepEmTM);
+    G4Positron::Positron()->SetTrackingManager(hepEmTM);
+  }
+
   // generic ion
   particle = G4GenericIon::GenericIon();
   G4ionIonisation* ionIoni = new G4ionIonisation();