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

Last change on this file since 830 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 17.8 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.22 2006/06/29 19:40:25 gunter Exp $
28// GEANT4 tag $Name: $
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 if (verboseLevel > 0)
113 {
114 G4cout << GetProcessName() << " is created " << G4endl
115 << "Energy range: "
116 << lowEnergyLimit / keV << " keV - "
117 << highEnergyLimit / GeV << " GeV"
118 << G4endl;
119 }
120}
121
122// destructor
123
124G4LowEnergyPolarizedCompton::~G4LowEnergyPolarizedCompton()
125{
126 delete meanFreePathTable;
127 delete crossSectionHandler;
128 delete scatterFunctionData;
129 delete rangeTest;
130}
131
132
133void G4LowEnergyPolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& )
134{
135
136 crossSectionHandler->Clear();
137 G4String crossSectionFile = "comp/ce-cs-";
138 crossSectionHandler->LoadData(crossSectionFile);
139 delete meanFreePathTable;
140 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
141}
142
143G4VParticleChange* G4LowEnergyPolarizedCompton::PostStepDoIt(const G4Track& aTrack,
144 const G4Step& aStep)
145{
146 // The scattered gamma energy is sampled according to Klein - Nishina formula.
147 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
148 // GEANT4 internal units
149 //
150 // Note : Effects due to binding of atomic electrons are negliged.
151
152 aParticleChange.Initialize(aTrack);
153
154 // Dynamic particle quantities
155 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
156 G4double gammaEnergy0 = incidentPhoton->GetKineticEnergy();
157 G4ThreeVector gammaPolarization0 = incidentPhoton->GetPolarization();
158
159 // gammaPolarization0 = gammaPolarization0.unit(); //
160
161 // Protection: a polarisation parallel to the
162 // direction causes problems;
163 // in that case find a random polarization
164
165 G4ThreeVector gammaDirection0 = incidentPhoton->GetMomentumDirection();
166 // ---- MGP ---- Next two lines commented out to remove compilation warnings
167 // G4double scalarproduct = gammaPolarization0.dot(gammaDirection0);
168 // G4double angle = gammaPolarization0.angle(gammaDirection0);
169
170 // Make sure that the polarization vector is perpendicular to the
171 // gamma direction. If not
172
173 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
174 { // only for testing now
175 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
176 }
177 else
178 {
179 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
180 {
181 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
182 }
183 }
184
185 // End of Protection
186
187 // Within energy limit?
188
189 if(gammaEnergy0 <= lowEnergyLimit)
190 {
191 aParticleChange.ProposeTrackStatus(fStopAndKill);
192 aParticleChange.ProposeEnergy(0.);
193 aParticleChange.ProposeLocalEnergyDeposit(gammaEnergy0);
194 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
195 }
196
197 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
198
199 // Select randomly one element in the current material
200
201 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
202 G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
203
204 // Sample the energy and the polarization of the scattered photon
205
206 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
207
208 G4double epsilon0 = 1./(1. + 2*E0_m);
209 G4double epsilon0Sq = epsilon0*epsilon0;
210 G4double alpha1 = - std::log(epsilon0);
211 G4double alpha2 = 0.5*(1.- epsilon0Sq);
212
213 G4double wlGamma = h_Planck*c_light/gammaEnergy0;
214 G4double gammaEnergy1;
215 G4ThreeVector gammaDirection1;
216
217 do {
218 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
219 {
220 epsilon = std::exp(-alpha1*G4UniformRand());
221 epsilonSq = epsilon*epsilon;
222 }
223 else
224 {
225 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
226 epsilon = std::sqrt(epsilonSq);
227 }
228
229 onecost = (1.- epsilon)/(epsilon*E0_m);
230 sinThetaSqr = onecost*(2.-onecost);
231
232 // Protection
233 if (sinThetaSqr > 1.)
234 {
235 if (verboseLevel>0) G4cout
236 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
237 << "sin(theta)**2 = "
238 << sinThetaSqr
239 << "; set to 1"
240 << G4endl;
241 sinThetaSqr = 1.;
242 }
243 if (sinThetaSqr < 0.)
244 {
245 if (verboseLevel>0) G4cout
246 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
247 << "sin(theta)**2 = "
248 << sinThetaSqr
249 << "; set to 0"
250 << G4endl;
251 sinThetaSqr = 0.;
252 }
253 // End protection
254
255 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
256 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
257 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
258
259 } while(greject < G4UniformRand()*Z);
260
261
262 // ****************************************************
263 // Phi determination
264 // ****************************************************
265
266 G4double phi = SetPhi(epsilon,sinThetaSqr);
267
268 //
269 // scattered gamma angles. ( Z - axis along the parent gamma)
270 //
271
272 G4double cosTheta = 1. - onecost;
273
274 // Protection
275
276 if (cosTheta > 1.)
277 {
278 if (verboseLevel>0) G4cout
279 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
280 << "cosTheta = "
281 << cosTheta
282 << "; set to 1"
283 << G4endl;
284 cosTheta = 1.;
285 }
286 if (cosTheta < -1.)
287 {
288 if (verboseLevel>0) G4cout
289 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
290 << "cosTheta = "
291 << cosTheta
292 << "; set to -1"
293 << G4endl;
294 cosTheta = -1.;
295 }
296 // End protection
297
298
299 G4double sinTheta = std::sqrt (sinThetaSqr);
300
301 // Protection
302 if (sinTheta > 1.)
303 {
304 if (verboseLevel>0) G4cout
305 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
306 << "sinTheta = "
307 << sinTheta
308 << "; set to 1"
309 << G4endl;
310 sinTheta = 1.;
311 }
312 if (sinTheta < -1.)
313 {
314 if (verboseLevel>0) G4cout
315 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
316 << "sinTheta = "
317 << sinTheta
318 << "; set to -1"
319 << G4endl;
320 sinTheta = -1.;
321 }
322 // End protection
323
324
325 G4double dirx = sinTheta*std::cos(phi);
326 G4double diry = sinTheta*std::sin(phi);
327 G4double dirz = cosTheta ;
328
329 //
330 // update G4VParticleChange for the scattered photon
331 //
332
333 gammaEnergy1 = epsilon*gammaEnergy0;
334
335 // New polarization
336
337 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
338 sinThetaSqr,
339 phi,
340 cosTheta);
341
342 // Set new direction
343 G4ThreeVector tmpDirection1( dirx,diry,dirz );
344 gammaDirection1 = tmpDirection1;
345
346 // Change reference frame.
347
348 SystemOfRefChange(gammaDirection0,gammaDirection1,
349 gammaPolarization0,gammaPolarization1);
350
351 if (gammaEnergy1 > 0.)
352 {
353 aParticleChange.ProposeEnergy( gammaEnergy1 ) ;
354 aParticleChange.ProposeMomentumDirection( gammaDirection1 );
355 aParticleChange.ProposePolarization( gammaPolarization1 );
356 }
357 else
358 {
359 aParticleChange.ProposeEnergy(0.) ;
360 aParticleChange.ProposeTrackStatus(fStopAndKill);
361 }
362
363 //
364 // kinematic of the scattered electron
365 //
366
367 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 ;
368
369 // Generate the electron only if with large enough range w.r.t. cuts and safety
370
371 G4double safety = aStep.GetPostStepPoint()->GetSafety();
372
373 if (rangeTest->Escape(G4Electron::Electron(),couple,ElecKineEnergy,safety))
374 {
375 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
376 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
377 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
378 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
379 aParticleChange.SetNumberOfSecondaries(1);
380 aParticleChange.AddSecondary(electron);
381 aParticleChange.ProposeLocalEnergyDeposit(0.);
382 }
383 else
384 {
385 aParticleChange.SetNumberOfSecondaries(0);
386 aParticleChange.ProposeLocalEnergyDeposit(ElecKineEnergy);
387 }
388
389 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
390
391}
392
393
394G4double G4LowEnergyPolarizedCompton::SetPhi(G4double energyRate,
395 G4double sinSqrTh)
396{
397 G4double rand1;
398 G4double rand2;
399 G4double phiProbability;
400 G4double phi;
401 G4double a, b;
402
403 do
404 {
405 rand1 = G4UniformRand();
406 rand2 = G4UniformRand();
407 phiProbability=0.;
408 phi = twopi*rand1;
409
410 a = 2*sinSqrTh;
411 b = energyRate + 1/energyRate;
412
413 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
414
415
416
417 }
418 while ( rand2 > phiProbability );
419 return phi;
420}
421
422
423G4ThreeVector G4LowEnergyPolarizedCompton::SetPerpendicularVector(G4ThreeVector& a)
424{
425 G4double dx = a.x();
426 G4double dy = a.y();
427 G4double dz = a.z();
428 G4double x = dx < 0.0 ? -dx : dx;
429 G4double y = dy < 0.0 ? -dy : dy;
430 G4double z = dz < 0.0 ? -dz : dz;
431 if (x < y) {
432 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
433 }else{
434 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
435 }
436}
437
438G4ThreeVector G4LowEnergyPolarizedCompton::GetRandomPolarization(G4ThreeVector& direction0)
439{
440 G4ThreeVector d0 = direction0.unit();
441 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
442 G4ThreeVector a0 = a1.unit(); // unit vector
443
444 G4double rand1 = G4UniformRand();
445
446 G4double angle = twopi*rand1; // random polar angle
447 G4ThreeVector b0 = d0.cross(a0); // cross product
448
449 G4ThreeVector c;
450
451 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
452 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
453 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
454
455 G4ThreeVector c0 = c.unit();
456
457 return c0;
458
459}
460
461
462G4ThreeVector G4LowEnergyPolarizedCompton::GetPerpendicularPolarization
463(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
464{
465
466 //
467 // The polarization of a photon is always perpendicular to its momentum direction.
468 // Therefore this function removes those vector component of gammaPolarization, which
469 // points in direction of gammaDirection
470 //
471 // Mathematically we search the projection of the vector a on the plane E, where n is the
472 // plains normal vector.
473 // The basic equation can be found in each geometry book (e.g. Bronstein):
474 // p = a - (a o n)/(n o n)*n
475
476 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
477}
478
479
480G4ThreeVector G4LowEnergyPolarizedCompton::SetNewPolarization(G4double epsilon,
481 G4double sinSqrTh,
482 G4double phi,
483 G4double costheta)
484{
485 G4double rand1;
486 G4double rand2;
487 G4double cosPhi = std::cos(phi);
488 G4double sinPhi = std::sin(phi);
489 G4double sinTheta = std::sqrt(sinSqrTh);
490 G4double cosSqrPhi = cosPhi*cosPhi;
491 // G4double cossqrth = 1.-sinSqrTh;
492 // G4double sinsqrphi = sinPhi*sinPhi;
493 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
494
495
496 // Determination of Theta
497
498 G4double thetaProbability;
499 G4double theta;
500 G4double a, b;
501 G4double cosTheta;
502
503 do
504 {
505 rand1 = G4UniformRand();
506 rand2 = G4UniformRand();
507 thetaProbability=0.;
508 theta = twopi*rand1;
509 a = 4*normalisation*normalisation;
510 b = (epsilon + 1/epsilon) - 2;
511 thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
512 cosTheta = std::cos(theta);
513 }
514 while ( rand2 > thetaProbability );
515
516 G4double cosBeta = cosTheta;
517 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
518
519 G4ThreeVector gammaPolarization1;
520
521 G4double xParallel = normalisation*cosBeta;
522 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
523 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
524 G4double xPerpendicular = 0.;
525 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
526 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
527
528 G4double xTotal = (xParallel + xPerpendicular);
529 G4double yTotal = (yParallel + yPerpendicular);
530 G4double zTotal = (zParallel + zPerpendicular);
531
532 gammaPolarization1.setX(xTotal);
533 gammaPolarization1.setY(yTotal);
534 gammaPolarization1.setZ(zTotal);
535
536 return gammaPolarization1;
537
538}
539
540
541void G4LowEnergyPolarizedCompton::SystemOfRefChange(G4ThreeVector& direction0,
542 G4ThreeVector& direction1,
543 G4ThreeVector& polarization0,
544 G4ThreeVector& polarization1)
545{
546 // direction0 is the original photon direction ---> z
547 // polarization0 is the original photon polarization ---> x
548 // need to specify y axis in the real reference frame ---> y
549 G4ThreeVector Axis_Z0 = direction0.unit();
550 G4ThreeVector Axis_X0 = polarization0.unit();
551 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
552
553 G4double direction_x = direction1.getX();
554 G4double direction_y = direction1.getY();
555 G4double direction_z = direction1.getZ();
556
557 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
558 G4double polarization_x = polarization1.getX();
559 G4double polarization_y = polarization1.getY();
560 G4double polarization_z = polarization1.getZ();
561
562 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
563
564}
565
566
567G4bool G4LowEnergyPolarizedCompton::IsApplicable(const G4ParticleDefinition& particle)
568{
569 return ( &particle == G4Gamma::Gamma() );
570}
571
572
573G4double G4LowEnergyPolarizedCompton::GetMeanFreePath(const G4Track& track,
574 G4double,
575 G4ForceCondition*)
576{
577 const G4DynamicParticle* photon = track.GetDynamicParticle();
578 G4double energy = photon->GetKineticEnergy();
579 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
580 size_t materialIndex = couple->GetIndex();
581 G4double meanFreePath;
582 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
583 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
584 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
585 return meanFreePath;
586}
587
588
589
590
591
592
593
594
595
596
597
598
599
600
Note: See TracBrowser for help on using the repository browser.