Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4CoherentPairProduction Class Reference

#include <G4CoherentPairProduction.hh>

+ Inheritance diagram for G4CoherentPairProduction:

Public Member Functions

 G4CoherentPairProduction (const G4String &processName="cpp", G4ProcessType aType=fElectromagnetic)
 
 ~G4CoherentPairProduction ()=default
 
G4VParticleChange * PostStepDoIt (const G4Track &, const G4Step &) override
 
G4bool IsApplicable (const G4ParticleDefinition &aPD) override
 
void ProcessDescription (std::ostream &) const override
 
void Input (const G4Material *crystal, const G4String &lattice)
 special functions
 
void Input (const G4Material *crystal, const G4String &lattice, const G4String &filePath)
 
void Input (const G4ChannelingFastSimCrystalData *crystalData)
 
void ActivateIncoherentScattering ()
 
G4ChannelingFastSimCrystalData * GetCrystalData ()
 
G4double ModelMinPrimaryEnergy ()
 get cuts
 
G4double GetHighAngleLimit ()
 
G4double GetPPKineticEnergyCut ()
 
G4int GetSamplingPairsNumber ()
 
G4double GetChargeParticleAngleFactor ()
 
G4double GetNTrajectorySteps ()
 get number of trajectory steps of a single particle (e- or e+)
 
G4double GetEffectiveLrad ()
 
G4String GetG4RegionName ()
 get the name of G4Region in which the model is applicable
 
void SetLowEnergyLimit (G4double energy)
 set cuts
 
void SetHighAngleLimit (G4double angle)
 
void SetPPKineticEnergyCut (G4double kineticEnergyCut)
 
void SetSamplingPairsNumber (G4int nPairs)
 
void SetChargeParticleAngleFactor (G4double chargeParticleAngleFactor)
 
void SetNTrajectorySteps (G4int nTrajectorySteps)
 set number of trajectory steps of a single particle (e- or e+)
 
void SetG4RegionName (const G4String &nameG4Region)
 set the name of G4Region in which the model is applicable
 
G4double GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *condition) override
 
- Public Member Functions inherited from G4VDiscreteProcess
 G4VDiscreteProcess (const G4String &aName, G4ProcessType aType=fNotDefined)
 
 G4VDiscreteProcess (G4VDiscreteProcess &)
 
virtual ~G4VDiscreteProcess ()
 
G4VDiscreteProcess & operator= (const G4VDiscreteProcess &)=delete
 
virtual G4double PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual G4double AlongStepGetPhysicalInteractionLength (const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
 
virtual G4double AtRestGetPhysicalInteractionLength (const G4Track &, G4ForceCondition *)
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &, const G4Step &)
 
virtual G4VParticleChange * AlongStepDoIt (const G4Track &, const G4Step &)
 
virtual G4double GetCrossSection (const G4double, const G4MaterialCutsCouple *)
 
virtual G4double MinPrimaryEnergy (const G4ParticleDefinition *, const G4Material *)
 
- Public Member Functions inherited from G4VProcess
 G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)
 
 G4VProcess (const G4VProcess &right)
 
virtual ~G4VProcess ()
 
G4VProcess & operator= (const G4VProcess &)=delete
 
G4bool operator== (const G4VProcess &right) const
 
G4bool operator!= (const G4VProcess &right) const
 
G4double GetCurrentInteractionLength () const
 
void SetPILfactor (G4double value)
 
G4double GetPILfactor () const
 
G4double AlongStepGPIL (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)
 
G4double AtRestGPIL (const G4Track &track, G4ForceCondition *condition)
 
G4double PostStepGPIL (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void PreparePhysicsTable (const G4ParticleDefinition &)
 
virtual G4bool StorePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
virtual G4bool RetrievePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
 
const G4String & GetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
 
const G4String & GetProcessName () const
 
G4ProcessType GetProcessType () const
 
void SetProcessType (G4ProcessType)
 
G4int GetProcessSubType () const
 
void SetProcessSubType (G4int)
 
virtual const G4VProcess * GetCreatorProcess () const
 
virtual void StartTracking (G4Track *)
 
virtual void EndTracking ()
 
virtual void SetProcessManager (const G4ProcessManager *)
 
virtual const G4ProcessManager * GetProcessManager ()
 
virtual void ResetNumberOfInteractionLengthLeft ()
 
G4double GetNumberOfInteractionLengthLeft () const
 
G4double GetTotalNumberOfInteractionLengthTraversed () const
 
G4bool isAtRestDoItIsEnabled () const
 
G4bool isAlongStepDoItIsEnabled () const
 
G4bool isPostStepDoItIsEnabled () const
 
virtual void DumpInfo () const
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
virtual void SetMasterProcess (G4VProcess *masterP)
 
const G4VProcess * GetMasterProcess () const
 
virtual void BuildWorkerPhysicsTable (const G4ParticleDefinition &part)
 
virtual void PrepareWorkerPhysicsTable (const G4ParticleDefinition &)
 

Additional Inherited Members

- Static Public Member Functions inherited from G4VProcess
static const G4String & GetProcessTypeName (G4ProcessType)
 
- Protected Member Functions inherited from G4VProcess
void SubtractNumberOfInteractionLengthLeft (G4double prevStepSize)
 
void ClearNumberOfInteractionLengthLeft ()
 
- Protected Attributes inherited from G4VProcess
const G4ProcessManager * aProcessManager = nullptr
 
G4VParticleChange * pParticleChange = nullptr
 
G4ParticleChange aParticleChange
 
G4double theNumberOfInteractionLengthLeft = -1.0
 
G4double currentInteractionLength = -1.0
 
G4double theInitialNumberOfInteractionLength = -1.0
 
G4String theProcessName
 
G4String thePhysicsTableFileName
 
G4ProcessType theProcessType = fNotDefined
 
G4int theProcessSubType = -1
 
G4double thePILfactor = 1.0
 
G4int verboseLevel = 0
 
G4bool enableAtRestDoIt = true
 
G4bool enableAlongStepDoIt = true
 
G4bool enablePostStepDoIt = true
 

Detailed Description

Definition at line 47 of file G4CoherentPairProduction.hh.

Constructor & Destructor Documentation

◆ G4CoherentPairProduction()

G4CoherentPairProduction::G4CoherentPairProduction ( const G4String & processName = "cpp",
G4ProcessType aType = fElectromagnetic )

Definition at line 50 of file G4CoherentPairProduction.cc.

51 :
52 G4VDiscreteProcess(aName)
53{
54 SetProcessSubType(fCoherentPairProduction);
55}
fCoherentPairProduction
@ fCoherentPairProduction
Definition G4EmProcessSubType.hh:65
G4VDiscreteProcess::G4VDiscreteProcess
G4VDiscreteProcess(const G4String &aName, G4ProcessType aType=fNotDefined)
Definition G4VDiscreteProcess.cc:51
G4VProcess::SetProcessSubType
void SetProcessSubType(G4int)
Definition G4VProcess.hh:410

◆ ~G4CoherentPairProduction()

G4CoherentPairProduction::~G4CoherentPairProduction ( )
default

Member Function Documentation

◆ ActivateIncoherentScattering()

void G4CoherentPairProduction::ActivateIncoherentScattering ( )
inline

activate incoherent scattering (standard gamma conversion should be switched off in physics list)

Definition at line 79 of file G4CoherentPairProduction.hh.

79{fIncoherentScattering = true;}

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

◆ GetChargeParticleAngleFactor()

G4double G4CoherentPairProduction::GetChargeParticleAngleFactor ( )
inline

get the number of particle angles 1/gamma in pair production defining the width of the angular distribution of pair sampling in the Baier-Katkov Integral

Definition at line 96 of file G4CoherentPairProduction.hh.

96{return fChargeParticleAngleFactor;}

◆ GetCrystalData()

G4ChannelingFastSimCrystalData * G4CoherentPairProduction::GetCrystalData ( )
inline

Definition at line 81 of file G4CoherentPairProduction.hh.

81{return fCrystalData;}

◆ GetEffectiveLrad()

G4double G4CoherentPairProduction::GetEffectiveLrad ( )
inline

get effective radiation length (due to coherent process of pair production) simulated for the current photon

Definition at line 104 of file G4CoherentPairProduction.hh.

104{return fEffectiveLrad;}

◆ GetG4RegionName()

G4String G4CoherentPairProduction::GetG4RegionName ( )
inline

get the name of G4Region in which the model is applicable

Definition at line 107 of file G4CoherentPairProduction.hh.

107{return fG4RegionName;}

◆ GetHighAngleLimit()

G4double G4CoherentPairProduction::GetHighAngleLimit ( )
inline

Definition at line 86 of file G4CoherentPairProduction.hh.

86{return fHighAngleLimit;}

Referenced by GetMeanFreePath().

◆ GetMeanFreePath()

G4double G4CoherentPairProduction::GetMeanFreePath ( const G4Track & aTrack,
G4double ,
G4ForceCondition * condition )
overridevirtual

Implements G4VDiscreteProcess.

Definition at line 59 of file G4CoherentPairProduction.cc.

62{
63 //current logical volume
64 G4LogicalVolume* crystallogic;
65
66 //momentum direction and coordinates (see comments below)
67 G4ThreeVector momentumDirectionGamma,xyzGamma0,xyzGamma;
68 //angle of the photon in the local reference system of the volume
69 G4double txGamma0 = 0, tyGamma0 = 0;
70
71 *condition = NotForced;
72
73 //model activation
74 G4bool modelTrigger = false;
75
76 //photon energy
77 G4double eGamma = aTrack.GetTotalEnergy();
78
79 //energy cut, at the beginning, to not check everything else
80 if(eGamma > ModelMinPrimaryEnergy())
81 {
82 //current logical volume
83 crystallogic = aTrack.GetVolume()->GetLogicalVolume();
84
85 //the model works only in the G4Region fG4RegionName
86 if(crystallogic->GetRegion()->GetName()==fG4RegionName)
87 {
88 fCrystalData->SetGeometryParameters(crystallogic);
89
90 //the momentum direction of the photon in the local reference system of the volume
91 momentumDirectionGamma =
92 (aTrack.GetTouchableHandle()->GetHistory()->
93 GetTopTransform().NetRotation().inverse())*aTrack.GetMomentumDirection();
94
95 //the coordinates of the photon in the local reference system of the volume
96 xyzGamma0 =
97 aTrack.GetTouchableHandle()->GetHistory()->
98 GetTopTransform().TransformPoint(aTrack.GetPosition());
99
100 // the coordinates of the photon in the co-rotating reference system within
101 //a channel (elementary periodic cell)
102 xyzGamma = fCrystalData->CoordinatesFromBoxToLattice(xyzGamma0);
103
104 //angle of the photon in the local reference system of the volume
105 //(!!! ONLY FORWARD DIRECTION, momentumDirectionGamma.getZ()>0,
106 txGamma0 = std::atan(momentumDirectionGamma.x()/momentumDirectionGamma.z());
107 tyGamma0 = std::atan(momentumDirectionGamma.y()/momentumDirectionGamma.z());
108
109 //recalculate angle into the lattice reference system
110 G4double angle = fCrystalData->AngleXFromBoxToLattice(txGamma0,xyzGamma.z());
111 if (fCrystalData->GetModel()==2)
112 {
113 angle = std::sqrt(angle*angle+tyGamma0*tyGamma0);
114 }
115
116 //Applies the parameterisation not at the last step, only forward local direction
117 //above low energy limit and below angular limit
118 modelTrigger = (momentumDirectionGamma.z()>0. &&
119 std::abs(angle) < GetHighAngleLimit());
120 }
121 }
122
123 if(modelTrigger)
124 {
125 //execute the model
126
127 G4double x=0.,y=0.,z=0.;// the coordinates of charged particles
128 //in the co-rotating reference system within
129 //a channel (elementary periodic cell)
130 G4double tx0=0.,ty0=0.; // the angles of charged particles
131 // in the local reference system of the volume
132 G4double txPreStep0=0.,tyPreStep0=0.; // the same as tx0, ty0 before the step
133 // in the co-rotating reference system within
134 //a channel (elementary periodic cell)
135
136 G4ThreeVector scatteringAnglesAndEnergyLoss;//output of scattering functions
137
138 //coordinates in Runge-Kutta calculations
139 G4double x1=0.,x2=0.,x3=0.,x4=0.,y1=0.,y2=0.,y3=0.,y4=0.;
140 //angles in Runge-Kutta calculations
141 G4double tx1=0.,tx2=0.,tx3=0.,tx4=0.,ty1=0.,ty2=0.,ty3=0.,ty4=0.;
142 //variables in Runge-Kutta calculations
143 G4double kvx1=0.,kvx2=0.,kvx3=0.,kvx4=0.,kvy1=0.,kvy2=0.,kvy3=0.,kvy4=0.;
144 //simulation step along z (internal step of the model) and its parts
145 G4double dz=0.,dzd3=0.,dzd8=0.;//dzd3 = dz/3; dzd8 = dz/8;
146 //simulation step along the momentum direction
147 G4double momentumDirectionStep;
148 //effective simulation step (taking into account nuclear density along the trajectory)
149 G4double effectiveStep=0.;
150
151 // Baier-Katkov variables
152 G4double dzMeV=0.; //step in MeV^-1
153 G4double axt=0.,ayt=0.; //charged particle accelerations
154 G4double vxin=0.,vyin=0.;//the angles vs the photon (with incoherent scattering)
155 G4double vxno=0.,vyno=0.;//the angles vs the photon (without incoherent scattering)
156
157 G4double dzmod=0.;
158 G4double fa1=0.,faseBefore=0.,faseBeforedz=0.,faseBeforedzd2=0.;
159 G4double faseAfter=0.,fa2dfaseBefore2=0.;
160
161 G4double skJ=0, skIx=0., skIy=0.;
162 G4double sinfa1=0.,cosfa1=0.;
163
164 //2-vector is needed for an initial parameter collection of 1 pair
165 //vector of 2-vectors is an initial parameter collection of all sampling pair
166
167 //collection of etotal for a single pair
168 CLHEP::Hep2Vector twoVectorEtotal(0.,0.);
169
170 //collection of x for a single pair
171 CLHEP::Hep2Vector twoVectorX(0.,0.);
172
173 //collection of y for a single pair
174 CLHEP::Hep2Vector twoVectorY(0.,0.);
175
176 //collection of tx for a single pair
177 CLHEP::Hep2Vector twoVectorTX(0.,0.);
178
179 //collection of tx for a single pair
180 CLHEP::Hep2Vector twoVectorTY(0.,0.);
181
182 fullVectorEtotal.clear();
183 fullVectorX.clear();
184 fullVectorY.clear();
185 fullVectorTX.clear();
186 fullVectorTY.clear();
187 fPairProductionCDFdz.clear();
188 fPairProductionCDFdz.push_back(0.);//0th element equal to 0
189
190 const G4double charge[2] = {-1.,1.}; //particle charge
191 const G4String particleName[2] = {"e-", "e+"};
192
193 // the coordinates of a charged particle in the reference system within
194 //a channel (elementary periodic cell)
195 G4ThreeVector xyzparticle = xyzGamma;//changed below
196
197 //the idea of pair production simulation is analogical to radiation in G4BaierKatkov
198 //since the matrix element of these processes is the same => we solve inverse problem
199 //to radiation: sample the pairs, calculate their trajectories and then calculate the
200 //probabilities using Baier-Katkov analogically to radiation
201
202 //cycle by sampling e+- pairs
203 for(G4int i=0; i<fNMCPairs;i++)
204 {
205 //pair energy uniform sampling
206 G4double etotal = fMass + fPPKineticEnergyCut +
207 G4UniformRand()*(eGamma-2*(fMass+fPPKineticEnergyCut));//particle
208 //total energy
209
210 G4double phi = CLHEP::twopi*G4UniformRand();//necessary for pair kinematics
211
212 //the probability of the production of the current pair (will be simulated)
213 //per distance
214 G4double probabilityPPdz = 0.;
215
216 //cycle e- and e+ within single pair
217 for(G4int j=0; j<2;j++)
218 {
219 if(j==1){etotal=eGamma-etotal;} //2nd particle energy
220 twoVectorEtotal[j]=etotal;
221
222 //Baier-Katkov input
223 //intermediate variables to reduce calculations (the same names as in G4BaierKatkov)
224 G4double e2 = etotal*etotal;
225 G4double gammaInverse2 = fMass*fMass/(etotal*etotal);// 1/gamma^2
226 //normalization coefficient
227 G4double coefNorm = CLHEP::fine_structure_const/(8*(CLHEP::pi2))/(2.*fNMCPairs);
228 //G4double phi = CLHEP::twopi*G4UniformRand();//necessary for pair kinematics
229 G4double om = eGamma;
230 G4double eprime=om-etotal; //E'=omega-E
231 G4double eprime2 = eprime*eprime;
232 G4double e2pluseprime2 =e2+eprime2;
233 G4double omprime=etotal*om/eprime;//om'=E*om/(om-E)
234 G4double omprimed2=omprime/2;
235
236 //difference vs G4BaierKatkov: om -> etotal
237 G4double coefNorme2deprime2 = coefNorm*e2/eprime2; //e2/om/om;//e2/eprime2;
238
239 G4double gammaInverse2om = gammaInverse2*om*om;
240
241 //initialize intermediate integrals with zeros
242 G4double fa=0.,ss=0.,sc=0.,ssx=0.,ssy=0.,scx=0.,scy=0.;
243
244 //End of Baier-Katkov input
245
246 G4bool fbreak = false;//flag of the trajectory cycle break
247
248 //set fCrystalData parameters depending on the particle parameters
249 fCrystalData->SetParticleProperties(etotal, fMass,
250 charge[j], particleName[j]);
251
252 //needed just to setup the correct value of channel No in the crystal
253 //since later it may be changed during the trajectory calculation
254 fCrystalData->CoordinatesFromBoxToLattice(xyzGamma0);
255
256 //coordinate sampling: random x and y due to coordinate uncertainty
257 //in the interaction point
258 if(j==0)
259 {
260 x = fCrystalData->GetChannelWidthX()*G4UniformRand();
261 y = fCrystalData->GetChannelWidthY()*G4UniformRand();
262 }
263 else
264 {
265 x=twoVectorX[0];
266 y=twoVectorY[0];
267 }
268 twoVectorX[j] = x;
269 twoVectorY[j] = y;
270 //definite z as a coordinate of the photon (uncertainty of the
271 //interaction point is taking into account later by simulation
272 //of the position of pair production)
273 z = xyzGamma.z();
274
275 //angles of the photon in the co-rotating reference system within a channel =>
276 //angular distribution center
277 G4double tx = fCrystalData->AngleXFromBoxToLattice(txGamma0,z);
278 G4double ty = tyGamma0;
279 G4double momentumDirectionZGamma = 1./
280 std::sqrt(1.+std::pow(std::tan(tx),2)+
281 std::pow(std::tan(ty),2));
282
283 //angle sampling: depends on angular range within a particle trajectory
284 //defined by the Lindhard angle and on the angle of radiation proportional
285 //to 1/gamma
286
287 //range of MC integration on angles
288 G4double paramParticleAngle = fChargeParticleAngleFactor*fMass/etotal;
289
290 G4double axangle=0.;
291 if (fCrystalData->GetModel()==1)//1D model (only angle vs plane matters)
292 {
293 axangle = std::abs(tx);
294 }
295 else if (fCrystalData->GetModel()==2)//2D model
296 {
297 axangle = std::sqrt(tx*tx+ty*ty);
298 }
299
300 if(axangle>fCrystalData->GetLindhardAngle()+DBL_EPSILON)
301 {
302 paramParticleAngle+=axangle
303 -std::sqrt(axangle*axangle
304 -fCrystalData->GetLindhardAngle()
305 *fCrystalData->GetLindhardAngle());
306 }
307 else
308 {
309 paramParticleAngle+=fCrystalData->GetLindhardAngle();
310 }
311
312
313 //ONLY forward direction
314 if (paramParticleAngle>CLHEP::halfpi-DBL_EPSILON){paramParticleAngle=CLHEP::halfpi;}
315
316 G4double rho=1.;
317 G4double rhocut=CLHEP::halfpi/paramParticleAngle;//radial angular cut of
318 //the distribution
319 G4double norm=std::atan(rhocut*rhocut)*
320 CLHEP::pi*paramParticleAngle*paramParticleAngle;
321
322
323 //distribution with long tails (useful to not exclude particle angles
324 //after a strong single scattering)
325 //at ellipsescale < 1 => half of statistics
326 do
327 {
328 rho = std::sqrt(std::tan(CLHEP::halfpi*G4UniformRand()));
329 }
330 while (rho>rhocut);
331
332 //normalization coefficient for intergration on angles of charged particles
333 G4double angleNormCoef = (1.+rho*rho*rho*rho)*norm;
334
335 tx+=charge[j]*paramParticleAngle*rho*std::cos(phi);
336 twoVectorTX[j] = tx;
337 ty+=charge[j]*paramParticleAngle*rho*std::sin(phi);
338 twoVectorTY[j] = ty;
339
340 G4double zalongGamma = 0;//necessary for renormalization of PP probability
341 //depending on the trajectory length along Gamma direction
342 //starting the trajectory
343 //here we don't care about the boundaries of the crystal volume
344 //the trajectory is very short and the pair production probability obtained
345 //in Baier-Katkov will be extrapolated to the real step inside the crystal volume
346 for(G4int k=0; k<fNTrajectorySteps;k++)
347 {
348 //back to the local reference system of the volume
349 txPreStep0 = fCrystalData->AngleXFromLatticeToBox(tx,z);
350 tyPreStep0 = ty;
351
352 dz = fCrystalData->GetSimulationStep(tx,ty);
353 dzd3=dz/3;
354 dzd8=dz/8;
355
356 //trajectory calculation:
357 //Runge-Cutt "3/8"
358 //fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ() is due to dependence
359 //of the radius on x; GetCurv gets 1/R for the central ("central plane/axis")
360
361 //first step
362 kvx1=fCrystalData->Ex(x,y);
363 x1=x+tx*dzd3;
364 tx1=tx+(kvx1-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dzd3;
365 if (fCrystalData->GetModel()==2)
366 {
367 kvy1=fCrystalData->Ey(x,y);
368 y1=y+ty*dzd3;
369 ty1=ty+kvy1*dzd3;
370 }
371
372 //second step
373 kvx2=fCrystalData->Ex(x1,y1);
374 x2=x-tx*dzd3+tx1*dz;
375 tx2=tx-(kvx1-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dzd3+
376 (kvx2-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dz;
377 if (fCrystalData->GetModel()==2)
378 {
379 kvy2=fCrystalData->Ey(x1,y1);
380 y2=y-ty*dzd3+ty1*dz;
381 ty2=ty-kvy1*dzd3+kvy2*dz;
382 }
383
384 //third step
385 kvx3=fCrystalData->Ex(x2,y2);
386 x3=x+(tx-tx1+tx2)*dz;
387 tx3=tx+(kvx1-kvx2+kvx3-
388 fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dz;
389 if (fCrystalData->GetModel()==2)
390 {
391 kvy3=fCrystalData->Ey(x2,y2);
392 y3=y+(ty-ty1+ty2)*dz;
393 ty3=ty+(kvy1-kvy2+kvy3)*dz;
394 }
395
396 //fourth step
397 kvx4=fCrystalData->Ex(x3,y3);
398 x4=x+(tx+3.*tx1+3.*tx2+tx3)*dzd8;
399 tx4=tx+(kvx1+3.*kvx2+3.*kvx3+kvx4)*dzd8-
400 fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ()*dz;
401 if (fCrystalData->GetModel()==2)
402 {
403 kvy4=fCrystalData->Ey(x3,y3);
404 y4=y+(ty+3.*ty1+3.*ty2+ty3)*dzd8;
405 ty4=ty+(kvy1+3.*kvy2+3.*kvy3+kvy4)*dzd8;
406 }
407 else
408 {
409 y4 =y+ty*dz;
410 ty4=ty;
411 }
412
413 x=x4;
414 tx=tx4;
415 y=y4;
416 ty=ty4;
417
418 z+=dz*fCrystalData->GetCorrectionZ();//motion along the z coordinate
419 //("central plane/axis", no current plane/axis)
420
421 xyzparticle = fCrystalData->ChannelChange(x,y,z);
422 x=xyzparticle.x();
423 y=xyzparticle.y();
424 z=xyzparticle.z();
425
426 momentumDirectionStep =
427 dz*std::sqrt(1+std::pow(std::tan(tx),2)+std::pow(std::tan(ty),2));
428 zalongGamma += dz/momentumDirectionZGamma;
429
430 //default scattering and energy loss 0
431 scatteringAnglesAndEnergyLoss.set(0.,0.,0.);
432
433 if(fIncoherentScattering)
434 {
435 //calculate separately for each element of the crystal
436 for (G4int ii = 0; ii < fCrystalData->GetNelements(); ii++)
437 {
438 //effective step taking into account nuclear density along the trajectory
439 effectiveStep = momentumDirectionStep*
440 fCrystalData->NuclearDensity(x,y,ii);
441 //Coulomb scattering on screened atomic potential
442 //(both multiple and single)
443 scatteringAnglesAndEnergyLoss +=
444 fCrystalData->CoulombAtomicScattering(effectiveStep,
445 momentumDirectionStep,
446 ii);
447 }
448 //electron scattering and coherent part of ionization energy losses
449 scatteringAnglesAndEnergyLoss += fCrystalData->CoulombElectronScattering(
450 fCrystalData->MinIonizationEnergy(x,y),
451 fCrystalData->ElectronDensity(x,y),
452 momentumDirectionStep);
453 tx += scatteringAnglesAndEnergyLoss.x();
454 ty += scatteringAnglesAndEnergyLoss.y();
455 }
456
457 //To avoid backward direction
458 if(std::abs(tx)>CLHEP::halfpi-DBL_EPSILON||
459 std::abs(ty)>CLHEP::halfpi-DBL_EPSILON)
460 {
461 G4cout << "Warning: particle angle is beyond +-pi/2 range => "
462 "skipping the calculation of its probability" << G4endl;
463 fbreak = true;
464 break;
465 }
466
467 //**********Baier-Katkov start
468
469 //back to the local reference system of the volume
470 tx0 = fCrystalData->AngleXFromLatticeToBox(tx,z);
471 ty0 = ty;
472
473 dzMeV=momentumDirectionStep/CLHEP::hbarc;// in MeV^-1
474
475 // accelerations
476 axt=(tx0-scatteringAnglesAndEnergyLoss.x()-txPreStep0)/dzMeV;
477 ayt=(ty0-scatteringAnglesAndEnergyLoss.y()-tyPreStep0)/dzMeV;
478
479 //the angles vs the photon (with incoherent scattering)
480 vxin = tx0-txGamma0;
481 vyin = ty0-tyGamma0;
482 //the angles vs the photon (without incoherent scattering)
483 vxno = vxin-scatteringAnglesAndEnergyLoss.x();
484 vyno = vyin-scatteringAnglesAndEnergyLoss.y();
485
486 //phase difference before scattering
487 faseBefore=omprimed2*(gammaInverse2+vxno*vxno+vyno*vyno);//phi' t<ti//MeV
488
489 faseBeforedz = faseBefore*dzMeV;
490 faseBeforedzd2 = faseBeforedz/2.;
491 fa+=faseBeforedz; //
492 fa1=fa-faseBeforedzd2;//
493 dzmod=2*std::sin(faseBeforedzd2)/faseBefore;//MeV^-1
494
495 //phi''/faseBefore^2
496 fa2dfaseBefore2 = omprime*(axt*vxno+ayt*vyno)/(faseBefore*faseBefore);
497
498 //phase difference after scattering
499 faseAfter=omprimed2*(gammaInverse2+vxin*vxin+vyin*vyin);//phi' ti+O//MeV
500
501 skJ=1/faseAfter-1/faseBefore-fa2dfaseBefore2*dzmod;//MeV^-1
502 skIx=vxin/faseAfter-vxno/faseBefore+dzmod*(axt/faseBefore-
503 vxno*fa2dfaseBefore2);
504 skIy=vyin/faseAfter-vyno/faseBefore+dzmod*(ayt/faseBefore-
505 vyno*fa2dfaseBefore2);
506
507 sinfa1 = std::sin(fa1);
508 cosfa1 = std::cos(fa1);
509
510 ss+=sinfa1*skJ;//sum sin integral J of BK
511 sc+=cosfa1*skJ;//sum cos integral J of BK
512 ssx+=sinfa1*skIx;// sum sin integral Ix of BK
513 ssy+=sinfa1*skIy;// sum sin integral Iy of BK
514 scx+=cosfa1*skIx;// sum cos integral Ix of BK
515 scy+=cosfa1*skIy;// sum cos integral Iy of BK
516 }
517
518 //only of the trajectory cycle was not broken
519 if(!fbreak)
520 {
521 G4double i2=ssx*ssx+scx*scx+ssy*ssy+scy*scy;//MeV^-2
522 G4double j2=ss*ss+sc*sc;//MeV^-2
523
524 probabilityPPdz += coefNorme2deprime2*angleNormCoef*
525 (i2*e2pluseprime2+j2*gammaInverse2om)/zalongGamma;
526 }
527 }
528
529 //filling the CDF of probabilities of the production of sampling pairs
530 fPairProductionCDFdz.push_back(fPairProductionCDFdz[i]+probabilityPPdz);
531 //**********Baier-Katkov end
532
533 //accumulation of initial parameters of sampling pairs
534 fullVectorEtotal.push_back(twoVectorEtotal);
535 fullVectorX.push_back(twoVectorX);
536 fullVectorY.push_back(twoVectorY);
537 fullVectorTX.push_back(twoVectorTX);
538 fullVectorTY.push_back(twoVectorTY);
539 }
540
541 //photon mean free path
542 //fPairProductionCDFdz.back() = full pair production probability
543 //simulated for the current photon along photon direction
544 G4double lMeanFreePath = 1/fPairProductionCDFdz.back();
545
546 fEffectiveLrad = 7.*lMeanFreePath/9.;//only for scoring purpose
547
548 return lMeanFreePath;
549 }
550 else
551 {
552 //dummy process, does not occur
553 return DBL_MAX;
554 }
555
556}
condition
G4double condition(const G4ErrorSymMatrix &m)
NotForced
@ NotForced
Definition G4ForceCondition.hh:47
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4LogicalVolume::GetRegion
G4Region * GetRegion() const
G4Region::GetName
const G4String & GetName() const
G4TouchableHistory::GetHistory
virtual const G4NavigationHistory * GetHistory() const
G4Track::GetVolume
G4VPhysicalVolume * GetVolume() const
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4Track::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const
G4Track::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4Track::GetTotalEnergy
G4double GetTotalEnergy() const
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
DBL_EPSILON
#define DBL_EPSILON
Definition templates.hh:66
DBL_MAX
#define DBL_MAX
Definition templates.hh:62

◆ GetNTrajectorySteps()

G4double G4CoherentPairProduction::GetNTrajectorySteps ( )
inline

get number of trajectory steps of a single particle (e- or e+)

Definition at line 99 of file G4CoherentPairProduction.hh.

99{return fNTrajectorySteps;}

◆ GetPPKineticEnergyCut()

G4double G4CoherentPairProduction::GetPPKineticEnergyCut ( )
inline

Definition at line 87 of file G4CoherentPairProduction.hh.

87{return fPPKineticEnergyCut;}

◆ GetSamplingPairsNumber()

G4int G4CoherentPairProduction::GetSamplingPairsNumber ( )
inline

get the number of pairs in sampling of Baier-Katkov Integral (MC integration by e+- energy and angles <=> e+- momentum)

Definition at line 91 of file G4CoherentPairProduction.hh.

91{return fNMCPairs;}

◆ Input() [1/3]

void G4CoherentPairProduction::Input ( const G4ChannelingFastSimCrystalData * crystalData)

Definition at line 674 of file G4CoherentPairProduction.cc.

675{
676 //setting the class with containing all
677 //the crystal material and crystal lattice data and
678 //Channeling scattering and ionization processes
679 //fCrystalData = new G4ChannelingFastSimCrystalData();
680
681 fCrystalData = const_cast<G4ChannelingFastSimCrystalData*>(crystalData);
682}

◆ Input() [2/3]

void G4CoherentPairProduction::Input ( const G4Material * crystal,
const G4String & lattice )
inline

special functions

Definition at line 66 of file G4CoherentPairProduction.hh.

68 {Input(crystal,lattice,"");}
G4CoherentPairProduction::Input
void Input(const G4Material *crystal, const G4String &lattice)
special functions
Definition G4CoherentPairProduction.hh:66

Referenced by G4CoherentPairProductionPhysics::ConstructProcess(), and Input().

◆ Input() [3/3]

void G4CoherentPairProduction::Input ( const G4Material * crystal,
const G4String & lattice,
const G4String & filePath )

Definition at line 660 of file G4CoherentPairProduction.cc.

663{
664 //initializing the class with containing all
665 //the crystal material and crystal lattice data and
666 //Channeling scattering and ionization processes
667 fCrystalData = new G4ChannelingFastSimCrystalData();
668 //setting all the crystal material and lattice data
669 fCrystalData->SetMaterialProperties(crystal,lattice,filePath);
670}

◆ IsApplicable()

G4bool G4CoherentPairProduction::IsApplicable ( const G4ParticleDefinition & aPD)
inlineoverridevirtual

Reimplemented from G4VProcess.

Definition at line 57 of file G4CoherentPairProduction.hh.

58 {
59 return(aPD.GetParticleName() == "gamma");
60 }
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96

◆ ModelMinPrimaryEnergy()

G4double G4CoherentPairProduction::ModelMinPrimaryEnergy ( )
inline

get cuts

Definition at line 85 of file G4CoherentPairProduction.hh.

85{ return fLowEnergyLimit;}

Referenced by GetMeanFreePath().

◆ PostStepDoIt()

G4VParticleChange * G4CoherentPairProduction::PostStepDoIt ( const G4Track & aTrack,
const G4Step & aStep )
overridevirtual

Reimplemented from G4VDiscreteProcess.

Definition at line 560 of file G4CoherentPairProduction.cc.

562{
563 //example with no physical sense
564 aParticleChange.Initialize(aTrack);
565 //G4LogicalVolume* aLV = aTrack.GetVolume()->GetLogicalVolume();
566
567 const G4ParticleDefinition* chargedParticleDefinition[2] =
568 {G4Electron::Electron(),G4Positron::Positron()};
569
570 // the coordinates of the photon in the local reference system of the volume
571 G4ThreeVector xyzGamma0 =
572 aTrack.GetTouchableHandle()->GetHistory()->
573 GetTopTransform().TransformPoint(aTrack.GetPosition());
574
575 // the coordinates of the photon in the co-rotating reference system within
576 //a channel (elementary periodic cell)
577 G4ThreeVector xyzGamma = fCrystalData->CoordinatesFromBoxToLattice(xyzGamma0);
578
579 //global time
580 G4double tGlobalGamma = aTrack.GetGlobalTime();
581
582 G4double ksi1 = G4UniformRand()*fPairProductionCDFdz.back();
583
584 //randomly choosing the pair to be produced from the sampling list
585 //according to the probabilities calculated in the Baier-Katkov integral
586 G4int ipair = FindVectorIndex(fPairProductionCDFdz,ksi1)-1;//index of
587 //a pair produced
588
589 // the coordinates of a charged particle in the reference system within
590 //a channel (elementary periodic cell)
591 G4ThreeVector xyzparticle;
592 //cycle e- and e+ within single pair
593 for(G4int j=0; j<2;j++)
594 {
595 xyzparticle.set(fullVectorX[ipair][j],fullVectorY[ipair][j],xyzGamma.z());
596
597 //in the local reference system of the volume
598 G4ThreeVector newParticleCoordinateXYZ =
599 fCrystalData->CoordinatesFromLatticeToBox(xyzparticle);
600 //the same in the global reference system
601 newParticleCoordinateXYZ =
602 aTrack.GetTouchableHandle()->GetHistory()->
603 GetTopTransform().Inverse().TransformPoint(newParticleCoordinateXYZ);
604
605 //back to the local reference system of the volume
606 G4double tx0 = fCrystalData->AngleXFromLatticeToBox(fullVectorTX[ipair][j],xyzGamma.z());
607 G4double ty0 = fullVectorTY[ipair][j];
608
609 G4double momentumDirectionZ = 1./
610 std::sqrt(1.+std::pow(std::tan(tx0),2)+
611 std::pow(std::tan(ty0),2));
612
613 //momentum direction vector of the charged particle produced
614 //in the local reference system of the volume
615 G4ThreeVector momentumDirectionParticle = G4ThreeVector(momentumDirectionZ*std::tan(tx0),
616 momentumDirectionZ*std::tan(ty0),
617 momentumDirectionZ);
618 //the same in the global reference system
619 momentumDirectionParticle =
620 (aTrack.GetTouchableHandle()->GetHistory()->GetTopTransform().NetRotation()) *
621 momentumDirectionParticle;
622
623 G4DynamicParticle* chargedParticle =
624 new G4DynamicParticle(chargedParticleDefinition[j],
625 momentumDirectionParticle,
626 fullVectorEtotal[ipair][j]-fMass);
627
628 // Create the track for the secondary particle
629 G4Track* secondaryTrack = new G4Track(chargedParticle,
630 tGlobalGamma,
631 newParticleCoordinateXYZ);
632 secondaryTrack->SetTouchableHandle(aStep.GetPostStepPoint()->GetTouchableHandle());
633 secondaryTrack->SetParentID(aTrack.GetTrackID());
634
635 //generation of a secondary charged particle
636 aParticleChange.AddSecondary(secondaryTrack);
637 }
638
639 //killing the photon
640 aParticleChange.ProposeTrackStatus(fStopAndKill);
641
642 return &aParticleChange;
643}
fStopAndKill
@ fStopAndKill
Definition G4TrackStatus.hh:46
G4Electron::Electron
static G4Electron * Electron()
Definition G4Electron.cc:91
G4Positron::Positron
static G4Positron * Positron()
Definition G4Positron.cc:90
G4StepPoint::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const
G4Step::GetPostStepPoint
G4StepPoint * GetPostStepPoint() const
G4Track::GetTrackID
G4int GetTrackID() const
G4Track::SetTouchableHandle
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4Track::GetGlobalTime
G4double GetGlobalTime() const
G4Track::SetParentID
void SetParentID(const G4int aValue)
G4VProcess::aParticleChange
G4ParticleChange aParticleChange
Definition G4VProcess.hh:331

◆ ProcessDescription()

void G4CoherentPairProduction::ProcessDescription ( std::ostream & out) const
overridevirtual

Reimplemented from G4VProcess.

Definition at line 686 of file G4CoherentPairProduction.cc.

687{
688 out << " Coherent pair production";
689 G4VDiscreteProcess::ProcessDescription(out);
690}
G4VProcess::ProcessDescription
virtual void ProcessDescription(std::ostream &outfile) const
Definition G4VProcess.cc:181

◆ SetChargeParticleAngleFactor()

void G4CoherentPairProduction::SetChargeParticleAngleFactor ( G4double chargeParticleAngleFactor)
inline

set the number of particle angles 1/gamma in pair production defining the width of the angular distribution of pair sampling in the Baier-Katkov Integral

Definition at line 121 of file G4CoherentPairProduction.hh.

122 {fChargeParticleAngleFactor = chargeParticleAngleFactor;}

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

◆ SetG4RegionName()

void G4CoherentPairProduction::SetG4RegionName ( const G4String & nameG4Region)
inline

set the name of G4Region in which the model is applicable

Definition at line 129 of file G4CoherentPairProduction.hh.

129{fG4RegionName=nameG4Region;}

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

◆ SetHighAngleLimit()

void G4CoherentPairProduction::SetHighAngleLimit ( G4double angle)
inline

Definition at line 111 of file G4CoherentPairProduction.hh.

111{fHighAngleLimit=angle;}

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

◆ SetLowEnergyLimit()

void G4CoherentPairProduction::SetLowEnergyLimit ( G4double energy)
inline

set cuts

Definition at line 110 of file G4CoherentPairProduction.hh.

110{fLowEnergyLimit=energy;}
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

◆ SetNTrajectorySteps()

void G4CoherentPairProduction::SetNTrajectorySteps ( G4int nTrajectorySteps)
inline

set number of trajectory steps of a single particle (e- or e+)

Definition at line 125 of file G4CoherentPairProduction.hh.

126 {fNTrajectorySteps = nTrajectorySteps;}

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

◆ SetPPKineticEnergyCut()

void G4CoherentPairProduction::SetPPKineticEnergyCut ( G4double kineticEnergyCut)
inline

Definition at line 112 of file G4CoherentPairProduction.hh.

112{fPPKineticEnergyCut=kineticEnergyCut;}

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

◆ SetSamplingPairsNumber()

void G4CoherentPairProduction::SetSamplingPairsNumber ( G4int nPairs)
inline

set the number of pairs in sampling of Baier-Katkov Integral (MC integration by e+- energy and angles <=> e+- momentum)

Definition at line 116 of file G4CoherentPairProduction.hh.

116{fNMCPairs = nPairs;}

Referenced by G4CoherentPairProductionPhysics::ConstructProcess().

The documentation for this class was generated from the following files: