| 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 //
27 //
28 // -------------------------------------------
29 // GEANT 4 class implementation file
30 // CERN Geneva Switzerland
31 //
32
33 // --------- G4LowEnergyPolarizedCompton class
34 //
35 // by G.Depaola & F.Longo (21 may 20
36 //
37 // 21 May 2001 - MGP Modified to inherit
38 // Applies same algorit
39 // if the incoming phot
40 // Temporary protection
41 // of polarisation || i
42 //
43 // 17 October 2001 - F.Longo - Revised accordi
44 //
45 // 21 February 2002 - F.Longo Revisions with A
46 // - better descrip
47 // - system of ref
48 // 22 January 2003 - V.Ivanchenko Cut per reg
49 // 24 April 2003 - V.Ivanchenko Cut per reg
50 //
51 //
52 // *******************************************
53 //
54 // Corrections by Rui Curado da Silva (2000)
55 // New Implementation by G.Depaola & F.Longo
56 //
57 // - sampling of phi
58 // - polarization of scattered photon
59 //
60 // -------------------------------------------
61
62 #include "G4LowEnergyPolarizedCompton.hh"
63 #include "Randomize.hh"
64 #include "G4PhysicalConstants.hh"
65 #include "G4SystemOfUnits.hh"
66 #include "G4ParticleDefinition.hh"
67 #include "G4Track.hh"
68 #include "G4Step.hh"
69 #include "G4ForceCondition.hh"
70 #include "G4Gamma.hh"
71 #include "G4Electron.hh"
72 #include "G4DynamicParticle.hh"
73 #include "G4VParticleChange.hh"
74 #include "G4ThreeVector.hh"
75 #include "G4RDVCrossSectionHandler.hh"
76 #include "G4RDCrossSectionHandler.hh"
77 #include "G4RDVEMDataSet.hh"
78 #include "G4RDCompositeEMDataSet.hh"
79 #include "G4RDVDataSetAlgorithm.hh"
80 #include "G4RDLogLogInterpolation.hh"
81 #include "G4RDVRangeTest.hh"
82 #include "G4RDRangeTest.hh"
83 #include "G4MaterialCutsCouple.hh"
84
85 // constructor
86
87 G4LowEnergyPolarizedCompton::G4LowEnergyPolari
88 : G4VDiscreteProcess(processName),
89 lowEnergyLimit (250*eV), // i
90 highEnergyLimit(100*GeV),
91 intrinsicLowEnergyLimit(10*eV),
92 intrinsicHighEnergyLimit(100*GeV)
93 {
94 if (lowEnergyLimit < intrinsicLowEnergyLimit
95 highEnergyLimit > intrinsicHighEnergyLim
96 {
97 G4Exception("G4LowEnergyPolarizedCompton
98 "OutOfRange", FatalException
99 "Energy outside intrinsic pr
100 }
101
102 crossSectionHandler = new G4RDCrossSectionHa
103
104
105 G4RDVDataSetAlgorithm* scatterInterpolation
106 G4String scatterFile = "comp/ce-sf-";
107 scatterFunctionData = new G4RDCompositeEMDat
108 scatterFunctionData->LoadData(scatterFile);
109
110 meanFreePathTable = 0;
111
112 rangeTest = new G4RDRangeTest;
113
114 // For Doppler broadening
115 shellData.SetOccupancyData();
116
117
118 if (verboseLevel > 0)
119 {
120 G4cout << GetProcessName() << " is crea
121 << "Energy range: "
122 << lowEnergyLimit / keV << " keV - "
123 << highEnergyLimit / GeV << " GeV"
124 << G4endl;
125 }
126 }
127
128 // destructor
129
130 G4LowEnergyPolarizedCompton::~G4LowEnergyPolar
131 {
132 delete meanFreePathTable;
133 delete crossSectionHandler;
134 delete scatterFunctionData;
135 delete rangeTest;
136 }
137
138
139 void G4LowEnergyPolarizedCompton::BuildPhysics
140 {
141
142 crossSectionHandler->Clear();
143 G4String crossSectionFile = "comp/ce-cs-";
144 crossSectionHandler->LoadData(crossSectionFi
145 delete meanFreePathTable;
146 meanFreePathTable = crossSectionHandler->Bui
147
148 // For Doppler broadening
149 G4String file = "/doppler/shell-doppler";
150 shellData.LoadData(file);
151
152 }
153
154 G4VParticleChange* G4LowEnergyPolarizedCompton
155 const G4Step& aStep)
156 {
157 // The scattered gamma energy is sampled acc
158 // The random number techniques of Butcher &
159 // GEANT4 internal units
160 //
161 // Note : Effects due to binding of atomic e
162
163 aParticleChange.Initialize(aTrack);
164
165 // Dynamic particle quantities
166 const G4DynamicParticle* incidentPhoton = aT
167 G4double gammaEnergy0 = incidentPhoton->GetK
168 G4ThreeVector gammaPolarization0 = incidentP
169
170 // gammaPolarization0 = gammaPolarization0.
171
172 // Protection: a polarisation parallel to th
173 // direction causes problems;
174 // in that case find a random polarization
175
176 G4ThreeVector gammaDirection0 = incidentPhot
177 // ---- MGP ---- Next two lines commented ou
178 // G4double scalarproduct = gammaPolarizatio
179 // G4double angle = gammaPolarization0.angle
180
181 // Make sure that the polarization vector is
182 // gamma direction. If not
183
184 if(!(gammaPolarization0.isOrthogonal(gammaDi
185 { // only for testing now
186 gammaPolarization0 = GetRandomPolarizati
187 }
188 else
189 {
190 if ( gammaPolarization0.howOrthogonal(ga
191 {
192 gammaPolarization0 = GetPerpendicularPolar
193 }
194 }
195
196 // End of Protection
197
198 // Within energy limit?
199
200 if(gammaEnergy0 <= lowEnergyLimit)
201 {
202 aParticleChange.ProposeTrackStatus(fStop
203 aParticleChange.ProposeEnergy(0.);
204 aParticleChange.ProposeLocalEnergyDeposi
205 return G4VDiscreteProcess::PostStepDoIt(
206 }
207
208 G4double E0_m = gammaEnergy0 / electron_mass
209
210 // Select randomly one element in the curren
211
212 const G4MaterialCutsCouple* couple = aTrack.
213 G4int Z = crossSectionHandler->SelectRandomA
214
215 // Sample the energy and the polarization of
216
217 G4double epsilon, epsilonSq, onecost, sinThe
218
219 G4double epsilon0 = 1./(1. + 2*E0_m);
220 G4double epsilon0Sq = epsilon0*epsilon0;
221 G4double alpha1 = - std::log(epsilon0);
222 G4double alpha2 = 0.5*(1.- epsilon0Sq);
223
224 G4double wlGamma = h_Planck*c_light/gammaEne
225 G4double gammaEnergy1;
226 G4ThreeVector gammaDirection1;
227
228 do {
229 if ( alpha1/(alpha1+alpha2) > G4UniformRan
230 {
231 epsilon = std::exp(-alpha1*G4UniformRand()
232 epsilonSq = epsilon*epsilon;
233 }
234 else
235 {
236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4
237 epsilon = std::sqrt(epsilonSq);
238 }
239
240 onecost = (1.- epsilon)/(epsilon*E0_m);
241 sinThetaSqr = onecost*(2.-onecost);
242
243 // Protection
244 if (sinThetaSqr > 1.)
245 {
246 if (verboseLevel>0) G4cout
247 << " -- Warning -- G4LowEnergyPolarizedCom
248 << "sin(theta)**2 = "
249 << sinThetaSqr
250 << "; set to 1"
251 << G4endl;
252 sinThetaSqr = 1.;
253 }
254 if (sinThetaSqr < 0.)
255 {
256 if (verboseLevel>0) G4cout
257 << " -- Warning -- G4LowEnergyPolarizedCom
258 << "sin(theta)**2 = "
259 << sinThetaSqr
260 << "; set to 0"
261 << G4endl;
262 sinThetaSqr = 0.;
263 }
264 // End protection
265
266 G4double x = std::sqrt(onecost/2.) / (wlG
267 G4double scatteringFunction = scatterFunct
268 greject = (1. - epsilon*sinThetaSqr/(1.+ e
269
270 } while(greject < G4UniformRand()*Z);
271
272
273 // *****************************************
274 // Phi determination
275 // *****************************************
276
277 G4double phi = SetPhi(epsilon,sinThetaSqr);
278
279 //
280 // scattered gamma angles. ( Z - axis along
281 //
282
283 G4double cosTheta = 1. - onecost;
284
285 // Protection
286
287 if (cosTheta > 1.)
288 {
289 if (verboseLevel>0) G4cout
290 << " -- Warning -- G4LowEnergyPolarizedCompt
291 << "cosTheta = "
292 << cosTheta
293 << "; set to 1"
294 << G4endl;
295 cosTheta = 1.;
296 }
297 if (cosTheta < -1.)
298 {
299 if (verboseLevel>0) G4cout
300 << " -- Warning -- G4LowEnergyPolarizedCompt
301 << "cosTheta = "
302 << cosTheta
303 << "; set to -1"
304 << G4endl;
305 cosTheta = -1.;
306 }
307 // End protection
308
309
310 G4double sinTheta = std::sqrt (sinThetaSqr);
311
312 // Protection
313 if (sinTheta > 1.)
314 {
315 if (verboseLevel>0) G4cout
316 << " -- Warning -- G4LowEnergyPolarizedCompt
317 << "sinTheta = "
318 << sinTheta
319 << "; set to 1"
320 << G4endl;
321 sinTheta = 1.;
322 }
323 if (sinTheta < -1.)
324 {
325 if (verboseLevel>0) G4cout
326 << " -- Warning -- G4LowEnergyPolarizedCompt
327 << "sinTheta = "
328 << sinTheta
329 << "; set to -1"
330 << G4endl;
331 sinTheta = -1.;
332 }
333 // End protection
334
335
336 G4double dirx = sinTheta*std::cos(phi);
337 G4double diry = sinTheta*std::sin(phi);
338 G4double dirz = cosTheta ;
339
340
341 // oneCosT , eom
342
343
344
345 // Doppler broadening - Method based on:
346 // Y. Namito, S. Ban and H. Hirayama,
347 // "Implementation of the Doppler Broadening
348 // NIM A 349, pp. 489-494, 1994
349
350 // Maximum number of sampling iterations
351
352 G4int maxDopplerIterations = 1000;
353 G4double bindingE = 0.;
354 G4double photonEoriginal = epsilon * gammaEn
355 G4double photonE = -1.;
356 G4int iteration = 0;
357 G4double eMax = gammaEnergy0;
358
359 do
360 {
361 iteration++;
362 // Select shell based on shell occupancy
363 G4int shell = shellData.SelectRandomShel
364 bindingE = shellData.BindingEnergy(Z,she
365
366 eMax = gammaEnergy0 - bindingE;
367
368 // Randomly sample bound electron moment
369 G4double pSample = profileData.RandomSel
370 // Rescale from atomic units
371 G4double pDoppler = pSample * fine_struc
372 G4double pDoppler2 = pDoppler * pDoppler
373 G4double var2 = 1. + onecost * E0_m;
374 G4double var3 = var2*var2 - pDoppler2;
375 G4double var4 = var2 - pDoppler2 * cosTh
376 G4double var = var4*var4 - var3 + pDoppl
377 if (var > 0.)
378 {
379 G4double varSqrt = std::sqrt(var);
380 G4double scale = gammaEnergy0 / var3;
381 // Random select either root
382 if (G4UniformRand() < 0.5) photonE = (var4
383 else photonE = (var4 + varSqrt) * scale;
384 }
385 else
386 {
387 photonE = -1.;
388 }
389 } while ( iteration <= maxDopplerIterations
390 (photonE < 0. || photonE > eMax || phot
391
392 // End of recalculation of photon energy wit
393 // Revert to original if maximum number of i
394 if (iteration >= maxDopplerIterations)
395 {
396 photonE = photonEoriginal;
397 bindingE = 0.;
398 }
399
400 gammaEnergy1 = photonE;
401
402 // G4cout << "--> PHOTONENERGY1 = " << photo
403
404
405 /// Doppler Broadeing
406
407
408
409
410 //
411 // update G4VParticleChange for the scattere
412 //
413
414 // gammaEnergy1 = epsilon*gammaEnergy0;
415
416
417 // New polarization
418
419 G4ThreeVector gammaPolarization1 = SetNewPol
420 sinThetaSqr,
421 phi,
422 cosTheta);
423
424 // Set new direction
425 G4ThreeVector tmpDirection1( dirx,diry,dirz
426 gammaDirection1 = tmpDirection1;
427
428 // Change reference frame.
429
430 SystemOfRefChange(gammaDirection0,gammaDirec
431 gammaPolarization0,gammaPolarization1)
432
433 if (gammaEnergy1 > 0.)
434 {
435 aParticleChange.ProposeEnergy( gammaEner
436 aParticleChange.ProposeMomentumDirection
437 aParticleChange.ProposePolarization( gam
438 }
439 else
440 {
441 aParticleChange.ProposeEnergy(0.) ;
442 aParticleChange.ProposeTrackStatus(fStop
443 }
444
445 //
446 // kinematic of the scattered electron
447 //
448
449 G4double ElecKineEnergy = gammaEnergy0 - gam
450
451
452 // Generate the electron only if with large
453
454 G4double safety = aStep.GetPostStepPoint()->
455
456
457 if (rangeTest->Escape(G4Electron::Electron()
458 {
459 G4double ElecMomentum = std::sqrt(ElecKi
460 G4ThreeVector ElecDirection((gammaEnergy
461 gammaEnergy1 * gammaDirection1) * (
462 G4DynamicParticle* electron = new G4Dyna
463 aParticleChange.SetNumberOfSecondaries(1
464 aParticleChange.AddSecondary(electron);
465 // aParticleChange.ProposeLocalEner
466 aParticleChange.ProposeLocalEnergyDeposi
467 }
468 else
469 {
470 aParticleChange.SetNumberOfSecondaries(0
471 aParticleChange.ProposeLocalEnergyDeposi
472 }
473
474 return G4VDiscreteProcess::PostStepDoIt( aTr
475
476 }
477
478
479 G4double G4LowEnergyPolarizedCompton::SetPhi(G
480 G4double sinSqrTh)
481 {
482 G4double rand1;
483 G4double rand2;
484 G4double phiProbability;
485 G4double phi;
486 G4double a, b;
487
488 do
489 {
490 rand1 = G4UniformRand();
491 rand2 = G4UniformRand();
492 phiProbability=0.;
493 phi = twopi*rand1;
494
495 a = 2*sinSqrTh;
496 b = energyRate + 1/energyRate;
497
498 phiProbability = 1 - (a/b)*(std::cos(phi
499
500
501
502 }
503 while ( rand2 > phiProbability );
504 return phi;
505 }
506
507
508 G4ThreeVector G4LowEnergyPolarizedCompton::Set
509 {
510 G4double dx = a.x();
511 G4double dy = a.y();
512 G4double dz = a.z();
513 G4double x = dx < 0.0 ? -dx : dx;
514 G4double y = dy < 0.0 ? -dy : dy;
515 G4double z = dz < 0.0 ? -dz : dz;
516 if (x < y) {
517 return x < z ? G4ThreeVector(-dy,dx,0) : G
518 }else{
519 return y < z ? G4ThreeVector(dz,0,-dx) : G
520 }
521 }
522
523 G4ThreeVector G4LowEnergyPolarizedCompton::Get
524 {
525 G4ThreeVector d0 = direction0.unit();
526 G4ThreeVector a1 = SetPerpendicularVector(d0
527 G4ThreeVector a0 = a1.unit(); // unit vector
528
529 G4double rand1 = G4UniformRand();
530
531 G4double angle = twopi*rand1; // random pola
532 G4ThreeVector b0 = d0.cross(a0); // cross pr
533
534 G4ThreeVector c;
535
536 c.setX(std::cos(angle)*(a0.x())+std::sin(ang
537 c.setY(std::cos(angle)*(a0.y())+std::sin(ang
538 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang
539
540 G4ThreeVector c0 = c.unit();
541
542 return c0;
543
544 }
545
546
547 G4ThreeVector G4LowEnergyPolarizedCompton::Get
548 (const G4ThreeVector& gammaDirection, const G4
549 {
550
551 //
552 // The polarization of a photon is always pe
553 // Therefore this function removes those vec
554 // points in direction of gammaDirection
555 //
556 // Mathematically we search the projection o
557 // plains normal vector.
558 // The basic equation can be found in each g
559 // p = a - (a o n)/(n o n)*n
560
561 return gammaPolarization - gammaPolarization
562 }
563
564
565 G4ThreeVector G4LowEnergyPolarizedCompton::Set
566 G4double sinSqrTh,
567 G4double phi,
568 G4double costheta)
569 {
570 G4double rand1;
571 G4double rand2;
572 G4double cosPhi = std::cos(phi);
573 G4double sinPhi = std::sin(phi);
574 G4double sinTheta = std::sqrt(sinSqrTh);
575 G4double cosSqrPhi = cosPhi*cosPhi;
576 // G4double cossqrth = 1.-sinSqrTh;
577 // G4double sinsqrphi = sinPhi*sinPhi;
578 G4double normalisation = std::sqrt(1. - cosS
579
580
581 // Determination of Theta
582
583 // ---- MGP ---- Commented out the following
584 // warnings (unused variables)
585 // G4double thetaProbability;
586 G4double theta;
587 // G4double a, b;
588 // G4double cosTheta;
589
590 /*
591
592 depaola method
593
594 do
595 {
596 rand1 = G4UniformRand();
597 rand2 = G4UniformRand();
598 thetaProbability=0.;
599 theta = twopi*rand1;
600 a = 4*normalisation*normalisation;
601 b = (epsilon + 1/epsilon) - 2;
602 thetaProbability = (b + a*std::cos(theta
603 cosTheta = std::cos(theta);
604 }
605 while ( rand2 > thetaProbability );
606
607 G4double cosBeta = cosTheta;
608
609 */
610
611
612 // Dan Xu method (IEEE TNS, 52, 1160 (2005))
613
614 rand1 = G4UniformRand();
615 rand2 = G4UniformRand();
616
617 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi
618 {
619 if (rand2<0.5)
620 theta = pi/2.0;
621 else
622 theta = 3.0*pi/2.0;
623 }
624 else
625 {
626 if (rand2<0.5)
627 theta = 0;
628 else
629 theta = pi;
630 }
631 G4double cosBeta = std::cos(theta);
632 G4double sinBeta = std::sqrt(1-cosBeta*cosBe
633
634 G4ThreeVector gammaPolarization1;
635
636 G4double xParallel = normalisation*cosBeta;
637 G4double yParallel = -(sinSqrTh*cosPhi*sinPh
638 G4double zParallel = -(costheta*sinTheta*cos
639 G4double xPerpendicular = 0.;
640 G4double yPerpendicular = (costheta)*sinBeta
641 G4double zPerpendicular = -(sinTheta*sinPhi)
642
643 G4double xTotal = (xParallel + xPerpendicula
644 G4double yTotal = (yParallel + yPerpendicula
645 G4double zTotal = (zParallel + zPerpendicula
646
647 gammaPolarization1.setX(xTotal);
648 gammaPolarization1.setY(yTotal);
649 gammaPolarization1.setZ(zTotal);
650
651 return gammaPolarization1;
652
653 }
654
655
656 void G4LowEnergyPolarizedCompton::SystemOfRefC
657 G4ThreeVector& direction1,
658 G4ThreeVector& polarization0,
659 G4ThreeVector& polarization1)
660 {
661 // direction0 is the original photon directi
662 // polarization0 is the original photon pola
663 // need to specify y axis in the real refere
664 G4ThreeVector Axis_Z0 = direction0.unit();
665 G4ThreeVector Axis_X0 = polarization0.unit()
666 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_
667
668 G4double direction_x = direction1.getX();
669 G4double direction_y = direction1.getY();
670 G4double direction_z = direction1.getZ();
671
672 direction1 = (direction_x*Axis_X0 + directio
673 G4double polarization_x = polarization1.getX
674 G4double polarization_y = polarization1.getY
675 G4double polarization_z = polarization1.getZ
676
677 polarization1 = (polarization_x*Axis_X0 + po
678
679 }
680
681
682 G4bool G4LowEnergyPolarizedCompton::IsApplicab
683 {
684 return ( &particle == G4Gamma::Gamma() );
685 }
686
687
688 G4double G4LowEnergyPolarizedCompton::GetMeanF
689 G4double,
690 G4ForceCondition*)
691 {
692 const G4DynamicParticle* photon = track.GetD
693 G4double energy = photon->GetKineticEnergy()
694 const G4MaterialCutsCouple* couple = track.G
695 size_t materialIndex = couple->GetIndex();
696 G4double meanFreePath;
697 if (energy > highEnergyLimit) meanFreePath =
698 else if (energy < lowEnergyLimit) meanFreePa
699 else meanFreePath = meanFreePathTable->FindV
700 return meanFreePath;
701 }
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716