| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Modified by Z. Francis, S. Incerti to handl 27 // Modified by Z. Francis, S. Incerti to handle HZE
28 // && inverse rudd function sampling 26-10-201 28 // && inverse rudd function sampling 26-10-2010
29 // <<
30 // Rewitten by V.Ivanchenko 21.05.2023 <<
31 // <<
32 29
33 #include "G4EmCorrections.hh" <<
34 #include "G4DNARuddIonisationExtendedModel.hh" 30 #include "G4DNARuddIonisationExtendedModel.hh"
35 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
37 #include "G4UAtomicDeexcitation.hh" 33 #include "G4UAtomicDeexcitation.hh"
38 #include "G4LossTableManager.hh" 34 #include "G4LossTableManager.hh"
39 #include "G4NistManager.hh" <<
40 #include "G4DNAChemistryManager.hh" 35 #include "G4DNAChemistryManager.hh"
41 #include "G4DNAMolecularMaterial.hh" 36 #include "G4DNAMolecularMaterial.hh"
42 37
43 #include "G4IonTable.hh" 38 #include "G4IonTable.hh"
44 #include "G4DNARuddAngle.hh" 39 #include "G4DNARuddAngle.hh"
45 #include "G4DeltaAngle.hh" 40 #include "G4DeltaAngle.hh"
46 #include "G4Exp.hh" 41 #include "G4Exp.hh"
47 #include "G4Log.hh" 42 #include "G4Log.hh"
48 #include "G4Pow.hh" 43 #include "G4Pow.hh"
49 #include "G4Alpha.hh" 44 #include "G4Alpha.hh"
50 #include "G4Proton.hh" << 45
>> 46 static G4Pow * gpow = G4Pow::GetInstance();
>> 47
51 48
52 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 50
54 G4DNACrossSectionDataSet* G4DNARuddIonisationE << 51 using namespace std;
55 G4DNACrossSectionDataSet* G4DNARuddIonisationE <<
56 G4DNACrossSectionDataSet* G4DNARuddIonisationE <<
57 const std::vector<G4double>* G4DNARuddIonisati <<
58 <<
59 namespace <<
60 { <<
61 const G4double scaleFactor = CLHEP::m*CLHEP: <<
62 const G4double tolerance = 1*CLHEP::eV; <<
63 const G4double Ry = 13.6*CLHEP::eV; <<
64 <<
65 // Following values provided by M. Dingfelde <<
66 const G4double Bj[5] = {12.60*CLHEP::eV, 14. <<
67 32.20*CLHEP::eV, 539 <<
68 } <<
69 52
70 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 54
72 G4DNARuddIonisationExtendedModel::G4DNARuddIon 55 G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel(const G4ParticleDefinition*,
73 56 const G4String& nam)
74 : G4VEmModel(nam) << 57 :G4VEmModel(nam),isInitialised(false)
75 { 58 {
76 fEmCorrections = G4LossTableManager::Instanc << 59 slaterEffectiveCharge[0]=0.;
77 fGpow = G4Pow::GetInstance(); << 60 slaterEffectiveCharge[1]=0.;
78 fLowestEnergy = 100*CLHEP::eV; << 61 slaterEffectiveCharge[2]=0.;
79 fLimitEnergy = 1*CLHEP::keV; << 62 sCoefficient[0]=0.;
>> 63 sCoefficient[1]=0.;
>> 64 sCoefficient[2]=0.;
>> 65
>> 66 lowEnergyLimitForA[1] = 0 * eV;
>> 67 lowEnergyLimitForA[2] = 0 * eV;
>> 68 lowEnergyLimitForA[3] = 0 * eV;
>> 69 lowEnergyLimitOfModelForA[1] = 100 * eV;
>> 70 lowEnergyLimitOfModelForA[4] = 1 * keV;
>> 71 lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
>> 72 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
>> 73 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
>> 74 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
>> 75
>> 76 verboseLevel= 0;
>> 77 // Verbosity scale:
>> 78 // 0 = nothing
>> 79 // 1 = warning for energy non-conservation
>> 80 // 2 = details of energy budget
>> 81 // 3 = calculation of cross sections, file openings, sampling of atoms
>> 82 // 4 = entering in methods
>> 83
>> 84 if( verboseLevel>0 )
>> 85 {
>> 86 G4cout << "Rudd ionisation model is constructed " << G4endl;
>> 87 }
>> 88
>> 89 // Define default angular generator
>> 90 SetAngularDistribution(new G4DNARuddAngle());
80 91
81 // Mark this model as "applicable" for atomi << 92 // Mark this model as "applicable" for atomic deexcitation
82 SetDeexcitationFlag(true); << 93 SetDeexcitationFlag(true);
83 94
84 // Define default angular generator << 95 // Selection of stationary mode
85 SetAngularDistribution(new G4DNARuddAngle()) <<
86 96
87 if (nullptr == xshelium) { LoadData(); } << 97 statCode = false;
88 } 98 }
89 99
90 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 101
92 G4DNARuddIonisationExtendedModel::~G4DNARuddIo 102 G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel()
93 { 103 {
94 if(isFirst) { << 104 if(isIon) {
95 for(auto & i : xsdata) { delete i; } << 105 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 106 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
>> 107 {
>> 108 G4DNACrossSectionDataSet* table = pos->second;
>> 109 delete table;
>> 110 }
>> 111 } else {
>> 112 delete mainTable;
96 } 113 }
97 } 114 }
98 115
99 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 117
101 void G4DNARuddIonisationExtendedModel::LoadDat << 118 void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* particle,
>> 119 const G4DataVector& /*cuts*/)
102 { 120 {
103 // initialisation of static data once << 121 if (isInitialised) { return; }
104 isFirst = true; << 122 if (verboseLevel > 3)
105 G4String filename("dna/sigma_ionisation_h_ru << 123 G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
106 xsdata[0] = new G4DNACrossSectionDataSet(new <<
107 xsdata[0]->LoadData(filename); <<
108 <<
109 filename = "dna/sigma_ionisation_p_rudd"; <<
110 xsdata[1] = new G4DNACrossSectionDataSet(new <<
111 xsdata[1]->LoadData(filename); <<
112 <<
113 filename = "dna/sigma_ionisation_alphapluspl <<
114 xsdata[2] = new G4DNACrossSectionDataSet(new <<
115 xsdata[2]->LoadData(filename); <<
116 <<
117 filename = "dna/sigma_ionisation_li_rudd"; <<
118 xsdata[3] = new G4DNACrossSectionDataSet(new <<
119 xsdata[3]->LoadData(filename); <<
120 <<
121 filename = "dna/sigma_ionisation_be_rudd"; <<
122 xsdata[4] = new G4DNACrossSectionDataSet(new <<
123 xsdata[4]->LoadData(filename); <<
124 <<
125 filename = "dna/sigma_ionisation_b_rudd"; <<
126 xsdata[5] = new G4DNACrossSectionDataSet(new <<
127 xsdata[5]->LoadData(filename); <<
128 <<
129 filename = "dna/sigma_ionisation_c_rudd"; <<
130 xsdata[6] = new G4DNACrossSectionDataSet(new <<
131 xsdata[6]->LoadData(filename); <<
132 <<
133 filename = "dna/sigma_ionisation_n_rudd"; <<
134 xsdata[7] = new G4DNACrossSectionDataSet(new <<
135 xsdata[7]->LoadData(filename); <<
136 <<
137 filename = "dna/sigma_ionisation_o_rudd"; <<
138 xsdata[8] = new G4DNACrossSectionDataSet(new <<
139 xsdata[8]->LoadData(filename); <<
140 <<
141 filename = "dna/sigma_ionisation_si_rudd"; <<
142 xsdata[14] = new G4DNACrossSectionDataSet(ne <<
143 xsdata[14]->LoadData(filename); <<
144 <<
145 filename = "dna/sigma_ionisation_fe_rudd"; <<
146 xsdata[26] = new G4DNACrossSectionDataSet(ne <<
147 xsdata[26]->LoadData(filename); <<
148 <<
149 filename = "dna/sigma_ionisation_alphaplus_r <<
150 xsalphaplus = new G4DNACrossSectionDataSet(n <<
151 xsalphaplus->LoadData(filename); <<
152 <<
153 filename = "dna/sigma_ionisation_he_rudd"; <<
154 xshelium = new G4DNACrossSectionDataSet(new <<
155 xshelium->LoadData(filename); <<
156 <<
157 // to avoid possible threading problem fill <<
158 auto water = G4NistManager::Instance()->Find <<
159 fpWaterDensity = <<
160 G4DNAMolecularMaterial::Instance()->GetNum <<
161 } <<
162 124
163 //....oooOO0OOooo........oooOO0OOooo........oo << 125 // Energy limits
>> 126 G4String fileProton("dna/sigma_ionisation_p_rudd");
>> 127 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
>> 128 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
>> 129 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
>> 130 G4String fileHelium("dna/sigma_ionisation_he_rudd");
>> 131 G4String fileCarbon("dna/sigma_ionisation_c_rudd");
>> 132 G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
>> 133 G4String fileOxygen("dna/sigma_ionisation_o_rudd");
>> 134 G4String fileSilicon("dna/sigma_ionisation_si_rudd");
>> 135 G4String fileIron("dna/sigma_ionisation_fe_rudd");
164 136
165 void G4DNARuddIonisationExtendedModel::Initial << 137 G4String pname = particle->GetParticleName();
166 <<
167 { <<
168 if (p != fParticle) { SetParticle(p); } <<
169 138
170 // particle change object may be externally << 139 G4DNAGenericIonsManager *instance;
171 if (nullptr == fParticleChangeForGamma) { << 140 instance = G4DNAGenericIonsManager::Instance();
172 fParticleChangeForGamma = GetParticleChang << 141 protonDef = G4Proton::ProtonDefinition();
173 } << 142 hydrogenDef = instance->GetIon("hydrogen");
>> 143 alphaPlusPlusDef = G4Alpha::Alpha();
>> 144 alphaPlusDef = instance->GetIon("alpha+");
>> 145 heliumDef = instance->GetIon("helium");
174 146
175 // initialisation once in each thread << 147 carbonDef = instance->GetIon("carbon");
176 if (!isInitialised) { << 148 nitrogenDef = instance->GetIon("nitrogen");
177 isInitialised = true; << 149 oxygenDef = instance->GetIon("oxygen");
178 const G4String& pname = fParticle->GetPart << 150 siliconDef = instance->GetIon("silicon");
179 if (pname == "proton") { << 151 ironDef = instance->GetIon("iron");
180 idx = 1; << 152
181 xscurrent = xsdata[1]; << 153 G4String carbon;
182 fElow = fLowestEnergy; << 154 G4String nitrogen;
183 } else if (pname == "hydrogen") { << 155 G4String oxygen;
184 idx = 0; << 156 G4String silicon;
185 xscurrent = xsdata[0]; << 157 G4String iron;
186 fElow = fLowestEnergy; << 158
187 } else if (pname == "alpha") { << 159 G4double scaleFactor = 1 * m*m;
188 idx = 1; << 160 massC12 = carbonDef->GetPDGMass();
189 xscurrent = xsdata[2]; << 161
190 isHelium = true; << 162 // LIMITS AND DATA
191 fElow = fLimitEnergy; << 163
192 } else if (pname == "alpha+") { << 164 // **********************************************************************************************
193 idx = 1; << 165
194 isHelium = true; << 166 if(pname == "proton") {
195 xscurrent = xsalphaplus; << 167 localMinEnergy = lowEnergyLimitForA[1];
196 fElow = fLimitEnergy; << 168
197 // The following values are provided by << 169 // Cross section
198 slaterEffectiveCharge[0]=2.0; << 170 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
199 slaterEffectiveCharge[1]=2.0; << 171 mainTable->LoadData(fileProton);
200 slaterEffectiveCharge[2]=2.0; << 172
201 sCoefficient[0]=0.7; << 173 // **********************************************************************************************
202 sCoefficient[1]=0.15; << 174
203 sCoefficient[2]=0.15; << 175 } else if(pname == "hydrogen") {
204 } else if (pname == "helium") { << 176
205 idx = 0; << 177 localMinEnergy = lowEnergyLimitForA[1];
206 isHelium = true; << 178
207 fElow = fLimitEnergy; << 179 // Cross section
208 xscurrent = xshelium; << 180 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
209 slaterEffectiveCharge[0]=1.7; << 181 mainTable->LoadData(fileHydrogen);
210 slaterEffectiveCharge[1]=1.15; << 182
211 slaterEffectiveCharge[2]=1.15; << 183 // **********************************************************************************************
212 sCoefficient[0]=0.5; << 184
213 sCoefficient[1]=0.25; << 185 } else if(pname == "alpha") {
214 sCoefficient[2]=0.25; << 186
215 } else { << 187 localMinEnergy = lowEnergyLimitForA[4];
>> 188
>> 189 // Cross section
>> 190 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
>> 191 mainTable->LoadData(fileAlphaPlusPlus);
>> 192
>> 193 // **********************************************************************************************
>> 194
>> 195 } else if(pname == "alpha+") {
>> 196
>> 197 localMinEnergy = lowEnergyLimitForA[4];
>> 198
>> 199 // Cross section
>> 200 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
>> 201 mainTable->LoadData(fileAlphaPlus);
>> 202
>> 203 // **********************************************************************************************
>> 204
>> 205 } else if(pname == "helium") {
>> 206
>> 207 localMinEnergy = lowEnergyLimitForA[4];
>> 208
>> 209 // Cross section
>> 210 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
>> 211 mainTable->LoadData(fileHelium);
>> 212
>> 213 // **********************************************************************************************
>> 214
>> 215 } else if(pname == "GenericIon") {
>> 216
216 isIon = true; 217 isIon = true;
217 idx = -1; << 218 carbon = carbonDef->GetParticleName();
218 xscurrent = xsdata[1]; << 219 localMinEnergy = lowEnergyLimitForA[5]*massC12/proton_mass_c2;
219 fElow = fLowestEnergy; << 220
220 } << 221 // Cross section
221 // defined stationary mode << 222 mainTable = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
222 statCode = G4EmParameters::Instance()->DNA << 223 mainTable->LoadData(fileCarbon);
223 << 224
224 // initialise atomic de-excitation << 225 tableData[carbon] = mainTable;
225 fAtomDeexcitation = G4LossTableManager::In << 226
226 << 227 // **********************************************************************************************
227 if (verbose > 0) { << 228
228 G4cout << "### G4DNARuddIonisationExtend << 229 oxygen = oxygenDef->GetParticleName();
229 << "/n idx=" << idx << " isIon=" << << 230 tableFile[oxygen] = fileOxygen;
230 << " isHelium=" << isHelium << G4endl; << 231
>> 232 // Cross section
>> 233 G4DNACrossSectionDataSet* tableOxygen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 234 eV, scaleFactor );
>> 235 tableOxygen->LoadData(fileOxygen);
>> 236 tableData[oxygen] = tableOxygen;
>> 237
>> 238 // **********************************************************************************************
>> 239
>> 240 nitrogen = nitrogenDef->GetParticleName();
>> 241 tableFile[nitrogen] = fileNitrogen;
>> 242
>> 243 // Cross section
>> 244 G4DNACrossSectionDataSet* tableNitrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 245 eV, scaleFactor );
>> 246 tableNitrogen->LoadData(fileNitrogen);
>> 247 tableData[nitrogen] = tableNitrogen;
>> 248
>> 249 // **********************************************************************************************
>> 250
>> 251 silicon = siliconDef->GetParticleName();
>> 252 tableFile[silicon] = fileSilicon;
>> 253
>> 254 // Cross section
>> 255 G4DNACrossSectionDataSet* tableSilicon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 256 eV, scaleFactor );
>> 257 tableSilicon->LoadData(fileSilicon);
>> 258 tableData[silicon] = tableSilicon;
>> 259
>> 260 // **********************************************************************************************
>> 261
>> 262 iron = ironDef->GetParticleName();
>> 263 tableFile[iron] = fileIron;
>> 264
>> 265 // Cross section
>> 266
>> 267 G4DNACrossSectionDataSet* tableIron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 268 eV, scaleFactor );
>> 269 tableIron->LoadData(fileIron);
>> 270 tableData[iron] = tableIron;
231 } 271 }
232 } << 272 // **********************************************************************************************
>> 273
>> 274 if( verboseLevel>0 )
>> 275 {
>> 276 G4cout << "Rudd ionisation model is initialized " << G4endl
>> 277 << "Energy range for model: "
>> 278 << LowEnergyLimit() / keV << " keV - "
>> 279 << HighEnergyLimit() / keV << " keV for "
>> 280 << pname
>> 281 << " internal low energy limit E(keV)=" << localMinEnergy / keV
>> 282 << G4endl;
>> 283 }
>> 284
>> 285 // Initialize water density pointer
>> 286 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
>> 287
>> 288 //
>> 289 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
>> 290
>> 291 fParticleChangeForGamma = GetParticleChangeForGamma();
>> 292 isInitialised = true;
233 } 293 }
234 294
235 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236 296
237 void G4DNARuddIonisationExtendedModel::SetPart << 297 G4double G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(const G4Material* material,
238 { << 298 const G4ParticleDefinition* partDef,
239 fParticle = p; << 299 G4double k,
240 fMass = p->GetPDGMass(); << 300 G4double,
241 fMassRate = (isIon) ? CLHEP::proton_mass_c2/ << 301 G4double)
>> 302 {
>> 303 if (verboseLevel > 3)
>> 304 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
>> 305
>> 306 currParticle = GetDNAIonParticleDefinition(partDef);
>> 307 currentScaledEnergy = k;
>> 308 G4double e = k;
>> 309 G4double q2 = 1.0;
>> 310 currentTable = mainTable;
>> 311
>> 312 if (isIon){
>> 313 if (currParticle == nullptr) {//not DNA particle
>> 314 currentScaledEnergy *= massC12/partDef->GetPDGMass();
>> 315 G4double q = partDef->GetPDGCharge()/(eplus*6);
>> 316 q2 *= q*q;
>> 317 e = currentScaledEnergy;
>> 318 currParticle = carbonDef;
>> 319 }
>> 320 G4String pname = currParticle->GetParticleName();
>> 321 auto goodTable = tableData.find(pname);
>> 322 currentTable = goodTable->second;
>> 323 }
>> 324 // below low the ion should be stopped
>> 325 if (currentScaledEnergy < localMinEnergy) { return DBL_MAX; }
>> 326
>> 327 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
>> 328 G4double sigma = 0.0;
>> 329 if (nullptr != currentTable) {
>> 330 sigma = currentTable->FindValue(e)*q2;
>> 331 } else {
>> 332 G4cout << "G4DNARuddIonisationExtendedModel - no data table for "
>> 333 << partDef->GetParticleName() << G4endl;
>> 334 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(...)","em0002",
>> 335 FatalException,"Data table is not available for the model.");
>> 336 }
>> 337 if (verboseLevel > 2)
>> 338 {
>> 339 G4cout << "__________________________________" << G4endl;
>> 340 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
>> 341 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << partDef->GetParticleName() << G4endl;
>> 342 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
>> 343 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
>> 344 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
>> 345 }
>> 346 return sigma*waterDensity;
242 } 347 }
243 348
244 //....oooOO0OOooo........oooOO0OOooo........oo 349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245 350
246 G4double << 351 void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
247 G4DNARuddIonisationExtendedModel::CrossSection << 352 const G4MaterialCutsCouple* couple,
248 << 353 const G4DynamicParticle* particle,
249 << 354 G4double,
250 << 355 G4double)
251 { << 356 {
252 // check if model is applicable for given ma << 357 if (verboseLevel > 3)
253 G4double density = (material->GetIndex() < f << 358 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
254 ? (*fpWaterDensity)[material->GetIndex()] << 359
255 if (0.0 == density) { return 0.0; } << 360 // stop ion with energy below low energy limit
256 << 361 G4double k = particle->GetKineticEnergy();
257 // ion may be different << 362 if (currentScaledEnergy < localMinEnergy) {
258 if (fParticle != part) { SetParticle(part); << 363 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
259 << 364 fParticleChangeForGamma->ProposeTrackStatus(fStopButAlive);
260 // ion shoud be stopped - check on kinetic e << 365 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
261 if (kinE < fLowestEnergy) { return DBL_MAX; << 366 return;
>> 367 }
262 368
263 G4double e = kinE*fMassRate; << 369 // sampling of final state
>> 370 G4ParticleDefinition* definition = particle->GetDefinition();
>> 371 const G4ThreeVector& primaryDirection = particle->GetMomentumDirection();
>> 372
>> 373 G4int ionizationShell = RandomSelect(currentScaledEnergy);
>> 374
>> 375 // sample deexcitation
>> 376 // here we assume that H_{2}O electronic levels are the same as Oxygen.
>> 377 // this can be considered true with a rough 10% error in energy on K-shell,
>> 378
>> 379 G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
>> 380
>> 381 //SI: additional protection if tcs interpolation method is modified
>> 382 if (k < bindingEnergy) return;
>> 383 //
>> 384
>> 385 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
>> 386
>> 387 // is ionisation possible?
>> 388 G4double scatteredEnergy = k - bindingEnergy - secondaryKinetic;
>> 389 if(scatteredEnergy < 0.0) { return; }
>> 390
>> 391 G4int Z = 8;
>> 392
>> 393 G4ThreeVector deltaDirection =
>> 394 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
>> 395 Z, ionizationShell,
>> 396 couple->GetMaterial());
>> 397
>> 398 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
>> 399 fvect->push_back(dp);
>> 400
>> 401 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
>> 402
>> 403 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
>> 404 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
>> 405
>> 406 // SI: only atomic deexcitation from K shell is considered
>> 407 if(fAtomDeexcitation != nullptr && ionizationShell == 4)
>> 408 {
>> 409 const G4AtomicShell* shell
>> 410 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
>> 411 secNumberInit = fvect->size();
>> 412 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
>> 413 secNumberFinal = fvect->size();
>> 414
>> 415 if(secNumberFinal > secNumberInit)
>> 416 {
>> 417 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
>> 418 {
>> 419 //Check if there is enough residual energy
>> 420 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
>> 421 {
>> 422 //Ok, this is a valid secondary: keep it
>> 423 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
>> 424 }
>> 425 else
>> 426 {
>> 427 //Invalid secondary: not enough energy to create it!
>> 428 //Keep its energy in the local deposit
>> 429 delete (*fvect)[i];
>> 430 (*fvect)[i]=0;
>> 431 }
>> 432 }
>> 433 }
>> 434 }
>> 435
>> 436 //bindingEnergy has been decreased
>> 437 //by the amount of energy taken away by deexc. products
>> 438 if (!statCode)
>> 439 {
>> 440 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
>> 441 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
>> 442 }
>> 443 else
>> 444 {
>> 445 fParticleChangeForGamma->SetProposedKineticEnergy(k);
>> 446 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
>> 447 }
>> 448
>> 449 // TEST //////////////////////////
>> 450 // if (secondaryKinetic<0) abort();
>> 451 // if (scatteredEnergy<0) abort();
>> 452 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
>> 453 // if (k-scatteredEnergy<0) abort();
>> 454 /////////////////////////////////
>> 455
>> 456 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
>> 457 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
>> 458 ionizationShell,
>> 459 theIncomingTrack);
>> 460 // SI - not useful since low energy of model is 0 eV
>> 461 }
>> 462
>> 463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 464
>> 465 G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
>> 466 G4double k,
>> 467 G4int shell)
>> 468 {
>> 469 //-- Fast sampling method -----
>> 470 G4double proposed_energy;
>> 471 G4double random1;
>> 472 G4double value_sampling;
>> 473 G4double max1;
>> 474
>> 475 do
>> 476 {
>> 477 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
>> 478
>> 479 max1=0.;
>> 480
>> 481 for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
>> 482 max1=RejectionFunction(particleDefinition, k, en, shell);
>> 483
>> 484 random1 = G4UniformRand()*max1;
>> 485
>> 486 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
>> 487
>> 488 } while(random1 > value_sampling);
>> 489
>> 490 return(proposed_energy);
>> 491 }
>> 492
>> 493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 494
>> 495 G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
>> 496 G4double k,
>> 497 G4double proposed_ws,
>> 498 G4int ionizationLevelIndex)
>> 499 {
>> 500 const G4int j=ionizationLevelIndex;
>> 501 G4double Bj_energy, alphaConst;
>> 502 G4double Ry = 13.6*eV;
>> 503 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
>> 504
>> 505 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
>> 506
>> 507 // Following values provided by M. Dingfelder (priv. comm)
>> 508 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
>> 509
>> 510 if (j == 4)
>> 511 {
>> 512 alphaConst = 0.66;
>> 513 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
>> 514 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
>> 515 //---
>> 516 }
>> 517 else
>> 518 {
>> 519 alphaConst = 0.64;
>> 520 Bj_energy = Bj[ionizationLevelIndex];
>> 521 }
264 522
265 G4double sigma = (e > fElow) ? xscurrent->Fi << 523 G4double energyTransfer = proposed_ws + Bj_energy;
266 : xscurrent->FindValue(fElow) * e / fElow; << 524 proposed_ws/=Bj_energy;
267 525
268 if (idx == -1) { << 526 G4double tau = 0.;
269 sigma *= fEmCorrections->EffectiveChargeSq << 527 G4double A_ion = 0.;
270 } << 528 tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
>> 529 A_ion = particleDefinition->GetAtomicMass();
>> 530
>> 531 G4double v2;
>> 532 G4double beta2;
>> 533
>> 534 if((tau/MeV)<5.447761194e-2)
>> 535 {
>> 536 v2 = tau / Bj_energy;
>> 537 beta2 = 2.*tau / electron_mass_c2;
>> 538 }
>> 539 // Relativistic
>> 540 else
>> 541 {
>> 542 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
>> 543 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
>> 544 }
271 545
272 sigma *= density; << 546 G4double v = std::sqrt(v2);
>> 547 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
>> 548 G4double rejection_term = 1.+G4Exp(alphaConst*(proposed_ws - wc) / v);
>> 549 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
>> 550 //* (S/Bj_energy) ; Not needed anymore
>> 551
>> 552 G4bool isHelium = false;
>> 553
>> 554 if ( particleDefinition == protonDef
>> 555 || particleDefinition == hydrogenDef
>> 556 )
>> 557 {
>> 558 return(rejection_term);
>> 559 }
273 560
274 if (verbose > 1) { << 561 else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
275 G4cout << "G4DNARuddIonisationExtendedMode << 562 {
276 << " Ekin(keV)=" << kinE/CLHEP::keV << 563 G4double Z = particleDefinition->GetAtomicNumber();
277 << " sigma(cm^2)=" << sigma/CLHEP:: << 564
278 } << 565 G4double x = 100.*std::sqrt(beta2)/gpow->powA(Z, 2./3.);
279 return sigma; << 566 G4double Zeffion = Z*(1.-G4Exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
280 } << 567 rejection_term*=Zeffion*Zeffion;
>> 568 }
281 569
282 //....oooOO0OOooo........oooOO0OOooo........oo << 570 else if (particleDefinition == alphaPlusPlusDef )
>> 571 {
>> 572 isHelium = true;
>> 573 slaterEffectiveCharge[0]=0.;
>> 574 slaterEffectiveCharge[1]=0.;
>> 575 slaterEffectiveCharge[2]=0.;
>> 576 sCoefficient[0]=0.;
>> 577 sCoefficient[1]=0.;
>> 578 sCoefficient[2]=0.;
>> 579 }
283 580
284 void << 581 else if (particleDefinition == alphaPlusDef )
285 G4DNARuddIonisationExtendedModel::SampleSecond << 582 {
286 << 583 isHelium = true;
287 << 584 slaterEffectiveCharge[0]=2.0;
288 << 585 // The following values are provided by M. Dingfelder (priv. comm)
289 { << 586 slaterEffectiveCharge[1]=2.0;
290 const G4ParticleDefinition* pd = dpart->GetD << 587 slaterEffectiveCharge[2]=2.0;
291 if (fParticle != pd) { SetParticle(pd); } << 588 //
292 << 589 sCoefficient[0]=0.7;
293 // stop ion with energy below low energy lim << 590 sCoefficient[1]=0.15;
294 G4double kinE = dpart->GetKineticEnergy(); << 591 sCoefficient[2]=0.15;
295 // ion shoud be stopped - check on kinetic e << 592 }
296 if (kinE <= fLowestEnergy) { <<
297 fParticleChangeForGamma->SetProposedKineti <<
298 fParticleChangeForGamma->ProposeTrackStatu <<
299 fParticleChangeForGamma->ProposeLocalEnerg <<
300 return; <<
301 } <<
302 593
303 G4int shell = SelectShell(kinE*fMassRate); << 594 else if (particleDefinition == heliumDef )
304 G4double bindingEnergy = (useDNAWaterStructu << 595 {
305 ? waterStructure.IonisationEnergy(shell) : << 596 isHelium = true;
306 << 597 slaterEffectiveCharge[0]=1.7;
307 //Si: additional protection if tcs interpola << 598 slaterEffectiveCharge[1]=1.15;
308 if (kinE < bindingEnergy) { return; } << 599 slaterEffectiveCharge[2]=1.15;
309 << 600 sCoefficient[0]=0.5;
310 G4double esec = SampleElectronEnergy(kinE, s << 601 sCoefficient[1]=0.25;
311 G4double esum = 0.0; << 602 sCoefficient[2]=0.25;
312 <<
313 // sample deexcitation <<
314 // here we assume that H2O electronic levels <<
315 // this can be considered true with a rough <<
316 G4int Z = 8; <<
317 G4ThreeVector deltaDir = <<
318 GetAngularDistribution()->SampleDirectionF <<
319 <<
320 // SI: only atomic deexcitation from K shell <<
321 if(fAtomDeexcitation != nullptr && shell == <<
322 auto as = G4AtomicShellEnumerator(0); <<
323 auto ashell = fAtomDeexcitation->GetAtomic <<
324 fAtomDeexcitation->GenerateParticles(fvect <<
325 <<
326 // compute energy sum from de-excitation <<
327 for (auto const & ptr : *fvect) { <<
328 esum += ptr->GetKineticEnergy(); <<
329 } 603 }
330 } <<
331 // check energy balance <<
332 // remaining excitation energy of water mole <<
333 G4double exc = bindingEnergy - esum; <<
334 <<
335 // remaining projectile energy <<
336 G4double scatteredEnergy = kinE - bindingEne <<
337 if(scatteredEnergy < -tolerance || exc < -to <<
338 G4cout << "G4DNARuddIonisationExtendedMode <<
339 << "negative final E(keV)=" << scat <<
340 << kinE/CLHEP::keV << " " << pd->G <<
341 << " Edelta(keV)=" << esec/CLHEP::k <<
342 << G4endl; <<
343 } <<
344 604
345 // projectile << 605 if (isHelium)
346 if (!statCode) { << 606 {
347 fParticleChangeForGamma->SetProposedKineti <<
348 fParticleChangeForGamma->ProposeLocalEnerg <<
349 } else { <<
350 fParticleChangeForGamma->SetProposedKineti <<
351 fParticleChangeForGamma->ProposeLocalEnerg <<
352 } <<
353 607
354 // delta-electron << 608 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
355 auto dp = new G4DynamicParticle(G4Electron: <<
356 fvect->push_back(dp); <<
357 609
358 // create radical << 610 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
359 const G4Track* theIncomingTrack = fParticleC << 611 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
360 G4DNAChemistryManager::Instance()->CreateWat << 612 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
361 theIncomingTrack); <<
362 } <<
363 613
364 //....oooOO0OOooo........oooOO0OOooo........oo << 614 rejection_term*= zEff * zEff;
>> 615 }
365 616
366 G4int G4DNARuddIonisationExtendedModel::Select << 617 return (rejection_term);
367 { <<
368 G4double sum = 0.0; <<
369 G4double xs; <<
370 for (G4int i=0; i<5; ++i) { <<
371 auto ptr = xscurrent->GetComponent(i); <<
372 xs = (e > fElow) ? ptr->FindValue(e) : ptr <<
373 sum += xs; <<
374 fTemp[i] = sum; <<
375 } <<
376 sum *= G4UniformRand(); <<
377 for (G4int i=0; i<5; ++i) { <<
378 if (sum <= fTemp[i]) { return i; } <<
379 } <<
380 return 0; <<
381 } 618 }
382 619
>> 620
383 //....oooOO0OOooo........oooOO0OOooo........oo 621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
384 622
385 G4double G4DNARuddIonisationExtendedModel::Max << 623
>> 624 G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
>> 625 G4double k,
>> 626 G4int ionizationLevelIndex)
386 { 627 {
387 // kinematic limit <<
388 G4double tau = kine/fMass; <<
389 G4double gam = 1.0 + tau; <<
390 G4double emax = 2.0*CLHEP::electron_mass_c2* <<
391 <<
392 // Initialisation of sampling <<
393 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2; <<
394 if (shell == 4) { <<
395 //Data For Liquid Water K SHELL from Dingf <<
396 A1 = 1.25; <<
397 B1 = 0.5; <<
398 C1 = 1.00; <<
399 D1 = 1.00; <<
400 E1 = 3.00; <<
401 A2 = 1.10; <<
402 B2 = 1.30; <<
403 C2 = 1.00; <<
404 D2 = 0.00; <<
405 alphaConst = 0.66; <<
406 } else { <<
407 //Data For Liquid Water from Dingfelder (P <<
408 A1 = 1.02; <<
409 B1 = 82.0; <<
410 C1 = 0.45; <<
411 D1 = -0.80; <<
412 E1 = 0.38; <<
413 A2 = 1.07; <<
414 // Value provided by M. Dingfelder (priv. <<
415 B2 = 11.6; <<
416 C2 = 0.60; <<
417 D2 = 0.04; <<
418 alphaConst = 0.64; <<
419 } <<
420 bEnergy = Bj[shell]; <<
421 G4double v2 = 0.25*emax/(bEnergy*gam*gam); <<
422 v = std::sqrt(v2); <<
423 u = Ry/bEnergy; <<
424 wc = 4.*v2 - 2.*v - 0.25*u; <<
425 <<
426 G4double L1 = (C1 * fGpow->powA(v, D1)) / (1 <<
427 G4double L2 = C2 * fGpow->powA(v, D2); <<
428 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + <<
429 G4double H2 = (A2 / v2) + (B2 / (v2 * v2)); <<
430 <<
431 F1 = L1 + H1; <<
432 F2 = (L2 * H2) / (L2 + H2); <<
433 return emax; <<
434 } <<
435 628
436 //....oooOO0OOooo........oooOO0OOooo........oo << 629 const G4int j=ionizationLevelIndex;
437 630
438 G4double G4DNARuddIonisationExtendedModel::Sam << 631 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
439 << 632 //G4double alphaConst ;
440 { << 633 G4double Bj_energy;
441 G4double emax = MaxEnergy(kine, shell); << 634
442 // compute cumulative probability function << 635 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
443 G4double step = 1*CLHEP::eV; << 636 // Following values provided by M. Dingfelder (priv. comm)
444 auto nn = (G4int)(emax/step); << 637
445 nn = std::min(std::max(nn, 10), 100); << 638 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
446 step = emax/(G4double)nn; << 639
447 << 640 if (j == 4)
448 // find max probability << 641 {
449 G4double pmax = ProbabilityFunction(kine, 0. << 642 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
450 //G4cout << "## E(keV)=" << kine/keV << " em << 643 A1 = 1.25;
451 // << " pmax(0)=" << pmax << " shell=" << 644 B1 = 0.5;
452 << 645 C1 = 1.00;
453 G4double e0 = 0.0; // energy with max probab << 646 D1 = 1.00;
454 // 2 areas after point with max probability << 647 E1 = 3.00;
455 G4double e1 = emax; << 648 A2 = 1.10;
456 G4double e2 = emax; << 649 B2 = 1.30;
457 G4double p1 = 0.0; << 650 C2 = 1.00;
458 G4double p2 = 0.0; << 651 D2 = 0.00;
459 const G4double f = 0.25; << 652 //alphaConst = 0.66;
460 << 653 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
461 // find max probability << 654 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
462 G4double e = 0.0; << 655 //---
463 G4double p = 0.0; <<
464 for (G4int i=0; i<nn; ++i) { <<
465 e += step; <<
466 p = ProbabilityFunction(kine, e, shell); <<
467 if (p > pmax) { <<
468 pmax = p; <<
469 e0 = e; <<
470 } else { <<
471 break; <<
472 } 656 }
473 } << 657 else
474 // increase step to be more effective << 658 {
475 step *= 2.0; << 659 //Data For Liquid Water from Dingfelder (Protons in Water)
476 // 2-nd area << 660 A1 = 1.02;
477 for (G4int i=0; i<nn; ++i) { << 661 B1 = 82.0;
478 e += step; << 662 C1 = 0.45;
479 if (std::abs(e - emax) < step) { << 663 D1 = -0.80;
480 e1 = emax; << 664 E1 = 0.38;
481 break; << 665 A2 = 1.07;
482 } << 666 //B2 = 14.6; From Ding Paper
483 p = ProbabilityFunction(kine, e, shell); << 667 // Value provided by M. Dingfelder (priv. comm)
484 if (p < f*pmax) { << 668 B2 = 11.6;
485 p1 = p; << 669 //
486 e1 = e; << 670 C2 = 0.60;
487 break; << 671 D2 = 0.04;
>> 672 //alphaConst = 0.64;
>> 673
>> 674 Bj_energy = Bj[ionizationLevelIndex];
488 } 675 }
489 } << 676
490 // 3-d area << 677 G4double tau = 0.;
491 if (e < emax) { << 678 G4double A_ion = 0.;
492 for (G4int i=0; i<nn; ++i) { << 679 tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
493 e += step; << 680
494 if (std::abs(e - emax) < step) { << 681 A_ion = particle->GetAtomicMass();
495 e2 = emax; << 682
496 break; << 683 G4double v2;
497 } << 684 G4double beta2;
498 p = ProbabilityFunction(kine, e, shell); << 685 if((tau/MeV)<5.447761194e-2)
499 if (p < f*p1) { << 686 {
500 p2 = p; << 687 v2 = tau / Bj_energy;
501 e2 = e; << 688 beta2 = 2.*tau / electron_mass_c2;
502 break; <<
503 } <<
504 } 689 }
505 } << 690 // Relativistic
506 pmax *= 1.05; << 691 else
507 // regression method with 3 regions << 692 {
508 G4double s0 = pmax*e1; << 693 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
509 G4double s1 = s0 + p1 * (e2 - e1); << 694 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
510 G4double s2 = s1 + p2 * (emax - e2); <<
511 s0 = (s0 == s1) ? 1.0 : s0 / s2; <<
512 s1 = (s1 == s2) ? 1.0 : s1 / s2; <<
513 <<
514 //G4cout << "pmax=" << pmax << " e1(keV)=" < <<
515 // << " p2=" << p2 << " s0=" << s0 << " s1 <<
516 <<
517 // sampling <<
518 G4int count = 0; <<
519 G4double ymax, y, deltae; <<
520 for (G4int i = 0; i<100000; ++i) { <<
521 G4double q = G4UniformRand(); <<
522 if (q <= s0) { <<
523 ymax = pmax; <<
524 deltae = e1 * q / s0; <<
525 } else if (q <= s1) { <<
526 ymax = p1; <<
527 deltae = e1 + (e2 - e1) * (q - s0) / (s1 <<
528 } else { <<
529 ymax = p2; <<
530 deltae = e2 + (emax - e2) * (q - s1) / ( <<
531 } 695 }
532 y = ProbabilityFunction(kine, deltae, shel << 696
533 //G4cout << " " << i << ". deltae=" << << 697 G4double v = std::sqrt(v2);
534 // << " y=" << y << " ymax=" << ymax << 698 //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
535 if (y > ymax && count < 10) { << 699 G4double L1 = (C1* gpow->powA(v,(D1))) / (1.+ E1*gpow->powA(v, (D1+4.)));
536 ++count; << 700 G4double L2 = C2*gpow->powA(v,(D2));
537 G4cout << "G4DNARuddIonisationExtendedMo << 701 G4double H1 = (A1*G4Log(1.+v2)) / (v2+(B1/v2));
538 << fParticle->GetParticleName() << " E( << 702 G4double H2 = (A2/v2) + (B2/(v2*v2));
539 << " Edelta(keV)=" << deltae/CLHEP::keV << 703 G4double F1 = L1+H1;
540 << " y=" << y << " ymax=" << ymax << " << 704 G4double F2 = (L2*H2)/(L2+H2);
>> 705
>> 706 // ZF. generalized & relativistic version
>> 707 G4double maximumEnergy;
>> 708
>> 709 //---- maximum kinetic energy , non relativistic ------
>> 710 if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
>> 711 {
>> 712 maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
541 } 713 }
542 if (ymax * G4UniformRand() <= y) { << 714 //---- relativistic -----------------------------------
543 return deltae; << 715 else
>> 716 {
>> 717 G4double gamma = 1./sqrt(1.-beta2);
>> 718 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
>> 719 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
544 } 720 }
545 } << 721
546 deltae = std::min(e0 + step, 0.5*emax); << 722 //either it is transfered energy or secondary electron energy ...
547 return deltae; << 723 //maximumEnergy-=Bj_energy;
>> 724
>> 725 //-----------------------------------------------------
>> 726 G4double wmax = maximumEnergy/Bj_energy;
>> 727 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
>> 728 c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c
>> 729 G4double randVal = G4UniformRand();
>> 730 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
>> 731 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
>> 732 // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
>> 733 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
>> 734 proposed_ws*=Bj_energy;
>> 735
>> 736 return(proposed_ws);
>> 737
548 } 738 }
549 739
550 //....oooOO0OOooo........oooOO0OOooo........oo 740 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
551 741
552 G4double G4DNARuddIonisationExtendedModel::Pro << 742 G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
553 << 743 G4double energyTransferred,
554 << 744 G4double slaterEffectiveChg,
555 { << 745 G4double shellNumber)
556 // Shells ids are 0 1 2 3 4 (4 is k shell) << 746 {
557 // !!Attention, "energyTransfer" here is the << 747 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
558 // that the secondary kinetic en << 748 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
559 // << 749
560 // ds S F1(nu) + << 750 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
561 // ---- = G(k) * ---- ----------------- << 751 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
562 // dw Bj (1+w)^3 * [1 + e <<
563 // <<
564 // w is the secondary electron kinetic Energ <<
565 // <<
566 // All the other parameters can be found in <<
567 // <<
568 // M.Eugene Rudd, 1988, User-Friendly model <<
569 // electrons from protons or electron collis <<
570 // <<
571 G4double w = deltae/bEnergy; <<
572 G4double x = alphaConst*(w - wc)/v; <<
573 G4double y = (x > -15.) ? 1.0 + G4Exp(x) : 1 <<
574 <<
575 G4double res = CorrectionFactor(kine, shell) <<
576 (fGpow->powN((1. + w)/u, 3) * y); <<
577 <<
578 if (isHelium) { <<
579 G4double energyTransfer = deltae + bEnergy <<
580 G4double Zeff = 2.0 - <<
581 (sCoefficient[0] * S_1s(kine, energyTran <<
582 sCoefficient[1] * S_2s(kine, energyTran <<
583 sCoefficient[2] * S_2p(kine, energyTran <<
584 752
585 res *= Zeff * Zeff; << 753 return value;
586 } <<
587 return std::max(res, 0.0); <<
588 } 754 }
589 755
590 //....oooOO0OOooo........oooOO0OOooo........oo 756 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
591 757
592 G4double G4DNARuddIonisationExtendedModel::Com << 758 G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
593 const G4ParticleDefinition* p, G4doub << 759 G4double energyTransferred,
>> 760 G4double slaterEffectiveChg,
>> 761 G4double shellNumber)
594 { 762 {
595 if (fParticle != p) { SetParticle(p); } << 763 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
596 MaxEnergy(e, shell); << 764 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
597 return ProbabilityFunction(e, deltae, shell) << 765
>> 766 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
>> 767 G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
>> 768
>> 769 return value;
>> 770
598 } 771 }
599 772
600 //....oooOO0OOooo........oooOO0OOooo........oo 773 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
601 774
602 G4double G4DNARuddIonisationExtendedModel::S_1 << 775 G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
603 << 776 G4double energyTransferred,
604 << 777 G4double slaterEffectiveChg,
605 778 G4double shellNumber)
606 { 779 {
607 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) << 780 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
608 // Dingfelder, in Chattanooga 2005 proceedin << 781 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
>> 782
>> 783 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
>> 784 G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
609 785
610 G4double r = Rh(kine, energyTransfer, slater << 786 return value;
611 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. <<
612 return value; <<
613 } 787 }
614 788
615 //....oooOO0OOooo........oooOO0OOooo........oo 789 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
616 790
617 G4double G4DNARuddIonisationExtendedModel::S_2 << 791 G4double G4DNARuddIonisationExtendedModel::R(G4double t,
618 << 792 G4double energyTransferred,
619 << 793 G4double slaterEffectiveChg,
620 << 794 G4double shellNumber)
621 { 795 {
622 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) << 796 // tElectron = m_electron / m_alpha * t
623 // Dingfelder, in Chattanooga 2005 proceedin << 797 // Dingfelder, in Chattanooga 2005 proceedings, p 4
624 798
625 G4double r = Rh(kine, energyTransfer, slater << 799 G4double tElectron = 0.511/3728. * t;
626 G4double value = << 800 // The following values are provided by M. Dingfelder (priv. comm)
627 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * << 801 G4double H = 2.*13.60569172 * eV;
>> 802 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
628 803
629 return value; << 804 return value;
630 } 805 }
631 806
632 //....oooOO0OOooo........oooOO0OOooo........oo 807 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
633 808
634 G4double G4DNARuddIonisationExtendedModel::S_2 << 809 G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
635 <<
636 <<
637 <<
638 { 810 {
639 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^ << 811 // ZF Shortened
640 // Dingfelder, in Chattanooga 2005 proceedin << 812
>> 813 if (particleDefinition == hydrogenDef && shell < 4)
>> 814 {
>> 815 G4double value = ((G4Log(k/eV)/gpow->logZ(10))-4.2)/0.5;
>> 816 // The following values are provided by M. Dingfelder (priv. comm)
>> 817 return((0.6/(1+G4Exp(value))) + 0.9);
>> 818 }
>> 819 else
>> 820 {
>> 821 return(1.);
>> 822 }
>> 823 }
>> 824
>> 825 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 826
>> 827 G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k)
>> 828 {
>> 829 G4int level = 0;
>> 830
>> 831 // Retrieve data table corresponding to the current particle type
>> 832
>> 833 G4DNACrossSectionDataSet* table = currentTable;
641 834
642 G4double r = Rh(kine, energyTransfer, slater << 835 if (table != nullptr)
643 G4double value = << 836 {
644 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3 << 837 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
645 838
646 return value; << 839 const G4int n = (G4int)table->NumberOfComponents();
>> 840 G4int i(n);
>> 841 G4double value = 0.;
>> 842
>> 843 while (i>0)
>> 844 {
>> 845 --i;
>> 846 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
>> 847
>> 848 value += valuesBuffer[i];
>> 849 }
>> 850
>> 851 value *= G4UniformRand();
>> 852
>> 853 i = n;
>> 854
>> 855 while (i > 0)
>> 856 {
>> 857 --i;
>> 858
>> 859 if (valuesBuffer[i] > value)
>> 860 {
>> 861 delete[] valuesBuffer;
>> 862 return i;
>> 863 }
>> 864 value -= valuesBuffer[i];
>> 865 }
>> 866
>> 867 if (valuesBuffer) delete[] valuesBuffer;
>> 868
>> 869 }
>> 870 return level;
647 } 871 }
648 872
649 //....oooOO0OOooo........oooOO0OOooo........oo 873 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
650 874
651 G4double G4DNARuddIonisationExtendedModel::Rh( << 875 G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
652 <<
653 { 876 {
654 // The following values are provided by M. D << 877 G4double sigma = 0.;
655 // Dingfelder, in Chattanooga 2005 proceedin << 878
>> 879 const G4DynamicParticle* particle = track.GetDynamicParticle();
>> 880 G4double k = particle->GetKineticEnergy();
656 881
657 G4double escaled = CLHEP::electron_mass_c2/f << 882 G4double lowLim = 0;
658 const G4double H = 13.60569172 * CLHEP::eV; << 883 G4double highLim = 0;
659 G4double value = 2.0*std::sqrt(escaled / H)* <<
660 884
661 return value; << 885 const G4String& particleName = particle->GetDefinition()->GetParticleName();
>> 886
>> 887 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
>> 888 pos1 = lowEnergyLimit.find(particleName);
>> 889
>> 890 if (pos1 != lowEnergyLimit.end())
>> 891 {
>> 892 lowLim = pos1->second;
>> 893 }
>> 894
>> 895 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
>> 896 pos2 = highEnergyLimit.find(particleName);
>> 897
>> 898 if (pos2 != highEnergyLimit.end())
>> 899 {
>> 900 highLim = pos2->second;
>> 901 }
>> 902
>> 903 if (k >= lowLim && k <= highLim)
>> 904 {
>> 905 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 906 pos = tableData.find(particleName);
>> 907
>> 908 if (pos != tableData.end())
>> 909 {
>> 910 G4DNACrossSectionDataSet* table = pos->second;
>> 911 if (table != 0)
>> 912 {
>> 913 sigma = table->FindValue(k);
>> 914 }
>> 915 }
>> 916 else
>> 917 {
>> 918 G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
>> 919 FatalException,"Model not applicable to particle type.");
>> 920 }
>> 921 }
>> 922
>> 923 return sigma;
662 } 924 }
663 925
664 //....oooOO0OOooo........oooOO0OOooo........oo 926 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
665 927
666 G4double << 928 G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
667 G4DNARuddIonisationExtendedModel::CorrectionFa <<
668 { 929 {
669 // ZF Shortened << 930 return 0;
670 G4double res = 1.0; <<
671 if (shell < 4 && 0 != idx) { <<
672 const G4double ln10 = fGpow->logZ(10); <<
673 G4double x = 2.0*((G4Log(kine/CLHEP::eV)/l <<
674 // The following values are provided by M. <<
675 res = 0.6/(1.0 + G4Exp(x)) + 0.9; <<
676 } <<
677 return res; <<
678 } 931 }
679 932
680 //....oooOO0OOooo........oooOO0OOooo........oo 933 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 934
>> 935 G4ParticleDefinition* G4DNARuddIonisationExtendedModel::GetDNAIonParticleDefinition(const G4ParticleDefinition* particleDefinition)
>> 936 {
>> 937 //for proton, hydrogen, alphas
>> 938 if(particleDefinition->GetAtomicMass() <= 4)
>> 939 {
>> 940 return const_cast<G4ParticleDefinition*>(particleDefinition);
>> 941 }
>> 942 else{
>> 943 auto instance = G4DNAGenericIonsManager::Instance();
>> 944
>> 945 auto PDGEncoding = particleDefinition->GetPDGEncoding();
>> 946 if(PDGEncoding == 1000140280){
>> 947 return instance->GetIon("silicon");
>> 948 }else if(PDGEncoding == 1000260560){
>> 949 return instance->GetIon("iron");
>> 950 }else if(PDGEncoding == 1000080160){
>> 951 return instance->GetIon("oxygen");
>> 952 }else if(PDGEncoding == 1000070140){
>> 953 return instance->GetIon("nitrogen");
>> 954 }else if(PDGEncoding == 1000060120){
>> 955 return instance->GetIon("carbon");
>> 956 }
>> 957 //if there is no DNA particle, get nullptr
>> 958 return nullptr;
>> 959 }
>> 960 }
681 961