Skip to content

Commit

Permalink
An additional paramater to control if the e+ correction to MSC angle …
Browse files Browse the repository at this point in the history 
…should be applied.

A global flag has been added to control if the e+ correction to the theta_0 MSC
angle should be applied or not. It's ON (true) by default, i.e. there is e+
correction, but should be set to OFF (false) in the ATLAS simulation (as they
don't use it: responsible for moving away from their experimental data).
  • Loading branch information
@mnovak42
mnovak42 committed Jan 28, 2025
1 parent f3d834e commit afef3f7
Show file tree
Hide file tree
Showing 8 changed files with 36 additions and 15 deletions.
7 changes: 5 additions & 2 deletions G4HepEm/G4HepEm/src/G4HepEmConfig.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,17 +208,20 @@ void G4HepEmConfig::Dump() {
<< std::setw(5) << std::right
<< fG4HepEmParameters->fElectronTrackingCut/CLHEP::keV
<< " [keV] " << std::endl;
std::cout << std::left << std::setw(30) << " Min loss table energies " << " : "
std::cout << std::left << std::setw(30) << " Min loss table energy " << " : "
<< std::setw(5) << std::right
<< fG4HepEmParameters->fMinLossTableEnergy/CLHEP::keV
<< " [keV] " << std::endl;
std::cout << std::left << std::setw(30) << " Max loss table energies " << " : "
std::cout << std::left << std::setw(30) << " Max loss table energy " << " : "
<< std::setw(5) << std::right
<< fG4HepEmParameters->fMaxLossTableEnergy/CLHEP::TeV
<< " [TeV] " << std::endl;
std::cout << std::left << std::setw(30) << " Number of loss table bins " << " : "
<< std::setw(5) << std::right
<< fG4HepEmParameters->fNumLossTableBins << std::endl;
std::cout << std::left << std::setw(30) << " Positron corr. in MSC theta0 " << " : "
<< std::setw(5) << std::right
<< fG4HepEmParameters->fIsMSCPositronCor << std::endl;
std::cout << std::left << std::setw(30) << " Linear loss limit " << " : "
<< std::setw(5) << std::right
<< fG4HepEmParameters->fParametersPerRegion[0].fLinELossLimit*100
Expand Down
3 changes: 3 additions & 0 deletions G4HepEm/G4HepEmData/include/G4HepEmParameters.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,9 @@ struct G4HepEmParameters {
* in case of \f$e^-/e^+\f$ primary particles.*/
double fElectronBremModelLim = 1000; // 1 GeV

/** Flag to indicate if the e+ correction should be used in the MSC \theta_0 angle. */
bool fIsMSCPositronCor = true;

/** Number of detector regions */
int fNumRegions = 0;
/** A `G4HepEmRegionParmeters` array for the individual detector regions. */
Expand Down
2 changes: 2 additions & 0 deletions G4HepEm/G4HepEmDataJsonIO/src/G4HepEmDataJsonIOImpl.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,6 +179,7 @@ namespace nlohmann
j["fMaxLossTableEnergy"] = d->fMaxLossTableEnergy;
j["fNumLossTableBins"] = d->fNumLossTableBins;
j["fElectronBremModelLim"] = d->fElectronBremModelLim;
j["fIsMSCPositronCor"] = d->fIsMSCPositronCor;
j["fNumRegions"] = d->fNumRegions;
j["fParametersPerRegion"] =
make_span(d->fNumRegions, d->fParametersPerRegion);
Expand All @@ -200,6 +201,7 @@ namespace nlohmann
d->fMaxLossTableEnergy = j.at("fMaxLossTableEnergy").get<double>();
d->fNumLossTableBins = j.at("fNumLossTableBins").get<int>();
d->fElectronBremModelLim = j.at("fElectronBremModelLim").get<double>();
d->fIsMSCPositronCor = j.at("fIsMSCPositronCor").get<bool>();
d->fNumRegions = j.at("fNumRegions").get<int>();

d->fParametersPerRegion = new G4HepEmRegionParmeters[d->fNumRegions];
Expand Down
9 changes: 9 additions & 0 deletions G4HepEm/G4HepEmInit/src/G4HepEmParametersInit.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include "G4HepEmParameters.hh"

// g4 include
#include "G4Version.hh"
#include "G4EmParameters.hh"
#include "G4RegionStore.hh"
#include "G4MscStepLimitType.hh"
Expand All @@ -25,6 +26,14 @@ void InitHepEmParameters(struct G4HepEmParameters* hepEmPars) {
// energy limit between the 2 models (Seltzer-Berger and RelBrem) used for e-/e+
hepEmPars->fElectronBremModelLim = 1.0*CLHEP::GeV;

// flag to indicate if the e+ correction to the MSC theta0 angle should be used
// note: em parameter for this available only in Geant4 version >= 11.1
#if G4VERSION_NUMBER >= 1110
hepEmPars->fIsMSCPositronCor = G4EmParameters::Instance()->MscPositronCorrection();
#else
hepEmPars->fIsMSCPositronCor = true;
#endif

// get the number of detector regions and allocate the per-region data array
int numRegions = G4RegionStore::GetInstance()->size();
hepEmPars->fNumRegions = numRegions;
Expand Down
6 changes: 3 additions & 3 deletions G4HepEm/G4HepEmRun/include/G4HepEmElectronInteractionUMSC.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ public:
G4HepEmHostDevice
static void SampleScattering(G4HepEmData* hepEmData, G4HepEmMSCTrackData* mscData, double pStepLength,
double preStepEkin, double preStepTr1mfp, double postStepEkin, double postStepTr1mfp,
int imat, bool isElectron, G4HepEmRandomEngine* rnge);
int imat, bool isElectron, bool isPosCor, G4HepEmRandomEngine* rnge);



Expand All @@ -43,7 +43,7 @@ public:
static double SampleCosineTheta(double pStepLengt, double preStepEkin, double preStepTr1mfp,
double postStepEkin, double postStepTr1mfp, double umscTlimitMin,
double radLength, double zeff, const double* umscTailCoeff, const double* umscThetaCoeff,
bool isElectron, G4HepEmRandomEngine* rnge);
bool isElectron, bool isPosCor, G4HepEmRandomEngine* rnge);

// auxilary method for sampling cos(theta) in a simplified way: using an arbitrary pdf with correct mean and stdev
// (used in the above `SampleCosineTheta`)
Expand All @@ -53,7 +53,7 @@ public:
// auxilary method for computing theta0 (used in the above `SampleCosineTheta`)
G4HepEmHostDevice
static double ComputeTheta0(double stepInRadLength, double postStepEkin, double preStepEkin,
double zeff, const double* umscThetaCoeff, bool isElectron);
double zeff, const double* umscThetaCoeff, bool isElectron, bool isPosCor);

// auxilary method for computing the e+ correction to theta0 (used in the above `ComputeTheta0` but only in case of e+)
G4HepEmHostDevice
Expand Down
15 changes: 8 additions & 7 deletions G4HepEm/G4HepEmRun/include/G4HepEmElectronInteractionUMSC.icc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,10 +128,10 @@ void G4HepEmElectronInteractionUMSC::StepLimit(G4HepEmData* hepEmData, G4HepEmPa

void G4HepEmElectronInteractionUMSC::SampleScattering(G4HepEmData* hepEmData, G4HepEmMSCTrackData* mscData,
double pStepLength, double preStepEkin, double preStepTr1mfp, double postStepEkin, double postStepTr1mfp,
int imat, bool isElectron, G4HepEmRandomEngine* rnge) {
int imat, bool isElectron, bool isPosCor, G4HepEmRandomEngine* rnge) {
const struct G4HepEmMatData& matData = hepEmData->fTheMaterialData->fMaterialData[imat];
const double cost = SampleCosineTheta(pStepLength, preStepEkin, preStepTr1mfp, postStepEkin, postStepTr1mfp, mscData->fTlimitMin,
matData.fRadiationLength, matData.fZeff, matData.fUMSCTailCoeff, matData.fUMSCThetaCoeff, isElectron, rnge);
matData.fRadiationLength, matData.fZeff, matData.fUMSCTailCoeff, matData.fUMSCThetaCoeff, isElectron, isPosCor, rnge);
// no scattering so no dispacement in case cost = 1.0
if (std::abs(cost) >= 1.0) {
mscData->fIsNoScatteringInMSC = true;
Expand All @@ -150,7 +150,7 @@ void G4HepEmElectronInteractionUMSC::SampleScattering(G4HepEmData* hepEmData, G4

double G4HepEmElectronInteractionUMSC::SampleCosineTheta(double pStepLength, double preStepEkin, double preStepTr1mfp,
double postStepEkin, double postStepTr1mfp, double umscTlimitMin, double radLength, double zeff,
const double* umscTailCoeff, const double* umscThetaCoeff, bool isElectron, G4HepEmRandomEngine* rnge) {
const double* umscTailCoeff, const double* umscThetaCoeff, bool isElectron, bool isPosCor, G4HepEmRandomEngine* rnge) {
// NOTE: since I already know the finalEnergy in the electron Manager and that is already set to the track,
// I could get the logFinalEnergy from the updated track which could save up one log call
// compute the 1rst transport mfp at the post-step point energy
Expand Down Expand Up @@ -199,8 +199,8 @@ double G4HepEmElectronInteractionUMSC::SampleCosineTheta(double pStepLength, dou
const double tsmall = G4HepEmMin(umscTlimitMin, 1.0);
const bool stpNotExSmall = pStepLength > tsmall;
const double theta0 = stpNotExSmall
? ComputeTheta0(pStepLength/radLength, postStepEkin, preStepEkin, zeff, umscThetaCoeff, isElectron)
: ComputeTheta0(tsmall/radLength, postStepEkin, preStepEkin, zeff, umscThetaCoeff, isElectron)*std::sqrt(pStepLength/tsmall);
? ComputeTheta0(pStepLength/radLength, postStepEkin, preStepEkin, zeff, umscThetaCoeff, isElectron, isPosCor)
: ComputeTheta0(tsmall/radLength, postStepEkin, preStepEkin, zeff, umscThetaCoeff, isElectron, isPosCor)*std::sqrt(pStepLength/tsmall);
//
if (theta0 > kPi*0.166666) {
return SimpleScattering(xmeanth, x2meanth, rnge);
Expand Down Expand Up @@ -293,14 +293,15 @@ double G4HepEmElectronInteractionUMSC::SimpleScattering(double xmeanth, double x


// totry: all these could probably computed in `float`
double G4HepEmElectronInteractionUMSC::ComputeTheta0(double stepInRadLength, double postStepEkin, double preStepEkin, double zeff, const double* umscThetaCoeff, bool isElectron) {
double G4HepEmElectronInteractionUMSC::ComputeTheta0(double stepInRadLength, double postStepEkin, double preStepEkin, double zeff,
const double* umscThetaCoeff, bool isElectron, bool isPosCor) {
// ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
const double kHighland = 13.6; // note:: assumed to be in MeV
const double postInvBetaPc = (postStepEkin + kElectronMassC2)/(postStepEkin*(postStepEkin + 2.*kElectronMassC2));
const double invBetaPc = preStepEkin != postStepEkin
? std::sqrt(postInvBetaPc*(preStepEkin + kElectronMassC2)/(preStepEkin*(preStepEkin + 2.*kElectronMassC2)))
: postInvBetaPc;
const double y = isElectron
const double y = (isElectron || !isPosCor)
? stepInRadLength
: stepInRadLength * Theta0PositronCorrection(preStepEkin*postStepEkin, zeff);
return kHighland*std::sqrt(y)*invBetaPc*(umscThetaCoeff[0] + umscThetaCoeff[1]*G4HepEmLog(y));
Expand Down
3 changes: 2 additions & 1 deletion G4HepEm/G4HepEmRun/include/G4HepEmElectronManager.icc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -301,8 +301,9 @@ void G4HepEmElectronManager::SampleMSC(struct G4HepEmData* hepEmData, struct G4H
const double postStepTr1mfp = GetTransportMFP(elData, theImat, postStepEkin, postStepLEkin);
// - sample scattering: including net angular deflection and lateral dispacement that will be
// written into mscData::fDirection and mscData::fDisplacement
const bool isPosCor = hepEmPars->fIsMSCPositronCor;
G4HepEmElectronInteractionUMSC::SampleScattering(hepEmData, mscData, pStepLength, preStepEkin, mscData->fLambtr1, postStepEkin, postStepTr1mfp,
theImat, isElectron, rnge);
theImat, isElectron, isPosCor, rnge);
// NOTE: displacement will be applied in the caller where we have access to the required Geant4 functionality
// (and if its length is longer than a small minimal length and we are not ended up on boundary)
//
Expand Down
6 changes: 4 additions & 2 deletions testing/TestUtils/G4HepEmDataComparison.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,10 +72,12 @@ bool operator==(const G4HepEmParameters& lhs, const G4HepEmParameters& rhs)

return std::tie(lhs.fElectronTrackingCut, lhs.fMinLossTableEnergy,
lhs.fMaxLossTableEnergy, lhs.fNumLossTableBins,
lhs.fElectronBremModelLim, lhs.fNumRegions) ==
lhs.fElectronBremModelLim, lhs.fIsMSCPositronCor,
lhs.fNumRegions) ==
std::tie(rhs.fElectronTrackingCut, rhs.fMinLossTableEnergy,
rhs.fMaxLossTableEnergy, rhs.fNumLossTableBins,
rhs.fElectronBremModelLim, rhs.fNumRegions);
rhs.fElectronBremModelLim, rhs.fNumRegions,
rhs.fIsMSCPositronCor);
}

bool operator!=(const G4HepEmParameters& lhs, const G4HepEmParameters& rhs)
Expand Down

0 comments on commit afef3f7

Please sign in to comment.