Skip to content

Commit

Permalink
Merge pull request #38879 from hahnjo/transportation-with-msc
Browse files Browse the repository at this point in the history 
Introduce option to use G4TransportationWithMsc
  • Loading branch information
@cmsbuild
cmsbuild authored Jul 31, 2022
2 parents 161e17c + 4208c8c commit 0dd62d1
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 0dd62d1

Please sign in to comment.