Geant4 11.1.1
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 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 *, 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 SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
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 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 *)
 
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 = 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
 
size_t currentCoupleIndex = 0
 
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:48
G4IonTable::GetIonTable
static G4IonTable * GetIonTable()
Definition: G4IonTable.cc:170
G4PairProductionRelModel::fTheElectron
G4ParticleDefinition * fTheElectron
Definition: G4PairProductionRelModel.hh:179
G4Positron::Definition
static G4Positron * Definition()
Definition: G4Positron.cc:48
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:753

◆ ~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.

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:490
G4EmParameters::GetConversionType
G4int GetConversionType() const
Definition: G4EmParameters.cc:1084
G4EmParameters::Verbose
G4int Verbose() const
Definition: G4EmParameters.cc:983
G4PairProductionRelModel::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: G4PairProductionRelModel.cc:147

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];
392 // START Sampling
393 do {
394
395 rndmEngine->flatArray(6, rndmv6);
396
397 //////////////////////////////////////////////////
398 // pdf pow(x,c) with c = 1.4
399 // integral y = pow(x,(c+1))/(c+1) @ x = 1 => y = 1 /(1+c)
400 // invCdf exp( log(y /* *( c + 1.0 )/ (c + 1.0 ) */ ) /( c + 1.0) )
401 //////////////////////////////////////////////////
402 const G4double X1 =
403 G4Exp(G4Log(rndmv6[0])/(correctionIndex + 1.0));
404
405 const G4double x0 = G4Exp(xl1 + (xu1 - xl1)*rndmv6[1]);
406 const G4double dum0 = 1./(1.+x0);
407 const G4double cosTheta = (x0-1.)*dum0;
408 const G4double sinTheta = std::sqrt(4.*x0)*dum0;
409
410 const G4double PairInvMass = PairInvMassMin*G4Exp(X1*X1*lnPairInvMassRange);
411
412 // G4double rndmv3[3];
413 // rndmEngine->flatArray(3, rndmv3);
414
415 // cos and sin theta-lepton
416 const G4double cosThetaLept = std::cos(pi*rndmv6[2]);
417 // sin(ThetaLept) is always in [0,+1] if ThetaLept is in [0,pi]
418 const G4double sinThetaLept = std::sqrt((1.-cosThetaLept)*(1.+cosThetaLept));
419 // cos and sin phi-lepton
420 const G4double cosPhiLept = std::cos(twoPi*rndmv6[3]-pi);
421 // sin(PhiLept) is in [-1,0] if PhiLept in [-pi,0) and
422 // is in [0,+1] if PhiLept in [0,+pi]
423 const G4double sinPhiLept = std::copysign(std::sqrt((1.-cosPhiLept)*(1.+cosPhiLept)),rndmv6[3]-0.5);
424 // cos and sin phi
425 const G4double cosPhi = std::cos(twoPi*rndmv6[4]-pi);
426 const G4double sinPhi = std::copysign(std::sqrt((1.-cosPhi)*(1.+cosPhi)),rndmv6[4]-0.5);
427
428 //////////////////////////////////////////////////
429 // frames:
430 // 3 : the laboratory Lorentz frame, Geant4 axes definition
431 // 0 : the laboratory Lorentz frame, axes along photon direction and polarisation
432 // 1 : the center-of-mass Lorentz frame
433 // 2 : the pair Lorentz frame
434 //////////////////////////////////////////////////
435
436 // in the center-of-mass frame
437
438 const G4double RecEnergyCMS = (sCMSPlusRM2-PairInvMass*PairInvMass)*isqrts2;
439 const G4double LeptonEnergy2 = PairInvMass*0.5;
440
441 // New way of calucaltion thePRecoil to avoid underflow
442 G4double abp = std::max((2.0*GammaEnergy*RecoilMass -
443 PairInvMass*PairInvMass + 2.0*PairInvMass*RecoilMass)*
444 (2.0*GammaEnergy*RecoilMass -
445 PairInvMass*PairInvMass - 2.0*PairInvMass*RecoilMass),0.0);
446
447 G4double thePRecoil = std::sqrt(abp) * isqrts2;
448
449 // back to the center-of-mass frame
450 Recoil.set( thePRecoil*sinTheta*cosPhi,
451 thePRecoil*sinTheta*sinPhi,
452 thePRecoil*cosTheta,
453 RecEnergyCMS);
454
455 // in the pair frame
456 const G4double thePLepton = std::sqrt( (LeptonEnergy2-LeptonMass)
457 *(LeptonEnergy2+LeptonMass));
458
459 LeptonPlus.set(thePLepton*sinThetaLept*cosPhiLept,
460 thePLepton*sinThetaLept*sinPhiLept,
461 thePLepton*cosThetaLept,
462 LeptonEnergy2);
463
464 LeptonMinus.set(-LeptonPlus.x(),
465 -LeptonPlus.y(),
466 -LeptonPlus.z(),
467 LeptonEnergy2);
468
469
470 // Normalisation of final state phase space:
471 // Section 47 of Particle Data Group, Chin. Phys. C, 40, 100001 (2016)
472 // const G4double Norme = Recoil1.vect().mag() * LeptonPlus2.vect().mag();
473 const G4double Norme = Recoil.vect().mag() * LeptonPlus.vect().mag();
474
475 // e+, e- to CMS frame from pair frame
476
477 // boost vector from Pair to CMS
478 const G4ThreeVector pair2cms =
479 G4LorentzVector( -Recoil.x(), -Recoil.y(), -Recoil.z(),
480 sqrts-RecEnergyCMS).boostVector();
481
482 LeptonPlus.boost(pair2cms);
483 LeptonMinus.boost(pair2cms);
484
485 // back to the laboratory frame (make use of the CMS(0,0,Eg,Eg+RM)) form
486
487 Recoil.boostZ(betaCMS);
488 LeptonPlus.boostZ(betaCMS);
489 LeptonMinus.boostZ(betaCMS);
490
491 // Jacobian factors
492 const G4double Jacob0 = x0*dum0*dum0;
493 const G4double Jacob1 = 2.*X1*lnPairInvMassRange*PairInvMass;
494 const G4double Jacob2 = std::abs(sinThetaLept);
495
496 const G4double EPlus = LeptonPlus.t();
497 const G4double PPlus = LeptonPlus.vect().mag();
498 const G4double sinThetaPlus = LeptonPlus.vect().perp()/PPlus;
499 const G4double cosThetaPlus = LeptonPlus.vect().cosTheta();
500
501 const G4double pPX = LeptonPlus.x();
502 const G4double pPY = LeptonPlus.y();
503 const G4double dum1 = 1./std::sqrt( pPX*pPX + pPY*pPY );
504 const G4double cosPhiPlus = pPX*dum1;
505 const G4double sinPhiPlus = pPY*dum1;
506
507 // denominators:
508 // the two cancelling leading terms for forward emission at high energy, removed
509 const G4double elMassCTP = LeptonMass*cosThetaPlus;
510 const G4double ePlusSTP = EPlus*sinThetaPlus;
511 const G4double DPlus = (elMassCTP*elMassCTP + ePlusSTP*ePlusSTP)
512 /(EPlus + PPlus*cosThetaPlus);
513
514 const G4double EMinus = LeptonMinus.t();
515 const G4double PMinus = LeptonMinus.vect().mag();
516 const G4double sinThetaMinus = LeptonMinus.vect().perp()/PMinus;
517 const G4double cosThetaMinus = LeptonMinus.vect().cosTheta();
518
519 const G4double ePX = LeptonMinus.x();
520 const G4double ePY = LeptonMinus.y();
521 const G4double dum2 = 1./std::sqrt( ePX*ePX + ePY*ePY );
522 const G4double cosPhiMinus = ePX*dum2;
523 const G4double sinPhiMinus = ePY*dum2;
524
525 const G4double elMassCTM = LeptonMass*cosThetaMinus;
526 const G4double eMinSTM = EMinus*sinThetaMinus;
527 const G4double DMinus = (elMassCTM*elMassCTM + eMinSTM*eMinSTM)
528 /(EMinus + PMinus*cosThetaMinus);
529
530 // cos(phiMinus-PhiPlus)
531 const G4double cosdPhi = cosPhiPlus*cosPhiMinus + sinPhiPlus*sinPhiMinus;
532 const G4double PRec = Recoil.vect().mag();
533 const G4double q2 = PRec*PRec;
534
535 const G4double BigPhi = -LeptonMass2 / (GammaEnergy*GammaEnergy2 * q2*q2);
536
537 G4double FormFactor = 1.;
538 if (!iraw) {
539 if (itriplet) {
540 const G4double qun = factor1*iZ13*iZ13;
541 const G4double nun = qun * PRec;
542 if (nun < 1.) {
543 FormFactor = (nun < 0.01) ? (13.8-55.4*std::sqrt(nun))*nun
544 : std::sqrt(1-(nun-1)*(nun-1));
545 } // else FormFactor = 1 by default
546 } else {
547 const G4double dum3 = 217.*PRec*iZ13;
548 const G4double AFF = 1./(1. + dum3*dum3);
549 FormFactor = (1.-AFF)*(1-AFF);
550 }
551 } // else FormFactor = 1 by default
552
553 G4double betheheitler;
554 if (GammaPolarizationMag==0.) {
555 const G4double pPlusSTP = PPlus*sinThetaPlus;
556 const G4double pMinusSTM = PMinus*sinThetaMinus;
557 const G4double pPlusSTPperDP = pPlusSTP/DPlus;
558 const G4double pMinusSTMperDM = pMinusSTM/DMinus;
559 const G4double dunpol = BigPhi*(
560 pPlusSTPperDP *pPlusSTPperDP *(4.*EMinus*EMinus-q2)
561 + pMinusSTMperDM*pMinusSTMperDM*(4.*EPlus*EPlus - q2)
562 + 2.*pPlusSTPperDP*pMinusSTMperDM*cosdPhi
563 *(4.*EPlus*EMinus + q2 - 2.*GammaEnergy2)
564 - 2.*GammaEnergy2*(pPlusSTP*pPlusSTP+pMinusSTM*pMinusSTM)/(DMinus*DPlus));
565 betheheitler = dunpol * factor;
566 } else {
567 const G4double pPlusSTP = PPlus*sinThetaPlus;
568 const G4double pMinusSTM = PMinus*sinThetaMinus;
569 const G4double pPlusSTPCPPperDP = pPlusSTP*cosPhiPlus/DPlus;
570 const G4double pMinusSTMCPMperDM = pMinusSTM*cosPhiMinus/DMinus;
571 const G4double caa = 2.*(EPlus*pMinusSTMCPMperDM+EMinus*pPlusSTPCPPperDP);
572 const G4double cbb = pMinusSTMCPMperDM-pPlusSTPCPPperDP;
573 const G4double ccc = (pPlusSTP*pPlusSTP + pMinusSTM*pMinusSTM
574 +2.*pPlusSTP*pMinusSTM*cosdPhi)/ (DMinus*DPlus);
575 const G4double dtot= 2.*BigPhi*( caa*caa - q2*cbb*cbb - GammaEnergy2*ccc);
576 betheheitler = dtot * factor;
577 }
578 //
579 const G4double cross = Norme * Jacob0 * Jacob1 * Jacob2 * betheheitler
580 * FormFactor * RecoilMass / sqrts;
581 pdf = cross * (xu1 - xl1) / G4Exp(correctionIndex*G4Log(X1)); // cond1;
582 } while ( pdf < ymax * rndmv6[5] );
583 // END of Sampling
584
585 if ( fVerbose > 2 ) {
586 G4double recul = std::sqrt(Recoil.x()*Recoil.x()+Recoil.y()*Recoil.y()
587 +Recoil.z()*Recoil.z());
588 G4cout << "BetheHeitler5DModel GammaEnergy= " << GammaEnergy
589 << " PDF= " << pdf << " ymax= " << ymax
590 << " recul= " << recul << G4endl;
591 }
592
593 // back to Geant4 system
594
595 if ( fVerbose > 4 ) {
596 G4cout << "BetheHeitler5DModel GammaDirection " << GammaDirection << G4endl;
597 G4cout << "BetheHeitler5DModel GammaPolarization " << GammaPolarization << G4endl;
598 G4cout << "BetheHeitler5DModel GammaEnergy " << GammaEnergy << G4endl;
599 G4cout << "BetheHeitler5DModel Conv "
600 << (itriplet ? "triplet" : "nucl") << G4endl;
601 }
602
603 if (GammaPolarizationMag == 0.0) {
604 // set polarization axis orthohonal to direction
605 GammaPolarization = GammaDirection.orthogonal().unit();
606 } else {
607 // GammaPolarization not a unit vector
608 GammaPolarization /= GammaPolarizationMag;
609 }
610
611 // The unit norm vector that is orthogonal to the two others
612 G4ThreeVector yGrec = GammaDirection.cross(GammaPolarization);
613
614 // rotation from gamma ref. sys. to World
615 G4RotationMatrix GtoW(GammaPolarization,yGrec,GammaDirection);
616
617 Recoil.transform(GtoW);
618 LeptonPlus.transform(GtoW);
619 LeptonMinus.transform(GtoW);
620
621 if ( fVerbose > 2 ) {
622 G4cout << "BetheHeitler5DModel Recoil " << Recoil.x() << " " << Recoil.y() << " " << Recoil.z()
623 << " " << Recoil.t() << " " << G4endl;
624 G4cout << "BetheHeitler5DModel LeptonPlus " << LeptonPlus.x() << " " << LeptonPlus.y() << " "
625 << LeptonPlus.z() << " " << LeptonPlus.t() << " " << G4endl;
626 G4cout << "BetheHeitler5DModel LeptonMinus " << LeptonMinus.x() << " " << LeptonMinus.y() << " "
627 << LeptonMinus.z() << " " << LeptonMinus.t() << " " << G4endl;
628 }
629
630 // Create secondaries
631 auto aParticle1 = new G4DynamicParticle(fLepton1,LeptonMinus);
632 auto aParticle2 = new G4DynamicParticle(fLepton2,LeptonPlus);
633
634 // create G4DynamicParticle object for the particle3 ( recoil )
635 G4ParticleDefinition* RecoilPart;
636 if (itriplet) {
637 // triplet
638 RecoilPart = fTheElectron;
639 } else{
640 RecoilPart = theIonTable->GetIon(Z, A, 0);
641 }
642 auto aParticle3 = new G4DynamicParticle(RecoilPart,Recoil);
643
644 // Fill output vector
645 fvect->push_back(aParticle1);
646 fvect->push_back(aParticle2);
647 fvect->push_back(aParticle3);
648
649 // kill incident photon
650 fParticleChange->SetProposedKineticEnergy(0.);
651 fParticleChange->ProposeTrackStatus(fStopAndKill);
652}
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
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
G4endl
#define G4endl
Definition: G4ios.hh:57
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:198
G4Element::GetZasInt
G4int GetZasInt() const
Definition: G4Element.hh:132
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:178
G4PairProductionRelModel::fParticleChange
G4ParticleChangeForGamma * fParticleChange
Definition: G4PairProductionRelModel.hh:181
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:106
G4VEmModel::SelectIsotopeNumber
G4int SelectIsotopeNumber(const G4Element *) const
Definition: G4VEmModel.cc:276
G4VEmModel::SelectTargetAtom
const G4Element * SelectTargetAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:577
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:62
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:109

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: