Geant4 10.7.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")
 
virtual ~G4BetheHeitler5DModel ()
 
virtual 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)
 
- Public Member Functions inherited from G4PairProductionRelModel
 G4PairProductionRelModel (const G4ParticleDefinition *p=nullptr, const G4String &nam="BetheHeitlerLPM")
 
virtual ~G4PairProductionRelModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel) override
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX) override
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double) override
 
void SetLPMflag (G4bool val)
 
G4bool LPMflag () const
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
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 ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
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 *, G4double &eloss, G4double &niel, G4double length)
 
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 SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
virtual 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)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
G4int SelectIsotopeNumber (const G4Element *)
 
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 LPMFlag () 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 SetLPMFlag (G4bool val)
 
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 *)
 
const G4Element * GetCurrentElement () const
 
const G4Isotope * GetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (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 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
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const G4Material * pBaseMaterial
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 
G4double inveplus
 
G4double pFactor
 
- 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()

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

Definition at line 137 of file G4BetheHeitler5DModel.cc.

139 : G4PairProductionRelModel(pd, nam),fVerbose(1),fConversionType(0),
140 iraw(false),
141 fLepton1(G4Electron::Definition()),fLepton2(G4Positron::Definition()),
142 fConvMode(kEPair),
143 fTheMuPlus(G4MuonPlus::Definition()),fTheMuMinus(G4MuonMinus::Definition())
144{
145 theIonTable = G4IonTable::GetIonTable();
146 //Q: Do we need this on Model
147 SetLowEnergyLimit(2*fTheElectron->GetPDGMass());
148}
kEPair
const G4int kEPair
Definition: G4BetheHeitler5DModel.cc:131
G4Electron::Definition
static G4Electron * Definition()
Definition: G4Electron.cc:48
G4IonTable::GetIonTable
static G4IonTable * GetIonTable()
Definition: G4IonTable.cc:170
G4MuonMinus::Definition
static G4MuonMinus * Definition()
Definition: G4MuonMinus.cc:51
G4MuonPlus::Definition
static G4MuonPlus * Definition()
Definition: G4MuonPlus.cc:51
G4PairProductionRelModel::fTheElectron
G4ParticleDefinition * fTheElectron
Definition: G4PairProductionRelModel.hh:186
G4Positron::Definition
static G4Positron * Definition()
Definition: G4Positron.cc:48
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:764

◆ ~G4BetheHeitler5DModel()

G4BetheHeitler5DModel::~G4BetheHeitler5DModel ( )
virtual

Definition at line 152 of file G4BetheHeitler5DModel.cc.

153{}

Member Function Documentation

◆ Initialise()

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

Reimplemented from G4PairProductionRelModel.

Reimplemented in G4LivermoreGammaConversion5DModel.

Definition at line 157 of file G4BetheHeitler5DModel.cc.

159{
160 G4PairProductionRelModel::Initialise(part, vec);
161
162 G4EmParameters* theManager = G4EmParameters::Instance();
163 // place to initialise model parameters
164 // Verbosity levels: ( Can redefine as needed, but some consideration )
165 // 0 = nothing
166 // > 2 print results
167 // > 3 print rejection warning from transformation (fix bug from gammaray .. )
168 // > 4 print photon direction & polarisation
169 fVerbose = theManager->Verbose();
170 fConversionType = theManager->GetConversionType();
171 //////////////////////////////////////////////////////////////
172 // iraw :
173 // true : isolated electron or nucleus.
174 // false : inside atom -> screening form factor
175 iraw = theManager->OnIsolated();
176 // G4cout << "BH5DModel::Initialise verbose " << fVerbose
177 // << " isolated " << iraw << " ctype "<< fConversionType << G4endl;
178
179 //Q: Do we need this on Model
180 // The Leptons defined via SetLeptonPair(..) method
181 SetLowEnergyLimit(2*CLHEP::electron_mass_c2);
182
183 if (fConvMode == kEPair) {
184 assert(fLepton1->GetPDGEncoding() == fTheElectron->GetPDGEncoding()) ;
185 if (fVerbose > 3)
186 G4cout << "BH5DModel::Initialise conversion to e+ e-" << G4endl;
187 }
188
189 if (fConvMode == kMuPair) {
190 assert(fLepton1->GetPDGEncoding() == fTheMuMinus->GetPDGEncoding()) ;
191 if (fVerbose > 3)
192 G4cout << "BH5DModel::Initialise conversion to mu+ mu-" << G4endl;
193 }
194}
kMuPair
const G4int kMuPair
Definition: G4BetheHeitler5DModel.cc:132
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition: G4EmParameters.cc:67
G4EmParameters::OnIsolated
G4bool OnIsolated() const
Definition: G4EmParameters.cc:491
G4EmParameters::GetConversionType
G4int GetConversionType() const
Definition: G4EmParameters.cc:1054
G4EmParameters::Verbose
G4int Verbose() const
Definition: G4EmParameters.cc:975
G4PairProductionRelModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4PairProductionRelModel.cc:146

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

◆ SampleSecondaries()

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

Reimplemented from G4PairProductionRelModel.

Definition at line 280 of file G4BetheHeitler5DModel.cc.

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

200{
201 // Lepton1 - nagative charged particle
202 if ( p1->GetPDGEncoding() < 0 ){
203 if ( p1->GetPDGEncoding() ==
204 G4Positron::Definition()->GetPDGEncoding() ) {
205 SetConversionMode(kEPair);
206 fLepton1 = p2;
207 fLepton2 = p1;
208 // if (fVerbose)
209 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to e+ e-"
210 << G4endl;
211 } else if ( p1->GetPDGEncoding() ==
212 G4MuonPlus::Definition()->GetPDGEncoding() ) {
213 SetConversionMode(kMuPair);
214 fLepton1 = p2;
215 fLepton2 = p1;
216 // if (fVerbose)
217 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to mu+ mu-"
218 << G4endl;
219 } else {
220 // Exception
221 G4ExceptionDescription ed;
222 ed << "Model not applicable to particle(s) "
223 << p1->GetParticleName() << ", "
224 << p2->GetParticleName();
225 G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0002",
226 FatalException, ed);
227 }
228 } else {
229 if ( p1->GetPDGEncoding() ==
230 G4Electron::Definition()->GetPDGEncoding() ) {
231 SetConversionMode(kEPair);
232 fLepton1 = p1;
233 fLepton2 = p2;
234 // if (fVerbose)
235 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to e+ e-"
236 << G4endl;
237 } else if ( p1->GetPDGEncoding() ==
238 G4MuonMinus::Definition()->GetPDGEncoding() ) {
239 SetConversionMode(kMuPair);
240 fLepton1 = p1;
241 fLepton2 = p2;
242 // if (fVerbose)
243 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to mu+ mu-"
244 << G4endl;
245 } else {
246 // Exception
247 G4ExceptionDescription ed;
248 ed << "Model not applicable to particle(s) "
249 << p1->GetParticleName() << ", "
250 << p2->GetParticleName();
251 G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0002",
252 FatalException, ed);
253 }
254 }
255 if ( fLepton1->GetPDGEncoding() != fLepton2->GetAntiPDGEncoding() ) {
256 G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0007",
257 FatalErrorInArgument, "pair must be particle, antiparticle ");
258 G4cerr << "BH5DModel::SetLeptonPair BAD paricle/anti particle pair"
259 << fLepton1->GetParticleName() << ", "
260 << fLepton2->GetParticleName() << G4endl;
261 }
262}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
FatalErrorInArgument
@ FatalErrorInArgument
Definition: G4ExceptionSeverity.hh:62
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4cerr
G4GLOB_DLL std::ostream G4cerr
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:108

Referenced by G4GammaConversionToMuons::BuildPhysicsTable().

◆ SetVerbose()

void G4BetheHeitler5DModel::SetVerbose ( G4int  val)
inline

Definition at line 82 of file G4BetheHeitler5DModel.hh.

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