| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Author: Alexei Sytov
27 // Co-author: Gianfranco Paterno (testing)
28 // Using the key points of G4BaierKatkov and d
29 // partially described in L. Bandiera et al. E
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........oo
49
50 G4CoherentPairProduction::G4CoherentPairProduc
51
52 G4VDiscreteProcess(aName)
53 {
54 SetProcessSubType(fCoherentPairProduction)
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 G4double G4CoherentPairProduction::GetMeanFree
60 G4double,
61 G4ForceCon
62 {
63 //current logical volume
64 G4LogicalVolume* crystallogic;
65
66 //momentum direction and coordinates (see
67 G4ThreeVector momentumDirectionGamma,xyzGa
68 //angle of the photon in the local referen
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 che
80 if(eGamma > ModelMinPrimaryEnergy())
81 {
82 //current logical volume
83 crystallogic = aTrack.GetVolume()->Get
84
85 //the model works only in the G4Region
86 if(crystallogic->GetRegion()->GetName(
87 {
88 fCrystalData->SetGeometryParameter
89
90 //the momentum direction of the ph
91 momentumDirectionGamma =
92 (aTrack.GetTouchableHandle()->
93 GetTopTransform().NetRotat
94
95 //the coordinates of the photon in
96 xyzGamma0 =
97 aTrack.GetTouchableHandle()->G
98 GetTopTransform().Transfor
99
100 // the coordinates of the photon i
101 //a channel (elementary periodic c
102 xyzGamma = fCrystalData->Coordinat
103
104 //angle of the photon in the local
105 //(!!! ONLY FORWARD DIRECTION, mom
106 txGamma0 = std::atan(momentumDirec
107 tyGamma0 = std::atan(momentumDirec
108
109 //recalculate angle into the latti
110 G4double angle = fCrystalData->Ang
111 if (fCrystalData->GetModel()==2)
112 {
113 angle = std::sqrt(angle*angle+
114 }
115
116 //Applies the parameterisation not
117 //above low energy limit and below
118 modelTrigger = (momentumDirectionG
119 std::abs(angle) <
120 }
121 }
122
123 if(modelTrigger)
124 {
125 //execute the model
126
127 G4double x=0.,y=0.,z=0.;// the coordin
128 //in the co-rotating reference sys
129 //a channel (elementary periodic c
130 G4double tx0=0.,ty0=0.; // the angles
131 // in the local reference system o
132 G4double txPreStep0=0.,tyPreStep0=0.;
133 // in the co-rotating reference sy
134 //a channel (elementary periodic c
135
136 G4ThreeVector scatteringAnglesAndEnerg
137
138 //coordinates in Runge-Kutta calculati
139 G4double x1=0.,x2=0.,x3=0.,x4=0.,y1=0.
140 //angles in Runge-Kutta calculations
141 G4double tx1=0.,tx2=0.,tx3=0.,tx4=0.,t
142 //variables in Runge-Kutta calculation
143 G4double kvx1=0.,kvx2=0.,kvx3=0.,kvx4=
144 //simulation step along z (internal st
145 G4double dz=0.,dzd3=0.,dzd8=0.;//dzd3
146 //simulation step along the momentum d
147 G4double momentumDirectionStep;
148 //effective simulation step (taking in
149 G4double effectiveStep=0.;
150
151 // Baier-Katkov variables
152 G4double dzMeV=0.; //step in MeV^-1
153 G4double axt=0.,ayt=0.; //charged part
154 G4double vxin=0.,vyin=0.;//the angles
155 G4double vxno=0.,vyno=0.;//the angles
156
157 G4double dzmod=0.;
158 G4double fa1=0.,faseBefore=0.,faseBefo
159 G4double faseAfter=0.,fa2dfaseBefore2=
160
161 G4double skJ=0, skIx=0., skIy=0.;
162 G4double sinfa1=0.,cosfa1=0.;
163
164 //2-vector is needed for an initial pa
165 //vector of 2-vectors is an initial pa
166
167 //collection of etotal for a single pa
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.);//0
189
190 const G4double charge[2] = {-1.,1.}; /
191 const G4String particleName[2] = {"e-"
192
193 // the coordinates of a charged partic
194 //a channel (elementary periodic cell)
195 G4ThreeVector xyzparticle = xyzGamma;/
196
197 //the idea of pair production simulati
198 //since the matrix element of these pr
199 //to radiation: sample the pairs, calc
200 //probabilities using Baier-Katkov ana
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 + fPPKinet
207 G4UniformRand()*
208
209
210 G4double phi = CLHEP::twopi*G4Unif
211
212 //the probability of the productio
213 //per distance
214 G4double probabilityPPdz = 0.;
215
216 //cycle e- and e+ within single pa
217 for(G4int j=0; j<2;j++)
218 {
219 if(j==1){etotal=eGamma-etotal;
220 twoVectorEtotal[j]=etotal;
221
222 //Baier-Katkov input
223 //intermediate variables to re
224 G4double e2 = etotal*etotal;
225 G4double gammaInverse2 = fMass
226 //normalization coefficient
227 G4double coefNorm = CLHEP::fin
228 //G4double phi = CLHEP::twopi*
229 G4double om = eGamma;
230 G4double eprime=om-etotal; //E
231 G4double eprime2 = eprime*epri
232 G4double e2pluseprime2 =e2+epr
233 G4double omprime=etotal*om/epr
234 G4double omprimed2=omprime/2;
235
236 //difference vs G4BaierKatkov:
237 G4double coefNorme2deprime2 =
238
239 G4double gammaInverse2om = gam
240
241 //initialize intermediate inte
242 G4double fa=0.,ss=0.,sc=0.,ssx
243
244 //End of Baier-Katkov input
245
246 G4bool fbreak = false;//flag o
247
248 //set fCrystalData parameters
249 fCrystalData->SetParticlePrope
250
251
252 //needed just to setup the cor
253 //since later it may be change
254 fCrystalData->CoordinatesFromB
255
256 //coordinate sampling: random
257 //in the interaction point
258 if(j==0)
259 {
260 x = fCrystalData->GetChann
261 y = fCrystalData->GetChann
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 o
271 //interaction point is taking
272 //of the position of pair prod
273 z = xyzGamma.z();
274
275 //angles of the photon in the
276 //angular distribution center
277 G4double tx = fCrystalData->An
278 G4double ty = tyGamma0;
279 G4double momentumDirectionZGam
280
281
282
283 //angle sampling: depends on a
284 //defined by the Lindhard angl
285 //to 1/gamma
286
287 //range of MC integration on a
288 G4double paramParticleAngle =
289
290 G4double axangle=0.;
291 if (fCrystalData->GetModel()==
292 {
293 axangle = std::abs(tx);
294 }
295 else if (fCrystalData->GetMod
296 {
297 axangle = std::sqrt(tx*tx+
298 }
299
300 if(axangle>fCrystalData->GetLi
301 {
302 paramParticleAngle+=axangl
303 -std
304
305
306 }
307 else
308 {
309 paramParticleAngle+=fCryst
310 }
311
312
313 //ONLY forward direction
314 if (paramParticleAngle>CLHEP::
315
316 G4double rho=1.;
317 G4double rhocut=CLHEP::halfpi/
318
319 G4double norm=std::atan(rhocut
320 CLHEP::pi*para
321
322
323 //distribution with long tails
324 //after a strong single scatte
325 //at ellipsescale < 1 => half
326 do
327 {
328 rho = std::sqrt(std::tan(C
329 }
330 while (rho>rhocut);
331
332 //normalization coefficient fo
333 G4double angleNormCoef = (1.+r
334
335 tx+=charge[j]*paramParticleAng
336 twoVectorTX[j] = tx;
337 ty+=charge[j]*paramParticleAng
338 twoVectorTY[j] = ty;
339
340 G4double zalongGamma = 0;//nec
341 //depending on the traject
342 //starting the trajectory
343 //here we don't care about the
344 //the trajectory is very short
345 //in Baier-Katkov will be extr
346 for(G4int k=0; k<fNTrajectoryS
347 {
348 //back to the local refere
349 txPreStep0 = fCrystalData-
350 tyPreStep0 = ty;
351
352 dz = fCrystalData->GetSimu
353 dzd3=dz/3;
354 dzd8=dz/8;
355
356 //trajectory calculation:
357 //Runge-Cutt "3/8"
358 //fCrystalData->GetCurv(z)
359 //of the radius on x; GetC
360
361 //first step
362 kvx1=fCrystalData->Ex(x,y)
363 x1=x+tx*dzd3;
364 tx1=tx+(kvx1-fCrystalData-
365 if (fCrystalData->GetModel
366 {
367 kvy1=fCrystalData->Ey(
368 y1=y+ty*dzd3;
369 ty1=ty+kvy1*dzd3;
370 }
371
372 //second step
373 kvx2=fCrystalData->Ex(x1,y
374 x2=x-tx*dzd3+tx1*dz;
375 tx2=tx-(kvx1-fCrystalData-
376 (kvx2-fCrystalData->
377 if (fCrystalData->GetModel
378 {
379 kvy2=fCrystalData->Ey(
380 y2=y-ty*dzd3+ty1*dz;
381 ty2=ty-kvy1*dzd3+kvy2*
382 }
383
384 //third step
385 kvx3=fCrystalData->Ex(x2,y
386 x3=x+(tx-tx1+tx2)*dz;
387 tx3=tx+(kvx1-kvx2+kvx3-
388 fCrystalData->
389 if (fCrystalData->GetModel
390 {
391 kvy3=fCrystalData->Ey(
392 y3=y+(ty-ty1+ty2)*dz;
393 ty3=ty+(kvy1-kvy2+kvy3
394 }
395
396 //fourth step
397 kvx4=fCrystalData->Ex(x3,y
398 x4=x+(tx+3.*tx1+3.*tx2+tx3
399 tx4=tx+(kvx1+3.*kvx2+3.*kv
400 fCrystalData->GetCur
401 if (fCrystalData->GetModel
402 {
403 kvy4=fCrystalData->Ey(
404 y4=y+(ty+3.*ty1+3.*ty2
405 ty4=ty+(kvy1+3.*kvy2+3
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->GetCor
419 //("central plane/axis
420
421 xyzparticle = fCrystalData
422 x=xyzparticle.x();
423 y=xyzparticle.y();
424 z=xyzparticle.z();
425
426 momentumDirectionStep =
427 dz*std::sqrt(1+std::po
428 zalongGamma += dz/momentum
429
430 //default scattering and e
431 scatteringAnglesAndEnergyL
432
433 if(fIncoherentScattering)
434 {
435 //calculate separately
436 for (G4int ii = 0; ii
437 {
438 //effective step ta
439 effectiveStep = mo
440 fC
441 //Coulomb scatteri
442 //(both multiple a
443 scatteringAnglesAn
444 fCrystalData->
445
446
447 }
448 //electron scattering
449 scatteringAnglesAndEne
450 fCrystalData->MinI
451 fCrystalData->Elec
452 momentumDirectionS
453 tx += scatteringAngles
454 ty += scatteringAngles
455 }
456
457 //To avoid backward direct
458 if(std::abs(tx)>CLHEP::hal
459 std::abs(ty)>CLHEP::ha
460 {
461 G4cout << "Warning: pa
462 "skipping th
463 fbreak = true;
464 break;
465 }
466
467 //**********Baier-Katkov s
468
469 //back to the local refere
470 tx0 = fCrystalData->AngleX
471 ty0 = ty;
472
473 dzMeV=momentumDirectionSte
474
475 // accelerations
476 axt=(tx0-scatteringAnglesA
477 ayt=(ty0-scatteringAnglesA
478
479 //the angles vs the photon
480 vxin = tx0-txGamma0;
481 vyin = ty0-tyGamma0;
482 //the angles vs the photon
483 vxno = vxin-scatteringAngl
484 vyno = vyin-scatteringAngl
485
486 //phase difference before
487 faseBefore=omprimed2*(gamm
488
489 faseBeforedz = faseBefore*
490 faseBeforedzd2 = faseBefor
491 fa+=faseBeforedz; //
492 fa1=fa-faseBeforedzd2;//
493 dzmod=2*std::sin(faseBefor
494
495 //phi''/faseBefore^2
496 fa2dfaseBefore2 = omprime*
497
498 //phase difference after s
499 faseAfter=omprimed2*(gamma
500
501 skJ=1/faseAfter-1/faseBefo
502 skIx=vxin/faseAfter-vxno/f
503
504 skIy=vyin/faseAfter-vyno/f
505
506
507 sinfa1 = std::sin(fa1);
508 cosfa1 = std::cos(fa1);
509
510 ss+=sinfa1*skJ;//sum sin i
511 sc+=cosfa1*skJ;//sum cos i
512 ssx+=sinfa1*skIx;// sum si
513 ssy+=sinfa1*skIy;// sum si
514 scx+=cosfa1*skIx;// sum co
515 scy+=cosfa1*skIy;// sum co
516 }
517
518 //only of the trajectory cycle
519 if(!fbreak)
520 {
521 G4double i2=ssx*ssx+scx*sc
522 G4double j2=ss*ss+sc*sc;//
523
524 probabilityPPdz += coefNor
525 (i2*e2p
526 }
527 }
528
529 //filling the CDF of probabilities
530 fPairProductionCDFdz.push_back(fPa
531 //**********Baier-Katkov end
532
533 //accumulation of initial paramete
534 fullVectorEtotal.push_back(twoVect
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 p
543 //simulated for the current photon alo
544 G4double lMeanFreePath = 1/fPairProduc
545
546 fEffectiveLrad = 7.*lMeanFreePath/9.;/
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........oo
559
560 G4VParticleChange* G4CoherentPairProduction::P
561 con
562 {
563 //example with no physical sense
564 aParticleChange.Initialize(aTrack);
565 //G4LogicalVolume* aLV = aTrack.GetVolume(
566
567 const G4ParticleDefinition* chargedParticl
568 {G4Electron::Electron(),G4Positron::Po
569
570 // the coordinates of the photon in the lo
571 G4ThreeVector xyzGamma0 =
572 aTrack.GetTouchableHandle()->GetHistor
573 GetTopTransform().TransformPoint(a
574
575 // the coordinates of the photon in the co
576 //a channel (elementary periodic cell)
577 G4ThreeVector xyzGamma = fCrystalData->Coo
578
579 //global time
580 G4double tGlobalGamma = aTrack.GetGlobalTi
581
582 G4double ksi1 = G4UniformRand()*fPairProdu
583
584 //randomly choosing the pair to be produce
585 //according to the probabilities calculate
586 G4int ipair = FindVectorIndex(fPairProduct
587 //a pair produced
588
589 // the coordinates of a charged particle i
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],
596
597 //in the local reference system of the
598 G4ThreeVector newParticleCoordinateXYZ
599 fCrystalData->CoordinatesFromL
600 //the same in the global reference sys
601 newParticleCoordinateXYZ =
602 aTrack.GetTouchableHandle()->GetHi
603 GetTopTransform().Inverse().Tr
604
605 //back to the local reference system o
606 G4double tx0 = fCrystalData->AngleXFro
607 G4double ty0 = fullVectorTY[ipair][j];
608
609 G4double momentumDirectionZ = 1./
610 std::sqr
611
612
613 //momentum direction vector of the cha
614 //in the local reference system of the
615 G4ThreeVector momentumDirectionParticl
616
617
618 //the same in the global reference sys
619 momentumDirectionParticle =
620 (aTrack.GetTouchableHandle()->GetH
621 momentumDirectionParticle;
622
623 G4DynamicParticle* chargedParticle =
624 new G4DynamicParticle(chargedParti
625 momentumDirectio
626 fullVectorEtotal
627
628 // Create the track for the secondary
629 G4Track* secondaryTrack = new G4Track(
630
631
632 secondaryTrack->SetTouchableHandle(aSt
633 secondaryTrack->SetParentID(aTrack.Get
634
635 //generation of a secondary charged pa
636 aParticleChange.AddSecondary(secondary
637 }
638
639 //killing the photon
640 aParticleChange.ProposeTrackStatus(fStopAn
641
642 return &aParticleChange;
643 }
644
645 //....oooOO0OOooo........oooOO0OOooo........oo
646
647 G4int G4CoherentPairProduction::FindVectorInde
648 {
649 auto iteratorbegin = myvector.begin();
650 auto iteratorend = myvector.end();
651
652 //vector index (for non precise values low
653 auto loweriterator = std::lower_bound(iter
654 //return the index of the vector element
655 return (G4int)std::distance(iteratorbegin,
656 }
657
658 //....oooOO0OOooo........oooOO0OOooo........oo
659
660 void G4CoherentPairProduction::Input(const G4M
661 const G4S
662 const G4S
663 {
664 //initializing the class with containing a
665 //the crystal material and crystal lattice
666 //Channeling scattering and ionization pro
667 fCrystalData = new G4ChannelingFastSimCrys
668 //setting all the crystal material and lat
669 fCrystalData->SetMaterialProperties(crysta
670 }
671
672 //....oooOO0OOooo........oooOO0OOooo........oo
673
674 void G4CoherentPairProduction::Input(const G4C
675 {
676 //setting the class with containing all
677 //the crystal material and crystal lattice
678 //Channeling scattering and ionization pro
679 //fCrystalData = new G4ChannelingFastSimCr
680
681 fCrystalData = const_cast<G4ChannelingFast
682 }
683
684 //....oooOO0OOooo........oooOO0OOooo........oo
685
686 void G4CoherentPairProduction::ProcessDescript
687 {
688 out << " Coherent pair production";
689 G4VDiscreteProcess::ProcessDescription(out
690 }
691
692 //....oooOO0OOooo........oooOO0OOooo........oo
693