Geant4 11.3.0
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)
 
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)
 
void SetMasterThread (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:738

Referenced by G4BetheHeitler5DModel(), G4LivermoreGammaConversion5DModel::G4LivermoreGammaConversion5DModel(), and operator=().

◆ ~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:504
G4EmParameters::GetConversionType
G4int GetConversionType() const
Definition G4EmParameters.cc:1153
G4EmParameters::Verbose
G4int Verbose() const
Definition G4EmParameters.cc:1030
G4PairProductionRelModel::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition G4PairProductionRelModel.cc:151

Referenced by 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 // there is no probability to have a lepton with zero momentum
509 // X and Y components of momentum may be zero, in that case SinPhi=1, cosPhi=0
510 const G4double EPlus = LeptonPlus.t();
511 const G4double PPlus = LeptonPlus.vect().mag();
512 const G4double pPX = LeptonPlus.x();
513 const G4double pPY = LeptonPlus.y();
514 const G4double pPZ = LeptonPlus.z();
515 G4double sinPhiPlus = 1.0;
516 G4double cosPhiPlus = 0.0;
517 G4double sinThetaPlus = 0.0;
518 G4double cosThetaPlus = pPZ/PPlus;
519 if (cosThetaPlus < 1.0 && cosThetaPlus > -1.0) {
520 sinThetaPlus = std::sqrt((1.0 - cosThetaPlus)*(1.0 + cosThetaPlus));
521 sinPhiPlus = pPY/(PPlus*sinThetaPlus);
522 cosPhiPlus = pPX/(PPlus*sinThetaPlus);
523 }
524
525 // denominators:
526 // the two cancelling leading terms for forward emission at high energy, removed
527 const G4double elMassCTP = LeptonMass*cosThetaPlus;
528 const G4double ePlusSTP = EPlus*sinThetaPlus;
529 const G4double DPlus = (elMassCTP*elMassCTP + ePlusSTP*ePlusSTP)
530 /(EPlus + PPlus*cosThetaPlus);
531
532 // there is no probability to have a lepton with zero momentum
533 // X and Y components of momentum may be zero, in that case SinPhi=0, cosPhi=1
534 const G4double EMinus = LeptonMinus.t();
535 const G4double PMinus = LeptonMinus.vect().mag();
536 const G4double ePX = LeptonMinus.x();
537 const G4double ePY = LeptonMinus.y();
538 const G4double ePZ = LeptonMinus.z();
539 G4double sinPhiMinus = 0.0;
540 G4double cosPhiMinus = 1.0;
541 G4double sinThetaMinus = 0.0;
542 G4double cosThetaMinus = ePZ/PMinus;
543 if (cosThetaMinus < 1.0 && cosThetaMinus > -1.0) {
544 sinThetaMinus = std::sqrt((1.0 - cosThetaMinus)*(1.0 + cosThetaMinus));
545 sinPhiMinus = ePY/(PMinus*sinThetaMinus);
546 cosPhiMinus = ePX/(PMinus*sinThetaMinus);
547 }
548
549 const G4double elMassCTM = LeptonMass*cosThetaMinus;
550 const G4double eMinSTM = EMinus*sinThetaMinus;
551 const G4double DMinus = (elMassCTM*elMassCTM + eMinSTM*eMinSTM)
552 /(EMinus + PMinus*cosThetaMinus);
553
554 // cos(phiMinus-PhiPlus)
555 const G4double cosdPhi = cosPhiPlus*cosPhiMinus + sinPhiPlus*sinPhiMinus;
556 const G4double PRec = Recoil.vect().mag();
557 const G4double q2 = PRec*PRec;
558
559 const G4double BigPhi = -LeptonMass2 / (GammaEnergy*GammaEnergy2 * q2*q2);
560
561 G4double FormFactor = 1.;
562 if (!iraw) {
563 if (itriplet) {
564 const G4double qun = factor1*iZ13*iZ13;
565 const G4double nun = qun * PRec;
566 if (nun < 1.) {
567 FormFactor = (nun < 0.01) ? (13.8-55.4*std::sqrt(nun))*nun
568 : std::sqrt(1-(nun-1)*(nun-1));
569 } // else FormFactor = 1 by default
570 } else {
571 const G4double dum3 = 217.*PRec*iZ13;
572 const G4double AFF = 1./(1. + dum3*dum3);
573 FormFactor = (1.-AFF)*(1-AFF);
574 }
575 } // else FormFactor = 1 by default
576
577 if (GammaPolarizationMag==0.) {
578 const G4double pPlusSTP = PPlus*sinThetaPlus;
579 const G4double pMinusSTM = PMinus*sinThetaMinus;
580 const G4double pPlusSTPperDP = pPlusSTP/DPlus;
581 const G4double pMinusSTMperDM = pMinusSTM/DMinus;
582 const G4double dunpol = BigPhi*(
583 pPlusSTPperDP *pPlusSTPperDP *(4.*EMinus*EMinus-q2)
584 + pMinusSTMperDM*pMinusSTMperDM*(4.*EPlus*EPlus - q2)
585 + 2.*pPlusSTPperDP*pMinusSTMperDM*cosdPhi
586 *(4.*EPlus*EMinus + q2 - 2.*GammaEnergy2)
587 - 2.*GammaEnergy2*(pPlusSTP*pPlusSTP+pMinusSTM*pMinusSTM)/(DMinus*DPlus));
588 betheheitler = dunpol * factor;
589 } else {
590 const G4double pPlusSTP = PPlus*sinThetaPlus;
591 const G4double pMinusSTM = PMinus*sinThetaMinus;
592 const G4double pPlusSTPCPPperDP = pPlusSTP*cosPhiPlus/DPlus;
593 const G4double pMinusSTMCPMperDM = pMinusSTM*cosPhiMinus/DMinus;
594 const G4double caa = 2.*(EPlus*pMinusSTMCPMperDM+EMinus*pPlusSTPCPPperDP);
595 const G4double cbb = pMinusSTMCPMperDM-pPlusSTPCPPperDP;
596 const G4double ccc = (pPlusSTP*pPlusSTP + pMinusSTM*pMinusSTM
597 +2.*pPlusSTP*pMinusSTM*cosdPhi)/ (DMinus*DPlus);
598 const G4double dtot= 2.*BigPhi*( caa*caa - q2*cbb*cbb - GammaEnergy2*ccc);
599 betheheitler = dtot * factor;
600 }
601 //
602 const G4double cross = Norme * Jacob0 * Jacob1 * Jacob2 * betheheitler
603 * FormFactor * RecoilMass / sqrts;
604 pdf = cross * (xu1 - xl1) / G4Exp(correctionIndex*G4Log(X1)); // cond1;
605 } while ( pdf < ymax * rndmv6[5] );
606 // END of Sampling
607
608 if ( fVerbose > 2 ) {
609 G4double recul = std::sqrt(Recoil.x()*Recoil.x()+Recoil.y()*Recoil.y()
610 +Recoil.z()*Recoil.z());
611 G4cout << "BetheHeitler5DModel GammaEnergy= " << GammaEnergy
612 << " PDF= " << pdf << " ymax= " << ymax
613 << " recul= " << recul << G4endl;
614 }
615
616 // back to Geant4 system
617
618 if ( fVerbose > 4 ) {
619 G4cout << "BetheHeitler5DModel GammaDirection " << GammaDirection << G4endl;
620 G4cout << "BetheHeitler5DModel GammaPolarization " << GammaPolarization << G4endl;
621 G4cout << "BetheHeitler5DModel GammaEnergy " << GammaEnergy << G4endl;
622 G4cout << "BetheHeitler5DModel Conv "
623 << (itriplet ? "triplet" : "nucl") << G4endl;
624 }
625
626 if (GammaPolarizationMag == 0.0) {
627 // set polarization axis orthohonal to direction
628 GammaPolarization = GammaDirection.orthogonal().unit();
629 } else {
630 // GammaPolarization not a unit vector
631 GammaPolarization /= GammaPolarizationMag;
632 }
633
634 // The unit norm vector that is orthogonal to the two others
635 G4ThreeVector yGrec = GammaDirection.cross(GammaPolarization);
636
637 // rotation from gamma ref. sys. to World
638 G4RotationMatrix GtoW(GammaPolarization,yGrec,GammaDirection);
639
640 Recoil.transform(GtoW);
641 LeptonPlus.transform(GtoW);
642 LeptonMinus.transform(GtoW);
643
644 if ( fVerbose > 2 ) {
645 G4cout << "BetheHeitler5DModel Recoil " << Recoil.x() << " " << Recoil.y() << " " << Recoil.z()
646 << " " << Recoil.t() << " " << G4endl;
647 G4cout << "BetheHeitler5DModel LeptonPlus " << LeptonPlus.x() << " " << LeptonPlus.y() << " "
648 << LeptonPlus.z() << " " << LeptonPlus.t() << " " << G4endl;
649 G4cout << "BetheHeitler5DModel LeptonMinus " << LeptonMinus.x() << " " << LeptonMinus.y() << " "
650 << LeptonMinus.z() << " " << LeptonMinus.t() << " " << G4endl;
651 }
652
653 // Create secondaries
654 auto aParticle1 = new G4DynamicParticle(fLepton1,LeptonMinus);
655 auto aParticle2 = new G4DynamicParticle(fLepton2,LeptonPlus);
656
657 // create G4DynamicParticle object for the particle3 ( recoil )
658 G4ParticleDefinition* RecoilPart;
659 if (itriplet) {
660 // triplet
661 RecoilPart = fTheElectron;
662 } else{
663 RecoilPart = theIonTable->GetIon(Z, A, 0);
664 }
665 auto aParticle3 = new G4DynamicParticle(RecoilPart,Recoil);
666
667 // Fill output vector
668 fvect->push_back(aParticle1);
669 fvect->push_back(aParticle2);
670 fvect->push_back(aParticle3);
671
672 // kill incident photon
673 fParticleChange->SetProposedKineticEnergy(0.);
674 fParticleChange->ProposeTrackStatus(fStopAndKill);
675}
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
G4ParticleMomentum
G4ThreeVector G4ParticleMomentum
Definition G4ParticleMomentum.hh:39
G4RotationMatrix
CLHEP::HepRotation G4RotationMatrix
Definition G4RotationMatrix.hh:37
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
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::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
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
G4VEmModel::SelectIsotopeNumber
G4int SelectIsotopeNumber(const G4Element *) const
Definition G4VEmModel.cc:274
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
G4INCL::Math::twoPi
const G4double twoPi
Definition G4INCLGlobals.hh:69

◆ 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

◆ 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: