Skip to content

Commit

Permalink
Introduce option to use G4TransportationWithMsc
Browse files Browse the repository at this point in the history 
This process combines G4Transportation and multiple scattering, which
allows to move directly from one discrete interaction point to the
next, with all intermediate steps due to MSC combined into one G4Step.
  • Loading branch information
@hahnjo
hahnjo committed Jul 28, 2022
1 parent ca284ed commit 4208c8c
Show file tree
Hide file tree
Showing 2 changed files with 173 additions and 26 deletions.
88 changes: 78 additions & 10 deletions SimG4Core/PhysicsLists/src/CMSEmStandardPhysics.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,12 @@
#include "G4GammaGeneralProcess.hh"
#include "G4LossTableManager.hh"

#include "G4Version.hh"
#if G4VERSION_NUMBER >= 1110
#include "G4ProcessManager.hh"
#include "G4TransportationWithMsc.hh"
#endif

#include "G4RegionStore.hh"
#include "G4Region.hh"
#include <string>
Expand Down Expand Up @@ -121,31 +127,64 @@ void CMSEmStandardPhysics::ConstructProcess() {

G4eIonisation* eioni = new G4eIonisation();

G4eMultipleScattering* msc = new G4eMultipleScattering;
G4UrbanMscModel* msc1 = new G4UrbanMscModel();
G4WentzelVIModel* msc2 = new G4WentzelVIModel();
msc1->SetHighEnergyLimit(highEnergyLimit);
msc2->SetLowEnergyLimit(highEnergyLimit);
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);

// e-/e+ msc for HCAL and HGCAL using the Urban model
G4UrbanMscModel* msc3 = nullptr;
if (nullptr != aRegion || nullptr != bRegion) {
G4UrbanMscModel* msc3 = new G4UrbanMscModel();
msc3 = new G4UrbanMscModel();
msc3->SetHighEnergyLimit(highEnergyLimit);
msc3->SetRangeFactor(fRangeFactor);
msc3->SetGeomFactor(fGeomFactor);
msc3->SetSafetyFactor(fSafetyFactor);
msc3->SetLambdaLimit(fLambdaLimit);
msc3->SetStepLimitType(fStepLimitType);
msc3->SetLocked(true);
}

#if G4VERSION_NUMBER >= 1110
G4TransportationWithMscType transportationWithMsc = G4EmParameters::Instance()->TransportationWithMsc();
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
G4ProcessManager* procManager = particle->GetProcessManager();
// Remove default G4Transportation and replace with G4TransportationWithMsc.
G4VProcess* removed = procManager->RemoveProcess(0);
if (removed->GetProcessName() != "Transportation") {
G4Exception("CMSEmStandardPhysics::ConstructProcess",
"em0050",
FatalException,
"replaced process is not G4Transportation!");
}
G4TransportationWithMsc* transportWithMsc =
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
transportWithMsc->SetMultipleSteps(true);
}
transportWithMsc->AddMscModel(msc1);
transportWithMsc->AddMscModel(msc2);
if (nullptr != aRegion) {
transportWithMsc->AddMscModel(msc3, -1, aRegion);
}
if (nullptr != bRegion) {
transportWithMsc->AddMscModel(msc3, -1, bRegion);
}
procManager->AddProcess(transportWithMsc, -1, 0, 0);
} else
#endif
{
// Register as a separate process.
G4eMultipleScattering* msc = new G4eMultipleScattering;
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);
if (nullptr != aRegion) {
msc->AddEmModel(-1, msc3, aRegion);
}
if (nullptr != bRegion) {
msc->AddEmModel(-1, msc3, bRegion);
}
ph->RegisterProcess(msc, particle);
}

// single scattering
Expand All @@ -156,7 +195,6 @@ void CMSEmStandardPhysics::ConstructProcess() {
ssm->SetLowEnergyLimit(highEnergyLimit);
ssm->SetActivationLowEnergyLimit(highEnergyLimit);

ph->RegisterProcess(msc, particle);
ph->RegisterProcess(eioni, particle);
ph->RegisterProcess(new G4eBremsstrahlung(), particle);
ph->RegisterProcess(ss, particle);
Expand All @@ -165,31 +203,62 @@ void CMSEmStandardPhysics::ConstructProcess() {
particle = G4Positron::Positron();
eioni = new G4eIonisation();

msc = new G4eMultipleScattering();
msc1 = new G4UrbanMscModel();
msc2 = new G4WentzelVIModel();
msc1->SetHighEnergyLimit(highEnergyLimit);
msc2->SetLowEnergyLimit(highEnergyLimit);
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);

// e-/e+ msc for HCAL and HGCAL using the Urban model
if (nullptr != aRegion || nullptr != bRegion) {
G4UrbanMscModel* msc3 = new G4UrbanMscModel();
msc3 = new G4UrbanMscModel();
msc3->SetHighEnergyLimit(highEnergyLimit);
msc3->SetRangeFactor(fRangeFactor);
msc3->SetGeomFactor(fGeomFactor);
msc3->SetSafetyFactor(fSafetyFactor);
msc3->SetLambdaLimit(fLambdaLimit);
msc3->SetStepLimitType(fStepLimitType);
msc3->SetLocked(true);
}

#if G4VERSION_NUMBER >= 1110
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
G4ProcessManager* procManager = particle->GetProcessManager();
// Remove default G4Transportation and replace with G4TransportationWithMsc.
G4VProcess* removed = procManager->RemoveProcess(0);
if (removed->GetProcessName() != "Transportation") {
G4Exception("CMSEmStandardPhysics::ConstructProcess",
"em0050",
FatalException,
"replaced process is not G4Transportation!");
}
G4TransportationWithMsc* transportWithMsc =
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
transportWithMsc->SetMultipleSteps(true);
}
transportWithMsc->AddMscModel(msc1);
transportWithMsc->AddMscModel(msc2);
if (nullptr != aRegion) {
transportWithMsc->AddMscModel(msc3, -1, aRegion);
}
if (nullptr != bRegion) {
transportWithMsc->AddMscModel(msc3, -1, bRegion);
}
procManager->AddProcess(transportWithMsc, -1, 0, 0);
} else
#endif
{
// Register as a separate process.
G4eMultipleScattering* msc = new G4eMultipleScattering;
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);
if (nullptr != aRegion) {
msc->AddEmModel(-1, msc3, aRegion);
}
if (nullptr != bRegion) {
msc->AddEmModel(-1, msc3, bRegion);
}
ph->RegisterProcess(msc, particle);
}

// single scattering
Expand All @@ -200,7 +269,6 @@ void CMSEmStandardPhysics::ConstructProcess() {
ssm->SetLowEnergyLimit(highEnergyLimit);
ssm->SetActivationLowEnergyLimit(highEnergyLimit);

ph->RegisterProcess(msc, particle);
ph->RegisterProcess(eioni, particle);
ph->RegisterProcess(new G4eBremsstrahlung(), particle);
ph->RegisterProcess(new G4eplusAnnihilation(), particle);
Expand Down
111 changes: 95 additions & 16 deletions SimG4Core/PhysicsLists/src/CMSEmStandardPhysicsXS.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,12 @@
#include "G4BuilderType.hh"
#include "G4GammaGeneralProcess.hh"

#include "G4Version.hh"
#if G4VERSION_NUMBER >= 1110
#include "G4ProcessManager.hh"
#include "G4TransportationWithMsc.hh"
#endif

#include "G4RegionStore.hh"
#include "G4Region.hh"
#include "G4GammaGeneralProcess.hh"
Expand Down Expand Up @@ -147,36 +153,75 @@ void CMSEmStandardPhysicsXS::ConstructProcess() {
particle = G4Electron::Electron();

// multiple scattering
G4eMultipleScattering* msc = new G4eMultipleScattering();
G4UrbanMscModel* msc1 = new G4UrbanMscModel();
G4WentzelVIModel* msc2 = new G4WentzelVIModel();
msc1->SetHighEnergyLimit(highEnergyLimit);
msc2->SetLowEnergyLimit(highEnergyLimit);
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);

// msc for HCAL using the Urban model
G4UrbanMscModel* msc4 = nullptr;
if (nullptr != aRegion) {
G4UrbanMscModel* msc4 = new G4UrbanMscModel();
msc4 = new G4UrbanMscModel();
msc4->SetHighEnergyLimit(highEnergyLimit);
msc4->SetRangeFactor(fRangeFactor);
msc4->SetGeomFactor(fGeomFactor);
msc4->SetSafetyFactor(fSafetyFactor);
msc4->SetLambdaLimit(fLambdaLimit);
msc4->SetStepLimitType(fStepLimitType);
msc4->SetLocked(true);
msc->AddEmModel(-1, msc4, aRegion);
}

// msc GS with Mott-correction
G4GoudsmitSaundersonMscModel* msc3 = nullptr;
if (nullptr != bRegion) {
G4GoudsmitSaundersonMscModel* msc3 = new G4GoudsmitSaundersonMscModel();
msc3 = new G4GoudsmitSaundersonMscModel();
msc3->SetHighEnergyLimit(highEnergyLimit);
msc3->SetRangeFactor(0.08);
msc3->SetSkin(3.);
msc3->SetStepLimitType(fUseSafetyPlus);
msc3->SetLocked(true);
msc->AddEmModel(-1, msc3, bRegion);
}

#if G4VERSION_NUMBER >= 1110
G4TransportationWithMscType transportationWithMsc = G4EmParameters::Instance()->TransportationWithMsc();
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
G4ProcessManager* procManager = particle->GetProcessManager();
// Remove default G4Transportation and replace with G4TransportationWithMsc.
G4VProcess* removed = procManager->RemoveProcess(0);
if (removed->GetProcessName() != "Transportation") {
G4Exception("CMSEmStandardPhysics::ConstructProcess",
"em0050",
FatalException,
"replaced process is not G4Transportation!");
}
G4TransportationWithMsc* transportWithMsc =
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
transportWithMsc->SetMultipleSteps(true);
}
transportWithMsc->AddMscModel(msc1);
transportWithMsc->AddMscModel(msc2);
if (msc4 != nullptr) {
transportWithMsc->AddMscModel(msc4, -1, aRegion);
}
if (msc3 != nullptr) {
transportWithMsc->AddMscModel(msc3, -1, bRegion);
}
procManager->AddProcess(transportWithMsc, -1, 0, 0);
} else
#endif
{
// Register as a separate process.
G4eMultipleScattering* msc = new G4eMultipleScattering;
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);
if (msc4 != nullptr) {
msc->AddEmModel(-1, msc4, aRegion);
}
if (msc3 != nullptr) {
msc->AddEmModel(-1, msc3, bRegion);
}
ph->RegisterProcess(msc, particle);
}

// single scattering
Expand All @@ -203,7 +248,6 @@ void CMSEmStandardPhysicsXS::ConstructProcess() {
G4ePairProduction* ee = new G4ePairProduction();

// register processes
ph->RegisterProcess(msc, particle);
ph->RegisterProcess(eioni, particle);
ph->RegisterProcess(brem, particle);
ph->RegisterProcess(ee, particle);
Expand All @@ -213,36 +257,72 @@ void CMSEmStandardPhysicsXS::ConstructProcess() {
particle = G4Positron::Positron();

// multiple scattering
msc = new G4eMultipleScattering();
msc1 = new G4UrbanMscModel();
msc2 = new G4WentzelVIModel();
msc1->SetHighEnergyLimit(highEnergyLimit);
msc2->SetLowEnergyLimit(highEnergyLimit);
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);

// msc for HCAL using the Urban model
if (nullptr != aRegion) {
G4UrbanMscModel* msc4 = new G4UrbanMscModel();
msc4 = new G4UrbanMscModel();
msc4->SetHighEnergyLimit(highEnergyLimit);
msc4->SetRangeFactor(fRangeFactor);
msc4->SetGeomFactor(fGeomFactor);
msc4->SetSafetyFactor(fSafetyFactor);
msc4->SetLambdaLimit(fLambdaLimit);
msc4->SetStepLimitType(fStepLimitType);
msc4->SetLocked(true);
msc->AddEmModel(-1, msc4, aRegion);
}

// msc GS with Mott-correction
if (nullptr != bRegion) {
G4GoudsmitSaundersonMscModel* msc3 = new G4GoudsmitSaundersonMscModel();
msc3 = new G4GoudsmitSaundersonMscModel();
msc3->SetHighEnergyLimit(highEnergyLimit);
msc3->SetRangeFactor(0.08);
msc3->SetSkin(3.);
msc3->SetStepLimitType(fUseSafetyPlus);
msc3->SetLocked(true);
msc->AddEmModel(-1, msc3, bRegion);
}

#if G4VERSION_NUMBER >= 1110
if (transportationWithMsc != G4TransportationWithMscType::fDisabled) {
G4ProcessManager* procManager = particle->GetProcessManager();
// Remove default G4Transportation and replace with G4TransportationWithMsc.
G4VProcess* removed = procManager->RemoveProcess(0);
if (removed->GetProcessName() != "Transportation") {
G4Exception("CMSEmStandardPhysics::ConstructProcess",
"em0050",
FatalException,
"replaced process is not G4Transportation!");
}
G4TransportationWithMsc* transportWithMsc =
new G4TransportationWithMsc(G4TransportationWithMsc::ScatteringType::MultipleScattering);
if (transportationWithMsc == G4TransportationWithMscType::fMultipleSteps) {
transportWithMsc->SetMultipleSteps(true);
}
transportWithMsc->AddMscModel(msc1);
transportWithMsc->AddMscModel(msc2);
if (msc4 != nullptr) {
transportWithMsc->AddMscModel(msc4, -1, aRegion);
}
if (msc3 != nullptr) {
transportWithMsc->AddMscModel(msc3, -1, bRegion);
}
procManager->AddProcess(transportWithMsc, -1, 0, 0);
} else
#endif
{
// Register as a separate process.
G4eMultipleScattering* msc = new G4eMultipleScattering;
msc->SetEmModel(msc1);
msc->SetEmModel(msc2);
if (msc4 != nullptr) {
msc->AddEmModel(-1, msc4, aRegion);
}
if (msc3 != nullptr) {
msc->AddEmModel(-1, msc3, bRegion);
}
ph->RegisterProcess(msc, particle);
}

// single scattering
Expand All @@ -267,7 +347,6 @@ void CMSEmStandardPhysicsXS::ConstructProcess() {
br1->SetHighEnergyLimit(CLHEP::GeV);

// register processes
ph->RegisterProcess(msc, particle);
ph->RegisterProcess(eioni, particle);
ph->RegisterProcess(brem, particle);
ph->RegisterProcess(ee, particle);
Expand Down

0 comments on commit 4208c8c

Please sign in to comment.