Skip to content

Commit 2bdfa18

Browse files
cmsbuildcmsbuild
authored
Merge pull request #47464 from civanch/updated_config_g4hepem
[FULLSIM] Updated configuration of G4HepEm physics
2 parents 8b1962a + 8941bef commit 2bdfa18

File tree

1 file changed

+189
-161
lines changed

1 file changed

+189
-161
lines changed

‎SimG4Core/PhysicsLists/src/CMSEmStandardPhysics.cc‎

Lines changed: 189 additions & 161 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11
#include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysics.h"
2-
#include "SimG4Core/PhysicsLists/interface/CMSHepEmTrackingManager.h"
32
#include "FWCore/MessageLogger/interface/MessageLogger.h"
43

54
#include <CLHEP/Units/SystemOfUnits.h>
@@ -44,6 +43,9 @@
4443
#include "G4RegionStore.hh"
4544
#include "G4Region.hh"
4645

46+
#include "G4HepEmTrackingManager.hh"
47+
#include "G4HepEmConfig.hh"
48+
4749
CMSEmStandardPhysics::CMSEmStandardPhysics(G4int ver, const edm::ParameterSet& p)
4850
: G4VPhysicsConstructor("CMSEmStandard_emm") {
4951
SetVerboseLevel(ver);
@@ -78,15 +80,6 @@ CMSEmStandardPhysics::CMSEmStandardPhysics(G4int ver, const edm::ParameterSet& p
7880
edm::LogVerbatim("PhysicsList") << "EMM -> EMH: Forcing usage of G4HepEm";
7981
fG4HepEmActive = true;
8082
}
81-
if (fG4HepEmActive) {
82-
// At the moment, G4HepEm supports only one configuration of MSC, so use
83-
// the most generic parameters everywhere.
84-
param->SetMscRangeFactor(fRangeFactor);
85-
param->SetMscGeomFactor(fGeomFactor);
86-
param->SetMscSafetyFactor(fSafetyFactor);
87-
param->SetMscLambdaLimit(fLambdaLimit);
88-
param->SetMscStepLimitType(fStepLimitType);
89-
}
9083
}
9184

9285
void CMSEmStandardPhysics::ConstructParticle() {
@@ -117,179 +110,214 @@ void CMSEmStandardPhysics::ConstructProcess() {
117110
const G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion", false);
118111
const G4Region* bRegion = G4RegionStore::GetInstance()->GetRegion("HGCalRegion", false);
119112

120-
// Add gamma EM Processes
121-
G4ParticleDefinition* particle = G4Gamma::Gamma();
122-
123-
G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
113+
if (!fG4HepEmActive) {
114+
// Add gamma EM Processes
115+
G4ParticleDefinition* particle = G4Gamma::Gamma();
124116

125-
if (param->GeneralProcessActive()) {
126-
G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
127-
sp->AddEmProcess(pee);
128-
sp->AddEmProcess(new G4ComptonScattering());
129-
sp->AddEmProcess(new G4GammaConversion());
130-
G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
131-
ph->RegisterProcess(sp, particle);
117+
G4PhotoElectricEffect* pee = new G4PhotoElectricEffect();
132118

133-
} else {
134-
ph->RegisterProcess(pee, particle);
135-
ph->RegisterProcess(new G4ComptonScattering(), particle);
136-
ph->RegisterProcess(new G4GammaConversion(), particle);
137-
}
119+
if (param->GeneralProcessActive()) {
120+
G4GammaGeneralProcess* sp = new G4GammaGeneralProcess();
121+
sp->AddEmProcess(pee);
122+
sp->AddEmProcess(new G4ComptonScattering());
123+
sp->AddEmProcess(new G4GammaConversion());
124+
G4LossTableManager::Instance()->SetGammaGeneralProcess(sp);
125+
ph->RegisterProcess(sp, particle);
138126

139-
// e-
140-
particle = G4Electron::Electron();
141-
142-
G4UrbanMscModel* msc1 = new G4UrbanMscModel();
143-
G4WentzelVIModel* msc2 = new G4WentzelVIModel();
144-
msc1->SetHighEnergyLimit(highEnergyLimit);
145-
msc2->SetLowEnergyLimit(highEnergyLimit);
146-
147-
// e-/e+ msc for HCAL and HGCAL using the Urban model
148-
G4UrbanMscModel* msc3 = nullptr;
149-
if (nullptr != aRegion || nullptr != bRegion) {
150-
msc3 = new G4UrbanMscModel();
151-
msc3->SetHighEnergyLimit(highEnergyLimit);
152-
msc3->SetRangeFactor(fRangeFactor);
153-
msc3->SetGeomFactor(fGeomFactor);
154-
msc3->SetSafetyFactor(fSafetyFactor);
155-
msc3->SetLambdaLimit(fLambdaLimit);
156-
msc3->SetStepLimitType(fStepLimitType);
157-
msc3->SetLocked(true);
158-
}
127+
} else {
128+
ph->RegisterProcess(pee, particle);
129+
ph->RegisterProcess(new G4ComptonScattering(), particle);
130+
ph->RegisterProcess(new G4GammaConversion(), particle);
131+
}
159132

160-
G4TransportationWithMscType transportationWithMsc = param->TransportationWithMsc();
161-
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
162-
// Remove default G4Transportation and replace with G4TransportationWithMsc.
163-
G4ProcessManager* procManager = particle->GetProcessManager();
164-
G4VProcess* removed = procManager->RemoveProcess(0);
165-
if (removed->GetProcessName() != "Transportation") {
166-
G4Exception("CMSEmStandardPhysics::ConstructProcess",
167-
"em0050",
168-
FatalException,
169-
"replaced process is not G4Transportation!");
133+
// e-
134+
particle = G4Electron::Electron();
135+
136+
G4UrbanMscModel* msc1 = new G4UrbanMscModel();
137+
G4WentzelVIModel* msc2 = new G4WentzelVIModel();
138+
msc1->SetHighEnergyLimit(highEnergyLimit);
139+
msc2->SetLowEnergyLimit(highEnergyLimit);
140+
141+
// e-/e+ msc for HCAL and HGCAL using the Urban model
142+
G4UrbanMscModel* msc3 = nullptr;
143+
if (nullptr != aRegion || nullptr != bRegion) {
144+
msc3 = new G4UrbanMscModel();
145+
msc3->SetHighEnergyLimit(highEnergyLimit);
146+
msc3->SetRangeFactor(fRangeFactor);
147+
msc3->SetGeomFactor(fGeomFactor);
148+
msc3->SetSafetyFactor(fSafetyFactor);
149+
msc3->SetLambdaLimit(fLambdaLimit);
150+
msc3->SetStepLimitType(fStepLimitType);
151+
msc3->SetLocked(true);
170152
}
171-
G4TransportationWithMsc* transportWithMsc =
172-
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
173-
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
174-
transportWithMsc->SetMultipleSteps(true);
153+
154+
G4TransportationWithMscType transportationWithMsc = param->TransportationWithMsc();
155+
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
156+
// Remove default G4Transportation and replace with G4TransportationWithMsc.
157+
G4ProcessManager* procManager = particle->GetProcessManager();
158+
G4VProcess* removed = procManager->RemoveProcess(0);
159+
if (removed->GetProcessName() != "Transportation") {
160+
G4Exception("CMSEmStandardPhysics::ConstructProcess",
161+
"em0050",
162+
FatalException,
163+
"replaced process is not G4Transportation!");
164+
}
165+
G4TransportationWithMsc* transportWithMsc =
166+
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
167+
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
168+
transportWithMsc->SetMultipleSteps(true);
169+
}
170+
transportWithMsc->AddMscModel(msc1);
171+
transportWithMsc->AddMscModel(msc2);
172+
if (nullptr != aRegion) {
173+
transportWithMsc->AddMscModel(msc3, -1, aRegion);
174+
}
175+
if (nullptr != bRegion) {
176+
transportWithMsc->AddMscModel(msc3, -1, bRegion);
177+
}
178+
procManager->AddProcess(transportWithMsc, -1, 0, 0);
179+
} else {
180+
// Multiple scattering is registered as a separate process
181+
G4eMultipleScattering* msc = new G4eMultipleScattering;
182+
msc->SetEmModel(msc1);
183+
msc->SetEmModel(msc2);
184+
if (nullptr != aRegion) {
185+
msc->AddEmModel(-1, msc3, aRegion);
186+
}
187+
if (nullptr != bRegion) {
188+
msc->AddEmModel(-1, msc3, bRegion);
189+
}
190+
ph->RegisterProcess(msc, particle);
175191
}
176-
transportWithMsc->AddMscModel(msc1);
177-
transportWithMsc->AddMscModel(msc2);
178-
if (nullptr != aRegion) {
179-
transportWithMsc->AddMscModel(msc3, -1, aRegion);
192+
193+
// single scattering
194+
G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
195+
G4CoulombScattering* ss = new G4CoulombScattering();
196+
ss->SetEmModel(ssm);
197+
ss->SetMinKinEnergy(highEnergyLimit);
198+
ssm->SetLowEnergyLimit(highEnergyLimit);
199+
ssm->SetActivationLowEnergyLimit(highEnergyLimit);
200+
201+
ph->RegisterProcess(new G4eIonisation(), particle);
202+
ph->RegisterProcess(new G4eBremsstrahlung(), particle);
203+
ph->RegisterProcess(ss, particle);
204+
205+
// e+
206+
particle = G4Positron::Positron();
207+
208+
msc1 = new G4UrbanMscModel();
209+
msc2 = new G4WentzelVIModel();
210+
msc1->SetHighEnergyLimit(highEnergyLimit);
211+
msc2->SetLowEnergyLimit(highEnergyLimit);
212+
213+
// e-/e+ msc for HCAL and HGCAL using the Urban model
214+
if (nullptr != aRegion || nullptr != bRegion) {
215+
msc3 = new G4UrbanMscModel();
216+
msc3->SetHighEnergyLimit(highEnergyLimit);
217+
msc3->SetRangeFactor(fRangeFactor);
218+
msc3->SetGeomFactor(fGeomFactor);
219+
msc3->SetSafetyFactor(fSafetyFactor);
220+
msc3->SetLambdaLimit(fLambdaLimit);
221+
msc3->SetStepLimitType(fStepLimitType);
222+
msc3->SetLocked(true);
180223
}
181-
if (nullptr != bRegion) {
182-
transportWithMsc->AddMscModel(msc3, -1, bRegion);
224+
225+
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
226+
G4ProcessManager* procManager = particle->GetProcessManager();
227+
// Remove default G4Transportation and replace with G4TransportationWithMsc.
228+
G4VProcess* removed = procManager->RemoveProcess(0);
229+
if (removed->GetProcessName() != "Transportation") {
230+
G4Exception("CMSEmStandardPhysics::ConstructProcess",
231+
"em0050",
232+
FatalException,
233+
"replaced process is not G4Transportation!");
234+
}
235+
G4TransportationWithMsc* transportWithMsc =
236+
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
237+
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
238+
transportWithMsc->SetMultipleSteps(true);
239+
}
240+
transportWithMsc->AddMscModel(msc1);
241+
transportWithMsc->AddMscModel(msc2);
242+
if (nullptr != aRegion) {
243+
transportWithMsc->AddMscModel(msc3, -1, aRegion);
244+
}
245+
if (nullptr != bRegion) {
246+
transportWithMsc->AddMscModel(msc3, -1, bRegion);
247+
}
248+
procManager->AddProcess(transportWithMsc, -1, 0, 0);
249+
} else {
250+
// Register as a separate process.
251+
G4eMultipleScattering* msc = new G4eMultipleScattering;
252+
msc->SetEmModel(msc1);
253+
msc->SetEmModel(msc2);
254+
if (nullptr != aRegion) {
255+
msc->AddEmModel(-1, msc3, aRegion);
256+
}
257+
if (nullptr != bRegion) {
258+
msc->AddEmModel(-1, msc3, bRegion);
259+
}
260+
ph->RegisterProcess(msc, particle);
183261
}
184-
procManager->AddProcess(transportWithMsc, -1, 0, 0);
262+
263+
// single scattering
264+
ssm = new G4eCoulombScatteringModel();
265+
ss = new G4CoulombScattering();
266+
ss->SetEmModel(ssm);
267+
ss->SetMinKinEnergy(highEnergyLimit);
268+
ssm->SetLowEnergyLimit(highEnergyLimit);
269+
ssm->SetActivationLowEnergyLimit(highEnergyLimit);
270+
271+
ph->RegisterProcess(new G4eIonisation(), particle);
272+
ph->RegisterProcess(new G4eBremsstrahlung(), particle);
273+
ph->RegisterProcess(new G4eplusAnnihilation(), particle);
274+
ph->RegisterProcess(ss, particle);
275+
185276
} else {
186-
// Multiple scattering is registered as a separate process
187-
G4eMultipleScattering* msc = new G4eMultipleScattering;
188-
msc->SetEmModel(msc1);
189-
msc->SetEmModel(msc2);
190-
if (nullptr != aRegion) {
191-
msc->AddEmModel(-1, msc3, aRegion);
192-
}
193-
if (nullptr != bRegion) {
194-
msc->AddEmModel(-1, msc3, bRegion);
277+
// G4HepEm is Active
278+
if (verboseLevel > 0) {
279+
edm::LogVerbatim("PhysicsList") << "G4HepEm is active, registering G4HepEmTrackingManager";
195280
}
196-
ph->RegisterProcess(msc, particle);
197-
}
281+
auto* hepEmTM = new G4HepEmTrackingManager(verboseLevel);
198282

199-
// single scattering
200-
G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel();
201-
G4CoulombScattering* ss = new G4CoulombScattering();
202-
ss->SetEmModel(ssm);
203-
ss->SetMinKinEnergy(highEnergyLimit);
204-
ssm->SetLowEnergyLimit(highEnergyLimit);
205-
ssm->SetActivationLowEnergyLimit(highEnergyLimit);
206-
207-
ph->RegisterProcess(new G4eIonisation(), particle);
208-
ph->RegisterProcess(new G4eBremsstrahlung(), particle);
209-
ph->RegisterProcess(ss, particle);
210-
211-
// e+
212-
particle = G4Positron::Positron();
213-
214-
msc1 = new G4UrbanMscModel();
215-
msc2 = new G4WentzelVIModel();
216-
msc1->SetHighEnergyLimit(highEnergyLimit);
217-
msc2->SetLowEnergyLimit(highEnergyLimit);
218-
219-
// e-/e+ msc for HCAL and HGCAL using the Urban model
220-
if (nullptr != aRegion || nullptr != bRegion) {
221-
msc3 = new G4UrbanMscModel();
222-
msc3->SetHighEnergyLimit(highEnergyLimit);
223-
msc3->SetRangeFactor(fRangeFactor);
224-
msc3->SetGeomFactor(fGeomFactor);
225-
msc3->SetSafetyFactor(fSafetyFactor);
226-
msc3->SetLambdaLimit(fLambdaLimit);
227-
msc3->SetStepLimitType(fStepLimitType);
228-
msc3->SetLocked(true);
229-
}
283+
// configuration
284+
G4HepEmConfig* config = hepEmTM->GetConfig();
285+
// First set global configuration parameters:
286+
// ------------------------------------------
287+
// The default MSC `RangeFactor`, `SafetyFactor`, `StepLimitType` parameters
288+
// as well as `SetApplyCuts` and `SetLowestElectronEnergy` are taken from
289+
// `G4EmParameters` so we set only the step function parameters here.
230290

231-
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
232-
G4ProcessManager* procManager = particle->GetProcessManager();
233-
// Remove default G4Transportation and replace with G4TransportationWithMsc.
234-
G4VProcess* removed = procManager->RemoveProcess(0);
235-
if (removed->GetProcessName() != "Transportation") {
236-
G4Exception("CMSEmStandardPhysics::ConstructProcess",
237-
"em0050",
238-
FatalException,
239-
"replaced process is not G4Transportation!");
240-
}
241-
G4TransportationWithMsc* transportWithMsc =
242-
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
243-
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
244-
transportWithMsc->SetMultipleSteps(true);
245-
}
246-
transportWithMsc->AddMscModel(msc1);
247-
transportWithMsc->AddMscModel(msc2);
248-
if (nullptr != aRegion) {
249-
transportWithMsc->AddMscModel(msc3, -1, aRegion);
250-
}
251-
if (nullptr != bRegion) {
252-
transportWithMsc->AddMscModel(msc3, -1, bRegion);
253-
}
254-
procManager->AddProcess(transportWithMsc, -1, 0, 0);
255-
} else {
256-
// Register as a separate process.
257-
G4eMultipleScattering* msc = new G4eMultipleScattering;
258-
msc->SetEmModel(msc1);
259-
msc->SetEmModel(msc2);
291+
config->SetEnergyLossStepLimitFunctionParameters(0.8, 1.0 * CLHEP::mm);
292+
293+
// Then set special configuration for some regions:
294+
// ------------------------------------------------
260295
if (nullptr != aRegion) {
261-
msc->AddEmModel(-1, msc3, aRegion);
296+
// HCal region
297+
const G4String& rname = aRegion->GetName();
298+
config->SetMinimalMSCStepLimit(fStepLimitType == fMinimal, rname);
299+
config->SetMSCRangeFactor(fRangeFactor, rname);
300+
config->SetMSCSafetyFactor(fSafetyFactor, rname);
262301
}
302+
263303
if (nullptr != bRegion) {
264-
msc->AddEmModel(-1, msc3, bRegion);
304+
// HGCal region
305+
const G4String& rname = bRegion->GetName();
306+
config->SetMinimalMSCStepLimit(fStepLimitType == fMinimal, rname);
307+
config->SetMSCRangeFactor(fRangeFactor, rname);
308+
config->SetMSCSafetyFactor(fSafetyFactor, rname);
309+
310+
config->SetWoodcockTrackingRegion(rname);
311+
config->SetWDTEnergyLimit(0.5 * CLHEP::MeV);
265312
}
266-
ph->RegisterProcess(msc, particle);
267-
}
268-
269-
// single scattering
270-
ssm = new G4eCoulombScatteringModel();
271-
ss = new G4CoulombScattering();
272-
ss->SetEmModel(ssm);
273-
ss->SetMinKinEnergy(highEnergyLimit);
274-
ssm->SetLowEnergyLimit(highEnergyLimit);
275-
ssm->SetActivationLowEnergyLimit(highEnergyLimit);
276313

277-
ph->RegisterProcess(new G4eIonisation(), particle);
278-
ph->RegisterProcess(new G4eBremsstrahlung(), particle);
279-
ph->RegisterProcess(new G4eplusAnnihilation(), particle);
280-
ph->RegisterProcess(ss, particle);
281-
282-
if (fG4HepEmActive) {
283-
if (verboseLevel > 0) {
284-
edm::LogVerbatim("PhysicsList") << "G4HepEm is active, registering G4HepEmTrackingManager";
285-
}
286-
auto* hepEmTM = new CMSHepEmTrackingManager(highEnergyLimit);
287314
G4Electron::Electron()->SetTrackingManager(hepEmTM);
288315
G4Positron::Positron()->SetTrackingManager(hepEmTM);
316+
G4Gamma::Gamma()->SetTrackingManager(hepEmTM);
289317
}
290318

291319
// generic ion
292-
particle = G4GenericIon::GenericIon();
320+
G4ParticleDefinition* particle = G4GenericIon::GenericIon();
293321
G4ionIonisation* ionIoni = new G4ionIonisation();
294322
ph->RegisterProcess(hmsc, particle);
295323
ph->RegisterProcess(ionIoni, particle);

0 commit comments

Comments
 (0)