source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.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.7 KB
RevLine 
[968]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//
[1007]26// $Id: G4LivermorePolarizedComptonModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-02 $
[968]28//
29
30#include "G4LivermorePolarizedComptonModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
39 const G4String& nam)
[1007]40:G4VEmModel(nam),isInitialised(false)
[968]41{
42 lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
43 highEnergyLimit = 100 * GeV;
44 SetLowEnergyLimit(lowEnergyLimit);
45 SetHighEnergyLimit(highEnergyLimit);
46
47 verboseLevel= 0;
48 // Verbosity scale:
49 // 0 = nothing
50 // 1 = warning for energy non-conservation
51 // 2 = details of energy budget
52 // 3 = calculation of cross sections, file openings, sampling of atoms
53 // 4 = entering in methods
54
55 G4cout << "Livermore Polarized Compton is constructed " << G4endl
56 << "Energy range: "
57 << lowEnergyLimit / keV << " keV - "
58 << highEnergyLimit / GeV << " GeV"
59 << G4endl;
60
61}
62
63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64
65G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
66{
[1007]67 delete meanFreePathTable;
68 delete crossSectionHandler;
69 delete scatterFunctionData;
[968]70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
75 const G4DataVector& cuts)
76{
77 if (verboseLevel > 3)
78 G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
79
[1007]80 InitialiseElementSelectors(particle,cuts);
[968]81
82 // Energy limits
83
84 if (LowEnergyLimit() < lowEnergyLimit)
85 {
86 G4cout << "G4LivermorePolarizedComptonModel: low energy limit increased from " <<
87 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
88 SetLowEnergyLimit(lowEnergyLimit);
89 }
90
91 if (HighEnergyLimit() > highEnergyLimit)
92 {
93 G4cout << "G4LivermorePolarizedComptonModel: high energy limit decreased from " <<
94 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
95 SetHighEnergyLimit(highEnergyLimit);
96 }
97
98 // Reading of data files - all materials are read
99
100 crossSectionHandler = new G4CrossSectionHandler;
101 crossSectionHandler->Clear();
102 G4String crossSectionFile = "comp/ce-cs-";
103 crossSectionHandler->LoadData(crossSectionFile);
104
105 meanFreePathTable = 0;
106 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
107
108 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
109 G4String scatterFile = "comp/ce-sf-";
110 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
111 scatterFunctionData->LoadData(scatterFile);
112
113 // For Doppler broadening
114 shellData.SetOccupancyData();
115 G4String file = "/doppler/shell-doppler";
116 shellData.LoadData(file);
117
118 //
119 if (verboseLevel > 2)
120 G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
121
122 G4cout << "Livermore Polarized Compton model is initialized " << G4endl
123 << "Energy range: "
124 << LowEnergyLimit() / keV << " keV - "
125 << HighEnergyLimit() / GeV << " GeV"
126 << G4endl;
127
128 //
129
130 if(isInitialised) return;
131
132 if(pParticleChange)
133 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
134 else
135 fParticleChange = new G4ParticleChangeForGamma();
136
137 isInitialised = true;
138}
139
140//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141
142G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
143 const G4ParticleDefinition*,
144 G4double GammaEnergy,
145 G4double Z, G4double,
146 G4double, G4double)
147{
148 if (verboseLevel > 3)
149 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
150
151 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
152 return cs;
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156
157void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
158 const G4MaterialCutsCouple* couple,
159 const G4DynamicParticle* aDynamicGamma,
160 G4double,
161 G4double)
162{
163 // The scattered gamma energy is sampled according to Klein - Nishina formula.
164 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
165 // GEANT4 internal units
166 //
167 // Note : Effects due to binding of atomic electrons are negliged.
168
169 if (verboseLevel > 3)
170 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
171
172 G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
173 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
174
175 // Protection: a polarisation parallel to the
176 // direction causes problems;
177 // in that case find a random polarization
178
179 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
180
181 // Make sure that the polarization vector is perpendicular to the
182 // gamma direction. If not
183
184 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
185 { // only for testing now
186 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
187 }
188 else
189 {
190 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
191 {
192 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
193 }
194 }
195
196 // End of Protection
197
198 // Within energy limit?
199
200 if(gammaEnergy0 <= lowEnergyLimit)
201 {
202 fParticleChange->ProposeTrackStatus(fStopAndKill);
203 fParticleChange->SetProposedKineticEnergy(0.);
204 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
[1007]205 // SI - IS THE FOLLOWING RETURN NECESSARY ?
[968]206 return;
207 }
208
209 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
210
211 // Select randomly one element in the current material
212
213 G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
214
215 // Sample the energy and the polarization of the scattered photon
216
217 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
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/gammaEnergy0;
225 G4double gammaEnergy1;
226 G4ThreeVector gammaDirection1;
227
228 do {
229 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
230 {
231 epsilon = std::exp(-alpha1*G4UniformRand());
232 epsilonSq = epsilon*epsilon;
233 }
234 else
235 {
236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
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 G4cout
247 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
248 << "sin(theta)**2 = "
249 << sinThetaSqr
250 << "; set to 1"
251 << G4endl;
252 sinThetaSqr = 1.;
253 }
254 if (sinThetaSqr < 0.)
255 {
256 G4cout
257 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
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.) / (wlGamma/cm);;
267 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
268 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
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 the parent gamma)
281 //
282
283 G4double cosTheta = 1. - onecost;
284
285 // Protection
286
287 if (cosTheta > 1.)
288 {
289 G4cout
290 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
291 << "cosTheta = "
292 << cosTheta
293 << "; set to 1"
294 << G4endl;
295 cosTheta = 1.;
296 }
297 if (cosTheta < -1.)
298 {
299 G4cout
300 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
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 G4cout
316 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
317 << "sinTheta = "
318 << sinTheta
319 << "; set to 1"
320 << G4endl;
321 sinTheta = 1.;
322 }
323 if (sinTheta < -1.)
324 {
325 G4cout
326 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
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 // 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 //
401 // update G4VParticleChange for the scattered photon
402 //
403
404 // gammaEnergy1 = epsilon*gammaEnergy0;
405
406
407 // New polarization
408
409 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
410 sinThetaSqr,
411 phi,
412 cosTheta);
413
414 // Set new direction
415 G4ThreeVector tmpDirection1( dirx,diry,dirz );
416 gammaDirection1 = tmpDirection1;
417
418 // Change reference frame.
419
420 SystemOfRefChange(gammaDirection0,gammaDirection1,
421 gammaPolarization0,gammaPolarization1);
422
423 if (gammaEnergy1 > 0.)
424 {
425 fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
426 fParticleChange->ProposeMomentumDirection( gammaDirection1 );
427 fParticleChange->ProposePolarization( gammaPolarization1 );
428 }
429 else
430 {
431 fParticleChange->SetProposedKineticEnergy(0.) ;
432 fParticleChange->ProposeTrackStatus(fStopAndKill);
433 }
434
435 //
436 // kinematic of the scattered electron
437 //
438
439 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
440
441 // SI - Removed range test
442
443 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
444
445 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
446 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
447
448 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
449
450 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
451 fvect->push_back(dp);
452
453}
454
455//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
456
457G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
458 G4double sinSqrTh)
459{
460 G4double rand1;
461 G4double rand2;
462 G4double phiProbability;
463 G4double phi;
464 G4double a, b;
465
466 do
467 {
468 rand1 = G4UniformRand();
469 rand2 = G4UniformRand();
470 phiProbability=0.;
471 phi = twopi*rand1;
472
473 a = 2*sinSqrTh;
474 b = energyRate + 1/energyRate;
475
476 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
477
478
479
480 }
481 while ( rand2 > phiProbability );
482 return phi;
483}
484
485
486//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
487
488G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
489{
490 G4double dx = a.x();
491 G4double dy = a.y();
492 G4double dz = a.z();
493 G4double x = dx < 0.0 ? -dx : dx;
494 G4double y = dy < 0.0 ? -dy : dy;
495 G4double z = dz < 0.0 ? -dz : dz;
496 if (x < y) {
497 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
498 }else{
499 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
500 }
501}
502
503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
504
505G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
506{
507 G4ThreeVector d0 = direction0.unit();
508 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
509 G4ThreeVector a0 = a1.unit(); // unit vector
510
511 G4double rand1 = G4UniformRand();
512
513 G4double angle = twopi*rand1; // random polar angle
514 G4ThreeVector b0 = d0.cross(a0); // cross product
515
516 G4ThreeVector c;
517
518 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
519 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
520 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
521
522 G4ThreeVector c0 = c.unit();
523
524 return c0;
525
526}
527
528//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
529
530G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
531(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
532{
533
534 //
535 // The polarization of a photon is always perpendicular to its momentum direction.
536 // Therefore this function removes those vector component of gammaPolarization, which
537 // points in direction of gammaDirection
538 //
539 // Mathematically we search the projection of the vector a on the plane E, where n is the
540 // plains normal vector.
541 // The basic equation can be found in each geometry book (e.g. Bronstein):
542 // p = a - (a o n)/(n o n)*n
543
544 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
545}
546
547//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
548
549G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
550 G4double sinSqrTh,
551 G4double phi,
552 G4double costheta)
553{
554 G4double rand1;
555 G4double rand2;
556 G4double cosPhi = std::cos(phi);
557 G4double sinPhi = std::sin(phi);
558 G4double sinTheta = std::sqrt(sinSqrTh);
559 G4double cosSqrPhi = cosPhi*cosPhi;
560 // G4double cossqrth = 1.-sinSqrTh;
561 // G4double sinsqrphi = sinPhi*sinPhi;
562 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
563
564
565 // Determination of Theta
566
567 // ---- MGP ---- Commented out the following 3 lines to avoid compilation
568 // warnings (unused variables)
569 // G4double thetaProbability;
570 G4double theta;
571 // G4double a, b;
572 // G4double cosTheta;
573
574 /*
575
576 depaola method
577
578 do
579 {
580 rand1 = G4UniformRand();
581 rand2 = G4UniformRand();
582 thetaProbability=0.;
583 theta = twopi*rand1;
584 a = 4*normalisation*normalisation;
585 b = (epsilon + 1/epsilon) - 2;
586 thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
587 cosTheta = std::cos(theta);
588 }
589 while ( rand2 > thetaProbability );
590
591 G4double cosBeta = cosTheta;
592
593 */
594
595
596 // Dan Xu method (IEEE TNS, 52, 1160 (2005))
597
598 rand1 = G4UniformRand();
599 rand2 = G4UniformRand();
600
601 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
602 {
603 if (rand2<0.5)
604 theta = pi/2.0;
605 else
606 theta = 3.0*pi/2.0;
607 }
608 else
609 {
610 if (rand2<0.5)
611 theta = 0;
612 else
613 theta = pi;
614 }
615 G4double cosBeta = std::cos(theta);
616 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
617
618 G4ThreeVector gammaPolarization1;
619
620 G4double xParallel = normalisation*cosBeta;
621 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
622 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
623 G4double xPerpendicular = 0.;
624 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
625 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
626
627 G4double xTotal = (xParallel + xPerpendicular);
628 G4double yTotal = (yParallel + yPerpendicular);
629 G4double zTotal = (zParallel + zPerpendicular);
630
631 gammaPolarization1.setX(xTotal);
632 gammaPolarization1.setY(yTotal);
633 gammaPolarization1.setZ(zTotal);
634
635 return gammaPolarization1;
636
637}
638
639//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
640
641void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
642 G4ThreeVector& direction1,
643 G4ThreeVector& polarization0,
644 G4ThreeVector& polarization1)
645{
646 // direction0 is the original photon direction ---> z
647 // polarization0 is the original photon polarization ---> x
648 // need to specify y axis in the real reference frame ---> y
649 G4ThreeVector Axis_Z0 = direction0.unit();
650 G4ThreeVector Axis_X0 = polarization0.unit();
651 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
652
653 G4double direction_x = direction1.getX();
654 G4double direction_y = direction1.getY();
655 G4double direction_z = direction1.getZ();
656
657 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
658 G4double polarization_x = polarization1.getX();
659 G4double polarization_y = polarization1.getY();
660 G4double polarization_z = polarization1.getZ();
661
662 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
663
664}
665
666//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
667
668G4double G4LivermorePolarizedComptonModel::GetMeanFreePath(const G4Track& track,
669 G4double,
670 G4ForceCondition*)
671{
672 const G4DynamicParticle* photon = track.GetDynamicParticle();
673 G4double energy = photon->GetKineticEnergy();
674 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
675 size_t materialIndex = couple->GetIndex();
676 G4double meanFreePath;
677 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
678 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
679 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
680 return meanFreePath;
681}
682
683
Note: See TracBrowser for help on using the repository browser.