Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-15 01:08:02

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