| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Author: Alexei Sytov
27 // Co-author: Gianfranco Paterno (testing)
28 // Using the key points of G4BaierKatkov and developments of V.V. Tikhomirov,
29 // partially described in L. Bandiera et al. Eur. Phys. J. C 82, 699 (2022)
30
31 #include "G4CoherentPairProduction.hh"
32
33 #include "Randomize.hh"
34 #include "G4TouchableHistory.hh"
35 #include "G4TouchableHandle.hh"
36 #include "G4SystemOfUnits.hh"
37
38 #include "G4Track.hh"
39 #include "G4Gamma.hh"
40 #include "G4Electron.hh"
41 #include "G4Positron.hh"
42
43 #include "G4ParticleDefinition.hh"
44 #include "G4ProcessManager.hh"
45 #include "G4EmProcessSubType.hh"
46 #include "G4TransportationManager.hh"
47
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50 G4CoherentPairProduction::G4CoherentPairProduction(const G4String& aName,
51 G4ProcessType):
52 G4VDiscreteProcess(aName)
53 {
54 SetProcessSubType(fCoherentPairProduction);
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58
59 G4double G4CoherentPairProduction::GetMeanFreePath(const G4Track& aTrack,
60 G4double,
61 G4ForceCondition* condition)
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 }
557
558 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
559
560 G4VParticleChange* G4CoherentPairProduction::PostStepDoIt(const G4Track& aTrack,
561 const G4Step& aStep)
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 }
644
645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646
647 G4int G4CoherentPairProduction::FindVectorIndex(std::vector<G4double> &myvector, G4double value)
648 {
649 auto iteratorbegin = myvector.begin();
650 auto iteratorend = myvector.end();
651
652 //vector index (for non precise values lower_bound gives upper value)
653 auto loweriterator = std::lower_bound(iteratorbegin, iteratorend, value);
654 //return the index of the vector element
655 return (G4int)std::distance(iteratorbegin, loweriterator);
656 }
657
658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
659
660 void G4CoherentPairProduction::Input(const G4Material *crystal,
661 const G4String &lattice,
662 const G4String &filePath)
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 }
671
672 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
673
674 void G4CoherentPairProduction::Input(const G4ChannelingFastSimCrystalData *crystalData)
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 }
683
684 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
685
686 void G4CoherentPairProduction::ProcessDescription(std::ostream& out) const
687 {
688 out << " Coherent pair production";
689 G4VDiscreteProcess::ProcessDescription(out);
690 }
691
692 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
693