source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.cc@ 980

Last change on this file since 980 was 968, checked in by garnier, 17 years ago

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