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