00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4DNASancheSolvatationModel.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DNAWaterExcitationStructure.hh"
00043 #include "G4ParticleChangeForGamma.hh"
00044 #include "G4Electron.hh"
00045 #include "G4NistManager.hh"
00046 #include "G4DNAChemistryManager.hh"
00047 #include "G4DNAMolecularMaterial.hh"
00048
00049 G4DNASancheSolvatationModel::G4DNASancheSolvatationModel(const G4ParticleDefinition*,
00050 const G4String& nam):
00051 G4VEmModel(nam),fIsInitialised(false)
00052 {
00053 fVerboseLevel = 0 ;
00054 SetLowEnergyLimit(0.);
00055 G4DNAWaterExcitationStructure exStructure ;
00056 SetHighEnergyLimit(exStructure.ExcitationEnergy(0));
00057 fParticleChangeForGamma = 0;
00058
00059 fpWaterDensity = 0;
00060 }
00061
00062
00063 G4DNASancheSolvatationModel::~G4DNASancheSolvatationModel()
00064 {}
00065
00066
00067 void G4DNASancheSolvatationModel::Initialise(const G4ParticleDefinition* particleDefinition,
00068 const G4DataVector&)
00069 {
00070 #ifdef G4VERBOSE
00071 if (fVerboseLevel)
00072 G4cout << "Calling G4SancheSolvatationModel::Initialise()" << G4endl;
00073 #endif
00074 if (particleDefinition != G4Electron::ElectronDefinition())
00075 {
00076 G4ExceptionDescription exceptionDescription ;
00077 exceptionDescription << "Attempting to calculate cross section for wrong particle";
00078 G4Exception("G4DNASancheSolvatationModel::CrossSectionPerVolume","G4DNASancheSolvatationModel001",
00079 FatalErrorInArgument,exceptionDescription);
00080 return;
00081 }
00082
00083 if(!fIsInitialised)
00084 {
00085 fIsInitialised = true;
00086 fParticleChangeForGamma = GetParticleChangeForGamma();
00087 }
00088
00089
00090 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00091 }
00092
00093
00094 G4double G4DNASancheSolvatationModel::CrossSectionPerVolume(const G4Material* material,
00095 const G4ParticleDefinition*,
00096 G4double ekin,
00097 G4double,
00098 G4double)
00099 {
00100 #ifdef G4VERBOSE
00101 if (fVerboseLevel > 1)
00102 G4cout << "Calling CrossSectionPerVolume() of G4SancheSolvatationModel" << G4endl;
00103 #endif
00104
00105 if(ekin > HighEnergyLimit())
00106 {
00107 return 0.0;
00108 }
00109
00110 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00111
00112 if(waterDensity!= 0.0)
00113
00114 {
00115 if (ekin <= HighEnergyLimit())
00116 {
00117 return DBL_MAX;
00118 }
00119 }
00120 return 0. ;
00121 }
00122
00123
00124 G4ThreeVector G4DNASancheSolvatationModel::RadialDistributionOfProducts(G4double expectationValue) const
00125 {
00126 G4double sigma = std::sqrt(1.57)/2*expectationValue;
00127
00128 G4double XValueForfMax = std::sqrt(2.*sigma*sigma);
00129 G4double fMaxValue = std::sqrt(2./3.14) * 1./(sigma*sigma*sigma) *
00130 (XValueForfMax*XValueForfMax)*
00131 std::exp(-1./2. * (XValueForfMax*XValueForfMax)/(sigma*sigma));
00132
00133 G4double R;
00134
00135 do
00136 {
00137 G4double aRandomfValue = fMaxValue*G4UniformRand();
00138
00139 G4double sign;
00140 if(G4UniformRand() > 0.5)
00141 {
00142 sign = +1.;
00143 }
00144 else
00145 {
00146 sign = -1;
00147 }
00148
00149 R = expectationValue + sign*3.*sigma* G4UniformRand();
00150 G4double f = std::sqrt(2./3.14) * 1/std::pow(sigma, 3) * R*R * std::exp(-1./2. * R*R/(sigma*sigma));
00151
00152 if(aRandomfValue < f)
00153 {
00154 break;
00155 }
00156 }
00157 while(1);
00158
00159 G4double costheta = (2.*G4UniformRand()-1.);
00160 G4double theta = std::acos (costheta);
00161 G4double phi = 2.*pi*G4UniformRand();
00162
00163 G4double xDirection = R*std::cos(phi)* std::sin(theta);
00164 G4double yDirection = R*std::sin(theta)*std::sin(phi);
00165 G4double zDirection = R*costheta;
00166 G4ThreeVector RandDirection = G4ThreeVector(xDirection, yDirection, zDirection);
00167
00168 return RandDirection;
00169 }
00170
00171
00172 void G4DNASancheSolvatationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
00173 const G4MaterialCutsCouple*,
00174 const G4DynamicParticle* particle,
00175 G4double,
00176 G4double)
00177 {
00178 #ifdef G4VERBOSE
00179 if (fVerboseLevel)
00180 G4cout << "Calling SampleSecondaries() of G4SancheSolvatationModel" << G4endl;
00181 #endif
00182 G4double k = particle->GetKineticEnergy();
00183
00184 if (k <= HighEnergyLimit())
00185 {
00186 G4double r_mean =
00187 (-0.003*std::pow(k/eV,6) + 0.0749*std::pow(k/eV,5) - 0.7197*std::pow(k/eV,4)
00188 + 3.1384*std::pow(k/eV,3) - 5.6926*std::pow(k/eV,2) + 5.6237*k/eV - 0.7883)*nanometer;
00189
00190 G4ThreeVector displacement = RadialDistributionOfProducts (r_mean);
00191
00192
00193 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00194 G4ThreeVector finalPosition(theIncomingTrack->GetPosition()+displacement);
00195
00196 G4DNAChemistryManager::Instance()->CreateSolvatedElectron(theIncomingTrack,&finalPosition);
00197
00198 fParticleChangeForGamma->SetProposedKineticEnergy(25.e-3*eV);
00199 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00200 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
00201 }
00202 }