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

Last change on this file since 924 was 819, checked in by garnier, 16 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.