File indexing completed on 2023-01-03 23:24:49
0001 #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysics.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "G4SystemOfUnits.hh"
0005 #include "G4ParticleDefinition.hh"
0006 #include "G4EmParameters.hh"
0007 #include "G4EmBuilder.hh"
0008
0009 #include "G4ComptonScattering.hh"
0010 #include "G4GammaConversion.hh"
0011 #include "G4PhotoElectricEffect.hh"
0012 #include "G4LivermorePhotoElectricModel.hh"
0013
0014 #include "G4MscStepLimitType.hh"
0015
0016 #include "G4eMultipleScattering.hh"
0017 #include "G4hMultipleScattering.hh"
0018 #include "G4eCoulombScatteringModel.hh"
0019 #include "G4CoulombScattering.hh"
0020 #include "G4WentzelVIModel.hh"
0021 #include "G4UrbanMscModel.hh"
0022
0023 #include "G4eIonisation.hh"
0024 #include "G4eBremsstrahlung.hh"
0025 #include "G4eplusAnnihilation.hh"
0026
0027 #include "G4hIonisation.hh"
0028 #include "G4ionIonisation.hh"
0029
0030 #include "G4ParticleTable.hh"
0031 #include "G4Gamma.hh"
0032 #include "G4Electron.hh"
0033 #include "G4Positron.hh"
0034 #include "G4GenericIon.hh"
0035
0036 #include "G4PhysicsListHelper.hh"
0037 #include "G4BuilderType.hh"
0038 #include "G4GammaGeneralProcess.hh"
0039 #include "G4LossTableManager.hh"
0040
0041 #include "G4Version.hh"
0042 #if G4VERSION_NUMBER >= 1110
0043 #include "G4ProcessManager.hh"
0044 #include "G4TransportationWithMsc.hh"
0045 #endif
0046
0047 #include "G4RegionStore.hh"
0048 #include "G4Region.hh"
0049 #include <string>
0050
0051 CMSEmStandardPhysics::CMSEmStandardPhysics(G4int ver, const edm::ParameterSet& p)
0052 : G4VPhysicsConstructor("CMSEmStandard_emm") {
0053 SetVerboseLevel(ver);
0054
0055 G4EmParameters* param = G4EmParameters::Instance();
0056 param->SetDefaults();
0057 param->SetVerbose(ver);
0058 param->SetApplyCuts(true);
0059 param->SetStepFunction(0.8, 1 * CLHEP::mm);
0060 param->SetMscRangeFactor(0.2);
0061 param->SetMscStepLimitType(fMinimal);
0062 param->SetFluo(false);
0063 SetPhysicsType(bElectromagnetic);
0064 fRangeFactor = p.getParameter<double>("G4MscRangeFactor");
0065 fGeomFactor = p.getParameter<double>("G4MscGeomFactor");
0066 fSafetyFactor = p.getParameter<double>("G4MscSafetyFactor");
0067 fLambdaLimit = p.getParameter<double>("G4MscLambdaLimit") * CLHEP::mm;
0068 std::string msc = p.getParameter<std::string>("G4MscStepLimit");
0069 fStepLimitType = fUseSafety;
0070 if (msc == "UseSafetyPlus") {
0071 fStepLimitType = fUseSafetyPlus;
0072 }
0073 if (msc == "Minimal") {
0074 fStepLimitType = fMinimal;
0075 }
0076 double tcut = p.getParameter<double>("G4TrackingCut") * CLHEP::MeV;
0077 param->SetLowestElectronEnergy(tcut);
0078 param->SetLowestMuHadEnergy(tcut);
0079 }
0080
0081 void CMSEmStandardPhysics::ConstructParticle() {
0082
0083 G4EmBuilder::ConstructMinimalEmSet();
0084 }
0085
0086 void CMSEmStandardPhysics::ConstructProcess() {
0087 if (verboseLevel > 0) {
0088 edm::LogVerbatim("PhysicsList") << "### " << GetPhysicsName() << " Construct EM Processes";
0089 }
0090
0091
0092
0093
0094
0095 G4EmBuilder::PrepareEMPhysics();
0096
0097 G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
0098
0099 G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc");
0100 G4NuclearStopping* pnuc(nullptr);
0101
0102
0103 auto param = G4EmParameters::Instance();
0104 G4double highEnergyLimit = param->MscEnergyLimit();
0105
0106 const G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
0107 const G4Region* bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
0108
0109
0110 G4ParticleDefinition* particle = G4Gamma::Gamma();
0111
0112 G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
0113
0114 if (param->GeneralProcessActive()) {
0115 G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
0116 sp->AddEmProcess(pee);
0117 sp->AddEmProcess(new G4ComptonScattering());
0118 sp->AddEmProcess(new G4GammaConversion());
0119 G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
0120 ph->RegisterProcess(sp, particle);
0121
0122 } else {
0123 ph->RegisterProcess(pee, particle);
0124 ph->RegisterProcess(new G4ComptonScattering(), particle);
0125 ph->RegisterProcess(new G4GammaConversion(), particle);
0126 }
0127
0128
0129 particle = G4Electron::Electron();
0130
0131 G4eIonisation* eioni = new G4eIonisation();
0132
0133 G4UrbanMscModel* msc1 = new G4UrbanMscModel();
0134 G4WentzelVIModel* msc2 = new G4WentzelVIModel();
0135 msc1->SetHighEnergyLimit(highEnergyLimit);
0136 msc2->SetLowEnergyLimit(highEnergyLimit);
0137
0138
0139 G4UrbanMscModel* msc3 = nullptr;
0140 if (nullptr != aRegion || nullptr != bRegion) {
0141 msc3 = new G4UrbanMscModel();
0142 msc3->SetHighEnergyLimit(highEnergyLimit);
0143 msc3->SetRangeFactor(fRangeFactor);
0144 msc3->SetGeomFactor(fGeomFactor);
0145 msc3->SetSafetyFactor(fSafetyFactor);
0146 msc3->SetLambdaLimit(fLambdaLimit);
0147 msc3->SetStepLimitType(fStepLimitType);
0148 msc3->SetLocked(true);
0149 }
0150
0151 #if G4VERSION_NUMBER >= 1110
0152 G4TransportationWithMscType transportationWithMsc = param->TransportationWithMsc();
0153 if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0154
0155 G4ProcessManager* procManager = particle->GetProcessManager();
0156 G4VProcess* removed = procManager->RemoveProcess(0);
0157 if (removed->GetProcessName() != "Transportation") {
0158 G4Exception("CMSEmStandardPhysics::ConstructProcess",
0159 "em0050",
0160 FatalException,
0161 "replaced process is not G4Transportation!");
0162 }
0163 G4TransportationWithMsc* transportWithMsc =
0164 new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0165 if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0166 transportWithMsc->SetMultipleSteps(true);
0167 }
0168 transportWithMsc->AddMscModel(msc1);
0169 transportWithMsc->AddMscModel(msc2);
0170 if (nullptr != aRegion) {
0171 transportWithMsc->AddMscModel(msc3, -1, aRegion);
0172 }
0173 if (nullptr != bRegion) {
0174 transportWithMsc->AddMscModel(msc3, -1, bRegion);
0175 }
0176 procManager->AddProcess(transportWithMsc, -1, 0, 0);
0177 } else
0178 #endif
0179 {
0180
0181 G4eMultipleScattering* msc = new G4eMultipleScattering;
0182 msc->SetEmModel(msc1);
0183 msc->SetEmModel(msc2);
0184 if (nullptr != aRegion) {
0185 msc->AddEmModel(-1, msc3, aRegion);
0186 }
0187 if (nullptr != bRegion) {
0188 msc->AddEmModel(-1, msc3, bRegion);
0189 }
0190 ph->RegisterProcess(msc, particle);
0191 }
0192
0193
0194 G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
0195 G4CoulombScattering* ss = new G4CoulombScattering();
0196 ss->SetEmModel(ssm);
0197 ss->SetMinKinEnergy(highEnergyLimit);
0198 ssm->SetLowEnergyLimit(highEnergyLimit);
0199 ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0200
0201 ph->RegisterProcess(eioni, particle);
0202 ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0203 ph->RegisterProcess(ss, particle);
0204
0205
0206 particle = G4Positron::Positron();
0207 eioni = new G4eIonisation();
0208
0209 msc1 = new G4UrbanMscModel();
0210 msc2 = new G4WentzelVIModel();
0211 msc1->SetHighEnergyLimit(highEnergyLimit);
0212 msc2->SetLowEnergyLimit(highEnergyLimit);
0213
0214
0215 if (nullptr != aRegion || nullptr != bRegion) {
0216 msc3 = new G4UrbanMscModel();
0217 msc3->SetHighEnergyLimit(highEnergyLimit);
0218 msc3->SetRangeFactor(fRangeFactor);
0219 msc3->SetGeomFactor(fGeomFactor);
0220 msc3->SetSafetyFactor(fSafetyFactor);
0221 msc3->SetLambdaLimit(fLambdaLimit);
0222 msc3->SetStepLimitType(fStepLimitType);
0223 msc3->SetLocked(true);
0224 }
0225
0226 #if G4VERSION_NUMBER >= 1110
0227 if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
0228 G4ProcessManager* procManager = particle->GetProcessManager();
0229
0230 G4VProcess* removed = procManager->RemoveProcess(0);
0231 if (removed->GetProcessName() != "Transportation") {
0232 G4Exception("CMSEmStandardPhysics::ConstructProcess",
0233 "em0050",
0234 FatalException,
0235 "replaced process is not G4Transportation!");
0236 }
0237 G4TransportationWithMsc* transportWithMsc =
0238 new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
0239 if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
0240 transportWithMsc->SetMultipleSteps(true);
0241 }
0242 transportWithMsc->AddMscModel(msc1);
0243 transportWithMsc->AddMscModel(msc2);
0244 if (nullptr != aRegion) {
0245 transportWithMsc->AddMscModel(msc3, -1, aRegion);
0246 }
0247 if (nullptr != bRegion) {
0248 transportWithMsc->AddMscModel(msc3, -1, bRegion);
0249 }
0250 procManager->AddProcess(transportWithMsc, -1, 0, 0);
0251 } else
0252 #endif
0253 {
0254
0255 G4eMultipleScattering* msc = new G4eMultipleScattering;
0256 msc->SetEmModel(msc1);
0257 msc->SetEmModel(msc2);
0258 if (nullptr != aRegion) {
0259 msc->AddEmModel(-1, msc3, aRegion);
0260 }
0261 if (nullptr != bRegion) {
0262 msc->AddEmModel(-1, msc3, bRegion);
0263 }
0264 ph->RegisterProcess(msc, particle);
0265 }
0266
0267
0268 ssm = new G4eCoulombScatteringModel();
0269 ss = new G4CoulombScattering();
0270 ss->SetEmModel(ssm);
0271 ss->SetMinKinEnergy(highEnergyLimit);
0272 ssm->SetLowEnergyLimit(highEnergyLimit);
0273 ssm->SetActivationLowEnergyLimit(highEnergyLimit);
0274
0275 ph->RegisterProcess(eioni, particle);
0276 ph->RegisterProcess(new G4eBremsstrahlung(), particle);
0277 ph->RegisterProcess(new G4eplusAnnihilation(), particle);
0278 ph->RegisterProcess(ss, particle);
0279
0280
0281 particle = G4GenericIon::GenericIon();
0282 G4ionIonisation* ionIoni = new G4ionIonisation();
0283 ph->RegisterProcess(hmsc, particle);
0284 ph->RegisterProcess(ionIoni, particle);
0285
0286
0287 G4EmBuilder::ConstructCharged(hmsc, pnuc);
0288 }