952 lines
34 KiB
C++
952 lines
34 KiB
C++
//
|
|
// ********************************************************************
|
|
// * DISCLAIMER *
|
|
// * *
|
|
// * The following disclaimer summarizes all the specific disclaimers *
|
|
// * of contributors to this software. The specific disclaimers,which *
|
|
// * govern, are listed with their locations in: *
|
|
// * http://cern.ch/geant4/license *
|
|
// * *
|
|
// * Neither the authors of this software system, nor their employing *
|
|
// * institutes,nor the agencies providing financial support for this *
|
|
// * work make any representation or warranty, express or implied, *
|
|
// * regarding this software system or assume any liability for its *
|
|
// * use. *
|
|
// * *
|
|
// * This code implementation is the intellectual property of the *
|
|
// * GEANT4 collaboration. *
|
|
// * By copying, distributing or modifying the Program (or any work *
|
|
// * based on the Program) you indicate your acceptance of this *
|
|
// * statement, and all its terms. *
|
|
// ********************************************************************
|
|
//
|
|
//
|
|
// $Id: G4MultipleScattering52.cc,v 1.2 2004/12/01 19:37:14 vnivanch Exp $
|
|
// GEANT4 tag $Name: geant4-07-00-cand-03 $
|
|
//
|
|
// -----------------------------------------------------------------------------
|
|
// 16/05/01 value of cparm changed , L.Urban
|
|
// 18/05/01 V.Ivanchenko Clean up against Linux ANSI compilation
|
|
// 07/08/01 new methods Store/Retrieve PhysicsTable (mma)
|
|
// 23-08-01 new angle and z distribution,energy dependence reduced,
|
|
// Store,Retrieve methods commented out temporarily, L.Urban
|
|
// 27-08-01 in BuildPhysicsTable:aParticleType.GetParticleName()=="mu+" (mma)
|
|
// 28-08-01 GetContinuousStepLimit and AlongStepDoIt moved from .icc file (mma)
|
|
// 03-09-01 value of data member factlim changed, L.Urban
|
|
// 10-09-01 small change in GetContinuousStepLimit, L.Urban
|
|
// 11-09-01 G4MultipleScatteringx put as default G4MultipleScattering
|
|
// store/retrieve physics table reactivated (mma)
|
|
// 13-09-01 corr. in ComputeTransportCrossSection, L.Urban
|
|
// 14-09-01 protection in GetContinuousStepLimit, L.Urban
|
|
// 17-09-01 migration of Materials to pure STL (mma)
|
|
// 27-09-01 value of data member factlim changed, L.Urban
|
|
// 31-10-01 big fixed in PostStepDoIt,L.Urban
|
|
// 24-04-02 some minor changes in boundary algorithm, L.Urban
|
|
// 06-05-02 bug fixed in GetContinuousStepLimit, L.Urban
|
|
// 24-05-02 changes in angle distribution and boundary algorithm, L.Urban
|
|
// 11-06-02 bug fixed in ComputeTransportCrossSection, L.Urban
|
|
// 12-08-02 bug fixed in PostStepDoIt (lateral displacement), L.Urban
|
|
// 15-08-02 new angle distribution, L.Urban
|
|
// 26-09-02 angle distribution + boundary algorithm modified, L.Urban
|
|
// 15-10-02 temporary fix for proton scattering
|
|
// 30-10-02 modified angle distribution,mods in boundary algorithm,
|
|
// changes in data members, L.Urban
|
|
// 30-10-02 rename variable cm - Ecm, V.Ivanchenko
|
|
// 11-12-02 precision problem in ComputeTransportCrossSection
|
|
// for small Tkin/for heavy particles cured, L.Urban
|
|
// 05-02-03 changes in data members, new sampling for geom.
|
|
// path length, step dependence reduced with new
|
|
// method
|
|
// 17-03-03 cut per region, V.Ivanchenko
|
|
// 13-04-03 add initialisation in GetContinuesStepLimit
|
|
// + change table size (V.Ivanchenko)
|
|
// 26-04-03 fix problems of retrieve tables (M.Asai)
|
|
// 23-05-03 important change in angle distribution for muons/hadrons
|
|
// the central part now is similar to the Highland parametrization +
|
|
// minor correction in angle sampling algorithm (for all particles)
|
|
// (L.Urban)
|
|
// 24-05-03 bug in nuclear size corr.computation fixed thanks to Vladimir(L.Urban)
|
|
// 30-05-03 misprint in PostStepDoIt corrected(L.Urban)
|
|
// 08-08-03 This class is frozen at the release 5.2 (V.Ivanchenko)
|
|
// 08-11-04 Remove Store/Retrieve tables (V.Ivantchenko)
|
|
// -----------------------------------------------------------------------------
|
|
//
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
#include "G4MultipleScattering52.hh"
|
|
#include "G4StepStatus.hh"
|
|
#include "G4Navigator.hh"
|
|
#include "G4TransportationManager.hh"
|
|
#include "Randomize.hh"
|
|
#include "G4ProductionCutsTable.hh"
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
using namespace std;
|
|
|
|
G4MultipleScattering52::G4MultipleScattering52(const G4String& processName)
|
|
: G4VContinuousDiscreteProcess(processName),
|
|
theTransportMeanFreePathTable(0),
|
|
taubig(8.0),tausmall(1.e-14),taulim(1.e-5),
|
|
LowestKineticEnergy(0.1*keV),
|
|
HighestKineticEnergy(100.*TeV),
|
|
TotBin(100),
|
|
materialIndex(0),
|
|
tLast (0.0),
|
|
zLast (0.0),
|
|
boundary(true),
|
|
facrange(0.199),tlimit(1.e10*mm),tlimitmin(1.e-7*mm),
|
|
cf(1.001),
|
|
stepno(0),stepnolastmsc(-1000000),nsmallstep(5),
|
|
laststep(0.),
|
|
valueGPILSelectionMSC(NotCandidateForSelection),
|
|
zmean(0.),samplez(true),
|
|
range(1.),T0(1.),T1(1.),lambda0(1.),lambda1(-1.),
|
|
Tlow(0.),alam(1.),blam(1.),dtrl(0.15),
|
|
lambdam(-1.),clam(1.),zm(1.),cthm(1.),
|
|
fLatDisplFlag(true),
|
|
NuclCorrPar (0.0615),
|
|
FactPar(0.40),
|
|
facxsi(1.)
|
|
{ }
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4MultipleScattering52::~G4MultipleScattering52()
|
|
{
|
|
if(theTransportMeanFreePathTable)
|
|
{
|
|
theTransportMeanFreePathTable->clearAndDestroy();
|
|
delete theTransportMeanFreePathTable;
|
|
}
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void G4MultipleScattering52::BuildPhysicsTable(
|
|
const G4ParticleDefinition& aParticleType)
|
|
{
|
|
// set values of some data members
|
|
if((aParticleType.GetParticleName() == "e-") ||
|
|
(aParticleType.GetParticleName() == "e+"))
|
|
{
|
|
// parameters for e+/e-
|
|
alfa1 = 1.45 ;
|
|
alfa2 = 0.60 ;
|
|
alfa3 = 0.30 ;
|
|
b = 1. ;
|
|
xsi = facxsi*2.22 ;
|
|
c0 = 2.30 ;
|
|
}
|
|
else
|
|
{
|
|
// parameters for heavy particles
|
|
alfa1 = 1.10 ;
|
|
alfa2 = 0.14 ;
|
|
alfa3 = 0.07 ;
|
|
b = 1. ;
|
|
xsi = facxsi*2.70 ;
|
|
c0 = 1.40 ;
|
|
}
|
|
|
|
// ..............................
|
|
Tlow = aParticleType.GetPDGMass();
|
|
|
|
// tables are built for MATERIALS
|
|
const G4double sigmafactor = twopi*classic_electr_radius*
|
|
classic_electr_radius;
|
|
G4double KineticEnergy,AtomicNumber,AtomicWeight,sigma,lambda;
|
|
G4double density;
|
|
|
|
// destroy old tables if any
|
|
if (theTransportMeanFreePathTable)
|
|
{
|
|
theTransportMeanFreePathTable->clearAndDestroy();
|
|
delete theTransportMeanFreePathTable;
|
|
}
|
|
|
|
// create table
|
|
const G4ProductionCutsTable* theCoupleTable=
|
|
G4ProductionCutsTable::GetProductionCutsTable();
|
|
size_t numOfCouples = theCoupleTable->GetTableSize();
|
|
|
|
theTransportMeanFreePathTable = new G4PhysicsTable(numOfCouples);
|
|
|
|
// loop for materials
|
|
for (size_t i=0; i<numOfCouples; i++)
|
|
{
|
|
|
|
// create physics vector and fill it
|
|
G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
|
|
LowestKineticEnergy,HighestKineticEnergy,TotBin);
|
|
|
|
// get elements in the material
|
|
const G4MaterialCutsCouple* couple = theCoupleTable->
|
|
GetMaterialCutsCouple(i);
|
|
const G4Material* material = couple->GetMaterial();
|
|
const G4ElementVector* theElementVector = material->GetElementVector();
|
|
const G4double* NbOfAtomsPerVolume =
|
|
material->GetVecNbOfAtomsPerVolume();
|
|
const G4int NumberOfElements = material->GetNumberOfElements();
|
|
density = material->GetDensity();
|
|
|
|
// loop for kinetic energy values
|
|
for (G4int i=0; i<TotBin; i++)
|
|
{
|
|
KineticEnergy = aVector->GetLowEdgeEnergy(i);
|
|
sigma = 0.;
|
|
|
|
// loop for element in the material
|
|
for (G4int iel=0; iel<NumberOfElements; iel++)
|
|
{
|
|
AtomicNumber = (*theElementVector)[iel]->GetZ();
|
|
AtomicWeight = (*theElementVector)[iel]->GetA();
|
|
sigma += NbOfAtomsPerVolume[iel]*
|
|
ComputeTransportCrossSection(aParticleType,KineticEnergy,
|
|
AtomicNumber,AtomicWeight);
|
|
}
|
|
sigma *= sigmafactor;
|
|
lambda = 1./sigma;
|
|
aVector->PutValue(i,lambda);
|
|
}
|
|
|
|
theTransportMeanFreePathTable->insert(aVector);
|
|
|
|
}
|
|
|
|
if((aParticleType.GetParticleName() == "e-" ) ||
|
|
(aParticleType.GetParticleName() == "mu+" ) ||
|
|
(aParticleType.GetParticleName() == "Mu" ) ||
|
|
(aParticleType.GetParticleName() == "proton") ) PrintInfoDefinition();
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4double G4MultipleScattering52::ComputeTransportCrossSection(
|
|
const G4ParticleDefinition& aParticleType,
|
|
G4double KineticEnergy,
|
|
G4double AtomicNumber,G4double AtomicWeight)
|
|
{
|
|
const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
|
|
Bohr_radius*Bohr_radius/(hbarc*hbarc);
|
|
const G4double epsmin = 1.e-4 , epsmax = 1.e10;
|
|
|
|
const G4double Zdat[15] = { 4., 6.,13.,20.,26.,29.,32.,38.,47.,
|
|
50.,56.,64.,74.,79.,82. };
|
|
|
|
const G4double Tdat[23] = {0.0001*MeV,0.0002*MeV,0.0004*MeV,0.0007*MeV,
|
|
0.001*MeV,0.002*MeV,0.004*MeV,0.007*MeV,
|
|
0.01*MeV,0.02*MeV,0.04*MeV,0.07*MeV,
|
|
0.1*MeV,0.2*MeV,0.4*MeV,0.7*MeV,
|
|
1.*MeV,2.*MeV,4.*MeV,7.*MeV,10.*MeV,20.*MeV,
|
|
10000.0*MeV};
|
|
|
|
// corr. factors for e-/e+ lambda
|
|
G4double celectron[15][23] =
|
|
{{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
|
|
1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
|
|
1.112,1.108,1.100,1.093,1.089,1.087,0.7235 },
|
|
{1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
|
|
1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
|
|
1.109,1.105,1.097,1.090,1.086,1.082,0.7925 },
|
|
{2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
|
|
1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
|
|
1.131,1.124,1.113,1.104,1.099,1.098,0.9147 },
|
|
{3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
|
|
1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
|
|
1.112,1.105,1.096,1.089,1.085,1.098,0.9700 },
|
|
{6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
|
|
1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
|
|
1.073,1.070,1.064,1.059,1.056,1.056,1.0022 },
|
|
{9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
|
|
1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
|
|
1.074,1.070,1.063,1.059,1.056,1.052,1.0158 },
|
|
{11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
|
|
1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
|
|
1.068,1.064,1.059,1.054,1.051,1.050,1.0284 },
|
|
{18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
|
|
1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
|
|
1.039,1.037,1.034,1.031,1.030,1.036,1.0515 },
|
|
{18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
|
|
1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
|
|
1.031,1.028,1.024,1.022,1.021,1.024,1.0834 },
|
|
{14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
|
|
1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
|
|
1.020,1.017,1.015,1.013,1.013,1.020,1.0937 },
|
|
{14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
|
|
1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
|
|
0.995,0.993,0.993,0.993,0.993,1.011,1.1140 },
|
|
{22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
|
|
1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
|
|
0.974,0.972,0.973,0.974,0.975,0.987,1.1410 },
|
|
{50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
|
|
1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
|
|
0.950,0.947,0.949,0.952,0.954,0.963,1.1750 },
|
|
{65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
|
|
1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
|
|
0.941,0.938,0.940,0.944,0.946,0.954,1.1922 },
|
|
// {45.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239, // paper.....
|
|
{55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
|
|
1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
|
|
0.933,0.930,0.933,0.936,0.939,0.949,1.2026 }};
|
|
G4double cpositron[15][23] = {
|
|
{2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
|
|
1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
|
|
1.131,1.126,1.117,1.108,1.103,1.100,0.7235 },
|
|
{3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
|
|
1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
|
|
1.138,1.132,1.122,1.113,1.108,1.102,0.7925 },
|
|
{7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
|
|
1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
|
|
1.203,1.190,1.173,1.159,1.151,1.145,0.9147 },
|
|
{9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
|
|
1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
|
|
1.225,1.210,1.191,1.175,1.166,1.174,0.9700 },
|
|
{17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
|
|
1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
|
|
1.217,1.203,1.184,1.169,1.160,1.151,1.0022 },
|
|
{24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
|
|
1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
|
|
1.237,1.222,1.201,1.184,1.174,1.159,1.0158 },
|
|
{23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
|
|
1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
|
|
1.252,1.234,1.212,1.194,1.183,1.170,1.0284 },
|
|
{22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
|
|
2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
|
|
1.254,1.237,1.214,1.195,1.185,1.179,1.0515 },
|
|
{33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
|
|
2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
|
|
1.312,1.288,1.258,1.235,1.221,1.205,1.0834 },
|
|
{32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
|
|
2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
|
|
1.320,1.294,1.264,1.240,1.226,1.214,1.0937 },
|
|
{29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
|
|
2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
|
|
1.328,1.302,1.270,1.245,1.231,1.233,1.1140 },
|
|
{38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
|
|
2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
|
|
1.361,1.330,1.294,1.267,1.251,1.239,1.1410 },
|
|
{49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
|
|
3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
|
|
1.409,1.372,1.330,1.298,1.280,1.258,1.1750 },
|
|
{59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
|
|
3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
|
|
1.442,1.400,1.354,1.319,1.299,1.272,1.1922 },
|
|
{56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
|
|
3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
|
|
1.456,1.412,1.364,1.328,1.307,1.282,1.2026 }};
|
|
|
|
G4double sigma;
|
|
G4double Z23 = 2.*log(AtomicNumber)/3.; Z23 = exp(Z23);
|
|
|
|
G4double ParticleMass = aParticleType.GetPDGMass();
|
|
G4double ParticleKineticEnergy = KineticEnergy ;
|
|
|
|
// correction if particle .ne. e-/e+
|
|
// compute equivalent kinetic energy
|
|
// lambda depends on p*beta ....
|
|
G4double Mass = ParticleMass ;
|
|
if((aParticleType.GetParticleName() != "e-") &&
|
|
(aParticleType.GetParticleName() != "e+") )
|
|
{
|
|
G4double TAU = KineticEnergy/Mass ;
|
|
G4double c = Mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
|
|
G4double w = c-2. ;
|
|
G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
|
|
KineticEnergy = electron_mass_c2*tau ;
|
|
Mass = electron_mass_c2 ;
|
|
}
|
|
|
|
G4double Charge = aParticleType.GetPDGCharge();
|
|
G4double ChargeSquare = Charge*Charge/(eplus*eplus);
|
|
|
|
G4double TotalEnergy = KineticEnergy + Mass ;
|
|
G4double beta2 = KineticEnergy*(TotalEnergy+Mass)
|
|
/(TotalEnergy*TotalEnergy);
|
|
G4double bg2 = KineticEnergy*(TotalEnergy+Mass)
|
|
/(Mass*Mass);
|
|
|
|
G4double eps = epsfactor*bg2/Z23;
|
|
|
|
if (eps<epsmin) sigma = 2.*eps*eps;
|
|
else if(eps<epsmax) sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
|
|
else sigma = log(2.*eps)-1.+1./eps;
|
|
|
|
sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
|
|
|
|
// nuclear size effect correction for high energy
|
|
// ( a simple approximation at present)
|
|
G4double corrnuclsize,a,x0,w1,w2,w;
|
|
|
|
x0 = 1. - NuclCorrPar*ParticleMass/(ParticleKineticEnergy*
|
|
exp(log(AtomicWeight/(g/mole))/3.));
|
|
if ( (x0 < -1.) || (ParticleKineticEnergy <= 10.*MeV))
|
|
{ x0 = -1.; corrnuclsize = 1.;}
|
|
else
|
|
{ a = 1.+1./eps;
|
|
if (eps > epsmax) w1=log(2.*eps)+1./eps-3./(8.*eps*eps);
|
|
else w1=log((a+1.)/(a-1.))-2./(a+1.);
|
|
w = 1./((1.-x0)*eps);
|
|
if (w < epsmin) w2=-log(w)-1.+2.*w-1.5*w*w;
|
|
else w2 = log((a-x0)/(a-1.))-(1.-x0)/(a-x0);
|
|
corrnuclsize = w1/w2;
|
|
corrnuclsize = exp(-FactPar*ParticleMass/ParticleKineticEnergy)*
|
|
(corrnuclsize-1.)+1.;
|
|
}
|
|
|
|
// interpolate in AtomicNumber and beta2
|
|
// get bin number in Z
|
|
G4int iZ = 14;
|
|
while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
|
|
if (iZ==14) iZ = 13;
|
|
if (iZ==-1) iZ = 0 ;
|
|
|
|
G4double Z1 = Zdat[iZ];
|
|
G4double Z2 = Zdat[iZ+1];
|
|
G4double ratZ = (AtomicNumber-Z1)/(Z2-Z1);
|
|
|
|
// get bin number in T (beta2)
|
|
G4int iT = 22;
|
|
while ((iT>=0)&&(Tdat[iT]>=KineticEnergy)) iT -= 1;
|
|
if(iT==22) iT = 21;
|
|
if(iT==-1) iT = 0 ;
|
|
|
|
// calculate betasquare values
|
|
G4double T = Tdat[iT], E = T + electron_mass_c2;
|
|
G4double b2small = T*(E+electron_mass_c2)/(E*E);
|
|
T = Tdat[iT+1]; E = T + electron_mass_c2;
|
|
G4double b2big = T*(E+electron_mass_c2)/(E*E);
|
|
G4double ratb2 = (beta2-b2small)/(b2big-b2small);
|
|
|
|
G4double c1,c2,cc1,cc2,corr;
|
|
|
|
if (Charge < 0.)
|
|
{
|
|
c1 = celectron[iZ][iT];
|
|
c2 = celectron[iZ+1][iT];
|
|
cc1 = c1+ratZ*(c2-c1);
|
|
|
|
c1 = celectron[iZ][iT+1];
|
|
c2 = celectron[iZ+1][iT+1];
|
|
cc2 = c1+ratZ*(c2-c1);
|
|
|
|
corr = cc1+ratb2*(cc2-cc1);
|
|
sigma /= corr;
|
|
}
|
|
|
|
if (Charge > 0.)
|
|
{
|
|
c1 = cpositron[iZ][iT];
|
|
c2 = cpositron[iZ+1][iT];
|
|
cc1 = c1+ratZ*(c2-c1);
|
|
|
|
c1 = cpositron[iZ][iT+1];
|
|
c2 = cpositron[iZ+1][iT+1];
|
|
cc2 = c1+ratZ*(c2-c1);
|
|
|
|
corr = cc1+ratb2*(cc2-cc1);
|
|
sigma /= corr;
|
|
}
|
|
|
|
// nucl. size correction for particles other than e+/e- only at present !!!!
|
|
if((aParticleType.GetParticleName() != "e-") &&
|
|
(aParticleType.GetParticleName() != "e+") )
|
|
sigma /= corrnuclsize;
|
|
|
|
return sigma;
|
|
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4double G4MultipleScattering52::GetContinuousStepLimit(
|
|
const G4Track& track,
|
|
G4double,
|
|
G4double currentMinimumStep,
|
|
G4double&)
|
|
{
|
|
G4double zPathLength,tPathLength;
|
|
const G4DynamicParticle* aParticle;
|
|
G4double tau,zt,cz,cz1,grej,grej0;
|
|
const G4double expmax = 100., ztmax = (2.*expmax+1.)/(2.*expmax+3.) ;
|
|
const G4double tmax = 1.e20*mm ;
|
|
|
|
G4bool isOut;
|
|
|
|
// this process is not a candidate for selection by default
|
|
valueGPILSelectionMSC = NotCandidateForSelection;
|
|
|
|
tPathLength = currentMinimumStep;
|
|
|
|
const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
|
|
materialIndex = couple->GetIndex();
|
|
|
|
aParticle = track.GetDynamicParticle();
|
|
T0 = aParticle->GetKineticEnergy();
|
|
|
|
lambda0 = (*theTransportMeanFreePathTable)
|
|
(materialIndex)->GetValue(T0,isOut);
|
|
|
|
|
|
range = G4EnergyLossTables::GetRange(aParticle->GetDefinition(),
|
|
T0,couple);
|
|
//VI Initialisation at the beginning of the step
|
|
cthm = 1.;
|
|
lambda1 = -1.;
|
|
lambdam = -1.;
|
|
alam = range;
|
|
blam = 1.+alam/lambda0 ;
|
|
zm = 1.;
|
|
|
|
// special treatment near boundaries ?
|
|
if (boundary && range >= currentMinimumStep)
|
|
{
|
|
// step limitation at boundary ?
|
|
stepno = track.GetCurrentStepNumber() ;
|
|
if(stepno == 1)
|
|
{
|
|
stepnolastmsc = -1000000 ;
|
|
tlimit = 1.e10 ;
|
|
}
|
|
|
|
if(stepno > 1)
|
|
{
|
|
if(track.GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
|
|
{
|
|
stepnolastmsc = stepno ;
|
|
// if : diff.treatment for small/not small Z
|
|
if(range > lambda0)
|
|
tlimit = facrange*range ;
|
|
else
|
|
tlimit = facrange*lambda0 ;
|
|
if(tlimit < tlimitmin) tlimit = tlimitmin ;
|
|
laststep = tlimit ;
|
|
if(tPathLength > tlimit)
|
|
{
|
|
tPathLength = tlimit ;
|
|
valueGPILSelectionMSC = CandidateForSelection;
|
|
}
|
|
}
|
|
else if(stepno > stepnolastmsc)
|
|
{
|
|
if((stepno - stepnolastmsc) < nsmallstep)
|
|
{
|
|
if(tPathLength > tlimit)
|
|
{
|
|
laststep *= cf ;
|
|
tPathLength = laststep ;
|
|
valueGPILSelectionMSC = CandidateForSelection;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// do the true -> geom transformation
|
|
zmean = tPathLength;
|
|
|
|
tau = tPathLength/lambda0 ;
|
|
|
|
if (tau < tausmall || range < currentMinimumStep) zPathLength = tPathLength;
|
|
else
|
|
{
|
|
if(tPathLength/range < dtrl) zmean = lambda0*(1.-exp(-tau));
|
|
else
|
|
{
|
|
T1 = G4EnergyLossTables::GetPreciseEnergyFromRange(
|
|
aParticle->GetDefinition(),range-tPathLength,couple);
|
|
lambda1 = (*theTransportMeanFreePathTable)
|
|
(materialIndex)->GetValue(T1,isOut);
|
|
if(T0 < Tlow)
|
|
alam = range ;
|
|
else
|
|
alam = lambda0*tPathLength/(lambda0-lambda1) ;
|
|
blam = 1.+alam/lambda0 ;
|
|
if(tPathLength/range < 2.*dtrl)
|
|
{
|
|
zmean = alam*(1.-exp(blam*log(1.-tPathLength/alam)))/blam ;
|
|
lambdam = -1. ;
|
|
}
|
|
else
|
|
{
|
|
G4double w = 1.-0.5*tPathLength/alam ;
|
|
lambdam = lambda0*w ;
|
|
clam = 1.+alam/lambdam ;
|
|
cthm = exp(alam*log(w)/lambda0) ;
|
|
zm = alam*(1.-exp(blam*log(w)))/blam ;
|
|
zmean = zm + alam*(1.-exp(clam*log(w)))*cthm/clam ;
|
|
}
|
|
}
|
|
|
|
// sample z
|
|
zt = zmean/tPathLength ;
|
|
if (samplez && (zt < ztmax) && (zt > 0.5))
|
|
{
|
|
cz = 0.5*(3.*zt-1.)/(1.-zt) ;
|
|
if(tPathLength < exp(log(tmax)/(2.*cz)))
|
|
{
|
|
cz1 = 1.+cz ;
|
|
grej0 = exp(cz1*log(cz*tPathLength/cz1))/cz ;
|
|
do
|
|
{
|
|
zPathLength = tPathLength*exp(log(G4UniformRand())/cz1) ;
|
|
grej = exp(cz*log(zPathLength))*(tPathLength-zPathLength)/grej0 ;
|
|
} while (grej < G4UniformRand()) ;
|
|
}
|
|
else zPathLength = zmean;
|
|
}
|
|
else zPathLength = zmean;
|
|
}
|
|
// protection against z > lambda
|
|
if(zPathLength > lambda0)
|
|
zPathLength = lambda0 ;
|
|
|
|
tLast = tPathLength;
|
|
zLast = zPathLength;
|
|
|
|
return zPathLength;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4VParticleChange* G4MultipleScattering52::AlongStepDoIt(
|
|
const G4Track& track,const G4Step& step)
|
|
{
|
|
// only a geom path->true path transformation is performed
|
|
|
|
fParticleChange.Initialize(track);
|
|
|
|
G4double geomPathLength = step.GetStepLength();
|
|
|
|
G4double truePathLength = 0. ;
|
|
|
|
//VI change order of if operators
|
|
if(geomPathLength == zLast) truePathLength = tLast;
|
|
else if(geomPathLength/lambda0 < tausmall) truePathLength = geomPathLength;
|
|
else
|
|
{
|
|
if(lambda1 < 0.) truePathLength = -lambda0*log(1.-geomPathLength/lambda0) ;
|
|
else if(lambdam < 0.)
|
|
{
|
|
if(blam*geomPathLength/alam < 1.)
|
|
truePathLength = alam*(1.-exp(log(1.-blam*geomPathLength/alam)/
|
|
blam)) ;
|
|
else
|
|
truePathLength = tLast;
|
|
}
|
|
else
|
|
{
|
|
if(geomPathLength <= zm)
|
|
{
|
|
if(blam*geomPathLength/alam < 1.)
|
|
truePathLength = alam*(1.-exp(log(1.-blam*geomPathLength/alam)/
|
|
blam)) ;
|
|
else
|
|
truePathLength = 0.5*tLast;
|
|
|
|
lambdam = -1. ;
|
|
}
|
|
else
|
|
{
|
|
if(clam*(geomPathLength-zm)/(alam*cthm) < 1.)
|
|
truePathLength = 0.5*tLast + alam*(1.-
|
|
exp(log(1.-clam*(geomPathLength-zm)/(alam*cthm)))/clam) ;
|
|
else
|
|
truePathLength = tLast ;
|
|
}
|
|
}
|
|
// protection ....
|
|
if(truePathLength > tLast)
|
|
truePathLength = tLast ;
|
|
}
|
|
|
|
//VI truePath length cannot be smaller than geomPathLength
|
|
if (truePathLength < geomPathLength) truePathLength = geomPathLength;
|
|
fParticleChange.ProposeTrueStepLength(truePathLength);
|
|
|
|
return &fParticleChange;
|
|
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4VParticleChange* G4MultipleScattering52::PostStepDoIt(
|
|
const G4Track& trackData,
|
|
const G4Step& stepData)
|
|
{
|
|
// angle distribution parameters
|
|
const G4double kappa = 2.5, kappapl1 = kappa+1., kappami1 = kappa-1. ;
|
|
|
|
fParticleChange.Initialize(trackData);
|
|
G4double truestep = stepData.GetStepLength();
|
|
|
|
const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
|
|
G4double KineticEnergy = aParticle->GetKineticEnergy();
|
|
G4double Mass = aParticle->GetDefinition()->GetPDGMass() ;
|
|
|
|
// do nothing for stopped particles !
|
|
if(KineticEnergy > 0.)
|
|
{
|
|
// change direction first ( scattering )
|
|
G4double cth = 1.0 ;
|
|
G4double tau = truestep/lambda0 ;
|
|
|
|
if (tau < tausmall) cth = 1.;
|
|
else if(tau > taubig) cth = -1.+2.*G4UniformRand();
|
|
else
|
|
{
|
|
if(lambda1 > 0.)
|
|
{
|
|
if(lambdam < 0.)
|
|
tau = -alam*log(1.-truestep/alam)/lambda0 ;
|
|
else
|
|
tau = -log(cthm)-alam*log(1.-(truestep-0.5*tLast)/alam)/lambdam ;
|
|
}
|
|
|
|
if(tau > taubig) cth = -1.+2.*G4UniformRand();
|
|
else
|
|
{
|
|
const G4double amax=25. ;
|
|
const G4double tau0 = 0.02 ;
|
|
const G4double c_highland = 13.6*MeV, corr_highland=0.038 ;
|
|
|
|
const G4double x1fac1 = exp(-xsi) ;
|
|
const G4double x1fac2 = (1.-(1.+xsi)*x1fac1)/(1.-x1fac1) ;
|
|
const G4double x1fac3 = 1.3 ; // x1fac3 >= 1. !!!!!!!!!
|
|
|
|
G4double a,x0,c,xmean1,xmean2,
|
|
xmeanth,prob,qprob ;
|
|
G4double ea,eaa,b1,bx,eb1,ebx,cnorm1,cnorm2,f1x0,f2x0,w ;
|
|
|
|
// for heavy particles take the width of the cetral part
|
|
// from the Highland formula
|
|
// (Particle Physics Booklet, July 2002, eq. 26.10)
|
|
if(Mass > electron_mass_c2) // + other conditions (beta, x/X0,...?)
|
|
{
|
|
G4double Q = fabs(aParticle->GetDefinition()->GetPDGCharge()) ;
|
|
G4double X0 = trackData.GetMaterialCutsCouple()->
|
|
GetMaterial()->GetRadlen() ;
|
|
G4double xx0 = truestep/X0 ;
|
|
G4double betacp = KineticEnergy*(KineticEnergy+2.*Mass)/
|
|
(KineticEnergy+Mass) ;
|
|
G4double theta0=c_highland*Q*sqrt(xx0)*
|
|
(1.+corr_highland*log(xx0))/betacp ;
|
|
if(theta0 > tausmall)
|
|
a = 0.5/(1.-cos(theta0)) ;
|
|
else
|
|
a = 1./(theta0*theta0) ;
|
|
}
|
|
else
|
|
{
|
|
w = log(tau/tau0) ;
|
|
if(tau < tau0)
|
|
a = (alfa1-alfa2*w)/tau ;
|
|
else
|
|
a = (alfa1+alfa3*w)/tau ;
|
|
}
|
|
|
|
xmeanth = exp(-tau) ;
|
|
|
|
x0 = 1.-xsi/a ;
|
|
if(x0 < -1.) x0 = -1. ;
|
|
|
|
if(x0 == -1.)
|
|
{
|
|
// 1 model fuction only
|
|
// in order to have xmean1 > xmeanth -> qprob < 1
|
|
if((1.-1./a) < xmeanth)
|
|
a = 1./(1.-xmeanth) ;
|
|
|
|
if(a*(1.-x0) < amax)
|
|
ea = exp(-a*(1.-x0)) ;
|
|
else
|
|
ea = 0. ;
|
|
eaa = 1.-ea ;
|
|
xmean1 = 1.-1./a+(1.-x0)*ea/eaa ;
|
|
|
|
c = 2. ;
|
|
b1 = b+1. ;
|
|
bx = b1 ;
|
|
eb1 = b1 ;
|
|
ebx = b1 ;
|
|
xmean2 = 0. ;
|
|
|
|
prob = 1. ;
|
|
qprob = xmeanth/xmean1 ;
|
|
}
|
|
else
|
|
{
|
|
// 2 model fuctions
|
|
// in order to have xmean1 > xmeanth
|
|
if((1.-x1fac2/a) < xmeanth)
|
|
{
|
|
a = x1fac3*x1fac2/(1.-xmeanth) ;
|
|
if(a*(1.-x0) < amax)
|
|
ea = exp(-a*(1.-x0)) ;
|
|
else
|
|
ea = 0. ;
|
|
eaa = 1.-ea ;
|
|
xmean1 = 1.-1./a+(1.-x0)*ea/eaa ;
|
|
}
|
|
else
|
|
{
|
|
ea = x1fac1 ;
|
|
eaa = 1.-x1fac1 ;
|
|
xmean1 = 1.-x1fac2/a ;
|
|
}
|
|
|
|
// from continuity of the 1st derivatives
|
|
c = a*(b-x0) ;
|
|
if(a*tau < c0)
|
|
c = c0*(b-x0)/tau ;
|
|
|
|
if(c == 1.) c=1.000001 ;
|
|
if(c == 2.) c=2.000001 ;
|
|
if(c == 3.) c=3.000001 ;
|
|
|
|
b1 = b+1. ;
|
|
bx=b-x0 ;
|
|
eb1=exp((c-1.)*log(b1)) ;
|
|
ebx=exp((c-1.)*log(bx)) ;
|
|
xmean2 = (x0*eb1+ebx+(eb1*bx-b1*ebx)/(2.-c))/(eb1-ebx) ;
|
|
|
|
cnorm1 = a/eaa ;
|
|
f1x0 = cnorm1*exp(-a*(1.-x0)) ;
|
|
cnorm2 = (c-1.)*eb1*ebx/(eb1-ebx) ;
|
|
f2x0 = cnorm2/exp(c*log(b-x0)) ;
|
|
|
|
// from continuity at x=x0
|
|
prob = f2x0/(f1x0+f2x0) ;
|
|
// from xmean = xmeanth
|
|
qprob = (f1x0+f2x0)*xmeanth/(f2x0*xmean1+f1x0*xmean2) ;
|
|
}
|
|
|
|
// protection against prob or qprob > 1 and
|
|
// prob or qprob < 0
|
|
// ***************************************************************
|
|
if((qprob > 1.) || (qprob < 0.) || (prob > 1.) || (prob < 0.))
|
|
{
|
|
// this print possibility has been left intentionally
|
|
// for debugging purposes ..........................
|
|
G4bool pr = false ;
|
|
// pr = true ;
|
|
if(pr)
|
|
{
|
|
const G4double prlim = 0.10 ;
|
|
if((fabs((xmeanth-xmean2)/(xmean1-xmean2)-prob)/prob > prlim) ||
|
|
((xmeanth-xmean2)/(xmean1-xmean2) > 1.) ||
|
|
((xmeanth-xmean2)/(xmean1-xmean2) < 0.) )
|
|
{
|
|
G4cout.precision(5) ;
|
|
G4cout << "\nparticle=" << aParticle->GetDefinition()->
|
|
GetParticleName() << " in material "
|
|
<< trackData.GetMaterialCutsCouple()->
|
|
GetMaterial()->GetName() << " with kinetic energy "
|
|
<< KineticEnergy << " MeV," << G4endl ;
|
|
G4cout << " step length="
|
|
<< truestep << " mm" << G4endl ;
|
|
G4cout << "p=" << prob << " q=" << qprob << " -----> "
|
|
<< "p=" << (xmeanth-xmean2)/(xmean1-xmean2)
|
|
<< " q=" << 1. << G4endl ;
|
|
}
|
|
}
|
|
qprob = 1. ;
|
|
prob = (xmeanth-xmean2)/(xmean1-xmean2) ;
|
|
}
|
|
// **************************************************************
|
|
|
|
// sampling of costheta
|
|
if(G4UniformRand() < qprob)
|
|
{
|
|
if(G4UniformRand() < prob)
|
|
cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
|
|
else
|
|
cth = b-b1*bx/exp(log(ebx-G4UniformRand()*(ebx-eb1))/(c-1.)) ;
|
|
}
|
|
else
|
|
cth = -1.+2.*G4UniformRand() ;
|
|
}
|
|
}
|
|
|
|
G4double sth = sqrt(1.-cth*cth);
|
|
G4double phi = twopi*G4UniformRand();
|
|
G4double dirx = sth*cos(phi), diry = sth*sin(phi), dirz = cth;
|
|
|
|
G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection();
|
|
|
|
G4ThreeVector newDirection(dirx,diry,dirz);
|
|
newDirection.rotateUz(ParticleDirection);
|
|
fParticleChange.ProposeMomentumDirection(newDirection.x(),
|
|
newDirection.y(),
|
|
newDirection.z());
|
|
|
|
if (fLatDisplFlag)
|
|
{
|
|
// compute mean lateral displacement, only for safety > tolerance !
|
|
G4double safetyminustolerance = stepData.GetPostStepPoint()->GetSafety();
|
|
G4double rmean, etau;
|
|
|
|
if (safetyminustolerance > 0.)
|
|
{
|
|
if (tau < tausmall) rmean = 0.;
|
|
else if(tau < taulim) rmean = kappa*tau*tau*tau*(1.-kappapl1*tau/4.)/6.;
|
|
else
|
|
{
|
|
if(tau<taubig) etau = exp(-tau);
|
|
else etau = 0.;
|
|
rmean = -kappa*tau;
|
|
rmean = -exp(rmean)/(kappa*kappami1);
|
|
rmean += tau-kappapl1/kappa+kappa*etau/kappami1;
|
|
}
|
|
|
|
if (rmean>0.) rmean = 2.*lambda0*sqrt(rmean/3.);
|
|
else rmean = 0.;
|
|
|
|
// for rmean > 0) only
|
|
if (rmean > 0.)
|
|
{
|
|
if (rmean>safetyminustolerance) rmean = safetyminustolerance;
|
|
|
|
// sample direction of lateral displacement
|
|
phi = twopi*G4UniformRand();
|
|
dirx = cos(phi); diry = sin(phi); dirz = 0.;
|
|
|
|
G4ThreeVector latDirection(dirx,diry,dirz);
|
|
latDirection.rotateUz(ParticleDirection);
|
|
|
|
// compute new endpoint of the Step
|
|
G4ThreeVector newPosition = stepData.GetPostStepPoint()->GetPosition()
|
|
+ rmean*latDirection;
|
|
|
|
G4Navigator* navigator =
|
|
G4TransportationManager::GetTransportationManager()
|
|
->GetNavigatorForTracking();
|
|
navigator->LocateGlobalPointWithinVolume(newPosition);
|
|
|
|
fParticleChange.ProposePosition(newPosition);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return &fParticleChange;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void G4MultipleScattering52::PrintInfoDefinition()
|
|
{
|
|
G4String comments = " Tables of transport mean free paths.";
|
|
comments += "\n New model of MSC , computes the lateral \n";
|
|
comments += " displacement of the particle , too.";
|
|
|
|
G4cout << G4endl << GetProcessName() << ": " << comments
|
|
<< "\n PhysicsTables from "
|
|
<< G4BestUnit(LowestKineticEnergy ,"Energy")
|
|
<< " to " << G4BestUnit(HighestKineticEnergy,"Energy")
|
|
<< " in " << TotBin << " bins. \n";
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|