physics/src/G4LowEnergyRayleigh.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // --------------------------------------------------------------------
0027 //
0028 //
0029 // Author: A. Forti
0030 //         Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
0031 //
0032 // History:
0033 // --------
0034 // Added Livermore data table construction methods A. Forti
0035 // Added BuildMeanFreePath A. Forti
0036 // Added PostStepDoIt A. Forti
0037 // Added SelectRandomAtom A. Forti
0038 // Added map of the elements  A.Forti
0039 // 24.04.01 V.Ivanchenko remove RogueWave
0040 // 11.08.2001 MGP - Major revision according to a design iteration
0041 // 06.10.2001 MGP - Added strategy to test range for secondary generation
0042 // 05.06.2002 F.Longo and G.Depaola  - bug fixed in angular distribution
0043 // 20.10.2002 G. Depaola   - Change sampling method of theta
0044 // 22.01.2003 V.Ivanchenko - Cut per region
0045 // 24.04.2003 V.Ivanchenko - Cut per region mfpt
0046 //
0047 // --------------------------------------------------------------------
0048
0049 #include "G4LowEnergyRayleigh.hh"
0050 #include "Randomize.hh"
0051 #include "G4PhysicalConstants.hh"
0052 #include "G4SystemOfUnits.hh"
0053 #include "G4ParticleDefinition.hh"
0054 #include "G4Track.hh"
0055 #include "G4Step.hh"
0056 #include "G4ForceCondition.hh"
0057 #include "G4Gamma.hh"
0058 #include "G4Electron.hh"
0059 #include "G4DynamicParticle.hh"
0060 #include "G4VParticleChange.hh"
0061 #include "G4ThreeVector.hh"
0062 #include "G4RDVCrossSectionHandler.hh"
0063 #include "G4RDCrossSectionHandler.hh"
0064 #include "G4RDVEMDataSet.hh"
0065 #include "G4RDCompositeEMDataSet.hh"
0066 #include "G4RDVDataSetAlgorithm.hh"
0067 #include "G4RDLogLogInterpolation.hh"
0068
0069 #include "G4MaterialCutsCouple.hh"
0070
0071 G4LowEnergyRayleigh::G4LowEnergyRayleigh(const G4String& processName)
0072   : G4VDiscreteProcess(processName),
0073     lowEnergyLimit(250*eV),
0074     highEnergyLimit(100*GeV),
0075     intrinsicLowEnergyLimit(10*eV),
0076     intrinsicHighEnergyLimit(100*GeV)
0077 {
0078   if (lowEnergyLimit < intrinsicLowEnergyLimit ||
0079       highEnergyLimit > intrinsicHighEnergyLimit)
0080     {
0081       G4Exception("G4LowEnergyRayleigh::G4LowEnergyRayleigh()",
0082                   "OutOfRange", FatalException,
0083                   "Energy limit outside intrinsic process validity range!");
0084     }
0085
0086   crossSectionHandler = new G4RDCrossSectionHandler();
0087
0088   G4RDVDataSetAlgorithm* ffInterpolation = new G4RDLogLogInterpolation;
0089   G4String formFactorFile = "rayl/re-ff-";
0090   formFactorData = new G4RDCompositeEMDataSet(ffInterpolation,1.,1.);
0091   formFactorData->LoadData(formFactorFile);
0092
0093   meanFreePathTable = 0;
0094
0095    if (verboseLevel > 0)
0096      {
0097        G4cout << GetProcessName() << " is created " << G4endl
0098           << "Energy range: "
0099           << lowEnergyLimit / keV << " keV - "
0100           << highEnergyLimit / GeV << " GeV"
0101           << G4endl;
0102      }
0103 }
0104
0105 G4LowEnergyRayleigh::~G4LowEnergyRayleigh()
0106 {
0107   delete meanFreePathTable;
0108   delete crossSectionHandler;
0109   delete formFactorData;
0110 }
0111
0112 void G4LowEnergyRayleigh::BuildPhysicsTable(const G4ParticleDefinition& )
0113 {
0114
0115   crossSectionHandler->Clear();
0116   G4String crossSectionFile = "rayl/re-cs-";
0117   crossSectionHandler->LoadData(crossSectionFile);
0118
0119   delete meanFreePathTable;
0120   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
0121 }
0122
0123 G4VParticleChange* G4LowEnergyRayleigh::PostStepDoIt(const G4Track& aTrack,
0124                              const G4Step& aStep)
0125 {
0126
0127   aParticleChange.Initialize(aTrack);
0128
0129   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
0130   G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
0131
0132   if (photonEnergy0 <= lowEnergyLimit)
0133     {
0134       aParticleChange.ProposeTrackStatus(fStopAndKill);
0135       aParticleChange.ProposeEnergy(0.);
0136       aParticleChange.ProposeLocalEnergyDeposit(photonEnergy0);
0137       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
0138     }
0139
0140   //  G4double e0m = photonEnergy0 / electron_mass_c2 ;
0141   G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
0142
0143   // Select randomly one element in the current material
0144   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
0145   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
0146
0147   // Sample the angle of the scattered photon
0148
0149   G4double wlPhoton = h_Planck*c_light/photonEnergy0;
0150
0151   G4double gReject,x,dataFormFactor;
0152   G4double randomFormFactor;
0153   G4double cosTheta;
0154   G4double sinTheta;
0155   G4double fcostheta;
0156
0157   do
0158     {
0159       do
0160       {
0161       cosTheta = 2. * G4UniformRand() - 1.;
0162       fcostheta = ( 1. + cosTheta*cosTheta)/2.;
0163       } while (fcostheta < G4UniformRand());
0164
0165       G4double sinThetaHalf = std::sqrt((1. - cosTheta) / 2.);
0166       x = sinThetaHalf / (wlPhoton/cm);
0167       if (x > 1.e+005)
0168          dataFormFactor = formFactorData->FindValue(x,Z-1);
0169       else
0170          dataFormFactor = formFactorData->FindValue(0.,Z-1);
0171       randomFormFactor = G4UniformRand() * Z * Z;
0172       sinTheta = std::sqrt(1. - cosTheta*cosTheta);
0173       gReject = dataFormFactor * dataFormFactor;
0174
0175     } while( gReject < randomFormFactor);
0176
0177   // Scattered photon angles. ( Z - axis along the parent photon)
0178   G4double phi = twopi * G4UniformRand() ;
0179   G4double dirX = sinTheta*std::cos(phi);
0180   G4double dirY = sinTheta*std::sin(phi);
0181   G4double dirZ = cosTheta;
0182
0183   // Update G4VParticleChange for the scattered photon
0184   G4ThreeVector photonDirection1(dirX, dirY, dirZ);
0185
0186   photonDirection1.rotateUz(photonDirection0);
0187   aParticleChange.ProposeEnergy(photonEnergy0);
0188   aParticleChange.ProposeMomentumDirection(photonDirection1);
0189
0190   aParticleChange.SetNumberOfSecondaries(0);
0191
0192   return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0193 }
0194
0195 G4bool G4LowEnergyRayleigh::IsApplicable(const G4ParticleDefinition& particle)
0196 {
0197   return ( &particle == G4Gamma::Gamma() );
0198 }
0199
0200 G4double G4LowEnergyRayleigh::GetMeanFreePath(const G4Track& track,
0201                           G4double, // previousStepSize
0202                           G4ForceCondition*)
0203 {
0204   const G4DynamicParticle* photon = track.GetDynamicParticle();
0205   G4double energy = photon->GetKineticEnergy();
0206   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
0207   size_t materialIndex = couple->GetIndex();
0208
0209   G4double meanFreePath;
0210   if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
0211   else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
0212   else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
0213   return meanFreePath;
0214 }