source: trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyPolarizedCompton.cc@ 1014

Last change on this file since 1014 was 1007, checked in by garnier, 17 years ago

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