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

Last change on this file since 1058 was 1055, checked in by garnier, 17 years ago

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