| 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 // CPA100 ionisation model class for electrons 26 // CPA100 ionisation model class for electrons
27 // 27 //
28 // Based on the work of M. Terrissol and M. C. 28 // Based on the work of M. Terrissol and M. C. Bordage
29 // 29 //
30 // Users are requested to cite the following p 30 // Users are requested to cite the following papers:
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Do 31 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terr 32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
33 // M. Bardies, N. Lampe, S. Incerti, Phys. M 33 // M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
34 // 34 //
35 // Authors of this class: 35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bor 36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
37 // 37 //
38 // 15.01.2014: creation 38 // 15.01.2014: creation
39 // 39 //
40 // Based on the study by S. Zein et. al. Nucl. <<
41 // 1/2/2023 : Hoang added modification <<
42 40
43 #include "G4DNACPA100IonisationModel.hh" 41 #include "G4DNACPA100IonisationModel.hh"
44 <<
45 #include "G4DNAChemistryManager.hh" <<
46 #include "G4DNAMaterialManager.hh" <<
47 #include "G4DNAMolecularMaterial.hh" <<
48 #include "G4EnvironmentUtils.hh" <<
49 #include "G4LossTableManager.hh" <<
50 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh"
52 #include "G4UAtomicDeexcitation.hh" 44 #include "G4UAtomicDeexcitation.hh"
53 << 45 #include "G4LossTableManager.hh"
54 #include <fstream> << 46 #include "G4DNAChemistryManager.hh"
>> 47 #include "G4DNAMolecularMaterial.hh"
55 48
56 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57 50
58 using namespace std; 51 using namespace std;
59 52
60 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 54
62 G4DNACPA100IonisationModel::G4DNACPA100Ionisat 55 G4DNACPA100IonisationModel::G4DNACPA100IonisationModel(const G4ParticleDefinition*,
63 56 const G4String& nam)
64 : G4VDNAModel(nam, "all") << 57 :G4VEmModel(nam),isInitialised(false)
65 { 58 {
66 fpGuanine = G4Material::GetMaterial("G4_GUAN << 59 verboseLevel= 0;
67 fpG4_WATER = G4Material::GetMaterial("G4_WAT << 60 // Verbosity scale:
68 fpDeoxyribose = G4Material::GetMaterial("G4_ << 61 // 0 = nothing
69 fpCytosine = G4Material::GetMaterial("G4_CYT << 62 // 1 = warning for energy non-conservation
70 fpThymine = G4Material::GetMaterial("G4_THYM << 63 // 2 = details of energy budget
71 fpAdenine = G4Material::GetMaterial("G4_ADEN << 64 // 3 = calculation of cross sections, file openings, sampling of atoms
72 fpPhosphate = G4Material::GetMaterial("G4_PH << 65 // 4 = entering in methods
73 fpParticle = G4Electron::ElectronDefinition( << 66
>> 67 if( verboseLevel>0 )
>> 68 {
>> 69 G4cout << "CPA100 ionisation model is constructed " << G4endl;
>> 70 }
>> 71
>> 72 SetLowEnergyLimit(11*eV);
>> 73 SetHighEnergyLimit(255955*eV);
>> 74
>> 75 // Mark this model as "applicable" for atomic deexcitation
>> 76 SetDeexcitationFlag(true);
>> 77 fAtomDeexcitation = 0;
>> 78 fParticleChangeForGamma = 0;
>> 79 fpMolWaterDensity = 0;
>> 80
>> 81 // Selection of computation method
>> 82
>> 83 // useDcs = true if usage of dcs for sampling of secondaries
>> 84 // useDcs = false if usage of composition sampling (DEFAULT)
>> 85
>> 86 useDcs = true;
>> 87
>> 88 // if useDcs is true, one has the following choice
>> 89 // fasterCode = true for usage of cumulated dcs (DEFAULT)
>> 90 // fasterCode = false for usage of non-cumulated dcs
>> 91
>> 92 fasterCode = true;
>> 93
>> 94 // Selection of stationary mode
>> 95
>> 96 statCode = false;
74 } 97 }
75 98
76 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 100
78 void G4DNACPA100IonisationModel::Initialise(co << 101 G4DNACPA100IonisationModel::~G4DNACPA100IonisationModel()
79 co <<
80 { 102 {
81 if (isInitialised) { << 103 // Cross section
82 return; <<
83 } <<
84 if (verboseLevel > 3) { <<
85 G4cout << "Calling G4DNACPA100IonisationMo <<
86 } <<
87 104
88 if (!G4DNAMaterialManager::Instance()->IsLoc << 105 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
89 if (p != fpParticle) { << 106 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
90 std::ostringstream oss; << 107 {
91 oss << " Model is not applied for this p << 108 G4DNACrossSectionDataSet* table = pos->second;
92 G4Exception("G4DNACPA100IonisationModel: << 109 delete table;
93 FatalException, oss.str().c_ <<
94 } <<
95 <<
96 const char* path = G4FindDataDir("G4LEDATA <<
97 <<
98 if (path == nullptr) { <<
99 G4Exception("G4DNACPA100IonisationModel: <<
100 "G4LEDATA environment variab <<
101 return; <<
102 } <<
103 <<
104 std::size_t index; <<
105 if (fpG4_WATER != nullptr) { <<
106 index = fpG4_WATER->GetIndex(); <<
107 G4String eFullFileName = ""; <<
108 fasterCode ? eFullFileName = "/dna/sigma <<
109 : eFullFileName = "/dna/sigma <<
110 AddCrossSectionData(index, p, "dna/sigma <<
111 1.e-20 * m * m); <<
112 SetLowELimit(index, p, 11 * eV); <<
113 SetHighELimit(index, p, 255955 * eV); <<
114 } <<
115 if (fpGuanine != nullptr) { <<
116 index = fpGuanine->GetIndex(); <<
117 G4String eFullFileName = ""; <<
118 if(useDcs) { <<
119 fasterCode ? eFullFileName = "/dna/sig <<
120 : eFullFileName = "/dna/sig <<
121 } <<
122 AddCrossSectionData(index, p, "dna/sigma <<
123 1. * cm * cm); <<
124 SetLowELimit(index, p, 11 * eV); <<
125 SetHighELimit(index, p, 1 * MeV); <<
126 } <<
127 if (fpDeoxyribose != nullptr) { <<
128 index = fpDeoxyribose->GetIndex(); <<
129 G4String eFullFileName = ""; <<
130 if(useDcs) { <<
131 eFullFileName = "/dna/sigmadiff_cumula <<
132 } <<
133 AddCrossSectionData(index, p, "dna/sigma <<
134 1. * cm * cm); <<
135 SetLowELimit(index, p, 11 * eV); <<
136 SetHighELimit(index, p, 1 * MeV); <<
137 } <<
138 if (fpCytosine != nullptr) { <<
139 index = fpCytosine->GetIndex(); <<
140 G4String eFullFileName = ""; <<
141 if(useDcs) { <<
142 fasterCode ? eFullFileName = "/dna/sig <<
143 : eFullFileName = "/dna/sig <<
144 } <<
145 AddCrossSectionData(index, p, "dna/sigma <<
146 1. * cm * cm); <<
147 SetLowELimit(index, p, 11 * eV); <<
148 SetHighELimit(index, p, 1 * MeV); <<
149 } <<
150 if (fpThymine != nullptr) { <<
151 index = fpThymine->GetIndex(); <<
152 G4String eFullFileName = ""; <<
153 if(useDcs) { <<
154 fasterCode ? eFullFileName = "/dna/sig <<
155 : eFullFileName = "/dna/sig <<
156 } <<
157 AddCrossSectionData(index, p, "dna/sigma <<
158 1. * cm * cm); <<
159 SetLowELimit(index, p, 11 * eV); <<
160 SetHighELimit(index, p, 1 * MeV); <<
161 } <<
162 if (fpAdenine != nullptr) { <<
163 index = fpAdenine->GetIndex(); <<
164 G4String eFullFileName = ""; <<
165 if(useDcs) { <<
166 fasterCode ? eFullFileName = "/dna/sig <<
167 : eFullFileName = "/dna/sig <<
168 } <<
169 AddCrossSectionData(index, p, "dna/sigma <<
170 1. * cm * cm); <<
171 SetLowELimit(index, p, 11 * eV); <<
172 SetHighELimit(index, p, 1 * MeV); <<
173 } <<
174 if (fpPhosphate != nullptr) { <<
175 index = fpPhosphate->GetIndex(); <<
176 G4String eFullFileName = ""; <<
177 if(useDcs) { <<
178 eFullFileName = "dna/sigmadiff_cumulat <<
179 } <<
180 AddCrossSectionData(index, p, "dna/sigma <<
181 1. * cm * cm); <<
182 SetLowELimit(index, p, 11 * eV); <<
183 SetHighELimit(index, p, 1 * MeV); <<
184 } <<
185 LoadCrossSectionData(p); <<
186 G4DNAMaterialManager::Instance()->SetMaste <<
187 fpModelData = this; <<
188 } <<
189 else { <<
190 auto dataModel = dynamic_cast<G4DNACPA100I <<
191 G4DNAMaterialManager::Instance()->GetMod <<
192 if (dataModel == nullptr) { <<
193 G4cout << "G4DNACPA100IonisationModel::C <<
194 throw; <<
195 } 110 }
196 fpModelData = dataModel; <<
197 } <<
198 111
199 fAtomDeexcitation = G4LossTableManager::Inst << 112 // Final state
>> 113
>> 114 eVecm.clear();
200 115
201 fParticleChangeForGamma = GetParticleChangeF <<
202 isInitialised = true; <<
203 } 116 }
204 117
205 G4double G4DNACPA100IonisationModel::CrossSect << 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206 << 119
207 << 120 void G4DNACPA100IonisationModel::Initialise(const G4ParticleDefinition* particle,
>> 121 const G4DataVector& /*cuts*/)
208 { 122 {
209 // initialise the cross section value (outpu <<
210 G4double sigma(0); <<
211 123
212 // Get the current particle name << 124 if (verboseLevel > 3)
213 const G4String& particleName = p->GetParticl << 125 G4cout << "Calling G4DNACPA100IonisationModel::Initialise()" << G4endl;
214 126
215 if (p != fpParticle) { << 127 // Energy limits
216 G4Exception("G4DNACPA100IonisationModel::C << 128
217 "No model is registered for th << 129 // The following file is proved by M. Terrissol et al. (sigion3)
218 } << 130
>> 131 G4String fileElectron("dna/sigma_ionisation_e_cpa100_form_rel");
>> 132
>> 133 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
>> 134
>> 135 G4String electron;
>> 136
>> 137 G4double scaleFactor = 1.e-20 * m*m;
>> 138
>> 139 char *path = getenv("G4LEDATA");
>> 140
>> 141 // *** ELECTRON
>> 142
>> 143 electron = electronDef->GetParticleName();
>> 144
>> 145 tableFile[electron] = fileElectron;
>> 146
>> 147 // Cross section
>> 148
>> 149 G4DNACrossSectionDataSet* tableE =
>> 150 new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
>> 151
>> 152 //G4DNACrossSectionDataSet* tableE =
>> 153 // new G4DNACrossSectionDataSet(new G4DNACPA100LogLogInterpolation, eV,scaleFactor );
>> 154
>> 155 tableE->LoadData(fileElectron);
219 156
220 auto matID = material->GetIndex(); << 157 tableData[electron] = tableE;
221 158
222 // Set the low and high energy limits << 159 // Final state
223 G4double lowLim = fpModelData->GetLowELimit( <<
224 G4double highLim = fpModelData->GetHighELimi <<
225 160
226 // Check that we are in the correct energy r << 161 // ******************************
227 if (ekin >= lowLim && ekin < highLim) { <<
228 // Get the map with all the model data tab <<
229 auto tableData = fpModelData->GetData(); <<
230 162
231 if ((*tableData)[matID][p] == nullptr) { << 163 if (useDcs)
232 G4Exception("G4DNACPA100IonisationModel: << 164 {
233 "No model is registered"); << 165
>> 166 std::ostringstream eFullFileName;
>> 167
>> 168 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat";
>> 169
>> 170 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_cpa100_rel.dat";
>> 171
>> 172 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
>> 173
>> 174 if (!eDiffCrossSection)
>> 175 {
>> 176 if (fasterCode) G4Exception("G4DNACPA100IonisationModel::Initialise","em0003",
>> 177 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat");
>> 178
>> 179 if (!fasterCode) G4Exception("G4DNACPA100IonisationModel::Initialise","em0003",
>> 180 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_cpa100_rel.dat");
234 } 181 }
235 else { << 182
236 sigma = (*tableData)[matID][p]->FindValu << 183 // Clear the arrays for re-initialization case (MT mode)
>> 184 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
>> 185
>> 186 eTdummyVec.clear();
>> 187 eVecm.clear();
>> 188 eProbaShellMap->clear();
>> 189 eDiffCrossSectionData->clear();
>> 190 eNrjTransfData->clear();
>> 191
>> 192 //
>> 193
>> 194 eTdummyVec.push_back(0.);
>> 195 while(!eDiffCrossSection.eof())
>> 196 {
>> 197 G4double tDummy;
>> 198 G4double eDummy;
>> 199 eDiffCrossSection>>tDummy>>eDummy;
>> 200 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
>> 201 for (G4int j=0; j<5; j++)
>> 202 {
>> 203 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
>> 204
>> 205 if (fasterCode)
>> 206 {
>> 207 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
>> 208 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
>> 209 }
>> 210
>> 211 // SI - only if eof is not reached
>> 212 if (!eDiffCrossSection.eof() && !fasterCode)
>> 213 eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
>> 214
>> 215 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
>> 216
>> 217 }
237 } 218 }
238 219
239 if (verboseLevel > 2) { << 220 //
240 auto MolDensity = << 221
241 (*G4DNAMolecularMaterial::Instance()-> << 222 } // end of if (useDcs)
>> 223
>> 224 // ******************************
>> 225
>> 226 //
>> 227
>> 228 if( verboseLevel>0 )
>> 229 {
>> 230 G4cout << "CPA100 ionisation model is initialized " << G4endl
>> 231 << "Energy range: "
>> 232 << LowEnergyLimit() / eV << " eV - "
>> 233 << HighEnergyLimit() / keV << " keV for "
>> 234 << particle->GetParticleName()
>> 235 << G4endl;
>> 236 }
>> 237
>> 238 // Initialize water density pointer
>> 239 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
>> 240
>> 241 // AD
>> 242 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
>> 243
>> 244 if (isInitialised) return;
>> 245 fParticleChangeForGamma = GetParticleChangeForGamma();
>> 246 isInitialised = true;
>> 247 }
>> 248
>> 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 250
>> 251 G4double G4DNACPA100IonisationModel::CrossSectionPerVolume( const G4Material* material,
>> 252 const G4ParticleDefinition* particleDefinition,
>> 253 G4double ekin,
>> 254 G4double,
>> 255 G4double)
>> 256 {
>> 257
>> 258 if (verboseLevel > 3)
>> 259 G4cout << "Calling CrossSectionPerVolume() of G4DNACPA100IonisationModel" << G4endl;
>> 260
>> 261 if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
>> 262
>> 263 // Calculate total cross section for model
>> 264
>> 265 G4double sigma=0;
>> 266
>> 267 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
>> 268
>> 269 const G4String& particleName = particleDefinition->GetParticleName();
>> 270
>> 271 if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
>> 272 {
>> 273 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 274 pos = tableData.find(particleName);
>> 275
>> 276 if (pos != tableData.end())
>> 277 {
>> 278 G4DNACrossSectionDataSet* table = pos->second;
>> 279 if (table != 0) sigma = table->FindValue(ekin);
>> 280 }
>> 281 else
>> 282 {
>> 283 G4Exception("G4DNACPA100IonisationModel::CrossSectionPerVolume","em0002",
>> 284 FatalException,"Model not applicable to particle type.");
>> 285 }
>> 286 }
>> 287
>> 288 if (verboseLevel > 2)
>> 289 {
242 G4cout << "_____________________________ 290 G4cout << "__________________________________" << G4endl;
243 G4cout << "°°° G4DNACPA100IonisationM << 291 G4cout << "G4DNACPA100IonisationModel - XS INFO START" << G4endl;
244 G4cout << "°°° Kinetic energy(eV)=" < << 292 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
245 G4cout << "°°° lowLim (eV) = " << low << 293 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
246 G4cout << "°°° Materials = " << (*G4M << 294 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
247 G4cout << "°°° Cross section per " << << 295 G4cout << "G4DNACPA100IonisationModel - XS INFO END" << G4endl;
248 << G4endl; <<
249 G4cout << "°°° Cross section per Phos <<
250 << sigma * MolDensity / (1. / cm) <<
251 G4cout << "°°° G4DNACPA100IonisationM <<
252 } 296 }
253 } <<
254 297
255 auto MolDensity = (*G4DNAMolecularMaterial:: << 298 return sigma*waterDensity;
256 return sigma * MolDensity; <<
257 } 299 }
258 300
259 //....oooOO0OOooo........oooOO0OOooo........oo 301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260 302
261 void G4DNACPA100IonisationModel::SampleSeconda << 303 void G4DNACPA100IonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
262 std::vector<G4DynamicParticle*>* fvect, << 304 const G4MaterialCutsCouple* ,//must be set!
263 const G4MaterialCutsCouple* couple, // must << 305 const G4DynamicParticle* particle,
264 const G4DynamicParticle* particle, G4double, << 306 G4double,
>> 307 G4double)
265 { 308 {
266 if (verboseLevel > 3) { << 309 if (verboseLevel > 3)
267 G4cout << "Calling SampleSecondaries() of 310 G4cout << "Calling SampleSecondaries() of G4DNACPA100IonisationModel" << G4endl;
268 } <<
269 auto k = particle->GetKineticEnergy(); <<
270 311
271 const G4Material* material = couple->GetMate << 312 G4double k = particle->GetKineticEnergy();
272 313
273 auto MatID = material->GetIndex(); << 314 const G4String& particleName = particle->GetDefinition()->GetParticleName();
274 315
275 auto p = particle->GetDefinition(); << 316 if (k >= LowEnergyLimit() && k <= HighEnergyLimit())
>> 317 {
>> 318 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
>> 319 G4double particleMass = particle->GetDefinition()->GetPDGMass();
>> 320 G4double totalEnergy = k + particleMass;
>> 321 G4double pSquare = k * (totalEnergy + particleMass);
>> 322 G4double totalMomentum = std::sqrt(pSquare);
>> 323
>> 324 G4int ionizationShell = -1;
>> 325
>> 326 ionizationShell = RandomSelect(k,particleName);
>> 327
>> 328 //SI: PROTECTION FOR G4LOGLOGINTERPOLATION ON UPPER VALUE
>> 329 if (k<waterStructure.IonisationEnergy(ionizationShell)) { return; }
>> 330
>> 331 G4double bindingEnergy = 0;
>> 332 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
>> 333
>> 334 G4double secondaryKinetic=-1000*eV;
>> 335
>> 336 if (useDcs && !fasterCode)
>> 337 secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
>> 338
>> 339 if (useDcs && fasterCode)
>> 340 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
>> 341
>> 342 if (!useDcs)
>> 343 secondaryKinetic = RandomizeEjectedElectronEnergyFromCompositionSampling(particle->GetDefinition(),k,ionizationShell);
>> 344
>> 345 // Quick test
>> 346 /*
>> 347 FILE* myFile;
>> 348 myFile=fopen("nrj.txt","a");
>> 349 fprintf(myFile,"%e\n", secondaryKinetic/eV );
>> 350 fclose(myFile);
>> 351 */
>> 352
>> 353 G4double cosTheta = 0.;
>> 354 G4double phi = 0.;
>> 355 RandomizeEjectedElectronDirection(particle->GetDefinition(), k,secondaryKinetic, cosTheta, phi);
>> 356
>> 357 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
>> 358 G4double dirX = sinTheta*std::cos(phi);
>> 359 G4double dirY = sinTheta*std::sin(phi);
>> 360 G4double dirZ = cosTheta;
>> 361 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
>> 362 deltaDirection.rotateUz(primaryDirection);
>> 363
>> 364 // SI - For atom. deexc. tagging - 23/05/2017
>> 365 if (secondaryKinetic>0)
>> 366 {
>> 367 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
>> 368 fvect->push_back(dp);
>> 369 }
>> 370 //
276 371
277 auto lowLim = fpModelData->GetLowELimit(MatI << 372 if (particle->GetDefinition() == G4Electron::ElectronDefinition())
278 auto highLim = fpModelData->GetHighELimit(Ma << 373 {
279 << 374 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
280 // Check if we are in the correct energy ran << 375
281 if (k >= lowLim && k < highLim) { << 376 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
282 const auto& primaryDirection = particle->G << 377 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
283 auto particleMass = particle->GetDefinitio << 378 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
284 auto totalEnergy = k + particleMass; << 379 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
285 auto pSquare = k * (totalEnergy + particle << 380 finalPx /= finalMomentum;
286 auto totalMomentum = std::sqrt(pSquare); << 381 finalPy /= finalMomentum;
287 G4int shell = -1; << 382 finalPz /= finalMomentum;
288 G4double bindingEnergy, secondaryKinetic; <<
289 shell = fpModelData->RandomSelectShell(k, <<
290 bindingEnergy = iStructure.IonisationEnerg <<
291 <<
292 if (k < bindingEnergy) { <<
293 return; <<
294 } <<
295 <<
296 auto info = std::make_tuple(MatID, k, shel <<
297 <<
298 secondaryKinetic = -1000 * eV; <<
299 if (fpG4_WATER->GetIndex() != MatID) {//fo <<
300 secondaryKinetic = fpModelData->Randomiz <<
301 }else if(fasterCode){ <<
302 secondaryKinetic = fpModelData->Random <<
303 }else{ <<
304 secondaryKinetic = fpModelData->Random <<
305 } <<
306 383
307 G4double cosTheta = 0.; << 384 G4ThreeVector direction;
308 G4double phi = 0.; << 385 direction.set(finalPx,finalPy,finalPz);
309 RandomizeEjectedElectronDirection(particle << 386
310 phi); << 387 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
311 << 388 }
312 G4double sinTheta = std::sqrt(1. - cosThet << 389
313 G4double dirX = sinTheta * std::cos(phi); << 390 else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
314 G4double dirY = sinTheta * std::sin(phi); << 391
315 G4double dirZ = cosTheta; << 392 // SI - For atom. deexc. tagging - 23/05/2017
316 G4ThreeVector deltaDirection(dirX, dirY, d << 393
317 deltaDirection.rotateUz(primaryDirection); << 394 // AM: sample deexcitation
318 << 395 // here we assume that H_{2}O electronic levels are the same of Oxigen.
319 // SI - For atom. deexc. tagging - 23/05/2 << 396 // this can be considered true with a rough 10% error in energy on K-shell,
320 if (secondaryKinetic > 0) { << 397
321 auto dp = new G4DynamicParticle(G4Electr << 398 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
322 fvect->push_back(dp); << 399 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
323 } << 400
324 << 401 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
325 if (particle->GetDefinition() != fpParticl << 402
326 fParticleChangeForGamma->ProposeMomentum << 403 // SI: only atomic deexcitation from K shell is considered
327 } << 404 if(fAtomDeexcitation && ionizationShell == 4)
328 else { << 405 {
329 G4double deltaTotalMomentum = << 406 G4int Z = 8;
330 std::sqrt(secondaryKinetic * (secondar << 407 const G4AtomicShell* shell =
331 G4double finalPx = << 408 fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
332 totalMomentum * primaryDirection.x() - << 409 secNumberInit = fvect->size();
333 G4double finalPy = << 410 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
334 totalMomentum * primaryDirection.y() - << 411 secNumberFinal = fvect->size();
335 G4double finalPz = << 412
336 totalMomentum * primaryDirection.z() - << 413 if(secNumberFinal > secNumberInit)
337 G4double finalMomentum = std::sqrt(final << 414 {
338 finalPx /= finalMomentum; << 415 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
339 finalPy /= finalMomentum; << 416 {
340 finalPz /= finalMomentum; << 417 //Check if there is enough residual energy
341 << 418 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
342 G4ThreeVector direction; << 419 {
343 direction.set(finalPx, finalPy, finalPz) << 420 //Ok, this is a valid secondary: keep it
344 << 421 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
345 fParticleChangeForGamma->ProposeMomentum << 422 }
346 } << 423 else
347 << 424 {
348 // SI - For atom. deexc. tagging - 23/05/2 << 425 //Invalid secondary: not enough energy to create it!
349 << 426 //Keep its energy in the local deposit
350 // AM: sample deexcitation << 427 delete (*fvect)[i];
351 // here we assume that H_{2}O electronic l << 428 (*fvect)[i]=0;
352 // this can be considered true with a roug << 429 }
353 << 430 }
354 G4double scatteredEnergy = k - bindingEner <<
355 <<
356 // SI: only atomic deexcitation from K she <<
357 // Hoang: only for water <<
358 if (material == G4Material::GetMaterial("G <<
359 std::size_t secNumberInit = 0; // need <<
360 std::size_t secNumberFinal = 0; // So I <<
361 if ((fAtomDeexcitation != nullptr) && sh <<
362 G4int Z = 8; <<
363 auto Kshell = fAtomDeexcitation->GetAt <<
364 secNumberInit = fvect->size(); <<
365 fAtomDeexcitation->GenerateParticles(f <<
366 secNumberFinal = fvect->size(); <<
367 if (secNumberFinal > secNumberInit) { <<
368 for (std::size_t i = secNumberInit; <<
369 // Check if there is enough residu <<
370 if (bindingEnergy >= ((*fvect)[i]) <<
371 // Ok, this is a valid secondary <<
372 bindingEnergy -= ((*fvect)[i])-> <<
373 } <<
374 else { <<
375 // Invalid secondary: not enough <<
376 // Keep its energy in the local <<
377 delete (*fvect)[i]; <<
378 (*fvect)[i] = nullptr; <<
379 } <<
380 } 431 }
>> 432
>> 433 }
>> 434
>> 435 //This should never happen
>> 436 if(bindingEnergy < 0.0)
>> 437 G4Exception("G4DNACPA100IonisatioModel1::SampleSecondaries()",
>> 438 "em2050",FatalException,"Negative local energy deposit");
>> 439
>> 440 //bindingEnergy has been decreased
>> 441 //by the amount of energy taken away by deexc. products
>> 442
>> 443 if (!statCode)
>> 444 {
>> 445 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
>> 446 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
>> 447 }
>> 448 else
>> 449 {
>> 450 fParticleChangeForGamma->SetProposedKineticEnergy(k);
>> 451 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
381 } 452 }
382 } <<
383 } <<
384 453
385 // This should never happen << 454 // TEST //////////////////////////
386 if (bindingEnergy < 0.0) { << 455 // if (secondaryKinetic<0) abort();
387 G4Exception("G4DNACPA100IonisatioModel1: << 456 // if (scatteredEnergy<0) abort();
388 "Negative local energy depos << 457 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
389 } << 458 // if (k-scatteredEnergy<0) abort();
390 if (!statCode) { << 459 /////////////////////////////////
391 fParticleChangeForGamma->SetProposedKine << 460
392 fParticleChangeForGamma->ProposeLocalEne << 461 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
393 } << 462 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
394 else { << 463 ionizationShell,
395 fParticleChangeForGamma->SetProposedKine << 464 theIncomingTrack);
396 fParticleChangeForGamma->ProposeLocalEne <<
397 } <<
398 <<
399 // only water for chemistry <<
400 if (fpG4_WATER != nullptr && material == G <<
401 const G4Track* theIncomingTrack = fParti <<
402 G4DNAChemistryManager::Instance()->Creat <<
403 <<
404 } 465 }
405 } << 466
406 } 467 }
407 468
408 //....oooOO0OOooo........oooOO0OOooo........oo 469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
409 470
410 G4double G4DNACPA100IonisationModel::Randomize << 471 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
>> 472 G4double k, G4int shell)
411 { 473 {
412 auto MatID = std::get<0>(info); << 474 // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
413 auto k = std::get<1>(info); << 475
414 auto shell = std::get<2>(info); << 476 if (particleDefinition == G4Electron::ElectronDefinition())
415 G4double maximumEnergyTransfer = 0.; << 477 {
416 auto IonLevel = iStructure.IonisationEnergy( << 478 G4double maximumEnergyTransfer=0.;
417 (k + IonLevel) / 2. > k ? maximumEnergyTrans << 479 if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
418 << 480 else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
419 G4double crossSectionMaximum = 0.; << 481
420 << 482 // SI : original method
421 G4double minEnergy = IonLevel; << 483 /*
422 G4double maxEnergy = maximumEnergyTransfer; << 484 G4double crossSectionMaximum = 0.;
423 << 485 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
424 // nEnergySteps can be optimized - 100 by de << 486 {
425 G4int nEnergySteps = 50; << 487 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
426 << 488 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
427 G4double value(minEnergy); << 489 }
428 G4double stpEnergy(std::pow(maxEnergy / valu << 490 */
429 G4int step(nEnergySteps); << 491
430 G4double differentialCrossSection = 0.; << 492 // SI : alternative method
431 while (step > 0) { << 493
432 step--; << 494 G4double crossSectionMaximum = 0.;
433 differentialCrossSection = DifferentialCro << 495
434 << 496 G4double minEnergy = waterStructure.IonisationEnergy(shell);
435 if (differentialCrossSection > 0) { << 497 G4double maxEnergy = maximumEnergyTransfer;
436 crossSectionMaximum = differentialCrossS << 498
437 break; << 499 // nEnergySteps can be optimized - 100 by default
438 } << 500 G4int nEnergySteps = 50;
439 value *= stpEnergy; << 501
440 } << 502 // *** METHOD 1
>> 503 // FOR SLOW COMPUTATION ONLY
>> 504 /*
>> 505 G4double value(minEnergy);
>> 506 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
>> 507 G4int step(nEnergySteps);
>> 508 while (step>0)
>> 509 {
>> 510 step--;
>> 511 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
>> 512 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
>> 513 value*=stpEnergy;
>> 514 }
>> 515 */
>> 516
>> 517 // *** METHOD 2 : Faster method for CPA100 only since DCS is monotonously decreasing
>> 518 // FOR SLOW COMPUTATION ONLY
>> 519
>> 520 G4double value(minEnergy);
>> 521 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
>> 522 G4int step(nEnergySteps);
>> 523 G4double differentialCrossSection = 0.;
>> 524 while (step>0)
>> 525 {
>> 526 step--;
>> 527 differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
>> 528 if(differentialCrossSection >0)
>> 529 {
>> 530 crossSectionMaximum=differentialCrossSection;
>> 531 break;
>> 532 }
>> 533 value*=stpEnergy;
>> 534 }
441 535
442 G4double secondaryElectronKineticEnergy = 0. << 536 //
443 do { <<
444 secondaryElectronKineticEnergy = G4Uniform <<
445 } while (G4UniformRand() * crossSectionMaxim <<
446 > DifferentialCrossSection(info, (s <<
447 537
448 return secondaryElectronKineticEnergy; << 538 G4double secondaryElectronKineticEnergy=0.;
>> 539 do
>> 540 {
>> 541 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
>> 542 } while(G4UniformRand()*crossSectionMaximum >
>> 543 DifferentialCrossSection(particleDefinition, k/eV,
>> 544 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
>> 545
>> 546 return secondaryElectronKineticEnergy;
>> 547
>> 548 }
>> 549
>> 550 return 0;
449 } 551 }
450 552
451 //....oooOO0OOooo........oooOO0OOooo........oo 553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
452 554
453 void G4DNACPA100IonisationModel::RandomizeEjec 555 void G4DNACPA100IonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition*,
454 << 556 G4double k,
455 << 557 G4double secKinetic,
456 << 558 G4double & cosTheta,
457 { << 559 G4double & phi )
458 phi = twopi * G4UniformRand(); << 560 {
459 G4double sin2O = (1. - secKinetic / k) / (1. << 561
460 cosTheta = std::sqrt(1. - sin2O); << 562 phi = twopi * G4UniformRand();
>> 563 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
>> 564 cosTheta = std::sqrt(1.-sin2O);
>> 565
461 } 566 }
462 567
463 //....oooOO0OOooo........oooOO0OOooo........oo 568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
464 569
465 G4double G4DNACPA100IonisationModel::Different << 570 G4double G4DNACPA100IonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
466 << 571 G4double k,
>> 572 G4double energyTransfer,
>> 573 G4int ionizationLevelIndex)
467 { 574 {
468 auto MatID = std::get<0>(info); << 575 G4double sigma = 0.;
469 auto k = std::get<1>(info) / eV; // in eV u << 576
470 auto shell = std::get<2>(info); << 577 if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
471 G4double sigma = 0.; << 578 {
472 G4double shellEnergy = iStructure.Ionisation << 579 G4double valueT1 = 0;
473 G4double kSE(energyTransfer - shellEnergy); << 580 G4double valueT2 = 0;
474 << 581 G4double valueE21 = 0;
475 if (energyTransfer >= shellEnergy) { << 582 G4double valueE22 = 0;
476 G4double valueT1 = 0; << 583 G4double valueE12 = 0;
477 G4double valueT2 = 0; << 584 G4double valueE11 = 0;
478 G4double valueE21 = 0; << 585
479 G4double valueE22 = 0; << 586 G4double xs11 = 0;
480 G4double valueE12 = 0; << 587 G4double xs12 = 0;
481 G4double valueE11 = 0; << 588 G4double xs21 = 0;
482 << 589 G4double xs22 = 0;
483 G4double xs11 = 0; << 590
484 G4double xs12 = 0; << 591 if (particleDefinition == G4Electron::ElectronDefinition())
485 G4double xs21 = 0; << 592 {
486 G4double xs22 = 0; << 593 // Protection against out of boundary access
487 << 594 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
488 auto t2 = std::upper_bound(fTMapWithVec[Ma << 595 //
489 fTMapWithVec[Ma << 596
490 auto t1 = t2 - 1; << 597 // k should be in eV and energy transfer eV also
491 << 598
492 if (kSE <= fEMapWithVector[MatID][fpPartic << 599 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
493 && kSE <= fEMapWithVector[MatID][fpPar << 600
494 { << 601 std::vector<G4double>::iterator t1 = t2-1;
495 auto e12 = std::upper_bound(fEMapWithVec << 602
496 fEMapWithVec << 603 // SI : the following condition avoids situations where energyTransfer >last vector element
497 auto e11 = e12 - 1; << 604
498 << 605 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
499 auto e22 = std::upper_bound(fEMapWithVec << 606 {
500 fEMapWithVec << 607 std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
501 auto e21 = e22 - 1; << 608 std::vector<G4double>::iterator e11 = e12-1;
502 << 609
503 valueT1 = *t1; << 610 std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
504 valueT2 = *t2; << 611 std::vector<G4double>::iterator e21 = e22-1;
505 valueE21 = *e21; << 612
506 valueE22 = *e22; << 613 valueT1 =*t1;
507 valueE12 = *e12; << 614 valueT2 =*t2;
508 valueE11 = *e11; << 615 valueE21 =*e21;
509 << 616 valueE22 =*e22;
510 xs11 = diffCrossSectionData[MatID][fpPar << 617 valueE12 =*e12;
511 xs12 = diffCrossSectionData[MatID][fpPar << 618 valueE11 =*e11;
512 xs21 = diffCrossSectionData[MatID][fpPar << 619
513 xs22 = diffCrossSectionData[MatID][fpPar << 620 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
514 } << 621 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
515 << 622 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
516 G4double xsProduct = xs11 * xs12 * xs21 * << 623 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
517 << 624
518 if (xsProduct != 0.) { << 625 }
519 sigma = QuadInterpolator(valueE11, value << 626
520 valueT1, valueT << 627 }
521 } <<
522 } <<
523 628
524 return sigma; << 629 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
>> 630 if (xsProduct != 0.)
>> 631 {
>> 632 sigma = QuadInterpolator( valueE11, valueE12,
>> 633 valueE21, valueE22,
>> 634 xs11, xs12,
>> 635 xs21, xs22,
>> 636 valueT1, valueT2,
>> 637 k, energyTransfer);
>> 638 }
>> 639
>> 640 }
>> 641 return sigma;
525 } 642 }
526 643
527 //....oooOO0OOooo........oooOO0OOooo........oo 644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
528 645
529 G4double G4DNACPA100IonisationModel::Interpola << 646 G4double G4DNACPA100IonisationModel::Interpolate( G4double e1,
530 << 647 G4double e2,
>> 648 G4double e,
>> 649 G4double xs1,
>> 650 G4double xs2)
531 { 651 {
532 G4double value = 0.; <<
533 652
534 // Log-log interpolation by default << 653 G4double value = 0.;
535 654
536 if (e1 != 0 && e2 != 0 && (std::log10(e2) - << 655 // Log-log interpolation by default
537 G4double a = (std::log10(xs2) - std::log10 <<
538 G4double b = std::log10(xs2) - a * std::lo <<
539 G4double sigma = a * std::log10(e) + b; <<
540 value = (std::pow(10., sigma)); <<
541 } <<
542 656
543 // Switch to lin-lin interpolation << 657 if (e1!=0 && e2!=0 && (std::log10(e2)-std::log10(e1)) !=0 && !fasterCode && useDcs)
544 /* << 658 {
545 if ((e2-e1)!=0) << 659 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
546 { << 660 G4double b = std::log10(xs2) - a*std::log10(e2);
547 G4double d1 = xs1; << 661 G4double sigma = a*std::log10(e) + b;
548 G4double d2 = xs2; << 662 value = (std::pow(10.,sigma));
549 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1 << 663 }
550 } <<
551 */ <<
552 664
553 // Switch to log-lin interpolation for faste << 665 // Switch to lin-lin interpolation
>> 666 /*
>> 667 if ((e2-e1)!=0)
>> 668 {
>> 669 G4double d1 = xs1;
>> 670 G4double d2 = xs2;
>> 671 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
>> 672 }
>> 673 */
554 674
555 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 & << 675 // Switch to log-lin interpolation for faster code
556 G4double d1 = std::log10(xs1); <<
557 G4double d2 = std::log10(xs2); <<
558 value = std::pow(10., (d1 + (d2 - d1) * (e <<
559 } <<
560 676
561 // Switch to lin-lin interpolation for faste << 677 if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0 && fasterCode && useDcs )
562 // in case one of xs1 or xs2 (=cum proba) va << 678 {
>> 679 G4double d1 = std::log10(xs1);
>> 680 G4double d2 = std::log10(xs2);
>> 681 value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
>> 682 }
563 683
564 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) << 684 // Switch to lin-lin interpolation for faster code
565 G4double d1 = xs1; << 685 // in case one of xs1 or xs2 (=cum proba) value is zero
566 G4double d2 = xs2; << 686
567 value = (d1 + (d2 - d1) * (e - e1) / (e2 - << 687 if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0) && fasterCode && useDcs )
568 } << 688 {
569 return value; << 689 G4double d1 = xs1;
>> 690 G4double d2 = xs2;
>> 691 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
>> 692 }
>> 693
>> 694 /*
>> 695 G4cout
>> 696 << e1 << " "
>> 697 << e2 << " "
>> 698 << e << " "
>> 699 << xs1 << " "
>> 700 << xs2 << " "
>> 701 << value
>> 702 << G4endl;
>> 703 */
>> 704
>> 705 return value;
570 } 706 }
571 707
572 //....oooOO0OOooo........oooOO0OOooo........oo 708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
573 709
574 G4double G4DNACPA100IonisationModel::QuadInter << 710 G4double G4DNACPA100IonisationModel::QuadInterpolator(G4double e11, G4double e12,
575 << 711 G4double e21, G4double e22,
576 << 712 G4double xs11, G4double xs12,
577 << 713 G4double xs21, G4double xs22,
578 { << 714 G4double t1, G4double t2,
579 G4double interpolatedvalue1 = Interpolate(e1 << 715 G4double t, G4double e)
580 G4double interpolatedvalue2 = Interpolate(e2 << 716 {
581 G4double value = Interpolate(t1, t2, t, inte << 717 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
582 << 718 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
583 return value; << 719 G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
584 } << 720
585 << 721 return value;
586 G4double <<
587 G4DNACPA100IonisationModel::RandomizeEjectedEl <<
588 { <<
589 auto MatID = std::get<0>(info); <<
590 auto shell = std::get<2>(info); <<
591 G4double secondaryElectronKineticEnergy = <<
592 RandomTransferedEnergy(info) * eV - iStruc <<
593 if (secondaryElectronKineticEnergy < 0.) { <<
594 return 0.; <<
595 } <<
596 return secondaryElectronKineticEnergy; <<
597 } 722 }
598 723
599 G4double G4DNACPA100IonisationModel::RandomTra << 724 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 725
>> 726 G4int G4DNACPA100IonisationModel::RandomSelect(G4double k, const G4String& particle )
600 { 727 {
601 auto materialID = std::get<0>(info); << 728 G4int level = 0;
602 auto k = std::get<1>(info) / eV; // data ta <<
603 auto shell = std::get<2>(info); <<
604 G4double ejectedElectronEnergy = 0.; <<
605 G4double valueK1 = 0; <<
606 G4double valueK2 = 0; <<
607 G4double valueCumulCS21 = 0; <<
608 G4double valueCumulCS22 = 0; <<
609 G4double valueCumulCS12 = 0; <<
610 G4double valueCumulCS11 = 0; <<
611 G4double secElecE11 = 0; <<
612 G4double secElecE12 = 0; <<
613 G4double secElecE21 = 0; <<
614 G4double secElecE22 = 0; <<
615 729
616 if (k == fTMapWithVec[materialID][fpParticle << 730 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
617 k = k * (1. - 1e-12); << 731 pos = tableData.find(particle);
618 } <<
619 732
620 G4double random = G4UniformRand(); << 733 if (pos != tableData.end())
621 auto k2 = std::upper_bound(fTMapWithVec[mate << 734 {
622 fTMapWithVec[mate << 735 G4DNACrossSectionDataSet* table = pos->second;
623 auto k1 = k2 - 1; <<
624 736
625 if (random <= fProbaShellMap[materialID][fpP << 737 if (table != 0)
626 && random <= fProbaShellMap[materialID][ << 738 {
627 { << 739 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
628 auto cumulCS12 = << 740 const size_t n(table->NumberOfComponents());
629 std::upper_bound(fProbaShellMap[material << 741 size_t i(n);
630 fProbaShellMap[material << 742 G4double value = 0.;
631 auto cumulCS11 = cumulCS12 - 1; << 743
632 // Second one. << 744 //Verification
633 auto cumulCS22 = << 745 /*
634 std::upper_bound(fProbaShellMap[material << 746 G4double tmp=200*keV;
635 fProbaShellMap[material << 747 G4cout << table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) << G4endl;
636 auto cumulCS21 = cumulCS22 - 1; << 748 G4cout << table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) << G4endl;
637 << 749 G4cout << table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) << G4endl;
638 valueK1 = *k1; << 750 G4cout << table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m) << G4endl;
639 valueK2 = *k2; << 751 G4cout << table->GetComponent(4)->FindValue(tmp)/(1e-20*m*m) << G4endl;
640 valueCumulCS11 = *cumulCS11; << 752 G4cout <<
641 valueCumulCS12 = *cumulCS12; << 753 table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) +
642 valueCumulCS21 = *cumulCS21; << 754 table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) +
643 valueCumulCS22 = *cumulCS22; << 755 table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) +
644 << 756 table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m)
645 secElecE11 = fEnergySecondaryData[material << 757 << G4endl;
646 secElecE12 = fEnergySecondaryData[material << 758 abort();
647 secElecE21 = fEnergySecondaryData[material << 759 */
648 secElecE22 = fEnergySecondaryData[material << 760 //
649 << 761 //Dump
650 if (valueCumulCS11 == 0. && valueCumulCS12 << 762 //
651 auto interpolatedvalue2 = << 763 /*
652 Interpolate(valueCumulCS21, valueCumul << 764 G4double minEnergy = 10.985 * eV;
653 G4double valueNrjTransf = Interpolate(va << 765 G4double maxEnergy = 255955. * eV;
654 return valueNrjTransf; << 766 G4int nEnergySteps = 1000;
>> 767 G4double energy(minEnergy);
>> 768 G4double stpEnergy(std::pow(maxEnergy/energy, 1./static_cast<G4double>(nEnergySteps-1)));
>> 769 G4int step(nEnergySteps);
>> 770 system ("rm -rf ionisation-cpa100.out");
>> 771 FILE* myFile=fopen("ionisation-cpa100.out","a");
>> 772 while (step>0)
>> 773 {
>> 774 step--;
>> 775 fprintf (myFile,"%16.9le %16.9le %16.9le %16.9le %16.9le %16.9le %16.9le \n",
>> 776 energy/eV,
>> 777 table->GetComponent(0)->FindValue(energy)/(1e-20*m*m),
>> 778 table->GetComponent(1)->FindValue(energy)/(1e-20*m*m),
>> 779 table->GetComponent(2)->FindValue(energy)/(1e-20*m*m),
>> 780 table->GetComponent(3)->FindValue(energy)/(1e-20*m*m),
>> 781 table->GetComponent(4)->FindValue(energy)/(1e-20*m*m),
>> 782 table->GetComponent(0)->FindValue(energy)/(1e-20*m*m)+
>> 783 table->GetComponent(1)->FindValue(energy)/(1e-20*m*m)+
>> 784 table->GetComponent(2)->FindValue(energy)/(1e-20*m*m)+
>> 785 table->GetComponent(3)->FindValue(energy)/(1e-20*m*m)+
>> 786 table->GetComponent(4)->FindValue(energy)/(1e-20*m*m)
>> 787 );
>> 788 energy*=stpEnergy;
>> 789 }
>> 790 fclose (myFile);
>> 791 abort();
>> 792 */
>> 793 //
>> 794 // end of dump
>> 795 //
>> 796 // Test of diff XS
>> 797 // G4double nrj1 = .26827E+04; // in eV
>> 798 // G4double nrj2 = .57991E+03; // in eV
>> 799 // Shells run from 0 to 4
>> 800 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 0)/(1e-20*m*m) << G4endl;
>> 801 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 1)/(1e-20*m*m) << G4endl;
>> 802 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 2)/(1e-20*m*m) << G4endl;
>> 803 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 3)/(1e-20*m*m) << G4endl;
>> 804 // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 4)/(1e-20*m*m) << G4endl;
>> 805 // abort();
>> 806 //
>> 807
>> 808 while (i>0)
>> 809 {
>> 810 i--;
>> 811 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
>> 812 value += valuesBuffer[i];
>> 813 }
>> 814
>> 815 value *= G4UniformRand();
>> 816
>> 817 i = n;
>> 818
>> 819 while (i > 0)
>> 820 {
>> 821 i--;
>> 822 if (valuesBuffer[i] > value)
>> 823 {
>> 824 delete[] valuesBuffer;
>> 825
>> 826 return i;
>> 827 }
>> 828 value -= valuesBuffer[i];
>> 829 }
>> 830
>> 831 if (valuesBuffer) delete[] valuesBuffer;
>> 832
>> 833 }
>> 834 }
>> 835 else
>> 836 {
>> 837 G4Exception("G4DNACPA100IonisationModel::RandomSelect","em0002",
>> 838 FatalException,"Model not applicable to particle type.");
655 } 839 }
656 } <<
657 840
658 if (random > fProbaShellMap[materialID][fpPa << 841 return level;
659 auto cumulCS22 = << 842 }
660 std::upper_bound(fProbaShellMap[material <<
661 fProbaShellMap[material <<
662 auto cumulCS21 = cumulCS22 - 1; <<
663 valueK1 = *k1; <<
664 valueK2 = *k2; <<
665 valueCumulCS21 = *cumulCS21; <<
666 valueCumulCS22 = *cumulCS22; <<
667 843
668 secElecE21 = fEnergySecondaryData[material << 844 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
669 secElecE22 = fEnergySecondaryData[material <<
670 845
671 G4double interpolatedvalue2 = << 846 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs
672 Interpolate(valueCumulCS21, valueCumulCS << 847 (G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
>> 848 {
>> 849 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
673 850
674 G4double value = Interpolate(valueK1, valu << 851 G4double secondaryElectronKineticEnergy = 0.;
675 return value; <<
676 } <<
677 G4double nrjTransfProduct = secElecE11 * sec <<
678 852
679 if (nrjTransfProduct != 0.) { << 853 secondaryElectronKineticEnergy=
680 ejectedElectronEnergy = << 854 RandomTransferedEnergy(particleDefinition, k/eV, shell)*eV-waterStructure.IonisationEnergy(shell);
681 QuadInterpolator(valueCumulCS11, valueCu << 855
682 secElecE12, secElecE21, << 856 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
683 } << 857 if (secondaryElectronKineticEnergy<0.) return 0.;
684 return ejectedElectronEnergy; << 858 //
>> 859
>> 860 return secondaryElectronKineticEnergy;
685 } 861 }
686 862
687 //....oooOO0OOooo........oooOO0OOooo........oo 863 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
688 864
689 G4double << 865 G4double G4DNACPA100IonisationModel::RandomTransferedEnergy
690 G4DNACPA100IonisationModel::RandomizeEjectedEl << 866 (G4ParticleDefinition* particleDefinition,G4double k, G4int ionizationLevelIndex)
691 { 867 {
692 auto MatID = std::get<0>(info); <<
693 auto tt = std::get<1>(info); <<
694 auto shell = std::get<2>(info); <<
695 // ***** METHOD by M. C. Bordage ***** (opti <<
696 // Composition sampling method based on eq <<
697 868
698 //// Defining constants << 869 G4double random = G4UniformRand();
699 G4double alfa = 1. / 137; // fine structure <<
700 G4double e_charge = 1.6e-19; // electron ch <<
701 G4double e_mass = 9.1e-31; // electron mass <<
702 G4double c = 3e8; // speed of light in vacu <<
703 G4double mc2 = e_mass * c * c / e_charge; / <<
704 870
705 G4double BB = iStructure.IonisationEnergy(sh << 871 G4double nrj = 0.;
706 872
707 if (tt <= BB) return 0.; << 873 G4double valueK1 = 0;
>> 874 G4double valueK2 = 0;
>> 875 G4double valuePROB21 = 0;
>> 876 G4double valuePROB22 = 0;
>> 877 G4double valuePROB12 = 0;
>> 878 G4double valuePROB11 = 0;
>> 879
>> 880 G4double nrjTransf11 = 0;
>> 881 G4double nrjTransf12 = 0;
>> 882 G4double nrjTransf21 = 0;
>> 883 G4double nrjTransf22 = 0;
708 884
709 G4double b_prime = BB / mc2; // binding ene << 885 if (particleDefinition == G4Electron::ElectronDefinition())
710 G4double beta_b2 = 1. - 1. / ((1 + b_prime) << 886 {
711 887
712 //// Indicent energy << 888 // Protection against out of boundary access
713 //// tt is the incident electron energy << 889 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
>> 890 //
714 891
715 G4double t_prime = tt / mc2; // incident en << 892 // k should be in eV
716 G4double t = tt / BB; // reduced incident e <<
717 893
718 G4double D = (1 + 2 * t_prime) / ((1 + t_pri << 894 std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
719 G4double F = b_prime * b_prime / ((1 + t_pri <<
720 895
721 G4double beta_t2 = 1 - 1 / ((1 + t_prime) * << 896 std::vector<G4double>::iterator k1 = k2-1;
722 897
723 G4double PHI_R = std::cos(std::sqrt(alfa * a << 898 /*
724 * std::log(beta_t2 << 899 G4cout << "----> k=" << k
725 G4double G_R = std::log(beta_t2 / (1 - beta_ << 900 << " " << *k1
>> 901 << " " << *k2
>> 902 << " " << random
>> 903 << " " << ionizationLevelIndex
>> 904 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
>> 905 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
>> 906 << G4endl;
>> 907 */
726 908
727 G4double tplus1 = t + 1; << 909 // SI : the following condition avoids situations where random >last vector element
728 G4double tminus1 = t - 1; <<
729 G4double tplus12 = tplus1 * tplus1; <<
730 G4double ZH1max = 1 + F - (PHI_R * D * (2 * <<
731 G4double ZH2max = 1 - PHI_R * D / 4; <<
732 910
733 G4double A1_p = ZH1max * tminus1 / tplus1; << 911 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
734 G4double A2_p = ZH2max * tminus1 / (t * tplu << 912 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
735 G4double A3_p = ((tplus12 - 4) / tplus12) * <<
736 913
737 G4double AAA = A1_p + A2_p + A3_p; << 914 {
738 915
739 G4double AA1_R = A1_p / AAA; << 916 std::vector<G4double>::iterator prob12 =
740 G4double AA2_R = (A1_p + A2_p) / AAA; << 917 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
>> 918 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
741 919
742 G4int FF = 0; << 920 std::vector<G4double>::iterator prob11 = prob12-1;
743 G4double fx = 0; <<
744 G4double gx = 0; <<
745 G4double gg = 0; <<
746 G4double wx = 0; <<
747 921
748 G4double r1 = 0; <<
749 G4double r2 = 0; <<
750 G4double r3 = 0; <<
751 922
752 // << 923 std::vector<G4double>::iterator prob22 =
>> 924 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
>> 925 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
753 926
754 do { << 927 std::vector<G4double>::iterator prob21 = prob22-1;
755 r1 = G4UniformRand(); <<
756 r2 = G4UniformRand(); <<
757 r3 = G4UniformRand(); <<
758 <<
759 if (r1 > AA2_R) <<
760 FF = 3; <<
761 else if ((r1 > AA1_R) && (r1 < AA2_R)) <<
762 FF = 2; <<
763 else <<
764 FF = 1; <<
765 928
766 switch (FF) { << 929 valueK1 =*k1;
767 case 1: { << 930 valueK2 =*k2;
768 fx = r2 * tminus1 / tplus1; << 931 valuePROB21 =*prob21;
769 wx = 1 / (1 - fx) - 1; << 932 valuePROB22 =*prob22;
770 gg = PHI_R * D * (wx + 1) / tplus1; << 933 valuePROB12 =*prob12;
771 gx = 1 - gg; << 934 valuePROB11 =*prob11;
772 gx = gx - gg * (wx + 1) / (2 * (t - wx << 935
773 gx = gx + F * (wx + 1) * (wx + 1); << 936
774 gx = gx / ZH1max; << 937 /*
775 break; << 938 G4cout << " " << random << " " << valuePROB11 << " "
776 } << 939 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
>> 940 */
>> 941
>> 942 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
>> 943 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
>> 944 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
>> 945 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
>> 946
>> 947 /*
>> 948 G4cout << " " << ionizationLevelIndex << " "
>> 949 << random << " " <<valueK1 << " " << valueK2 << G4endl;
>> 950
>> 951 G4cout << " " << random << " " << nrjTransf11 << " "
>> 952 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
>> 953 */
777 954
778 case 2: { <<
779 fx = tplus1 + r2 * tminus1; <<
780 wx = t * tminus1 * r2 / fx; <<
781 gx = 1 - (PHI_R * D * (t - wx) / (2 * <<
782 gx = gx / ZH2max; <<
783 break; <<
784 } 955 }
785 956
786 case 3: { << 957
787 fx = 1 - r2 * (tplus12 - 4) / tplus12; << 958 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
788 wx = std::sqrt(1 / fx) - 1; << 959
789 gg = (wx + 1) / (t - wx); << 960 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
790 gx = (1 + gg * gg * gg) / 2; << 961 {
791 break; << 962
>> 963 std::vector<G4double>::iterator prob22 =
>> 964
>> 965 std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
>> 966 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
>> 967
>> 968 std::vector<G4double>::iterator prob21 = prob22-1;
>> 969
>> 970 valueK1 =*k1;
>> 971 valueK2 =*k2;
>> 972 valuePROB21 =*prob21;
>> 973 valuePROB22 =*prob22;
>> 974
>> 975 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
>> 976
>> 977 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
>> 978 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
>> 979
>> 980 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
>> 981
>> 982 // zero is explicitely set
>> 983
>> 984 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
>> 985
>> 986 /*
>> 987 G4cout << " " << ionizationLevelIndex << " "
>> 988 << random << " " <<valueK1 << " " << valueK2 << G4endl;
>> 989
>> 990 G4cout << " " << random << " " << nrjTransf11 << " "
>> 991 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
>> 992
>> 993 G4cout << "ici" << " " << value << G4endl;
>> 994 */
>> 995
>> 996 return value;
792 } 997 }
793 } // switch <<
794 998
795 } while (r3 > gx); << 999 }
>> 1000
>> 1001 // End electron case
796 1002
797 return wx * BB; << 1003 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
>> 1004
>> 1005 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
>> 1006
>> 1007 if (nrjTransfProduct != 0.)
>> 1008 {
>> 1009 nrj = QuadInterpolator( valuePROB11, valuePROB12,
>> 1010 valuePROB21, valuePROB22,
>> 1011 nrjTransf11, nrjTransf12,
>> 1012 nrjTransf21, nrjTransf22,
>> 1013 valueK1, valueK2,
>> 1014 k, random);
>> 1015 }
>> 1016
>> 1017 //G4cout << nrj << endl;
>> 1018
>> 1019 return nrj ;
798 } 1020 }
799 1021
800 //....oooOO0OOooo........oooOO0OOooo........oo << 1022 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1023
>> 1024 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCompositionSampling
>> 1025 (G4ParticleDefinition*, G4double tt, G4int shell)
>> 1026 {
>> 1027 //G4cout << "*** Rejection method for " << " " << particleDefinition->GetParticleName() << G4endl;
>> 1028
>> 1029 // ***** METHOD 1 ***** (sequential)
>> 1030 /*
>> 1031
>> 1032 // ww is KINETIC ENERGY OF SECONDARY ELECTRON
>> 1033 G4double un=1.;
>> 1034 G4double deux=2.;
>> 1035
>> 1036 G4double bb = waterStructure.IonisationEnergy(shell);
>> 1037 G4double uu = waterStructure.UEnergy(shell);
>> 1038
>> 1039 if (tt<=bb) return 0.;
>> 1040
>> 1041 G4double t = tt/bb;
>> 1042 G4double u = uu/bb;
>> 1043 G4double tp1 = t + un;
>> 1044 G4double tu1 = t + u + un;
>> 1045 G4double tm1 = t - un;
>> 1046 G4double tp12 = tp1 * tp1;
>> 1047 G4double dlt = std::log(t);
>> 1048
>> 1049 G4double a1 = t * tm1 / tu1 / tp12;
>> 1050 G4double a2 = tm1 / tu1 / t / tp1 / deux;
>> 1051 G4double a3 = dlt * (tp12 - deux * deux ) / tu1 / tp12;
>> 1052 G4double ato = a1 + a2 + a3;
>> 1053
>> 1054 // 15
>> 1055
>> 1056 G4double r1 =G4UniformRand();
>> 1057 G4double r2 =G4UniformRand();
>> 1058 G4double r3 =G4UniformRand();
>> 1059
>> 1060 while (r1<=a1/ato)
>> 1061 {
>> 1062 G4double fx1=r2*tm1/tp1;
>> 1063 G4double wx1=un/(un-fx1)-un;
>> 1064 G4double gx1=(t-wx1)/t;
>> 1065 if(r3 <= gx1) return wx1*bb;
>> 1066
>> 1067 r1 =G4UniformRand();
>> 1068 r2 =G4UniformRand();
>> 1069 r3 =G4UniformRand();
>> 1070
>> 1071 }
>> 1072
>> 1073 // 20
>> 1074
>> 1075 while (r1<=(a1+a2)/ato)
>> 1076 {
>> 1077 G4double fx2=tp1+r2*tm1;
>> 1078 G4double wx2=t-t*tp1/fx2;
>> 1079 G4double gx2=deux*(un-(t-wx2)/tp1);
>> 1080 if(r3 <= gx2) return wx2*bb;
>> 1081
>> 1082 // REPEAT 15
>> 1083 r1 =G4UniformRand();
>> 1084 r2 =G4UniformRand();
>> 1085 r3 =G4UniformRand();
>> 1086
>> 1087 while (r1<=a1/ato)
>> 1088 {
>> 1089 G4double fx1=r2*tm1/tp1;
>> 1090 G4double wx1=un/(un-fx1)-un;
>> 1091 G4double gx1=(t-wx1)/t;
>> 1092 if(r3 <= gx1) return wx1*bb;
>> 1093 r1 =G4UniformRand();
>> 1094 r2 =G4UniformRand();
>> 1095 r3 =G4UniformRand();
>> 1096 }
>> 1097 // END 15
>> 1098
>> 1099 }
>> 1100
>> 1101 // 30
>> 1102
>> 1103 G4double wx3=std::sqrt(un/(un-r2*(tp12-deux*deux)/tp12))-un;
>> 1104 G4double gg3=(wx3+un)/(t-wx3);
>> 1105 G4double gx3=(un+gg3*gg3*gg3)/deux;
>> 1106
>> 1107 while (r3>gx3)
>> 1108 {
>> 1109
>> 1110 // 15
>> 1111
>> 1112 r1 =G4UniformRand();
>> 1113 r2 =G4UniformRand();
>> 1114 r3 =G4UniformRand();
>> 1115
>> 1116 while (r1<=a1/ato)
>> 1117 {
>> 1118 G4double fx1=r2*tm1/tp1;
>> 1119 G4double wx1=un/(un-fx1)-un;
>> 1120 G4double gx1=(t-wx1)/t;
>> 1121 if(r3 <= gx1) return wx1*bb;
>> 1122
>> 1123 r1 =G4UniformRand();
>> 1124 r2 =G4UniformRand();
>> 1125 r3 =G4UniformRand();
801 1126
802 void G4DNACPA100IonisationModel::ReadDiffCSFil <<
803 <<
804 <<
805 { <<
806 const char* path = G4FindDataDir("G4LEDATA") <<
807 if (path == nullptr) { <<
808 G4Exception("G4DNACPA100IonisationModel::R <<
809 "G4LEDATA environment variable <<
810 return; <<
811 } 1127 }
812 1128
813 std::ostringstream fullFileName; << 1129 // 20
814 fullFileName << path << "/" << file << ".dat << 1130
>> 1131 while (r1<=(a1+a2)/ato)
>> 1132 {
>> 1133 G4double fx2=tp1+r2*tm1;
>> 1134 G4double wx2=t-t*tp1/fx2;
>> 1135 G4double gx2=deux*(un-(t-wx2)/tp1);
>> 1136 if(r3 <= gx2)return wx2*bb;
>> 1137
>> 1138 // REPEAT 15
>> 1139 r1 =G4UniformRand();
>> 1140 r2 =G4UniformRand();
>> 1141 r3 =G4UniformRand();
>> 1142
>> 1143 while (r1<=a1/ato)
>> 1144 {
>> 1145 G4double fx1=r2*tm1/tp1;
>> 1146 G4double wx1=un/(un-fx1)-un;
>> 1147 G4double gx1=(t-wx1)/t;
>> 1148 if(r3 <= gx1) return wx1*bb;
>> 1149
>> 1150 r1 =G4UniformRand();
>> 1151 r2 =G4UniformRand();
>> 1152 r3 =G4UniformRand();
>> 1153 }
>> 1154 //
815 1155
816 std::ifstream diffCrossSection(fullFileName. <<
817 std::stringstream endPath; <<
818 if (!diffCrossSection) { <<
819 endPath << "Missing data file: " << file; <<
820 G4Exception("G4DNACPA100IonisationModel::I <<
821 endPath.str().c_str()); <<
822 } 1156 }
823 1157
824 // load data from the file << 1158 wx3=std::sqrt(un/(un-r2*(tp12-deux*deux)/tp12))-un;
825 fTMapWithVec[materialID][p].push_back(0.); << 1159 gg3=(wx3+un)/(t-wx3);
>> 1160 gx3=(un+gg3*gg3*gg3)/deux;
826 1161
827 G4String line; << 1162 }
828 1163
829 while (!diffCrossSection.eof()) { << 1164 //
830 G4double T, E; <<
831 diffCrossSection >> T >> E; <<
832 1165
833 if (T != fTMapWithVec[materialID][p].back( << 1166 return wx3*bb;
834 fTMapWithVec[materialID][p].push_back(T) << 1167 */
835 } <<
836 1168
837 // T is incident energy, E is the energy t << 1169 // ***** METHOD by M. C. Bordage ***** (optimized)
838 if (T != fTMapWithVec[materialID][p].back( <<
839 fTMapWithVec[materialID][p].push_back(T) <<
840 } <<
841 1170
842 auto eshell = (G4int)iStructure.NumberOfLe << 1171 G4double un=1.;
843 for (G4int shell = 0; shell < eshell; ++sh << 1172 G4double deux=2.;
844 diffCrossSection >> diffCrossSectionData << 1173
845 if (fasterCode) { << 1174 G4double bb = waterStructure.IonisationEnergy(shell);
846 fEnergySecondaryData[materialID][p][sh << 1175 G4double uu = waterStructure.UEnergy(shell);
847 [diffCrossSectionD << 1176
>> 1177 if (tt<=bb) return 0.;
>> 1178
>> 1179 G4double t = tt/bb;
>> 1180 G4double u = uu/bb;
>> 1181 G4double tp1 = t + un;
>> 1182 G4double tu1 = t + u + un;
>> 1183 G4double tm1 = t - un;
>> 1184 G4double tp12 = tp1 * tp1;
>> 1185 G4double dlt = std::log(t);
>> 1186
>> 1187 G4double a1 = t * tm1 / tu1 / tp12;
>> 1188 G4double a2 = tm1 / tu1 / t / tp1 / deux;
>> 1189 G4double a3 = dlt * (tp12 - deux * deux ) / tu1 / tp12;
>> 1190 G4double ato = a1 + a2 + a3;
>> 1191
>> 1192 G4double A1 = a1/ato;
>> 1193 G4double A2 = (a1+a2)/ato;
>> 1194 G4int F = 0;
>> 1195 G4double fx=0;
>> 1196 G4double gx=0;
>> 1197 G4double gg=0;
>> 1198 G4double wx=0;
>> 1199
>> 1200 G4double r1=0;
>> 1201 G4double r2=0;
>> 1202 G4double r3=0;
>> 1203
>> 1204 //
>> 1205
>> 1206 do
>> 1207 {
>> 1208 r1 =G4UniformRand();
>> 1209 r2 =G4UniformRand();
>> 1210 r3 =G4UniformRand();
>> 1211
>> 1212 if (r1>A2)
>> 1213 F=3;
>> 1214 else if ((r1>A1) && (r1< A2))
>> 1215 F=2;
>> 1216 else
>> 1217 F=1;
>> 1218
>> 1219 switch (F)
>> 1220 {
>> 1221 case 1:
>> 1222 {
>> 1223 fx=r2*tm1/tp1;
>> 1224 wx=un/(un-fx)-un;
>> 1225 gx=(t-wx)/t;
>> 1226 break;
>> 1227 }
>> 1228
>> 1229 case 2:
>> 1230 {
>> 1231 fx=tp1+r2*tm1;
>> 1232 wx=t-t*tp1/fx;
>> 1233 gx=deux*(un-(t-wx)/tp1);
>> 1234 break;
>> 1235 }
>> 1236
>> 1237 case 3:
>> 1238 {
>> 1239 fx=un-r2*(tp12-deux*deux)/tp12;
>> 1240 wx=sqrt(un/fx)-un;
>> 1241 gg=(wx+un)/(t-wx);
>> 1242 gx=(un+gg*gg*gg)/deux;
>> 1243 break;
>> 1244 }
>> 1245 } // switch
>> 1246
>> 1247 } while (r3>gx);
>> 1248
>> 1249 return wx*bb;
848 1250
849 fProbaShellMap[materialID][p][shell][T <<
850 diffCrossSectionData[materialID][p][ <<
851 } <<
852 else { <<
853 diffCrossSectionData[materialID][p][sh <<
854 fEMapWithVector[materialID][p][T].push <<
855 } <<
856 } <<
857 } <<
858 } 1251 }
859 1252