source: trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyRayleigh.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File size: 8.3 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//
28// $Id: G4LowEnergyRayleigh.cc,v 1.40 2009/06/11 15:47:08 mantero Exp $
29// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
30//
31// Author: A. Forti
32//         Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
33//
34// History:
35// --------
36// Added Livermore data table construction methods A. Forti
37// Added BuildMeanFreePath A. Forti
38// Added PostStepDoIt A. Forti
39// Added SelectRandomAtom A. Forti
40// Added map of the elements  A.Forti
41// 24.04.01 V.Ivanchenko remove RogueWave
42// 11.08.2001 MGP - Major revision according to a design iteration
43// 06.10.2001 MGP - Added strategy to test range for secondary generation
44// 05.06.2002 F.Longo and G.Depaola  - bug fixed in angular distribution
45// 20.10.2002 G. Depaola   - Change sampling method of theta
46// 22.01.2003 V.Ivanchenko - Cut per region
47// 24.04.2003 V.Ivanchenko - Cut per region mfpt
48//
49// --------------------------------------------------------------------
50
51#include "G4LowEnergyRayleigh.hh"
52#include "Randomize.hh"
53#include "G4ParticleDefinition.hh"
54#include "G4Track.hh"
55#include "G4Step.hh"
56#include "G4ForceCondition.hh"
57#include "G4Gamma.hh"
58#include "G4Electron.hh"
59#include "G4DynamicParticle.hh"
60#include "G4VParticleChange.hh"
61#include "G4ThreeVector.hh"
62#include "G4VCrossSectionHandler.hh"
63#include "G4CrossSectionHandler.hh"
64#include "G4VEMDataSet.hh"
65#include "G4CompositeEMDataSet.hh"
66#include "G4VDataSetAlgorithm.hh"
67#include "G4LogLogInterpolation.hh"
68
69#include "G4MaterialCutsCouple.hh"
70
71G4LowEnergyRayleigh::G4LowEnergyRayleigh(const G4String& processName)
72  : G4VDiscreteProcess(processName),
73    lowEnergyLimit(250*eV),
74    highEnergyLimit(100*GeV),
75    intrinsicLowEnergyLimit(10*eV),
76    intrinsicHighEnergyLimit(100*GeV)
77{
78  if (lowEnergyLimit < intrinsicLowEnergyLimit ||
79      highEnergyLimit > intrinsicHighEnergyLimit)
80    {
81      G4Exception("G4LowEnergyRayleigh::G4LowEnergyRayleigh - energy limit outside intrinsic process validity range");
82    }
83
84  crossSectionHandler = new G4CrossSectionHandler();
85
86  G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
87  G4String formFactorFile = "rayl/re-ff-";
88  formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
89  formFactorData->LoadData(formFactorFile);
90
91  meanFreePathTable = 0;
92
93   if (verboseLevel > 0)
94     {
95       G4cout << GetProcessName() << " is created " << G4endl
96              << "Energy range: "
97              << lowEnergyLimit / keV << " keV - "
98              << highEnergyLimit / GeV << " GeV"
99              << G4endl;
100     }
101
102   G4cout << G4endl;
103   G4cout << "*******************************************************************************" << G4endl;
104   G4cout << "*******************************************************************************" << G4endl;
105   G4cout << "   The class G4LowEnergyRayleigh is NOT SUPPORTED ANYMORE. " << G4endl;
106   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
107   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
108   G4cout << "*******************************************************************************" << G4endl;
109   G4cout << "*******************************************************************************" << G4endl;
110   G4cout << G4endl;
111}
112
113G4LowEnergyRayleigh::~G4LowEnergyRayleigh()
114{
115  delete meanFreePathTable;
116  delete crossSectionHandler;
117  delete formFactorData;
118}
119
120void G4LowEnergyRayleigh::BuildPhysicsTable(const G4ParticleDefinition& )
121{
122
123  crossSectionHandler->Clear();
124  G4String crossSectionFile = "rayl/re-cs-";
125  crossSectionHandler->LoadData(crossSectionFile);
126
127  delete meanFreePathTable;
128  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
129}
130
131G4VParticleChange* G4LowEnergyRayleigh::PostStepDoIt(const G4Track& aTrack,
132                                                     const G4Step& aStep)
133{
134
135  aParticleChange.Initialize(aTrack);
136
137  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
138  G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
139
140  if (photonEnergy0 <= lowEnergyLimit)
141    {
142      aParticleChange.ProposeTrackStatus(fStopAndKill);
143      aParticleChange.ProposeEnergy(0.);
144      aParticleChange.ProposeLocalEnergyDeposit(photonEnergy0);
145      return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
146    }
147
148  //  G4double e0m = photonEnergy0 / electron_mass_c2 ;
149  G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
150
151  // Select randomly one element in the current material
152  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
153  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
154
155  // Sample the angle of the scattered photon
156
157  G4double wlPhoton = h_Planck*c_light/photonEnergy0;
158
159  G4double gReject,x,dataFormFactor;
160  G4double randomFormFactor;
161  G4double cosTheta;
162  G4double sinTheta;
163  G4double fcostheta;
164
165  do
166    {
167      do
168      {
169      cosTheta = 2. * G4UniformRand() - 1.;
170      fcostheta = ( 1. + cosTheta*cosTheta)/2.;
171      } while (fcostheta < G4UniformRand());
172
173      G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
174      x = sinThetaHalf / (wlPhoton/cm);
175      if (x > 1.e+005)
176         dataFormFactor = formFactorData->FindValue(x,Z-1);
177      else
178         dataFormFactor = formFactorData->FindValue(0.,Z-1);
179      randomFormFactor = G4UniformRand() * Z * Z;
180      sinTheta = std::sqrt(1. - cosTheta*cosTheta);
181      gReject = dataFormFactor * dataFormFactor;
182
183    } while( gReject < randomFormFactor);
184
185  // Scattered photon angles. ( Z - axis along the parent photon)
186  G4double phi = twopi * G4UniformRand() ;
187  G4double dirX = sinTheta*std::cos(phi);
188  G4double dirY = sinTheta*std::sin(phi);
189  G4double dirZ = cosTheta;
190
191  // Update G4VParticleChange for the scattered photon
192  G4ThreeVector photonDirection1(dirX, dirY, dirZ);
193
194  photonDirection1.rotateUz(photonDirection0);
195  aParticleChange.ProposeEnergy(photonEnergy0);
196  aParticleChange.ProposeMomentumDirection(photonDirection1);
197
198  aParticleChange.SetNumberOfSecondaries(0);
199
200  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
201}
202
203G4bool G4LowEnergyRayleigh::IsApplicable(const G4ParticleDefinition& particle)
204{
205  return ( &particle == G4Gamma::Gamma() ); 
206}
207
208G4double G4LowEnergyRayleigh::GetMeanFreePath(const G4Track& track, 
209                                              G4double, // previousStepSize
210                                              G4ForceCondition*)
211{
212  const G4DynamicParticle* photon = track.GetDynamicParticle();
213  G4double energy = photon->GetKineticEnergy();
214  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
215  size_t materialIndex = couple->GetIndex();
216
217  G4double meanFreePath;
218  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
219  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
220  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
221  return meanFreePath;
222}
Note: See TracBrowser for help on using the repository browser.