PhysicsLists/src/CMSEmStandardPhysicsEMZ.cc

0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsEMZ.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include <CLHEP/Units/SystemOfUnits.h>
0005 #include "G4ParticleDefinition.hh"
0006 #include "G4LossTableManager.hh"
0007 #include "G4EmParameters.hh"
0008 #include "G4EmBuilder.hh"
0009
0010 #include "G4ComptonScattering.hh"
0011 #include "G4GammaConversion.hh"
0012 #include "G4PhotoElectricEffect.hh"
0013 #include "G4RayleighScattering.hh"
0014 #include "G4PEEffectFluoModel.hh"
0015 #include "G4KleinNishinaModel.hh"
0016 #include "G4LowEPComptonModel.hh"
0017 #include "G4BetheHeitler5DModel.hh"
0018 #include "G4LivermorePhotoElectricModel.hh"
0019
0020 #include "G4eMultipleScattering.hh"
0021 #include "G4hMultipleScattering.hh"
0022 #include "G4MscStepLimitType.hh"
0023 #include "G4UrbanMscModel.hh"
0024 #include "G4GoudsmitSaundersonMscModel.hh"
0025 #include "G4DummyModel.hh"
0026 #include "G4WentzelVIModel.hh"
0027 #include "G4CoulombScattering.hh"
0028 #include "G4eCoulombScatteringModel.hh"
0029
0030 #include "G4eIonisation.hh"
0031 #include "G4eBremsstrahlung.hh"
0032 #include "G4Generator2BS.hh"
0033 #include "G4SeltzerBergerModel.hh"
0034 #include "G4ePairProduction.hh"
0035 #include "G4UniversalFluctuation.hh"
0036
0037 #include "G4eplusAnnihilation.hh"
0038
0039 #include "G4hIonisation.hh"
0040 #include "G4ionIonisation.hh"
0041
0042 #include "G4Gamma.hh"
0043 #include "G4Electron.hh"
0044 #include "G4Positron.hh"
0045 #include "G4GenericIon.hh"
0046
0047 #include "G4IonParametrisedLossModel.hh"
0048 #include "G4NuclearStopping.hh"
0049 #include "G4ePairProduction.hh"
0050 #include "G4LivermoreIonisationModel.hh"
0051 #include "G4PenelopeIonisationModel.hh"
0052
0053 #include "G4PhysicsListHelper.hh"
0054 #include "G4BuilderType.hh"
0055 #include "G4GammaGeneralProcess.hh"
0056
0057 #include "G4RegionStore.hh"
0058 #include "G4Region.hh"
0059 #include "G4GammaGeneralProcess.hh"
0060
0061 #include <CLHEP/Units/SystemOfUnits.h>
0062
0063 CMSEmStandardPhysicsEMZ::CMSEmStandardPhysicsEMZ(G4int ver, const edm::ParameterSet& p)
0064     : G4VPhysicsConstructor("CMSEmStandard_emz") {
0065   SetVerboseLevel(ver);
0066   // EM parameters specific for this EM physics configuration
0067   G4EmParameters* param = G4EmParameters::Instance();
0068   param->SetDefaults();
0069   param->SetVerbose(ver);
0070   param->SetMinEnergy(100 * CLHEP::eV);
0071   param->SetLowestElectronEnergy(100 * CLHEP::eV);
0072   param->SetNumberOfBinsPerDecade(20);
0073   param->ActivateAngularGeneratorForIonisation(true);
0074   param->SetStepFunction(0.2, 10 * CLHEP::um);
0075   param->SetStepFunctionMuHad(0.1, 50 * CLHEP::um);
0076   param->SetStepFunctionLightIons(0.1, 20 * CLHEP::um);
0077   param->SetStepFunctionIons(0.1, 1 * CLHEP::um);
0078   param->SetUseMottCorrection(true);           // use Mott-correction for e-/e+ msc gs
0079   param->SetMscStepLimitType(fUseSafetyPlus);  // for e-/e+ msc gs
0080   param->SetMscSkin(3);                        // error-free stepping for e-/e+ msc gs
0081   param->SetMscRangeFactor(0.08);              // error-free stepping for e-/e+ msc gs
0082   param->SetMuHadLateralDisplacement(true);
0083   param->SetFluo(true);
0084   param->SetUseICRU90Data(true);
0085   param->SetMaxNIELEnergy(1 * CLHEP::MeV);
0086   double tcut = p.getParameter<double>("G4TrackingCut") * CLHEP::MeV;
0087   param->SetLowestElectronEnergy(tcut);
0088   param->SetLowestMuHadEnergy(tcut);
0089   SetPhysicsType(bElectromagnetic);
0090 }
0091
0092 void CMSEmStandardPhysicsEMZ::ConstructParticle() {
0093   // minimal set of particles for EM physics
0094   G4EmBuilder::ConstructMinimalEmSet();
0095 }
0096
0097 void CMSEmStandardPhysicsEMZ::ConstructProcess() {
0098   if (verboseLevel > 0) {
0099     edm::LogVerbatim("PhysicsList") << "### " << GetPhysicsName() << " Construct Processes";
0100   }
0101
0102   // This EM builder takes default models of Geant4 10 EMV.
0103   // Multiple scattering by Urban for all particles
0104   // except e+e- below 100 MeV for which the Urban model is used
0105   G4EmBuilder::PrepareEMPhysics();
0106
0107   G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0108
0109   // processes used by several particles
0110   G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc");
0111   G4double nielEnergyLimit = G4EmParameters::Instance()->MaxNIELEnergy();
0112   G4NuclearStopping* pnuc = nullptr;
0113   if (nielEnergyLimit > 0.0) {
0114     pnuc = new G4NuclearStopping();
0115     pnuc->SetMaxKinEnergy(nielEnergyLimit);
0116   }
0117
0118   // high energy limit for e+- scattering models
0119   G4double highEnergyLimit = G4EmParameters::Instance()->MscEnergyLimit();
0120
0121   // Add gamma EM Processes
0122   G4ParticleDefinition* particle = G4Gamma::Gamma();
0123
0124   // Photoelectric
0125   G4PhotoElectricEffect* pe = new G4PhotoElectricEffect();
0126   G4VEmModel* theLivermorePEModel = new G4LivermorePhotoElectricModel();
0127   pe->SetEmModel(theLivermorePEModel);
0128
0129   // Compton scattering
0130   G4ComptonScattering* cs = new G4ComptonScattering;
0131   cs->SetEmModel(new G4KleinNishinaModel());
0132   G4VEmModel* theLowEPComptonModel = new G4LowEPComptonModel();
0133   theLowEPComptonModel->SetHighEnergyLimit(20 * CLHEP::MeV);
0134   cs->AddEmModel(0, theLowEPComptonModel);
0135
0136   // Gamma conversion
0137   G4GammaConversion* gc = new G4GammaConversion();
0138   G4VEmModel* conv = new G4BetheHeitler5DModel();
0139   gc->SetEmModel(conv);
0140
0141   if (G4EmParameters::Instance()->GeneralProcessActive()) {
0142     G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
0143     sp->AddEmProcess(pe);
0144     sp->AddEmProcess(cs);
0145     sp->AddEmProcess(gc);
0146     sp->AddEmProcess(new G4RayleighScattering());
0147     G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
0148     ph->RegisterProcess(sp, particle);
0149   } else {
0150     ph->RegisterProcess(pe, particle);
0151     ph->RegisterProcess(cs, particle);
0152     ph->RegisterProcess(gc, particle);
0153     ph->RegisterProcess(new G4RayleighScattering(), particle);
0154   }
0155
0156   // e-
0157   particle = G4Electron::Electron();
0158
0159   // multiple scattering
0160   G4eMultipleScattering* msc = new G4eMultipleScattering();
0161   // e-/e+ msc gs with Mott-correction
0162   // (Mott-correction is set through G4EmParameters)
0163   G4GoudsmitSaundersonMscModel* msc1 = new G4GoudsmitSaundersonMscModel();
0164   G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0165   msc1->SetHighEnergyLimit(highEnergyLimit);
0166   msc2->SetLowEnergyLimit(highEnergyLimit);
0167   msc->SetEmModel(msc1);
0168   msc->SetEmModel(msc2);
0169
0170   G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0171   G4CoulombScattering* ss = new G4CoulombScattering();
0172   ss->SetEmModel(ssm);
0173   ss->SetMinKinEnergy(highEnergyLimit);
0174   ssm->SetLowEnergyLimit(highEnergyLimit);
0175   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0176
0177   // single scattering
0178   ssm = new G4eCoulombScatteringModel();
0179   ss = new G4CoulombScattering();
0180   ss->SetEmModel(ssm);
0181   ss->SetMinKinEnergy(highEnergyLimit);
0182   ssm->SetLowEnergyLimit(highEnergyLimit);
0183   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0184
0185   // ionisation
0186   G4eIonisation* eioni = new G4eIonisation();
0187   G4VEmModel* theIoniLiv = new G4LivermoreIonisationModel();
0188   theIoniLiv->SetHighEnergyLimit(0.1 * CLHEP::MeV);
0189   eioni->AddEmModel(0, theIoniLiv, new G4UniversalFluctuation());
0190
0191   // bremsstrahlung
0192   G4eBremsstrahlung* brem = new G4eBremsstrahlung();
0193   G4SeltzerBergerModel* br1 = new G4SeltzerBergerModel();
0194   G4eBremsstrahlungRelModel* br2 = new G4eBremsstrahlungRelModel();
0195   br1->SetAngularDistribution(new G4Generator2BS());
0196   br2->SetAngularDistribution(new G4Generator2BS());
0197   brem->SetEmModel(br1);
0198   brem->SetEmModel(br2);
0199   br1->SetHighEnergyLimit(CLHEP::GeV);
0200
0201   G4ePairProduction* ee = new G4ePairProduction();
0202
0203   // register processes
0204   ph->RegisterProcess(msc, particle);
0205   ph->RegisterProcess(eioni, particle);
0206   ph->RegisterProcess(brem, particle);
0207   ph->RegisterProcess(ee, particle);
0208   ph->RegisterProcess(ss, particle);
0209
0210   // e+
0211   particle = G4Positron::Positron();
0212
0213   // multiple scattering
0214   msc = new G4eMultipleScattering();
0215   // e-/e+ msc gs with Mott-correction
0216   // (Mott-correction is set through G4EmParameters)
0217   msc1 = new G4GoudsmitSaundersonMscModel();
0218   msc2 = new G4WentzelVIModel();
0219   msc1->SetHighEnergyLimit(highEnergyLimit);
0220   msc2->SetLowEnergyLimit(highEnergyLimit);
0221   msc->SetEmModel(msc1);
0222   msc->SetEmModel(msc2);
0223
0224   // single scattering
0225   ssm = new G4eCoulombScatteringModel();
0226   ss = new G4CoulombScattering();
0227   ss->SetEmModel(ssm);
0228   ss->SetMinKinEnergy(highEnergyLimit);
0229   ssm->SetLowEnergyLimit(highEnergyLimit);
0230   ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0231
0232   // ionisation
0233   eioni = new G4eIonisation();
0234   /*
0235   G4VEmModel* pen = new G4PenelopeIonisationModel();
0236   pen->SetHighEnergyLimit(0.1*CLHEP::MeV);
0237   eioni->AddEmModel(0, pen, new G4UniversalFluctuation());
0238   */
0239   // bremsstrahlung
0240   brem = new G4eBremsstrahlung();
0241   br1 = new G4SeltzerBergerModel();
0242   br2 = new G4eBremsstrahlungRelModel();
0243   br1->SetAngularDistribution(new G4Generator2BS());
0244   br2->SetAngularDistribution(new G4Generator2BS());
0245   brem->SetEmModel(br1);
0246   brem->SetEmModel(br2);
0247   br1->SetHighEnergyLimit(CLHEP::GeV);
0248
0249   // register processes
0250   ph->RegisterProcess(msc, particle);
0251   ph->RegisterProcess(eioni, particle);
0252   ph->RegisterProcess(brem, particle);
0253   ph->RegisterProcess(ee, particle);
0254   ph->RegisterProcess(new G4eplusAnnihilation(), particle);
0255   ph->RegisterProcess(ss, particle);
0256
0257   // generic ion
0258   particle = G4GenericIon::GenericIon();
0259   G4ionIonisation* ionIoni = new G4ionIonisation();
0260   ionIoni->SetEmModel(new G4IonParametrisedLossModel());
0261   ph->RegisterProcess(hmsc, particle);
0262   ph->RegisterProcess(ionIoni, particle);
0263   if (nullptr != pnuc)
0264     ph->RegisterProcess(pnuc, particle);
0265
0266   // muons, hadrons, ions
0267   G4EmBuilder::ConstructCharged(hmsc, pnuc);
0268 }