Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4BetheHeitler5DModel Class Reference

#include <G4BetheHeitler5DModel.hh>

+ Inheritance diagram for G4BetheHeitler5DModel:

Public Member Functions

 G4BetheHeitler5DModel (const G4ParticleDefinition *p=nullptr, const G4String &nam="BetheHeitler5D")
 
 ~G4BetheHeitler5DModel () override
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
void SampleSecondaries (std::vector< G4DynamicParticle * > *fvect, const G4MaterialCutsCouple *couple, const G4DynamicParticle *aDynamicGamma, G4double, G4double) override
 
void SetVerbose (G4int val)
 
void SetLeptonPair (const G4ParticleDefinition *p1, const G4ParticleDefinition *p2)
 
G4BetheHeitler5DModel & operator= (const G4BetheHeitler5DModel &right)=delete
 
 G4BetheHeitler5DModel (const G4BetheHeitler5DModel &)=delete
 
- Public Member Functions inherited from G4PairProductionRelModel
 G4PairProductionRelModel (const G4ParticleDefinition *p=nullptr, const G4String &nam="BetheHeitlerLPM")
 
 ~G4PairProductionRelModel () override
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel) override
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX) override
 
void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
 
void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double) override
 
void SetLPMflag (G4bool val)
 
G4bool LPMflag () const
 
G4PairProductionRelModel & operator= (const G4PairProductionRelModel &right)=delete
 
 G4PairProductionRelModel (const G4PairProductionRelModel &)=delete
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, const G4double &length, G4double &eloss)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * GetCurrentElement (const G4Material *mat=nullptr) const
 
G4int SelectRandomAtomNumber (const G4Material *) const
 
const G4Isotope * GetCurrentIsotope (const G4Element *elm=nullptr) const
 
G4int SelectIsotopeNumber (const G4Element *) const
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
G4VEmModel * GetTripletModel ()
 
void SetTripletModel (G4VEmModel *)
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy) const
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
void SetUseBaseMaterials (G4bool val)
 
G4bool UseBaseMaterials () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 
void SetLPMFlag (G4bool)
 
G4VEmModel & operator= (const G4VEmModel &right)=delete
 
 G4VEmModel (const G4VEmModel &)=delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4PairProductionRelModel
void ComputePhi12 (const G4double delta, G4double &phi1, G4double &phi2)
 
G4double ScreenFunction1 (const G4double delta)
 
G4double ScreenFunction2 (const G4double delta)
 
void ScreenFunction12 (const G4double delta, G4double &f1, G4double &f2)
 
G4double ComputeParametrizedXSectionPerAtom (G4double gammaEnergy, G4double Z)
 
G4double ComputeXSectionPerAtom (G4double gammaEnergy, G4double Z)
 
G4double ComputeDXSectionPerAtom (G4double eplusEnergy, G4double gammaEnergy, G4double Z)
 
G4double ComputeRelDXSectionPerAtom (G4double eplusEnergy, G4double gammaEnergy, G4double Z)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Protected Attributes inherited from G4PairProductionRelModel
G4bool isFirstInstance {false}
 
G4bool fIsUseLPMCorrection
 
G4bool fIsUseCompleteScreening
 
G4double fLPMEnergy
 
G4double fParametrizedXSectionThreshold
 
G4double fCoulombCorrectionThreshold
 
G4Pow * fG4Calc
 
G4ParticleDefinition * fTheGamma
 
G4ParticleDefinition * fTheElectron
 
G4ParticleDefinition * fThePositron
 
G4ParticleChangeForGamma * fParticleChange
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData = nullptr
 
G4VParticleChange * pParticleChange = nullptr
 
G4PhysicsTable * xSectionTable = nullptr
 
const G4Material * pBaseMaterial = nullptr
 
const std::vector< G4double > * theDensityFactor = nullptr
 
const std::vector< G4int > * theDensityIdx = nullptr
 
G4double inveplus
 
G4double pFactor = 1.0
 
std::size_t currentCoupleIndex = 0
 
std::size_t basedCoupleIndex = 0
 
G4bool lossFlucFlag = true
 
- Static Protected Attributes inherited from G4PairProductionRelModel
static const G4int gMaxZet = 120
 
static const G4double gLPMconstant
 
static const G4double gXGL [8]
 
static const G4double gWGL [8]
 
static const G4double gFelLowZet [8]
 
static const G4double gFinelLowZet [8]
 
static const G4double gXSecFactor
 
static const G4double gEgLPMActivation = 100.*CLHEP::GeV
 
static std::vector< ElementData * > gElementData
 
static LPMFuncs gLPMFuncs
 

Detailed Description

Definition at line 64 of file G4BetheHeitler5DModel.hh.

Constructor & Destructor Documentation

◆ G4BetheHeitler5DModel() [1/2]

G4BetheHeitler5DModel::G4BetheHeitler5DModel ( const G4ParticleDefinition * p = nullptr,
const G4String & nam = "BetheHeitler5D" )
explicit

Definition at line 129 of file G4BetheHeitler5DModel.cc.

131 : G4PairProductionRelModel(pd, nam),
132 fLepton1(G4Electron::Definition()),fLepton2(G4Positron::Definition()),
133 fTheMuPlus(nullptr),fTheMuMinus(nullptr),
134 fVerbose(1),
135 fConversionType(0),
136 fConvMode(kEPair),
137 iraw(false)
138{
139 theIonTable = G4IonTable::GetIonTable();
140 //Q: Do we need this on Model
141 SetLowEnergyLimit(2*fTheElectron->GetPDGMass());
142}
kEPair
const G4int kEPair
Definition G4BetheHeitler5DModel.cc:123
G4Electron::Definition
static G4Electron * Definition()
Definition G4Electron.cc:45
G4IonTable::GetIonTable
static G4IonTable * GetIonTable()
Definition G4IonTable.cc:155
G4PairProductionRelModel::G4PairProductionRelModel
G4PairProductionRelModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="BetheHeitlerLPM")
Definition G4PairProductionRelModel.cc:120
G4PairProductionRelModel::fTheElectron
G4ParticleDefinition * fTheElectron
Definition G4PairProductionRelModel.hh:180
G4Positron::Definition
static G4Positron * Definition()
Definition G4Positron.cc:45
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition G4VEmModel.hh:745

◆ ~G4BetheHeitler5DModel()

G4BetheHeitler5DModel::~G4BetheHeitler5DModel ( )
overridedefault

◆ G4BetheHeitler5DModel() [2/2]

G4BetheHeitler5DModel::G4BetheHeitler5DModel ( const G4BetheHeitler5DModel & )
delete

Member Function Documentation

◆ Initialise()

void G4BetheHeitler5DModel::Initialise ( const G4ParticleDefinition * part,
const G4DataVector & vec )
overridevirtual

Implements G4VEmModel.

Reimplemented in G4LivermoreGammaConversion5DModel.

Definition at line 150 of file G4BetheHeitler5DModel.cc.

152{
153 G4PairProductionRelModel::Initialise(part, vec);
154
155 G4EmParameters* theManager = G4EmParameters::Instance();
156 // place to initialise model parameters
157 // Verbosity levels: ( Can redefine as needed, but some consideration )
158 // 0 = nothing
159 // > 2 print results
160 // > 3 print rejection warning from transformation (fix bug from gammaray .. )
161 // > 4 print photon direction & polarisation
162 fVerbose = theManager->Verbose();
163 fConversionType = theManager->GetConversionType();
164 //////////////////////////////////////////////////////////////
165 // iraw :
166 // true : isolated electron or nucleus.
167 // false : inside atom -> screening form factor
168 iraw = theManager->OnIsolated();
169 // G4cout << "BH5DModel::Initialise verbose " << fVerbose
170 // << " isolated " << iraw << " ctype "<< fConversionType << G4endl;
171
172 //Q: Do we need this on Model
173 // The Leptons defined via SetLeptonPair(..) method
174 SetLowEnergyLimit(2*CLHEP::electron_mass_c2);
175}
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition G4EmParameters.cc:70
G4EmParameters::OnIsolated
G4bool OnIsolated() const
Definition G4EmParameters.cc:493
G4EmParameters::GetConversionType
G4int GetConversionType() const
Definition G4EmParameters.cc:1098
G4EmParameters::Verbose
G4int Verbose() const
Definition G4EmParameters.cc:986
G4PairProductionRelModel::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition G4PairProductionRelModel.cc:151

Referenced by G4GammaConversionToMuons::BuildPhysicsTable(), and G4LivermoreGammaConversion5DModel::Initialise().

◆ operator=()

G4BetheHeitler5DModel & G4BetheHeitler5DModel::operator= ( const G4BetheHeitler5DModel & right)
delete

◆ SampleSecondaries()

void G4BetheHeitler5DModel::SampleSecondaries ( std::vector< G4DynamicParticle * > * fvect,
const G4MaterialCutsCouple * couple,
const G4DynamicParticle * aDynamicGamma,
G4double ,
G4double  )
overridevirtual

Implements G4VEmModel.

Definition at line 238 of file G4BetheHeitler5DModel.cc.

242{
243 // MeV
244 static const G4double ElectronMass = CLHEP::electron_mass_c2;
245
246 const G4double LeptonMass = fLepton1->GetPDGMass();
247 const G4double LeptonMass2 = LeptonMass*LeptonMass;
248
249 static const G4double alpha0 = CLHEP::fine_structure_const;
250 // mm
251 static const G4double r0 = CLHEP::classic_electr_radius;
252 // mbarn
253 static const G4double r02 = r0*r0*1.e+25;
254 static const G4double twoPi = CLHEP::twopi;
255 static const G4double factor = alpha0 * r02 / (twoPi*twoPi);
256 // static const G4double factor1 = pow((6.0 * pi),(1.0/3.0))/(8.*alpha0*ElectronMass);
257 static const G4double factor1 = 2.66134007899/(8.*alpha0*ElectronMass);
258 //
259 G4double PairInvMassMin = 2.*LeptonMass;
260 G4double TrThreshold = 2.0 * ( (LeptonMass2)/ElectronMass + LeptonMass);
261
262 //
263 static const G4double nu[2][10] = {
264 //electron
265 { 0.0227436, 0.0582046, 3.0322675, 2.8275065, -0.0034004,
266 1.1212766, 1.8989468, 68.3492750, 0.0211186, 14.4},
267 //muon
268 {0.67810E-06, 0.86037E+05, 2.0008395, 1.6739719, -0.0057279,
269 1.4222, 0.0, 263230.0, 0.0521, 51.1338}
270 };
271 static const G4double tr[2][10] = {
272 //electron
273 { 0.0332350, 4.3942537, 2.8515925, 2.6351695, -0.0031510,
274 1.5737305, 1.8104647, 20.6434021, -0.0272586, 28.9},
275 //muon
276 {0.10382E-03, 0.14408E+17, 4.1368679, 3.2662121, -0.0163091,
277 0.0000, 0.0, 0.0, 0.0000, 1.0000}
278 };
279 //
280 static const G4double para[2][3][2] = {
281 //electron
282 { {11., -16.},{-1.17, -2.95},{-2., -0.5} },
283 //muon
284 { {17.5, 1.},{-1.17, -2.95},{2., 6.} }
285 };
286 //
287 static const G4double correctionIndex = 1.4;
288 //
289 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
290 // Protection, Will not be true tot cross section = 0
291 if ( GammaEnergy <= PairInvMassMin) { return; }
292
293 const G4double GammaEnergy2 = GammaEnergy*GammaEnergy;
294
295 //////////////////////////////////////////////////////////////
296 const G4ParticleMomentum GammaDirection =
297 aDynamicGamma->GetMomentumDirection();
298 G4ThreeVector GammaPolarization = aDynamicGamma->GetPolarization();
299
300 // The protection polarization perpendicular to the direction vector,
301 // as it done in G4LivermorePolarizedGammaConversionModel,
302 // assuming Direction is unitary vector
303 // (projection to plane) p_proj = p - (p o d)/(d o d) x d
304 if ( GammaPolarization.howOrthogonal(GammaDirection) != 0) {
305 GammaPolarization -= GammaPolarization.dot(GammaDirection) * GammaDirection;
306 }
307 // End of Protection
308 //
309 const G4double GammaPolarizationMag = GammaPolarization.mag();
310
311 //////////////////////////////////////////////////////////////
312 // target element
313 // select randomly one element constituting the material
314 const G4Element* anElement = SelectTargetAtom(couple, fTheGamma, GammaEnergy,
315 aDynamicGamma->GetLogKineticEnergy() );
316 // Atomic number
317 const G4int Z = anElement->GetZasInt();
318 const G4int A = SelectIsotopeNumber(anElement);
319 const G4double iZ13 = 1./anElement->GetIonisation()->GetZ3();
320 const G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z);
321
322 const G4double NuThreshold = 2.0 * ( (LeptonMass2)/targetMass + LeptonMass);
323 // No conversion possible below nuclear threshold
324 if ( GammaEnergy <= NuThreshold) { return; }
325
326 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
327
328 // itriplet : true -- triplet, false -- nuclear.
329 G4bool itriplet = false;
330 if (fConversionType == 1) {
331 itriplet = false;
332 } else if (fConversionType == 2) {
333 itriplet = true;
334 if ( GammaEnergy <= TrThreshold ) return;
335 } else if ( GammaEnergy > TrThreshold ) {
336 // choose triplet or nuclear from a triplet/nuclear=1/Z
337 // total cross section ratio.
338 // approximate at low energies !
339 if(rndmEngine->flat()*(Z+1) < 1.) {
340 itriplet = true;
341 }
342 }
343
344 //
345 const G4double RecoilMass = itriplet ? ElectronMass : targetMass;
346 const G4double RecoilMass2 = RecoilMass*RecoilMass;
347 const G4double sCMS = 2.*RecoilMass*GammaEnergy + RecoilMass2;
348 const G4double sCMSPlusRM2 = sCMS + RecoilMass2;
349 const G4double sqrts = std::sqrt(sCMS);
350 const G4double isqrts2 = 1./(2.*sqrts);
351 //
352 const G4double PairInvMassMax = sqrts-RecoilMass;
353 const G4double PairInvMassRange = PairInvMassMax/PairInvMassMin;
354 const G4double lnPairInvMassRange = G4Log(PairInvMassRange);
355
356 // initial state. Defines z axis of "0" frame as along photon propagation.
357 // Since CMS(0., 0., GammaEnergy, GammaEnergy+RecoilMass) set some constants
358 const G4double betaCMS = G4LorentzVector(0.0,0.0,GammaEnergy,GammaEnergy+RecoilMass).beta();
359
360 // maximum value of pdf
361 const G4double EffectiveZ = iraw ? 0.5 : Z;
362 const G4double Threshold = itriplet ? TrThreshold : NuThreshold;
363 const G4double AvailableEnergy = GammaEnergy - Threshold;
364 const G4double LogAvailableEnergy = G4Log(AvailableEnergy);
365 //
366 const G4double MaxDiffCross = itriplet
367 ? MaxDiffCrossSection(tr[fConvMode],
368 EffectiveZ, AvailableEnergy, LogAvailableEnergy)
369 : MaxDiffCrossSection(nu[fConvMode],
370 EffectiveZ, AvailableEnergy, LogAvailableEnergy);
371 //
372 // 50% safety marging factor
373 const G4double ymax = 1.5 * MaxDiffCross;
374 // x1 bounds
375 const G4double xu1 = (LogAvailableEnergy > para[fConvMode][2][0])
376 ? para[fConvMode][0][0] +
377 para[fConvMode][1][0]*LogAvailableEnergy
378 : para[fConvMode][0][0] +
379 para[fConvMode][2][0]*para[fConvMode][1][0];
380 const G4double xl1 = (LogAvailableEnergy > para[fConvMode][2][1])
381 ? para[fConvMode][0][1] +
382 para[fConvMode][1][1]*LogAvailableEnergy
383 : para[fConvMode][0][1] +
384 para[fConvMode][2][1]*para[fConvMode][1][1];
385 //
386 G4LorentzVector Recoil;
387 G4LorentzVector LeptonPlus;
388 G4LorentzVector LeptonMinus;
389 G4double pdf = 0.;
390
391 G4double rndmv6[6] = {0.0};
392 const G4double corrFac = 1.0/(correctionIndex + 1.0);
393 const G4double expLowLim = -20.;
394 const G4double logLowLim = G4Exp(expLowLim/corrFac);
395 G4double z0, z1, z2, x0, x1;
396 G4double betheheitler, sinTheta, cosTheta, dum0;
397 // START Sampling
398 do {
399
400 rndmEngine->flatArray(6, rndmv6);
401
402 //////////////////////////////////////////////////
403 // pdf pow(x,c) with c = 1.4
404 // integral y = pow(x,(c+1))/(c+1) @ x = 1 => y = 1 /(1+c)
405 // invCdf exp( log(y /* *( c + 1.0 )/ (c + 1.0 ) */ ) /( c + 1.0) )
406 //////////////////////////////////////////////////
407
408 z0 = (rndmv6[0] > logLowLim) ? G4Log(rndmv6[0])*corrFac : expLowLim;
409 G4double X1 = (z0 > expLowLim) ? G4Exp(z0) : 0.0;
410 z1 = xl1 + (xu1 - xl1)*rndmv6[1];
411 if (z1 > expLowLim) {
412 x0 = G4Exp(z1);
413 dum0 = 1.0/(1.0 + x0);
414 x1 = dum0*x0;
415 cosTheta = -1.0 + 2.0*x1;
416 sinTheta = 2*std::sqrt(x1*(1.0 - x1));
417 } else {
418 x0 = 0.0;
419 dum0 = 1.0;
420 cosTheta = -1.0;
421 sinTheta = 0.0;
422 }
423
424 z2 = X1*X1*lnPairInvMassRange;
425 const G4double PairInvMass = PairInvMassMin*((z2 > 1.e-3) ? G4Exp(z2) : 1 + z2 + 0.5*z2*z2);
426
427 // cos and sin theta-lepton
428 const G4double cosThetaLept = std::cos(pi*rndmv6[2]);
429 // sin(ThetaLept) is always in [0,+1] if ThetaLept is in [0,pi]
430 const G4double sinThetaLept = std::sqrt((1.-cosThetaLept)*(1.+cosThetaLept));
431 // cos and sin phi-lepton
432 const G4double cosPhiLept = std::cos(twoPi*rndmv6[3]-pi);
433 // sin(PhiLept) is in [-1,0] if PhiLept in [-pi,0) and
434 // is in [0,+1] if PhiLept in [0,+pi]
435 const G4double sinPhiLept = std::copysign(std::sqrt((1.-cosPhiLept)*(1.+cosPhiLept)),rndmv6[3]-0.5);
436 // cos and sin phi
437 const G4double cosPhi = std::cos(twoPi*rndmv6[4]-pi);
438 const G4double sinPhi = std::copysign(std::sqrt((1.-cosPhi)*(1.+cosPhi)),rndmv6[4]-0.5);
439
440 //////////////////////////////////////////////////
441 // frames:
442 // 3 : the laboratory Lorentz frame, Geant4 axes definition
443 // 0 : the laboratory Lorentz frame, axes along photon direction and polarisation
444 // 1 : the center-of-mass Lorentz frame
445 // 2 : the pair Lorentz frame
446 //////////////////////////////////////////////////
447
448 // in the center-of-mass frame
449
450 const G4double RecEnergyCMS = (sCMSPlusRM2-PairInvMass*PairInvMass)*isqrts2;
451 const G4double LeptonEnergy2 = PairInvMass*0.5;
452
453 // New way of calucaltion thePRecoil to avoid underflow
454 G4double abp = std::max((2.0*GammaEnergy*RecoilMass -
455 PairInvMass*PairInvMass + 2.0*PairInvMass*RecoilMass)*
456 (2.0*GammaEnergy*RecoilMass -
457 PairInvMass*PairInvMass - 2.0*PairInvMass*RecoilMass),0.0);
458
459 G4double thePRecoil = std::sqrt(abp) * isqrts2;
460
461 // back to the center-of-mass frame
462 Recoil.set( thePRecoil*sinTheta*cosPhi,
463 thePRecoil*sinTheta*sinPhi,
464 thePRecoil*cosTheta,
465 RecEnergyCMS);
466
467 // in the pair frame
468 const G4double thePLepton = std::sqrt( (LeptonEnergy2-LeptonMass)
469 *(LeptonEnergy2+LeptonMass));
470
471 LeptonPlus.set(thePLepton*sinThetaLept*cosPhiLept,
472 thePLepton*sinThetaLept*sinPhiLept,
473 thePLepton*cosThetaLept,
474 LeptonEnergy2);
475
476 LeptonMinus.set(-LeptonPlus.x(),
477 -LeptonPlus.y(),
478 -LeptonPlus.z(),
479 LeptonEnergy2);
480
481
482 // Normalisation of final state phase space:
483 // Section 47 of Particle Data Group, Chin. Phys. C, 40, 100001 (2016)
484 // const G4double Norme = Recoil1.vect().mag() * LeptonPlus2.vect().mag();
485 const G4double Norme = Recoil.vect().mag() * LeptonPlus.vect().mag();
486
487 // e+, e- to CMS frame from pair frame
488
489 // boost vector from Pair to CMS
490 const G4ThreeVector pair2cms =
491 G4LorentzVector( -Recoil.x(), -Recoil.y(), -Recoil.z(),
492 sqrts-RecEnergyCMS).boostVector();
493
494 LeptonPlus.boost(pair2cms);
495 LeptonMinus.boost(pair2cms);
496
497 // back to the laboratory frame (make use of the CMS(0,0,Eg,Eg+RM)) form
498
499 Recoil.boostZ(betaCMS);
500 LeptonPlus.boostZ(betaCMS);
501 LeptonMinus.boostZ(betaCMS);
502
503 // Jacobian factors
504 const G4double Jacob0 = x0*dum0*dum0;
505 const G4double Jacob1 = 2.*X1*lnPairInvMassRange*PairInvMass;
506 const G4double Jacob2 = std::abs(sinThetaLept);
507
508 const G4double EPlus = LeptonPlus.t();
509 const G4double PPlus = LeptonPlus.vect().mag();
510 const G4double sinThetaPlus = LeptonPlus.vect().perp()/PPlus;
511 const G4double cosThetaPlus = LeptonPlus.vect().cosTheta();
512
513 const G4double pPX = LeptonPlus.x();
514 const G4double pPY = LeptonPlus.y();
515 const G4double dum1 = 1./std::sqrt( pPX*pPX + pPY*pPY );
516 const G4double cosPhiPlus = pPX*dum1;
517 const G4double sinPhiPlus = pPY*dum1;
518
519 // denominators:
520 // the two cancelling leading terms for forward emission at high energy, removed
521 const G4double elMassCTP = LeptonMass*cosThetaPlus;
522 const G4double ePlusSTP = EPlus*sinThetaPlus;
523 const G4double DPlus = (elMassCTP*elMassCTP + ePlusSTP*ePlusSTP)
524 /(EPlus + PPlus*cosThetaPlus);
525
526 const G4double EMinus = LeptonMinus.t();
527 const G4double PMinus = LeptonMinus.vect().mag();
528 const G4double sinThetaMinus = LeptonMinus.vect().perp()/PMinus;
529 const G4double cosThetaMinus = LeptonMinus.vect().cosTheta();
530
531 const G4double ePX = LeptonMinus.x();
532 const G4double ePY = LeptonMinus.y();
533 const G4double dum2 = 1./std::sqrt( ePX*ePX + ePY*ePY );
534 const G4double cosPhiMinus = ePX*dum2;
535 const G4double sinPhiMinus = ePY*dum2;
536
537 const G4double elMassCTM = LeptonMass*cosThetaMinus;
538 const G4double eMinSTM = EMinus*sinThetaMinus;
539 const G4double DMinus = (elMassCTM*elMassCTM + eMinSTM*eMinSTM)
540 /(EMinus + PMinus*cosThetaMinus);
541
542 // cos(phiMinus-PhiPlus)
543 const G4double cosdPhi = cosPhiPlus*cosPhiMinus + sinPhiPlus*sinPhiMinus;
544 const G4double PRec = Recoil.vect().mag();
545 const G4double q2 = PRec*PRec;
546
547 const G4double BigPhi = -LeptonMass2 / (GammaEnergy*GammaEnergy2 * q2*q2);
548
549 G4double FormFactor = 1.;
550 if (!iraw) {
551 if (itriplet) {
552 const G4double qun = factor1*iZ13*iZ13;
553 const G4double nun = qun * PRec;
554 if (nun < 1.) {
555 FormFactor = (nun < 0.01) ? (13.8-55.4*std::sqrt(nun))*nun
556 : std::sqrt(1-(nun-1)*(nun-1));
557 } // else FormFactor = 1 by default
558 } else {
559 const G4double dum3 = 217.*PRec*iZ13;
560 const G4double AFF = 1./(1. + dum3*dum3);
561 FormFactor = (1.-AFF)*(1-AFF);
562 }
563 } // else FormFactor = 1 by default
564
565 if (GammaPolarizationMag==0.) {
566 const G4double pPlusSTP = PPlus*sinThetaPlus;
567 const G4double pMinusSTM = PMinus*sinThetaMinus;
568 const G4double pPlusSTPperDP = pPlusSTP/DPlus;
569 const G4double pMinusSTMperDM = pMinusSTM/DMinus;
570 const G4double dunpol = BigPhi*(
571 pPlusSTPperDP *pPlusSTPperDP *(4.*EMinus*EMinus-q2)
572 + pMinusSTMperDM*pMinusSTMperDM*(4.*EPlus*EPlus - q2)
573 + 2.*pPlusSTPperDP*pMinusSTMperDM*cosdPhi
574 *(4.*EPlus*EMinus + q2 - 2.*GammaEnergy2)
575 - 2.*GammaEnergy2*(pPlusSTP*pPlusSTP+pMinusSTM*pMinusSTM)/(DMinus*DPlus));
576 betheheitler = dunpol * factor;
577 } else {
578 const G4double pPlusSTP = PPlus*sinThetaPlus;
579 const G4double pMinusSTM = PMinus*sinThetaMinus;
580 const G4double pPlusSTPCPPperDP = pPlusSTP*cosPhiPlus/DPlus;
581 const G4double pMinusSTMCPMperDM = pMinusSTM*cosPhiMinus/DMinus;
582 const G4double caa = 2.*(EPlus*pMinusSTMCPMperDM+EMinus*pPlusSTPCPPperDP);
583 const G4double cbb = pMinusSTMCPMperDM-pPlusSTPCPPperDP;
584 const G4double ccc = (pPlusSTP*pPlusSTP + pMinusSTM*pMinusSTM
585 +2.*pPlusSTP*pMinusSTM*cosdPhi)/ (DMinus*DPlus);
586 const G4double dtot= 2.*BigPhi*( caa*caa - q2*cbb*cbb - GammaEnergy2*ccc);
587 betheheitler = dtot * factor;
588 }
589 //
590 const G4double cross = Norme * Jacob0 * Jacob1 * Jacob2 * betheheitler
591 * FormFactor * RecoilMass / sqrts;
592 pdf = cross * (xu1 - xl1) / G4Exp(correctionIndex*G4Log(X1)); // cond1;
593 } while ( pdf < ymax * rndmv6[5] );
594 // END of Sampling
595
596 if ( fVerbose > 2 ) {
597 G4double recul = std::sqrt(Recoil.x()*Recoil.x()+Recoil.y()*Recoil.y()
598 +Recoil.z()*Recoil.z());
599 G4cout << "BetheHeitler5DModel GammaEnergy= " << GammaEnergy
600 << " PDF= " << pdf << " ymax= " << ymax
601 << " recul= " << recul << G4endl;
602 }
603
604 // back to Geant4 system
605
606 if ( fVerbose > 4 ) {
607 G4cout << "BetheHeitler5DModel GammaDirection " << GammaDirection << G4endl;
608 G4cout << "BetheHeitler5DModel GammaPolarization " << GammaPolarization << G4endl;
609 G4cout << "BetheHeitler5DModel GammaEnergy " << GammaEnergy << G4endl;
610 G4cout << "BetheHeitler5DModel Conv "
611 << (itriplet ? "triplet" : "nucl") << G4endl;
612 }
613
614 if (GammaPolarizationMag == 0.0) {
615 // set polarization axis orthohonal to direction
616 GammaPolarization = GammaDirection.orthogonal().unit();
617 } else {
618 // GammaPolarization not a unit vector
619 GammaPolarization /= GammaPolarizationMag;
620 }
621
622 // The unit norm vector that is orthogonal to the two others
623 G4ThreeVector yGrec = GammaDirection.cross(GammaPolarization);
624
625 // rotation from gamma ref. sys. to World
626 G4RotationMatrix GtoW(GammaPolarization,yGrec,GammaDirection);
627
628 Recoil.transform(GtoW);
629 LeptonPlus.transform(GtoW);
630 LeptonMinus.transform(GtoW);
631
632 if ( fVerbose > 2 ) {
633 G4cout << "BetheHeitler5DModel Recoil " << Recoil.x() << " " << Recoil.y() << " " << Recoil.z()
634 << " " << Recoil.t() << " " << G4endl;
635 G4cout << "BetheHeitler5DModel LeptonPlus " << LeptonPlus.x() << " " << LeptonPlus.y() << " "
636 << LeptonPlus.z() << " " << LeptonPlus.t() << " " << G4endl;
637 G4cout << "BetheHeitler5DModel LeptonMinus " << LeptonMinus.x() << " " << LeptonMinus.y() << " "
638 << LeptonMinus.z() << " " << LeptonMinus.t() << " " << G4endl;
639 }
640
641 // Create secondaries
642 auto aParticle1 = new G4DynamicParticle(fLepton1,LeptonMinus);
643 auto aParticle2 = new G4DynamicParticle(fLepton2,LeptonPlus);
644
645 // create G4DynamicParticle object for the particle3 ( recoil )
646 G4ParticleDefinition* RecoilPart;
647 if (itriplet) {
648 // triplet
649 RecoilPart = fTheElectron;
650 } else{
651 RecoilPart = theIonTable->GetIon(Z, A, 0);
652 }
653 auto aParticle3 = new G4DynamicParticle(RecoilPart,Recoil);
654
655 // Fill output vector
656 fvect->push_back(aParticle1);
657 fvect->push_back(aParticle2);
658 fvect->push_back(aParticle3);
659
660 // kill incident photon
661 fParticleChange->SetProposedKineticEnergy(0.);
662 fParticleChange->ProposeTrackStatus(fStopAndKill);
663}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:180
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:227
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition G4LorentzVector.hh:33
fStopAndKill
@ fStopAndKill
Definition G4TrackStatus.hh:46
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::orthogonal
Hep3Vector orthogonal() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::Hep3Vector::mag
double mag() const
CLHEP::Hep3Vector::howOrthogonal
double howOrthogonal(const Hep3Vector &v) const
Definition SpaceVector.cc:215
CLHEP::Hep3Vector::cosTheta
double cosTheta() const
CLHEP::Hep3Vector::perp
double perp() const
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition LorentzVector.cc:173
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition LorentzVector.cc:55
CLHEP::HepLorentzVector::boostZ
HepLorentzVector & boostZ(double beta)
Definition LorentzVector.cc:223
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
CLHEP::HepLorentzVector::set
void set(double x, double y, double z, double t)
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
CLHEP::HepRandomEngine::flat
virtual double flat()=0
CLHEP::HepRandomEngine::flatArray
virtual void flatArray(const int size, double *vect)=0
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetLogKineticEnergy
G4double GetLogKineticEnergy() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetPolarization
const G4ThreeVector & GetPolarization() const
G4Element::GetIonisation
G4IonisParamElm * GetIonisation() const
Definition G4Element.hh:171
G4Element::GetZasInt
G4int GetZasInt() const
Definition G4Element.hh:120
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition G4IonTable.cc:445
G4IonisParamElm::GetZ3
G4double GetZ3() const
Definition G4IonisParamElm.hh:54
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition G4NucleiProperties.cc:49
G4PairProductionRelModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition G4PairProductionRelModel.hh:182
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition G4ParticleChangeForGamma.hh:106
G4VEmModel::SelectIsotopeNumber
G4int SelectIsotopeNumber(const G4Element *) const
Definition G4VEmModel.cc:273
G4VEmModel::SelectTargetAtom
const G4Element * SelectTargetAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition G4VEmModel.hh:576
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4INCL::Math::twoPi
const G4double twoPi
Definition G4INCLGlobals.hh:69

Referenced by G4GammaConversionToMuons::PostStepDoIt().

◆ SetLeptonPair()

void G4BetheHeitler5DModel::SetLeptonPair ( const G4ParticleDefinition * p1,
const G4ParticleDefinition * p2 )

Definition at line 179 of file G4BetheHeitler5DModel.cc.

181{
182 G4int pdg1 = p1->GetPDGEncoding();
183 G4int pdg2 = p2->GetPDGEncoding();
184 G4int pdg = std::abs(pdg1);
185 if ( pdg1 != -pdg2 || (pdg != 11 && pdg != 13) ) {
186 G4ExceptionDescription ed;
187 ed << " Wrong pair of leptons: " << p1->GetParticleName()
188 << " and " << p1->GetParticleName();
189 G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0007",
190 FatalErrorInArgument, ed, "");
191 } else {
192 if ( pdg == 11 ) {
193 SetConversionMode(kEPair);
194 if( pdg1 == 11 ) {
195 fLepton1 = p1;
196 fLepton2 = p2;
197 } else {
198 fLepton1 = p2;
199 fLepton2 = p1;
200 }
201 if (fVerbose > 0)
202 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to e+ e-"
203 << G4endl;
204 } else {
205 SetConversionMode(kMuPair);
206 if( pdg1 == 13 ) {
207 fLepton1 = p1;
208 fLepton2 = p2;
209 } else {
210 fLepton1 = p2;
211 fLepton2 = p1;
212 }
213 fTheMuPlus = fLepton2;
214 fTheMuMinus= fLepton1;
215 if (fVerbose > 0)
216 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to mu+ mu-"
217 << G4endl;
218 }
219 }
220}
kMuPair
const G4int kMuPair
Definition G4BetheHeitler5DModel.cc:124
FatalErrorInArgument
@ FatalErrorInArgument
Definition G4ExceptionSeverity.hh:68
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition G4Exception.hh:40
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96

Referenced by G4GammaConversionToMuons::BuildPhysicsTable().

◆ SetVerbose()

void G4BetheHeitler5DModel::SetVerbose ( G4int val)
inline

Definition at line 81 of file G4BetheHeitler5DModel.hh.

81{ fVerbose = val; }
The documentation for this class was generated from the following files: