File indexing completed on 2023-01-03 23:24:50
0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.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 "G4PhysicsListHelper.hh"
0048 #include "G4BuilderType.hh"
0049 #include "G4GammaGeneralProcess.hh"
0050
0051 #include "G4Version.hh"
0052 #if G4VERSION_NUMBER >= 1110
0053 #include "G4ProcessManager.hh"
0054 #include "G4TransportationWithMsc.hh"
0055 #endif
0056
0057 #include "G4RegionStore.hh"
0058 #include "G4Region.hh"
0059 #include "G4GammaGeneralProcess.hh"
0060
0061 #include "G4SystemOfUnits.hh"
0062
0063 CMSEmStandardPhysicsXS::CMSEmStandardPhysicsXS(G4int ver, const edm::ParameterSet& p)
0064 : G4VPhysicsConstructor("CMSEmStandard_emn") {
0065 SetVerboseLevel(ver);
0066
0067 G4EmParameters* param = G4EmParameters::Instance();
0068 param->SetDefaults();
0069 param->SetVerbose(ver);
0070 param->SetApplyCuts(true);
0071 param->SetMinEnergy(100 * CLHEP::eV);
0072 param->SetNumberOfBinsPerDecade(20);
0073 param->SetStepFunction(0.8, 1 * CLHEP::mm);
0074 param->SetMscRangeFactor(0.2);
0075 param->SetMscStepLimitType(fMinimal);
0076 param->SetFluo(true);
0077 param->SetUseMottCorrection(true);
0078 SetPhysicsType(bElectromagnetic);
0079 fRangeFactor = p.getParameter<double>("G4MscRangeFactor");
0080 fGeomFactor = p.getParameter<double>("G4MscGeomFactor");
0081 fSafetyFactor = p.getParameter<double>("G4MscSafetyFactor");
0082 fLambdaLimit = p.getParameter<double>("G4MscLambdaLimit") * CLHEP::mm;
0083 std::string msc = p.getParameter<std::string>("G4MscStepLimit");
0084 fStepLimitType = fUseSafety;
0085 if (msc == "UseSafetyPlus") {
0086 fStepLimitType = fUseSafetyPlus;
0087 }
0088 if (msc == "Minimal") {
0089 fStepLimitType = fMinimal;
0090 }
0091 double tcut = p.getParameter<double>("G4TrackingCut") * CLHEP::MeV;
0092 param->SetLowestElectronEnergy(tcut);
0093 param->SetLowestMuHadEnergy(tcut);
0094 }
0095
0096 void CMSEmStandardPhysicsXS::ConstructParticle() {
0097
0098 G4EmBuilder::ConstructMinimalEmSet();
0099 }
0100
0101 void CMSEmStandardPhysicsXS::ConstructProcess() {
0102 if (verboseLevel > 0) {
0103 edm::LogVerbatim("PhysicsList") << "### " << GetPhysicsName() << " Construct Processes";
0104 }
0105
0106
0107
0108
0109 G4EmBuilder::PrepareEMPhysics();
0110
0111 G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0112
0113
0114 G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc");
0115 G4NuclearStopping* pnuc(nullptr);
0116
0117
0118 G4double highEnergyLimit = G4EmParameters::Instance()->MscEnergyLimit();
0119 const G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
0120 const G4Region* bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
0121
0122
0123 G4ParticleDefinition* particle = G4Gamma::Gamma();
0124
0125
0126 G4PhotoElectricEffect* pe = new G4PhotoElectricEffect();
0127 G4VEmModel* theLivermorePEModel = new G4LivermorePhotoElectricModel();
0128 pe->SetEmModel(theLivermorePEModel);
0129
0130
0131 G4ComptonScattering* cs = new G4ComptonScattering;
0132 cs->SetEmModel(new G4KleinNishinaModel());
0133
0134
0135 G4GammaConversion* gc = new G4GammaConversion();
0136 G4VEmModel* conv = new G4BetheHeitler5DModel();
0137 gc->SetEmModel(conv);
0138
0139 if (G4EmParameters::Instance()->GeneralProcessActive()) {
0140 G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
0141 sp->AddEmProcess(pe);
0142 sp->AddEmProcess(cs);
0143 sp->AddEmProcess(gc);
0144 sp->AddEmProcess(new G4RayleighScattering());
0145 G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
0146 ph->RegisterProcess(sp, particle);
0147 } else {
0148 ph->RegisterProcess(pe, particle);
0149 ph->RegisterProcess(cs, particle);
0150 ph->RegisterProcess(gc, particle);
0151 ph->RegisterProcess(new G4RayleighScattering(), particle);
0152 }
0153
0154
0155 particle = G4Electron::Electron();
0156
0157
0158 G4UrbanMscModel* msc1 = new G4UrbanMscModel();
0159 G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0160 msc1->SetHighEnergyLimit(highEnergyLimit);
0161 msc2->SetLowEnergyLimit(highEnergyLimit);
0162
0163
0164 G4UrbanMscModel* msc4 = nullptr;
0165 if (nullptr != aRegion) {
0166 msc4 = new G4UrbanMscModel();
0167 msc4->SetHighEnergyLimit(highEnergyLimit);
0168 msc4->SetRangeFactor(fRangeFactor);
0169 msc4->SetGeomFactor(fGeomFactor);
0170 msc4->SetSafetyFactor(fSafetyFactor);
0171 msc4->SetLambdaLimit(fLambdaLimit);
0172 msc4->SetStepLimitType(fStepLimitType);
0173 msc4->SetLocked(true);
0174 }
0175
0176
0177 G4GoudsmitSaundersonMscModel* msc3 = nullptr;
0178 if (nullptr != bRegion) {
0179 msc3 = new G4GoudsmitSaundersonMscModel();
0180 msc3->SetHighEnergyLimit(highEnergyLimit);
0181 msc3->SetRangeFactor(0.08);
0182 msc3->SetSkin(3.);
0183 msc3->SetStepLimitType(fUseSafetyPlus);
0184 msc3->SetLocked(true);
0185 }
0186
0187 #if G4VERSION_NUMBER >= 1110
0188 G4TransportationWithMscType transportationWithMsc = G4EmParameters::Instance()->TransportationWithMsc();
0189 if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0190 G4ProcessManager* procManager = particle->GetProcessManager();
0191
0192 G4VProcess* removed = procManager->RemoveProcess(0);
0193 if (removed->GetProcessName() != "Transportation") {
0194 G4Exception("CMSEmStandardPhysics::ConstructProcess",
0195 "em0050",
0196 FatalException,
0197 "replaced process is not G4Transportation!");
0198 }
0199 G4TransportationWithMsc* transportWithMsc =
0200 new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0201 if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0202 transportWithMsc->SetMultipleSteps(true);
0203 }
0204 transportWithMsc->AddMscModel(msc1);
0205 transportWithMsc->AddMscModel(msc2);
0206 if (msc4 != nullptr) {
0207 transportWithMsc->AddMscModel(msc4, -1, aRegion);
0208 }
0209 if (msc3 != nullptr) {
0210 transportWithMsc->AddMscModel(msc3, -1, bRegion);
0211 }
0212 procManager->AddProcess(transportWithMsc, -1, 0, 0);
0213 } else
0214 #endif
0215 {
0216
0217 G4eMultipleScattering* msc = new G4eMultipleScattering;
0218 msc->SetEmModel(msc1);
0219 msc->SetEmModel(msc2);
0220 if (msc4 != nullptr) {
0221 msc->AddEmModel(-1, msc4, aRegion);
0222 }
0223 if (msc3 != nullptr) {
0224 msc->AddEmModel(-1, msc3, bRegion);
0225 }
0226 ph->RegisterProcess(msc, particle);
0227 }
0228
0229
0230 G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0231 G4CoulombScattering* ss = new G4CoulombScattering();
0232 ss->SetEmModel(ssm);
0233 ss->SetMinKinEnergy(highEnergyLimit);
0234 ssm->SetLowEnergyLimit(highEnergyLimit);
0235 ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0236
0237
0238 G4eIonisation* eioni = new G4eIonisation();
0239
0240
0241 G4eBremsstrahlung* brem = new G4eBremsstrahlung();
0242 G4SeltzerBergerModel* br1 = new G4SeltzerBergerModel();
0243 G4eBremsstrahlungRelModel* 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 G4ePairProduction* ee = new G4ePairProduction();
0251
0252
0253 ph->RegisterProcess(eioni, particle);
0254 ph->RegisterProcess(brem, particle);
0255 ph->RegisterProcess(ee, particle);
0256 ph->RegisterProcess(ss, particle);
0257
0258
0259 particle = G4Positron::Positron();
0260
0261
0262 msc1 = new G4UrbanMscModel();
0263 msc2 = new G4WentzelVIModel();
0264 msc1->SetHighEnergyLimit(highEnergyLimit);
0265 msc2->SetLowEnergyLimit(highEnergyLimit);
0266
0267
0268 if (nullptr != aRegion) {
0269 msc4 = new G4UrbanMscModel();
0270 msc4->SetHighEnergyLimit(highEnergyLimit);
0271 msc4->SetRangeFactor(fRangeFactor);
0272 msc4->SetGeomFactor(fGeomFactor);
0273 msc4->SetSafetyFactor(fSafetyFactor);
0274 msc4->SetLambdaLimit(fLambdaLimit);
0275 msc4->SetStepLimitType(fStepLimitType);
0276 msc4->SetLocked(true);
0277 }
0278
0279
0280 if (nullptr != bRegion) {
0281 msc3 = new G4GoudsmitSaundersonMscModel();
0282 msc3->SetHighEnergyLimit(highEnergyLimit);
0283 msc3->SetRangeFactor(0.08);
0284 msc3->SetSkin(3.);
0285 msc3->SetStepLimitType(fUseSafetyPlus);
0286 msc3->SetLocked(true);
0287 }
0288
0289 #if G4VERSION_NUMBER >= 1110
0290 if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0291 G4ProcessManager* procManager = particle->GetProcessManager();
0292
0293 G4VProcess* removed = procManager->RemoveProcess(0);
0294 if (removed->GetProcessName() != "Transportation") {
0295 G4Exception("CMSEmStandardPhysics::ConstructProcess",
0296 "em0050",
0297 FatalException,
0298 "replaced process is not G4Transportation!");
0299 }
0300 G4TransportationWithMsc* transportWithMsc =
0301 new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0302 if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0303 transportWithMsc->SetMultipleSteps(true);
0304 }
0305 transportWithMsc->AddMscModel(msc1);
0306 transportWithMsc->AddMscModel(msc2);
0307 if (msc4 != nullptr) {
0308 transportWithMsc->AddMscModel(msc4, -1, aRegion);
0309 }
0310 if (msc3 != nullptr) {
0311 transportWithMsc->AddMscModel(msc3, -1, bRegion);
0312 }
0313 procManager->AddProcess(transportWithMsc, -1, 0, 0);
0314 } else
0315 #endif
0316 {
0317
0318 G4eMultipleScattering* msc = new G4eMultipleScattering;
0319 msc->SetEmModel(msc1);
0320 msc->SetEmModel(msc2);
0321 if (msc4 != nullptr) {
0322 msc->AddEmModel(-1, msc4, aRegion);
0323 }
0324 if (msc3 != nullptr) {
0325 msc->AddEmModel(-1, msc3, bRegion);
0326 }
0327 ph->RegisterProcess(msc, particle);
0328 }
0329
0330
0331 ssm = new G4eCoulombScatteringModel();
0332 ss = new G4CoulombScattering();
0333 ss->SetEmModel(ssm);
0334 ss->SetMinKinEnergy(highEnergyLimit);
0335 ssm->SetLowEnergyLimit(highEnergyLimit);
0336 ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0337
0338
0339 eioni = new G4eIonisation();
0340
0341
0342 brem = new G4eBremsstrahlung();
0343 br1 = new G4SeltzerBergerModel();
0344 br2 = new G4eBremsstrahlungRelModel();
0345 br1->SetAngularDistribution(new G4Generator2BS());
0346 br2->SetAngularDistribution(new G4Generator2BS());
0347 brem->SetEmModel(br1);
0348 brem->SetEmModel(br2);
0349 br1->SetHighEnergyLimit(CLHEP::GeV);
0350
0351
0352 ph->RegisterProcess(eioni, particle);
0353 ph->RegisterProcess(brem, particle);
0354 ph->RegisterProcess(ee, particle);
0355 ph->RegisterProcess(new G4eplusAnnihilation(), particle);
0356 ph->RegisterProcess(ss, particle);
0357
0358
0359 particle = G4GenericIon::GenericIon();
0360 G4ionIonisation* ionIoni = new G4ionIonisation();
0361 ph->RegisterProcess(hmsc, particle);
0362 ph->RegisterProcess(ionIoni, particle);
0363
0364
0365 G4EmBuilder::ConstructCharged(hmsc, pnuc);
0366 }