Verbose settings.
This commit is contained in:
paraiso
@ -1,951 +0,0 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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......
|
||||
@ -1,264 +0,0 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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. *
|
||||
// ********************************************************************
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
// GEANT 4 class header file
|
||||
//
|
||||
// History:
|
||||
// 17 August 2004 P.Gumplinger and T.MacPhail
|
||||
// samples Michel spectrum including 1st order
|
||||
// radiative corrections
|
||||
// Reference: Florian Scheck "Muon Physics", in Physics Reports
|
||||
// (Review Section of Physics Letters) 44, No. 4 (1978)
|
||||
// 187-248. North-Holland Publishing Company, Amsterdam
|
||||
// at page 210 cc.
|
||||
//
|
||||
// W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
//
|
||||
#include "G4MuonDecayChannelWithSpin.hh"
|
||||
|
||||
#include "Randomize.hh"
|
||||
#include "G4DecayProducts.hh"
|
||||
#include "G4LorentzVector.hh"
|
||||
|
||||
G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4String& theParentName,
|
||||
G4double theBR)
|
||||
: G4MuonDecayChannel(theParentName,theBR)
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin()
|
||||
{
|
||||
}
|
||||
|
||||
G4DecayProducts *G4MuonDecayChannelWithSpin::DecayIt(G4double)
|
||||
{
|
||||
// This version assumes V-A coupling with 1st order radiative correctons,
|
||||
// the standard model Michel parameter values, but
|
||||
// gives incorrect energy spectrum for neutrinos
|
||||
|
||||
#ifdef G4VERBOSE
|
||||
if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
|
||||
#endif
|
||||
|
||||
if (parent == 0) FillParent();
|
||||
if (daughters == 0) FillDaughters();
|
||||
|
||||
// parent mass
|
||||
G4double parentmass = parent->GetPDGMass();
|
||||
|
||||
EMMU = parentmass;
|
||||
|
||||
//daughters'mass
|
||||
G4double daughtermass[3];
|
||||
G4double sumofdaughtermass = 0.0;
|
||||
for (G4int index=0; index<3; index++){
|
||||
daughtermass[index] = daughters[index]->GetPDGMass();
|
||||
sumofdaughtermass += daughtermass[index];
|
||||
}
|
||||
|
||||
EMASS = daughtermass[0];
|
||||
|
||||
//create parent G4DynamicParticle at rest
|
||||
G4ThreeVector dummy;
|
||||
G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
|
||||
//create G4Decayproducts
|
||||
G4DecayProducts *products = new G4DecayProducts(*parentparticle);
|
||||
delete parentparticle;
|
||||
|
||||
// calcurate electron energy
|
||||
|
||||
G4double michel_rho = 0.75; //Standard Model Michel rho
|
||||
G4double michel_delta = 0.75; //Standard Model Michel delta
|
||||
G4double michel_xsi = 1.00; //Standard Model Michel xsi
|
||||
G4double michel_eta = 0.00; //Standard Model eta
|
||||
|
||||
G4double rndm, x, ctheta;
|
||||
|
||||
G4double FG;
|
||||
G4double FG_max = 2.00;
|
||||
|
||||
G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
|
||||
G4double x0 = EMASS/W_mue;
|
||||
|
||||
G4double x0_squared = x0*x0;
|
||||
|
||||
// ***************************************************
|
||||
// x0 <= x <= 1. and -1 <= y <= 1
|
||||
//
|
||||
// F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y
|
||||
// ***************************************************
|
||||
|
||||
// ***** sampling F(x,y) directly (brute force) *****
|
||||
|
||||
do{
|
||||
|
||||
// Sample the positron energy by sampling from F
|
||||
|
||||
rndm = G4UniformRand();
|
||||
|
||||
x = x0 + rndm*(1.-x0);
|
||||
|
||||
G4double x_squared = x*x;
|
||||
|
||||
G4double F_IS, F_AS, G_IS, G_AS;
|
||||
|
||||
F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
|
||||
F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
|
||||
|
||||
G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
|
||||
G_IS = G_IS + michel_eta*(1.-x)*x0;
|
||||
|
||||
G_AS = 3.*(michel_xsi-1.)*(1.-x);
|
||||
G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
|
||||
G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
|
||||
|
||||
F_IS = F_IS + G_IS;
|
||||
F_AS = F_AS + G_AS;
|
||||
|
||||
// *** Radiative Corrections ***
|
||||
|
||||
G4double R_IS = F_c(x,x0);
|
||||
|
||||
G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
|
||||
|
||||
// *** Radiative Corrections ***
|
||||
|
||||
G4double R_AS = F_theta(x,x0);
|
||||
|
||||
rndm = G4UniformRand();
|
||||
|
||||
ctheta = 2.*rndm-1.;
|
||||
|
||||
G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
|
||||
|
||||
FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
|
||||
|
||||
if(FG>FG_max){
|
||||
G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl;
|
||||
FG_max = FG;
|
||||
}
|
||||
|
||||
rndm = G4UniformRand();
|
||||
|
||||
}while(FG<rndm*FG_max);
|
||||
|
||||
G4double energy = x * W_mue;
|
||||
|
||||
rndm = G4UniformRand();
|
||||
|
||||
G4double phi = twopi * rndm;
|
||||
|
||||
if(energy < EMASS) energy = EMASS;
|
||||
|
||||
// calculate daughter momentum
|
||||
G4double daughtermomentum[3];
|
||||
|
||||
daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
|
||||
|
||||
G4double stheta = std::sqrt(1.-ctheta*ctheta);
|
||||
G4double cphi = std::cos(phi);
|
||||
G4double sphi = std::sin(phi);
|
||||
|
||||
//Coordinates of the decay positron with respect to the muon spin
|
||||
|
||||
G4double px = stheta*cphi;
|
||||
G4double py = stheta*sphi;
|
||||
G4double pz = ctheta;
|
||||
|
||||
G4ThreeVector direction0(px,py,pz);
|
||||
|
||||
direction0.rotateUz(parent_polarization);
|
||||
|
||||
G4DynamicParticle * daughterparticle0
|
||||
= new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
|
||||
|
||||
products->PushProducts(daughterparticle0);
|
||||
|
||||
|
||||
// daughter 1 ,2 (neutrinos)
|
||||
// create neutrinos in the C.M frame of two neutrinos
|
||||
G4double energy2 = parentmass*(1.0 - x/2.0);
|
||||
G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
|
||||
G4double beta = -1.0*daughtermomentum[0]/energy2;
|
||||
G4double costhetan = 2.*G4UniformRand()-1.0;
|
||||
G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
|
||||
G4double phin = twopi*G4UniformRand()*rad;
|
||||
G4double sinphin = std::sin(phin);
|
||||
G4double cosphin = std::cos(phin);
|
||||
|
||||
G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
|
||||
G4DynamicParticle * daughterparticle1
|
||||
= new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
|
||||
G4DynamicParticle * daughterparticle2
|
||||
= new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
|
||||
|
||||
// boost to the muon rest frame
|
||||
G4LorentzVector p4;
|
||||
p4 = daughterparticle1->Get4Momentum();
|
||||
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
|
||||
daughterparticle1->Set4Momentum(p4);
|
||||
p4 = daughterparticle2->Get4Momentum();
|
||||
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
|
||||
daughterparticle2->Set4Momentum(p4);
|
||||
products->PushProducts(daughterparticle1);
|
||||
products->PushProducts(daughterparticle2);
|
||||
daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
|
||||
daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
|
||||
|
||||
// output message
|
||||
#ifdef G4VERBOSE
|
||||
if (GetVerboseLevel()>1) {
|
||||
G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
|
||||
G4cout << " create decay products in rest frame " <<G4endl;
|
||||
products->DumpInfo();
|
||||
}
|
||||
#endif
|
||||
return products;
|
||||
}
|
||||
|
||||
G4double G4MuonDecayChannelWithSpin::R_c(G4double x){
|
||||
|
||||
G4int n_max = (int)(100.*x);
|
||||
|
||||
if(n_max<10)n_max=10;
|
||||
|
||||
G4double L2 = 0.0;
|
||||
|
||||
for(G4int n=1; n<=n_max; n++){
|
||||
L2 += std::pow(x,n)/(n*n);
|
||||
}
|
||||
|
||||
G4double omega = std::log(EMMU/EMASS);
|
||||
|
||||
G4double r_c;
|
||||
|
||||
r_c = 2.*L2-(pi*pi/3.)-2.;
|
||||
r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x));
|
||||
r_c = r_c - std::log(x)*(2.*std::log(x)-1.);
|
||||
r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x);
|
||||
|
||||
return r_c;
|
||||
}
|
||||
@ -1,122 +0,0 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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. *
|
||||
// ********************************************************************
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
// GEANT 4 class header file
|
||||
//
|
||||
// History:
|
||||
// 17 August 2004 P.Gumplinger and T.MacPhail
|
||||
// samples Michel spectrum including 1st order
|
||||
// radiative corrections
|
||||
// Reference: Florian Scheck "Muon Physics", in Physics Reports
|
||||
// (Review Section of Physics Letters) 44, No. 4 (1978)
|
||||
// 187-248. North-Holland Publishing Company, Amsterdam
|
||||
// at page 210 cc.
|
||||
//
|
||||
// W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
#ifndef G4MuonDecayChannelWithSpin_hh
|
||||
#define G4MuonDecayChannelWithSpin_hh 1
|
||||
|
||||
#include "globals.hh"
|
||||
#include "G4ThreeVector.hh"
|
||||
#include "G4MuonDecayChannel.hh"
|
||||
|
||||
class G4MuonDecayChannelWithSpin : public G4MuonDecayChannel
|
||||
{
|
||||
// Class Decription
|
||||
// This class describes muon decay kinemtics.
|
||||
// This version assumes V-A coupling with 1st order radiative correctons,
|
||||
// the standard model Michel parameter values, but
|
||||
// gives incorrect energy spectrum for neutrinos
|
||||
|
||||
public: // With Description
|
||||
|
||||
//Constructors
|
||||
G4MuonDecayChannelWithSpin(const G4String& theParentName,
|
||||
G4double theBR);
|
||||
// Destructor
|
||||
virtual ~G4MuonDecayChannelWithSpin();
|
||||
|
||||
public: // With Description
|
||||
|
||||
virtual G4DecayProducts *DecayIt(G4double);
|
||||
|
||||
void SetPolarization(G4ThreeVector);
|
||||
const G4ThreeVector& GetPolarization() const;
|
||||
|
||||
private:
|
||||
|
||||
G4ThreeVector parent_polarization;
|
||||
|
||||
// Radiative Correction Factors
|
||||
|
||||
G4double F_c(G4double x, G4double x0);
|
||||
G4double F_theta(G4double x, G4double x0);
|
||||
G4double R_c(G4double x);
|
||||
|
||||
G4double EMMU;
|
||||
G4double EMASS;
|
||||
|
||||
};
|
||||
|
||||
inline void G4MuonDecayChannelWithSpin::SetPolarization(G4ThreeVector polar)
|
||||
{
|
||||
parent_polarization = polar;
|
||||
}
|
||||
|
||||
inline const G4ThreeVector& G4MuonDecayChannelWithSpin::GetPolarization() const
|
||||
{
|
||||
return parent_polarization;
|
||||
}
|
||||
|
||||
inline G4double G4MuonDecayChannelWithSpin::F_c(G4double x, G4double x0)
|
||||
{
|
||||
G4double omega = std::log(EMMU/EMASS);
|
||||
|
||||
G4double f_c;
|
||||
|
||||
f_c = (5.+17.*x-34.*x*x)*(omega+std::log(x))-22.*x+34.*x*x;
|
||||
f_c = (1.-x)/(3.*x*x)*f_c;
|
||||
f_c = (6.-4.*x)*R_c(x)+(6.-6.*x)*std::log(x) + f_c;
|
||||
f_c = (fine_structure_const/twopi) * (x*x-x0*x0) * f_c;
|
||||
|
||||
return f_c;
|
||||
}
|
||||
|
||||
inline G4double G4MuonDecayChannelWithSpin::F_theta(G4double x, G4double x0)
|
||||
{
|
||||
G4double omega = std::log(EMMU/EMASS);
|
||||
|
||||
G4double f_theta;
|
||||
|
||||
f_theta = (1.+x+34*x*x)*(omega+std::log(x))+3.-7.*x-32.*x*x;
|
||||
f_theta = f_theta + ((4.*(1.-x)*(1.-x))/x)*std::log(1.-x);
|
||||
f_theta = (1.-x)/(3.*x*x) * f_theta;
|
||||
f_theta = (2.-4.*x)*R_c(x)+(2.-6.*x)*std::log(x)-f_theta;
|
||||
f_theta = (fine_structure_const/twopi) * (x*x-x0*x0) * f_theta;
|
||||
|
||||
return f_theta;
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -1,119 +0,0 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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: G4Muonium.cc,v 1.8 2003/06/16 16:57:55 gunter Exp $
|
||||
// GEANT4 tag $Name: geant4-07-00-cand-01 $
|
||||
//
|
||||
//
|
||||
// ----------------------------------------------------------------------
|
||||
// GEANT 4 class implementation file
|
||||
//
|
||||
// History: first implementation, based on object model of
|
||||
// 4th April 1996, G.Cosmo
|
||||
// **********************************************************************
|
||||
// Added particle definitions, H.Kurashige, 19 April 1996
|
||||
// Added SetCuts implementation, L.Urban, 10 May 1996
|
||||
// Operators (+=, *=, ++, -> etc.) correctly used, P. Urban, 26/6/96
|
||||
// Add MuonPlusDefinition(), H.Kurashige 4 July 1996
|
||||
// Add MuoniumDefinition(),T.K.Paraiso 7 April 2005
|
||||
// ----------------------------------------------------------------------
|
||||
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
|
||||
#include "G4Muonium.hh"
|
||||
|
||||
#include "G4MuonDecayChannel.hh"
|
||||
#include "G4DecayTable.hh"
|
||||
#include "Randomize.hh"
|
||||
|
||||
// ######################################################################
|
||||
// ### MUONPLUS ###
|
||||
// ######################################################################
|
||||
|
||||
G4Muonium::G4Muonium(
|
||||
const G4String& aName, G4double mass,
|
||||
G4double width, G4double charge,
|
||||
G4int iSpin, G4int iParity,
|
||||
G4int iConjugation, G4int iIsospin,
|
||||
G4int iIsospin3, G4int gParity,
|
||||
const G4String& pType, G4int lepton,
|
||||
G4int baryon, G4int encoding,
|
||||
G4bool stable, G4double lifetime,
|
||||
G4DecayTable *decaytable )
|
||||
: G4ParticleDefinition( aName,mass,width,charge,iSpin,iParity,
|
||||
iConjugation,iIsospin,iIsospin3,gParity,pType,
|
||||
lepton,baryon,encoding,stable,lifetime,decaytable )
|
||||
{
|
||||
SetParticleSubType("mu");
|
||||
SetPDGStable(false);
|
||||
SetMuoniumGamma();
|
||||
|
||||
//create Decay Table
|
||||
G4DecayTable* table = GetDecayTable();
|
||||
if (table!=NULL) delete table;
|
||||
table = new G4DecayTable();
|
||||
|
||||
SetDecayTable(table);
|
||||
|
||||
|
||||
}
|
||||
|
||||
// ......................................................................
|
||||
// ... static member definitions ...
|
||||
// ......................................................................
|
||||
//
|
||||
// Arguments for constructor are as follows
|
||||
// name mass width charge
|
||||
// 2*spin parity C-conjugation
|
||||
// 2*Isospin 2*Isospin3 G-parity
|
||||
// type lepton number baryon number PDG encoding
|
||||
// stable lifetime decay table
|
||||
G4Muonium G4Muonium::theMuonium(
|
||||
"Mu", 0.1056584*GeV, 2.99591e-16*MeV, 0.*eplus,
|
||||
1, 0, 0,
|
||||
0, 0, 0,
|
||||
"lepton", -1, 0, -1313,
|
||||
true, 2197.03*ns, NULL
|
||||
);
|
||||
|
||||
G4Muonium* G4Muonium::MuoniumDefinition() {return &theMuonium;}
|
||||
|
||||
G4Muonium* G4Muonium::Muonium()
|
||||
{
|
||||
return &theMuonium;
|
||||
}
|
||||
|
||||
|
||||
void G4Muonium::SetMuoniumGamma()
|
||||
{
|
||||
G4double gamma;
|
||||
gamma = 0.5*((1.*eplus)/(0.1056584*GeV/(c_light*c_light))-(1.*eplus)/(0.51099906*MeV/(c_light*c_light)));
|
||||
|
||||
// SetGyromagneticRatio(gamma);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,90 +0,0 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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: G4Muonium.hh,v 1.7 2001/10/16 08:16:15 kurasige Exp $
|
||||
// GEANT4 tag $Name: geant4-07-00-cand-01 $
|
||||
//
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
// GEANT 4 class header file
|
||||
//
|
||||
// History: first implementation, based on object model of
|
||||
// 4-th April 1996, G.Cosmo
|
||||
// ****************************************************************
|
||||
// Added particle definitions, H.Kurashige, 19 April 1996
|
||||
// Added SetCuts, L.Urban, 12 June 1996
|
||||
// Added not static GetEnergyCuts() and GetLengthCuts(), G.Cosmo, 11 July 1996
|
||||
// Muonium defined as a Mu, T.K.Paraiso 07-04-05
|
||||
// ----------------------------------------------------------------
|
||||
|
||||
// Each class inheriting from G4VLepton
|
||||
// corresponds to a particle type; one and only one
|
||||
// instance for each class is guaranteed.
|
||||
|
||||
#ifndef G4Muonium_h
|
||||
#define G4Muonium_h 1
|
||||
|
||||
#include "globals.hh"
|
||||
#include "G4ios.hh"
|
||||
#include "G4ParticleDefinition.hh"
|
||||
|
||||
// ######################################################################
|
||||
// ### MUONPLUS ###
|
||||
// ######################################################################
|
||||
|
||||
class G4Muonium : public G4ParticleDefinition
|
||||
{
|
||||
private: // constructors are hide as private
|
||||
G4Muonium(
|
||||
const G4String& aName, G4double mass,
|
||||
G4double width, G4double charge,
|
||||
G4int iSpin, G4int iParity,
|
||||
G4int iConjugation, G4int iIsospin,
|
||||
G4int iIsospin3, G4int gParity,
|
||||
const G4String& pType, G4int lepton,
|
||||
G4int baryon, G4int encoding,
|
||||
G4bool stable, G4double lifetime,
|
||||
G4DecayTable *decaytable
|
||||
);
|
||||
|
||||
public:
|
||||
virtual ~G4Muonium(){;};
|
||||
|
||||
void SetMuoniumGamma();
|
||||
|
||||
static G4Muonium* MuoniumDefinition();
|
||||
static G4Muonium* Muonium();
|
||||
|
||||
private:
|
||||
static G4Muonium theMuonium;
|
||||
|
||||
};
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,460 +0,0 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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: G4VDecayChannel.cc,v 1.16 2004/12/02 08:09:00 kurasige Exp $
|
||||
// GEANT4 tag $Name: geant4-07-00-cand-03 $
|
||||
//
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
// GEANT 4 class header file
|
||||
//
|
||||
// History: first implementation, based on object model of
|
||||
// 27 July 1996 H.Kurashige
|
||||
// 30 May 1997 H.Kurashige
|
||||
// 23 Mar. 2000 H.Weber : add GetAngularMomentum
|
||||
// ------------------------------------------------------------
|
||||
|
||||
#include "G4ParticleDefinition.hh"
|
||||
#include "G4ParticleTable.hh"
|
||||
#include "G4DecayTable.hh"
|
||||
#include "G4DecayProducts.hh"
|
||||
#include "G4VDecayChannel.hh"
|
||||
|
||||
const G4String G4VDecayChannel::noName = " ";
|
||||
|
||||
G4VDecayChannel::G4VDecayChannel(const G4String &aName, G4int Verbose)
|
||||
:kinematics_name(aName),
|
||||
rbranch(0.0),
|
||||
numberOfDaughters(0),
|
||||
parent_name(0), daughters_name(0),
|
||||
particletable(0),
|
||||
parent(0), daughters(0),
|
||||
parent_mass(0.0), daughters_mass(0),
|
||||
verboseLevel(Verbose)
|
||||
{
|
||||
// set pointer to G4ParticleTable (static and singleton object)
|
||||
particletable = G4ParticleTable::GetParticleTable();
|
||||
}
|
||||
|
||||
G4VDecayChannel::G4VDecayChannel(const G4String &aName,
|
||||
const G4String& theParentName,
|
||||
G4double theBR,
|
||||
G4int theNumberOfDaughters,
|
||||
const G4String& theDaughterName1,
|
||||
const G4String& theDaughterName2,
|
||||
const G4String& theDaughterName3,
|
||||
const G4String& theDaughterName4 )
|
||||
:kinematics_name(aName),
|
||||
rbranch(theBR),
|
||||
numberOfDaughters(theNumberOfDaughters),
|
||||
parent_name(0), daughters_name(0),
|
||||
particletable(0),
|
||||
parent(0), daughters(0),
|
||||
parent_mass(0.0), daughters_mass(0),
|
||||
verboseLevel(1)
|
||||
{
|
||||
// set pointer to G4ParticleTable (static and singleton object)
|
||||
particletable = G4ParticleTable::GetParticleTable();
|
||||
|
||||
// parent name
|
||||
parent_name = new G4String(theParentName);
|
||||
|
||||
// cleate array
|
||||
daughters_name = new G4String*[numberOfDaughters];
|
||||
for (G4int index=0;index<numberOfDaughters;index++) daughters_name[index]=0;
|
||||
|
||||
// daughters' name
|
||||
if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
|
||||
if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
|
||||
if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
|
||||
if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
|
||||
}
|
||||
|
||||
|
||||
|
||||
G4VDecayChannel::G4VDecayChannel(const G4VDecayChannel &right)
|
||||
{
|
||||
kinematics_name = right.kinematics_name;
|
||||
verboseLevel = right.verboseLevel;
|
||||
rbranch = right.rbranch;
|
||||
|
||||
// copy parent name
|
||||
parent_name = new G4String(*right.parent_name);
|
||||
parent = 0;
|
||||
parent_mass = 0.0;
|
||||
|
||||
//create array
|
||||
numberOfDaughters = right.numberOfDaughters;
|
||||
|
||||
if ( numberOfDaughters >0 ) {
|
||||
daughters_name = new G4String*[numberOfDaughters];
|
||||
//copy daughters name
|
||||
for (G4int index=0; index < numberOfDaughters; index++)
|
||||
{
|
||||
daughters_name[index] = new G4String(*right.daughters_name[index]);
|
||||
}
|
||||
}
|
||||
|
||||
//
|
||||
daughters_mass = 0;
|
||||
daughters = 0;
|
||||
|
||||
// particle table
|
||||
particletable = G4ParticleTable::GetParticleTable();
|
||||
}
|
||||
|
||||
G4VDecayChannel & G4VDecayChannel::operator=(const G4VDecayChannel &right)
|
||||
{
|
||||
if (this != &right) {
|
||||
kinematics_name = right.kinematics_name;
|
||||
verboseLevel = right.verboseLevel;
|
||||
rbranch = right.rbranch;
|
||||
|
||||
// copy parent name
|
||||
parent_name = new G4String(*right.parent_name);
|
||||
|
||||
// clear daughters_name array
|
||||
ClearDaughtersName();
|
||||
|
||||
// recreate array
|
||||
numberOfDaughters = right.numberOfDaughters;
|
||||
if ( numberOfDaughters >0 ) {
|
||||
daughters_name = new G4String*[numberOfDaughters];
|
||||
//copy daughters name
|
||||
for (G4int index=0; index < numberOfDaughters; index++) {
|
||||
daughters_name[index] = new G4String(*right.daughters_name[index]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//
|
||||
parent = 0;
|
||||
daughters = 0;
|
||||
parent_mass = 0.0;
|
||||
daughters_mass = 0;
|
||||
|
||||
// particle table
|
||||
particletable = G4ParticleTable::GetParticleTable();
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
G4VDecayChannel::~G4VDecayChannel()
|
||||
{
|
||||
if (parent_name != 0) delete parent_name;
|
||||
ClearDaughtersName();
|
||||
if (daughters_mass != 0) delete [] daughters_mass;
|
||||
}
|
||||
|
||||
void G4VDecayChannel::ClearDaughtersName()
|
||||
{
|
||||
if ( daughters_name != 0) {
|
||||
if (numberOfDaughters>0) {
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>1) {
|
||||
G4cout << "G4VDecayChannel::ClearDaughtersName ";
|
||||
G4cout << "clear all daughters " << G4endl;
|
||||
}
|
||||
#endif
|
||||
for (G4int index=0; index < numberOfDaughters; index++) {
|
||||
if (daughters_name[index] != 0) delete daughters_name[index];
|
||||
}
|
||||
}
|
||||
delete [] daughters_name;
|
||||
daughters_name = 0;
|
||||
}
|
||||
//
|
||||
if (daughters != 0) delete [] daughters;
|
||||
if (daughters_mass != 0) delete [] daughters_mass;
|
||||
daughters = 0;
|
||||
daughters_mass = 0;
|
||||
|
||||
numberOfDaughters = 0;
|
||||
}
|
||||
|
||||
void G4VDecayChannel::SetNumberOfDaughters(G4int size)
|
||||
{
|
||||
if (size >0) {
|
||||
// remove old contents
|
||||
ClearDaughtersName();
|
||||
// cleate array
|
||||
daughters_name = new G4String*[size];
|
||||
for (G4int index=0;index<size;index++) daughters_name[index]=0;
|
||||
numberOfDaughters = size;
|
||||
}
|
||||
}
|
||||
|
||||
void G4VDecayChannel::SetDaughter(G4int anIndex,
|
||||
const G4String &particle_name)
|
||||
{
|
||||
// check numberOfDaughters is positive
|
||||
if (numberOfDaughters<=0) {
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>0) {
|
||||
G4cout << "G4VDecayChannel::SetDaughter: ";
|
||||
G4cout << "Number of daughters is not defined" << G4endl;
|
||||
}
|
||||
#endif
|
||||
return;
|
||||
}
|
||||
|
||||
// check existence of daughters_name array
|
||||
if (daughters_name == 0) {
|
||||
// cleate array
|
||||
daughters_name = new G4String*[numberOfDaughters];
|
||||
for (G4int index=0;index<numberOfDaughters;index++) {
|
||||
daughters_name[index]=0;
|
||||
}
|
||||
}
|
||||
|
||||
// check an index
|
||||
if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>0) {
|
||||
G4cout << "G4VDecayChannel::SetDaughter";
|
||||
G4cout << "index out of range " << anIndex << G4endl;
|
||||
}
|
||||
#endif
|
||||
} else {
|
||||
// delete the old name if it exists
|
||||
if (daughters_name[anIndex]!=0) delete daughters_name[anIndex];
|
||||
// fill the name
|
||||
daughters_name[anIndex] = new G4String(particle_name);
|
||||
// refill the array of daughters[] if it exists
|
||||
if (daughters != 0) FillDaughters();
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>1) {
|
||||
G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
|
||||
G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
|
||||
{
|
||||
if (parent_type != 0) SetDaughter(anIndex, parent_type->GetParticleName());
|
||||
}
|
||||
|
||||
void G4VDecayChannel::FillDaughters()
|
||||
{
|
||||
G4int index;
|
||||
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
|
||||
#endif
|
||||
if (daughters != 0) delete [] daughters;
|
||||
|
||||
// parent mass
|
||||
if (parent == 0) FillParent();
|
||||
G4double parentmass = parent->GetPDGMass();
|
||||
|
||||
//
|
||||
G4double sumofdaughtermass = 0.0;
|
||||
if ((numberOfDaughters <=0) || (daughters_name == 0) ){
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>0) {
|
||||
G4cout << "G4VDecayChannel::FillDaughters ";
|
||||
G4cout << "[ " << parent->GetParticleName() << " ]";
|
||||
G4cout << "numberOfDaughters is not defined yet";
|
||||
}
|
||||
#endif
|
||||
daughters = 0;
|
||||
G4Exception("G4VDecayChannel::FillDaughters");
|
||||
}
|
||||
|
||||
//create and set the array of pointers to daughter particles
|
||||
daughters = new G4ParticleDefinition*[numberOfDaughters];
|
||||
if (daughters_mass != 0) delete [] daughters_mass;
|
||||
daughters_mass = new G4double[numberOfDaughters];
|
||||
// loop over all daughters
|
||||
for (index=0; index < numberOfDaughters; index++) {
|
||||
if (daughters_name[index] == 0) {
|
||||
// daughter name is not defined
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>0) {
|
||||
G4cout << "G4VDecayChannel::FillDaughters ";
|
||||
G4cout << "[ " << parent->GetParticleName() << " ]";
|
||||
G4cout << index << "-th daughter is not defined yet" << G4endl;
|
||||
}
|
||||
#endif
|
||||
daughters[index] = 0;
|
||||
G4Exception("G4VDecayChannel::FillDaughters");
|
||||
}
|
||||
//search daughter particles in the particle table
|
||||
daughters[index] = particletable->FindParticle(*daughters_name[index]);
|
||||
if (daughters[index] == 0) {
|
||||
// can not find the daughter particle
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>0) {
|
||||
G4cout << "G4VDecayChannel::FillDaughters ";
|
||||
G4cout << "[ " << parent->GetParticleName() << " ]";
|
||||
G4cout << index << ":" << *daughters_name[index];
|
||||
G4cout << " is not defined !!" << G4endl;
|
||||
G4cout << " The BR of this decay mode is set to zero " << G4endl;
|
||||
}
|
||||
#endif
|
||||
SetBR(0.0);
|
||||
return;
|
||||
}
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>1) {
|
||||
G4cout << index << ":" << *daughters_name[index];
|
||||
G4cout << ":" << daughters[index] << G4endl;
|
||||
}
|
||||
#endif
|
||||
daughters_mass[index] = daughters[index]->GetPDGMass();
|
||||
sumofdaughtermass += daughters[index]->GetPDGMass();
|
||||
} // end loop over all daughters
|
||||
|
||||
// check sum of daghter mass
|
||||
G4double widthMass = parent->GetPDGWidth();
|
||||
if ( (parent->GetParticleType() != "nucleus") &&
|
||||
(sumofdaughtermass > parentmass + 5*widthMass) ){
|
||||
// !!! illegal mass !!!
|
||||
#ifdef G4VERBOSE
|
||||
if (GetVerboseLevel()>0) {
|
||||
G4cout << "G4VDecayChannel::FillDaughters ";
|
||||
G4cout << "[ " << parent->GetParticleName() << " ]";
|
||||
G4cout << " Energy/Momentum conserevation breaks " <<G4endl;
|
||||
if (GetVerboseLevel()>1) {
|
||||
G4cout << " parent:" << *parent_name;
|
||||
G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
|
||||
for (index=0; index < numberOfDaughters; index++){
|
||||
G4cout << " daughter " << index << ":" << *daughters_name[index];
|
||||
G4cout << " mass:" << daughters[index]->GetPDGMass()/GeV;
|
||||
G4cout << "[GeV/c/c]" <<G4endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void G4VDecayChannel::FillParent()
|
||||
{
|
||||
if (parent_name == 0) {
|
||||
// parent name is not defined
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>0) {
|
||||
G4cout << "G4VDecayChannel::FillParent ";
|
||||
G4cout << ": parent name is not defined !!" << G4endl;
|
||||
}
|
||||
#endif
|
||||
parent = 0;
|
||||
G4Exception("G4VDecayChannel::FillParent");
|
||||
}
|
||||
// search parent particle in the particle table
|
||||
parent = particletable->FindParticle(*parent_name);
|
||||
if (parent == 0) {
|
||||
// parent particle does not exist
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>0) {
|
||||
G4cout << "G4VDecayChannel::FillParent ";
|
||||
G4cout << *parent_name << " does not exist !!" << G4endl;
|
||||
}
|
||||
#endif
|
||||
G4Exception("G4VDecayChannel::FillParent");
|
||||
}
|
||||
parent_mass = parent->GetPDGMass();
|
||||
}
|
||||
|
||||
void G4VDecayChannel::SetParent(const G4ParticleDefinition * parent_type)
|
||||
{
|
||||
if (parent_type != 0) SetParent(parent_type->GetParticleName());
|
||||
}
|
||||
|
||||
G4int G4VDecayChannel::GetAngularMomentum()
|
||||
{
|
||||
// determine angular momentum
|
||||
|
||||
// fill pointers to daughter particles if not yet set
|
||||
if (daughters == 0) FillDaughters();
|
||||
|
||||
const G4int PiSpin = parent->GetPDGiSpin();
|
||||
const G4int PParity = parent->GetPDGiParity();
|
||||
if (2==numberOfDaughters) { // up to now we can only handle two particle decays
|
||||
const G4int D1iSpin = daughters[0]->GetPDGiSpin();
|
||||
const G4int D1Parity = daughters[0]->GetPDGiParity();
|
||||
const G4int D2iSpin = daughters[1]->GetPDGiSpin();
|
||||
const G4int D2Parity = daughters[1]->GetPDGiParity();
|
||||
const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
|
||||
const G4int MaxiSpin = D1iSpin + D2iSpin;
|
||||
const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int
|
||||
G4int lMin;
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>1) {
|
||||
G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
|
||||
G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
|
||||
}
|
||||
#endif
|
||||
for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings
|
||||
lMin = std::abs(PiSpin-j)/2;
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>1)
|
||||
G4cout << "-> checking 2*j=" << j << G4endl;
|
||||
#endif
|
||||
for (G4int l=lMin; l<=lMax; l++) {
|
||||
#ifdef G4VERBOSE
|
||||
if (verboseLevel>1)
|
||||
G4cout << " checking l=" << l << G4endl;
|
||||
#endif
|
||||
if (l%2==0) {
|
||||
if (PParity == D1Parity*D2Parity) { // check parity for this l
|
||||
return l;
|
||||
}
|
||||
} else {
|
||||
if (PParity == -1*D1Parity*D2Parity) { // check parity for this l
|
||||
return l;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
G4Exception ("G4VDecayChannel::GetAngularMomentum: Sorry, can't handle 3 particle decays (up to now)");
|
||||
}
|
||||
G4Exception ("G4VDecayChannel::GetAngularMomentum: Can't find angular momentum for this decay!");
|
||||
return 0;
|
||||
}
|
||||
|
||||
void G4VDecayChannel::DumpInfo()
|
||||
{
|
||||
G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
|
||||
G4cout << " : " ;
|
||||
for (G4int index=0; index < numberOfDaughters; index++)
|
||||
{
|
||||
if(daughters_name[index] != 0) {
|
||||
G4cout << " " << *(daughters_name[index]);
|
||||
} else {
|
||||
G4cout << " not defined ";
|
||||
}
|
||||
}
|
||||
G4cout << G4endl;
|
||||
}
|
||||
|
||||
const G4String& G4VDecayChannel::GetNoName() const
|
||||
{
|
||||
return noName;
|
||||
}
|
||||
@ -1,288 +0,0 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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: G4VDecayChannel.hh,v 1.11 2004/12/02 08:08:58 kurasige Exp $
|
||||
// GEANT4 tag $Name: geant4-07-00-cand-03 $
|
||||
//
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
// GEANT 4 class header file
|
||||
//
|
||||
// History: first implementation, based on object model of
|
||||
// 27 July 1996 H.Kurashige
|
||||
// 30 May 1997 H.Kurashige
|
||||
// 23 Mar. 2000 H.Weber : add GetAngularMomentum()
|
||||
// ------------------------------------------------------------
|
||||
#ifndef G4VDecayChannel_h
|
||||
#define G4VDecayChannel_h 1
|
||||
|
||||
#include "G4ios.hh"
|
||||
#include "globals.hh"
|
||||
#include <cmath>
|
||||
#include "G4ThreeVector.hh"
|
||||
|
||||
class G4ParticleDefinition;
|
||||
class G4DecayProducts;
|
||||
class G4ParticleTable;
|
||||
|
||||
class G4VDecayChannel
|
||||
{
|
||||
// Class Description
|
||||
// This class is a abstract class to describe decay kinematics
|
||||
//
|
||||
|
||||
public:
|
||||
//Constructors
|
||||
G4VDecayChannel(const G4String &aName, G4int Verbose = 1);
|
||||
G4VDecayChannel(const G4String &aName,
|
||||
const G4String& theParentName,
|
||||
G4double theBR,
|
||||
G4int theNumberOfDaughters,
|
||||
const G4String& theDaughterName1,
|
||||
const G4String& theDaughterName2 = "",
|
||||
const G4String& theDaughterName3 = "",
|
||||
const G4String& theDaughterName4 = "" );
|
||||
|
||||
// Destructor
|
||||
virtual ~G4VDecayChannel();
|
||||
|
||||
private:
|
||||
// copy constructor and assignment operatotr
|
||||
G4VDecayChannel(const G4VDecayChannel &);
|
||||
G4VDecayChannel & operator=(const G4VDecayChannel &);
|
||||
|
||||
public:
|
||||
// equality operators
|
||||
G4int operator==(const G4VDecayChannel &right) const {return (this == &right);}
|
||||
G4int operator!=(const G4VDecayChannel &right) const {return (this != &right);}
|
||||
|
||||
// less-than operator is defined for G4DecayTable
|
||||
G4int operator<(const G4VDecayChannel &right) const;
|
||||
|
||||
public: // With Description
|
||||
virtual G4DecayProducts* DecayIt(G4double parentMass = -1.0) = 0;
|
||||
|
||||
public: // With Description
|
||||
//get kinematics name
|
||||
G4String GetKinematicsName() const;
|
||||
//get branching ratio
|
||||
G4double GetBR() const;
|
||||
//get number of daughter particles
|
||||
G4int GetNumberOfDaughters() const;
|
||||
|
||||
//get the pointer to the parent particle
|
||||
G4ParticleDefinition * GetParent();
|
||||
//get the pointer to a daughter particle
|
||||
G4ParticleDefinition * GetDaughter(G4int anIndex);
|
||||
|
||||
//get the angular momentum of the decay
|
||||
G4int GetAngularMomentum();
|
||||
//get the name of the parent particle
|
||||
const G4String& GetParentName() const;
|
||||
//get the name of a daughter particle
|
||||
const G4String& GetDaughterName(G4int anIndex) const;
|
||||
|
||||
// get mass of parent
|
||||
G4double GetParentMass() const;
|
||||
G4double GetDaughterMass(G4int anIndex) const;
|
||||
G4ThreeVector GetParentPolarization();
|
||||
//set the parent particle (by name or by pointer)
|
||||
void SetParent(const G4ParticleDefinition * particle_type);
|
||||
void SetParent(const G4String &particle_name);
|
||||
//set branching ratio
|
||||
void SetBR(G4double value);
|
||||
//set number of daughter particles
|
||||
void SetNumberOfDaughters(G4int value);
|
||||
//set a daughter particle (by name or by pointer)
|
||||
void SetDaughter(G4int anIndex,
|
||||
const G4ParticleDefinition * particle_type);
|
||||
void SetDaughter(G4int anIndex,
|
||||
const G4String &particle_name);
|
||||
void SetParentPolarization(G4ThreeVector polar);
|
||||
protected:
|
||||
// kinematics name
|
||||
G4String kinematics_name;
|
||||
// branching ratio [0.0 - 1.0]
|
||||
G4double rbranch;
|
||||
// number of daughters
|
||||
G4int numberOfDaughters;
|
||||
// parent particle
|
||||
G4String* parent_name;
|
||||
//daughter particles
|
||||
G4String** daughters_name;
|
||||
|
||||
protected: // With Description
|
||||
// celar daughters array
|
||||
void ClearDaughtersName();
|
||||
|
||||
protected:
|
||||
// pointer to particle table
|
||||
G4ParticleTable* particletable;
|
||||
|
||||
// temporary buffers of pointers to G4ParticleDefinition
|
||||
G4ParticleDefinition* parent;
|
||||
G4ParticleDefinition** daughters;
|
||||
|
||||
// parent mass
|
||||
G4double parent_mass;
|
||||
G4double* daughters_mass;
|
||||
|
||||
//parent theParentPolarization
|
||||
G4ThreeVector theParentPolarization;
|
||||
|
||||
// fill daughters array
|
||||
void FillDaughters();
|
||||
// fill parent
|
||||
void FillParent();
|
||||
|
||||
public: // With Description
|
||||
void SetVerboseLevel(G4int value);
|
||||
G4int GetVerboseLevel() const;
|
||||
void DumpInfo();
|
||||
|
||||
private:
|
||||
const G4String& GetNoName() const;
|
||||
|
||||
private:
|
||||
// controle flag for output message
|
||||
G4int verboseLevel;
|
||||
// 0: Silent
|
||||
// 1: Warning message
|
||||
// 2: More
|
||||
|
||||
static const G4String noName;
|
||||
};
|
||||
|
||||
inline
|
||||
G4int G4VDecayChannel::operator<(const G4VDecayChannel &right) const
|
||||
{
|
||||
return (this->rbranch < right.rbranch);
|
||||
}
|
||||
|
||||
inline
|
||||
G4ParticleDefinition* G4VDecayChannel::GetDaughter(G4int anIndex)
|
||||
{
|
||||
//pointers to daughter particles are filled, if they are not set yet
|
||||
if (daughters == 0) FillDaughters();
|
||||
|
||||
//get the pointer to a daughter particle
|
||||
if ( (anIndex>=0) && (anIndex<numberOfDaughters) ) {
|
||||
return daughters[anIndex];
|
||||
} else {
|
||||
if (verboseLevel>0)
|
||||
G4cout << "G4VDecayChannel::GetDaughter index out of range "<<anIndex<<G4endl;
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
inline
|
||||
const G4String& G4VDecayChannel::GetDaughterName(G4int anIndex) const
|
||||
{
|
||||
if ( (anIndex>=0) && (anIndex<numberOfDaughters) ) {
|
||||
return *daughters_name[anIndex];
|
||||
} else {
|
||||
if (verboseLevel>0){
|
||||
G4cout << "G4VDecayChannel::GetDaughterName ";
|
||||
G4cout << "index out of range " << anIndex << G4endl;
|
||||
}
|
||||
return GetNoName();
|
||||
}
|
||||
}
|
||||
|
||||
inline
|
||||
G4double G4VDecayChannel::GetDaughterMass(G4int anIndex) const
|
||||
{
|
||||
if ( (anIndex>=0) && (anIndex<numberOfDaughters) ) {
|
||||
return daughters_mass[anIndex];
|
||||
} else {
|
||||
if (verboseLevel>0){
|
||||
G4cout << "G4VDecayChannel::GetDaughterMass ";
|
||||
G4cout << "index out of range " << anIndex << G4endl;
|
||||
}
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
inline
|
||||
G4ParticleDefinition* G4VDecayChannel::GetParent()
|
||||
{
|
||||
//the pointer to the parent particle is filled, if it is not set yet
|
||||
if (parent == 0) FillParent();
|
||||
//get the pointer to the parent particle
|
||||
return parent;
|
||||
}
|
||||
|
||||
inline
|
||||
const G4String& G4VDecayChannel::GetParentName() const
|
||||
{
|
||||
return *parent_name;
|
||||
}
|
||||
|
||||
inline
|
||||
G4double G4VDecayChannel::GetParentMass() const
|
||||
{
|
||||
return parent_mass;
|
||||
}
|
||||
|
||||
|
||||
inline
|
||||
void G4VDecayChannel::SetParent(const G4String &particle_name)
|
||||
{
|
||||
if (parent_name != 0) delete parent_name;
|
||||
parent_name = new G4String(particle_name);
|
||||
parent = 0;
|
||||
}
|
||||
|
||||
inline
|
||||
G4int G4VDecayChannel::GetNumberOfDaughters() const
|
||||
{
|
||||
return numberOfDaughters;
|
||||
}
|
||||
|
||||
inline
|
||||
G4String G4VDecayChannel::GetKinematicsName() const { return kinematics_name; }
|
||||
|
||||
inline
|
||||
void G4VDecayChannel::SetBR(G4double value){ rbranch = value; }
|
||||
|
||||
inline
|
||||
G4double G4VDecayChannel::GetBR() const { return rbranch; }
|
||||
|
||||
inline
|
||||
void G4VDecayChannel::SetVerboseLevel(G4int value){ verboseLevel = value; }
|
||||
|
||||
inline
|
||||
G4int G4VDecayChannel::GetVerboseLevel() const { return verboseLevel; }
|
||||
|
||||
|
||||
inline G4ThreeVector G4VDecayChannel::GetParentPolarization(){return theParentPolarization;}
|
||||
inline void G4VDecayChannel::SetParentPolarization(G4ThreeVector polar){theParentPolarization=polar;}
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//*
|
||||
// LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION Geant4 SIMULATION
|
||||
// ID : LEMuSR.cc , v 1.0
|
||||
// ID : LEMuSR.cc , v 1.2
|
||||
// AUTHOR: Taofiq PARAISO
|
||||
// DATE : 2004-06-24 09:57
|
||||
|
||||
@ -28,7 +28,7 @@
|
||||
#include "G4UIterminal.hh"
|
||||
#include "G4UItcsh.hh"
|
||||
#include "G4ios.hh"
|
||||
#include <iomanip.h>
|
||||
#include <iomanip>
|
||||
|
||||
// II LEMuSR CLASSES
|
||||
// a_ Mandatory Classes
|
||||
@ -81,16 +81,18 @@ int main(int argc,char** argv)//argc:: defines the user interface
|
||||
|
||||
LEMuSRPhysicsList* lemuPhysicsList = new LEMuSRPhysicsList();
|
||||
|
||||
//! 2.2 LEMuSR Action class
|
||||
LEMuSRPrimaryGeneratorAction* lemuPGA = new LEMuSRPrimaryGeneratorAction();
|
||||
|
||||
|
||||
//! 2.3 Setting the mandatory Initialization classes
|
||||
//! 2.2 Setting the mandatory Initialization classes
|
||||
runManager ->SetUserInitialization( lemuDetector );
|
||||
|
||||
runManager ->SetUserInitialization( lemuPhysicsList );
|
||||
|
||||
|
||||
//! 2.3 LEMuSR Action class
|
||||
LEMuSRPrimaryGeneratorAction* lemuPGA = new LEMuSRPrimaryGeneratorAction();
|
||||
|
||||
|
||||
|
||||
//! 2.4 Setting the mandatory Action class
|
||||
runManager ->SetUserAction( lemuPGA );
|
||||
|
||||
@ -143,17 +145,16 @@ int main(int argc,char** argv)//argc:: defines the user interface
|
||||
session = new G4UIterminal();
|
||||
#endif
|
||||
|
||||
// UI->ApplyCommand("/control/execute visual.mac");
|
||||
UI->ApplyCommand("/run/verbose 2");
|
||||
UI->ApplyCommand("/vis/sceneHandler/create");
|
||||
UI->ApplyCommand("/control/verbose 1");
|
||||
UI->ApplyCommand("/run/verbose 1");
|
||||
|
||||
// UI->ApplyCommand("/Detector/MagneticField on");
|
||||
// UI->ApplyCommand("/run/beamOn 1");
|
||||
#if defined LEMU_TEST_ASYM
|
||||
// UI->ApplyCommand("/Detector/ElectricField off");
|
||||
UI->ApplyCommand("/Detector/MagneticField 50");
|
||||
UI->ApplyCommand("/Detector/AsymCheck on");
|
||||
UI->ApplyCommand("/lemuGun/gunPosition 0 0 0");
|
||||
UI->ApplyCommand("/lemuGun/energy/defined 0");
|
||||
UI->ApplyCommand("/lemuGun/energy/offset 0");
|
||||
G4cout<<"\n READY TO TEST ASYMETRY! ";
|
||||
#endif
|
||||
|
||||
|
||||
@ -24,7 +24,7 @@
|
||||
|
||||
#include "G4UserSteppingAction.hh"
|
||||
#include "globals.hh"
|
||||
#include <fstream.h>
|
||||
#include <fstream>
|
||||
#include "LEMuSRVisManager.hh"
|
||||
#include "LEMuSRDetectorConstruction.hh"
|
||||
|
||||
|
||||
@ -24,7 +24,7 @@
|
||||
|
||||
#include "G4UserSteppingAction.hh"
|
||||
#include "globals.hh"
|
||||
#include <fstream.h>
|
||||
#include <fstream>
|
||||
#include "LEMuSRVisManager.hh"
|
||||
#include "LEMuSRDetectorConstruction.hh"
|
||||
|
||||
|
||||
@ -131,7 +131,9 @@ G4VParticleChange* LEMuSRDecay::DecayIt(const G4Track& aTrack, const G4Step& aSt
|
||||
testb++;
|
||||
}
|
||||
#ifdef G4VERBOSE
|
||||
if (GetVerboseLevel()>1) {
|
||||
G4cout<<"Decay Channel LEMuSR "<<testa<<" Decay Channel Geant "<<testb<<G4endl;
|
||||
}
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user