| 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 // Authors: S. Meylan and C. Villagrasa (IRSN, 26 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Models come from 27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459- 28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
29 // 29 //
30 30
31 #include "G4DNAPTBElasticModel.hh" 31 #include "G4DNAPTBElasticModel.hh"
32 <<
33 #include "G4DNAChampionElasticModel.hh" 32 #include "G4DNAChampionElasticModel.hh"
34 #include "G4DNAMaterialManager.hh" << 33 #include "G4PhysicalConstants.hh"
>> 34 #include "G4SystemOfUnits.hh"
35 #include "G4DNAMolecularMaterial.hh" 35 #include "G4DNAMolecularMaterial.hh"
36 #include "G4Proton.hh" 36 #include "G4Proton.hh"
37 #include "G4SystemOfUnits.hh" <<
38 37
39 G4DNAPTBElasticModel::G4DNAPTBElasticModel(con << 38 G4DNAPTBElasticModel::G4DNAPTBElasticModel(const G4String& applyToMaterial, const G4ParticleDefinition*,
40 con << 39 const G4String& nam)
41 : G4VDNAModel(nam, applyToMaterial) << 40 : G4VDNAModel(nam, applyToMaterial)
42 { << 41 {
43 if (verboseLevel > 0) { << 42 fKillBelowEnergy = 10*eV; // will be override by the limits defined for each material
44 G4cout << "PTB Elastic model is constructe << 43
45 } << 44 verboseLevel= 0;
46 fpTHF = G4Material::GetMaterial("THF", false << 45 // Verbosity scale:
47 fpPY = G4Material::GetMaterial("PY", false); << 46 // 0 = nothing
48 fpPU = G4Material::GetMaterial("PU", false); << 47 // 1 = warning for energy non-conservation
49 fpTMP = G4Material::GetMaterial("TMP", false << 48 // 2 = details of energy budget
50 fpG4_WATER = G4Material::GetMaterial("G4_WAT << 49 // 3 = calculation of cross sections, file openings, sampling of atoms
51 fpBackbone_THF = G4Material::GetMaterial("ba << 50 // 4 = entering in methods
52 fpCytosine_PY = G4Material::GetMaterial("cyt << 51
53 fpThymine_PY = G4Material::GetMaterial("thym << 52 if( verboseLevel>0 )
54 fpAdenine_PU = G4Material::GetMaterial("aden << 53 {
55 fpBackbone_TMP = G4Material::GetMaterial("ba << 54 G4cout << "PTB Elastic model is constructed " << G4endl;
56 fpGuanine_PU = G4Material::GetMaterial("guan << 55 }
57 fpN2 = G4Material::GetMaterial("N2", false); <<
58 } 56 }
59 57
60 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 59
62 void G4DNAPTBElasticModel::Initialise(const G4 << 60 G4DNAPTBElasticModel::~G4DNAPTBElasticModel()
63 const G4 <<
64 { 61 {
65 if (isInitialised) { <<
66 return; <<
67 } <<
68 if (verboseLevel > 3) <<
69 { <<
70 G4cout << "Calling G4DNAPTBElasticModel::I <<
71 } <<
72 if (particle != G4Electron::ElectronDefiniti <<
73 std::ostringstream oss; <<
74 oss << " Model is not applied for this par <<
75 G4Exception("G4DNAPTBElasticModel::G4DNAPT <<
76 oss.str().c_str()); <<
77 } <<
78 G4double scaleFactor = 1e-16 * cm * cm; <<
79 //****************************************** <<
80 // Cross section data <<
81 //****************************************** <<
82 <<
83 std::size_t index; <<
84 // MPietrzak, adding paths for N2 <<
85 if (fpN2 != nullptr) { <<
86 index = fpN2->GetIndex(); <<
87 AddCrossSectionData(index, particle, "dna/ <<
88 "dna/sigmadiff_cumulat <<
89 SetLowELimit(index, particle, 10 * eV); <<
90 SetHighELimit(index, particle, 1.02 * MeV) <<
91 } <<
92 // MPietrzak <<
93 <<
94 if (fpTHF != nullptr) { <<
95 index = fpTHF->GetIndex(); <<
96 AddCrossSectionData(index, particle, "dna/ <<
97 "dna/sigmadiff_cumulat <<
98 SetLowELimit(index, particle, 10 * eV); <<
99 SetHighELimit(index, particle, 1 * keV); <<
100 } <<
101 <<
102 if (fpPY != nullptr) { <<
103 index = fpPY->GetIndex(); <<
104 AddCrossSectionData(index, particle, "dna/ <<
105 "dna/sigmadiff_cumulat <<
106 SetLowELimit(index, particle, 10 * eV); <<
107 SetHighELimit(index, particle, 1 * keV); <<
108 } <<
109 <<
110 if (fpPU != nullptr) { <<
111 index = fpPU->GetIndex(); <<
112 AddCrossSectionData(index, particle, "dna/ <<
113 "dna/sigmadiff_cumulat <<
114 SetLowELimit(index, particle, 10 * eV); <<
115 SetHighELimit(index, particle, 1 * keV); <<
116 } <<
117 <<
118 if (fpTMP != nullptr) { <<
119 index = fpTMP->GetIndex(); <<
120 AddCrossSectionData(index, particle, "dna/ <<
121 "dna/sigmadiff_cumulat <<
122 SetLowELimit(index, particle, 10 * eV); <<
123 SetHighELimit(index, particle, 1 * keV); <<
124 } <<
125 //???? <<
126 if (fpG4_WATER != nullptr) { <<
127 index = fpG4_WATER->GetIndex(); <<
128 AddCrossSectionData(index, particle, "dna/ <<
129 "dna/sigmadiff_cumulat <<
130 SetLowELimit(index, particle, 10 * eV); <<
131 SetHighELimit(index, particle, 1 * keV); <<
132 } <<
133 // DNA materials <<
134 // <<
135 if (fpBackbone_THF != nullptr) { <<
136 index = fpBackbone_THF->GetIndex(); <<
137 AddCrossSectionData(index, particle, "dna/ <<
138 "dna/sigmadiff_cumulat <<
139 SetLowELimit(index, particle, 10 * eV); <<
140 SetHighELimit(index, particle, 1 * keV); <<
141 } <<
142 <<
143 if (fpCytosine_PY != nullptr) { <<
144 index = fpCytosine_PY->GetIndex(); <<
145 AddCrossSectionData(index, particle, "dna/ <<
146 "dna/sigmadiff_cumulat <<
147 SetLowELimit(index, particle, 10 * eV); <<
148 SetHighELimit(index, particle, 1 * keV); <<
149 } <<
150 <<
151 if (fpThymine_PY != nullptr) { <<
152 index = fpThymine_PY->GetIndex(); <<
153 AddCrossSectionData(index, particle, "dna/ <<
154 "dna/sigmadiff_cumulat <<
155 SetLowELimit(index, particle, 10 * eV); <<
156 SetHighELimit(index, particle, 1 * keV); <<
157 } <<
158 <<
159 if (fpAdenine_PU != nullptr) { <<
160 index = fpAdenine_PU->GetIndex(); <<
161 AddCrossSectionData(index, particle, "dna/ <<
162 "dna/sigmadiff_cumulat <<
163 SetLowELimit(index, particle, 10 * eV); <<
164 SetHighELimit(index, particle, 1 * keV); <<
165 } <<
166 if (fpGuanine_PU != nullptr) { <<
167 index = fpGuanine_PU->GetIndex(); <<
168 AddCrossSectionData(index, particle, "dna/ <<
169 "dna/sigmadiff_cumulat <<
170 SetLowELimit(index, particle, 10 * eV); <<
171 SetHighELimit(index, particle, 1 * keV); <<
172 } <<
173 <<
174 if (fpBackbone_TMP != nullptr) { <<
175 index = fpBackbone_TMP->GetIndex(); <<
176 AddCrossSectionData(index, particle, "dna/ <<
177 "dna/sigmadiff_cumulat <<
178 SetLowELimit(index, particle, 10 * eV); <<
179 SetHighELimit(index, particle, 1 * keV); <<
180 } <<
181 62
182 if (!G4DNAMaterialManager::Instance()->IsLoc <<
183 // Load the data <<
184 LoadCrossSectionData(particle); <<
185 G4DNAMaterialManager::Instance()->SetMaste <<
186 fpModelData = this; <<
187 } <<
188 else { <<
189 auto dataModel = dynamic_cast<G4DNAPTBElas <<
190 G4DNAMaterialManager::Instance()->GetMod <<
191 if (dataModel == nullptr) { <<
192 G4cout << "G4DNAPTBElasticModel::Initial <<
193 G4Exception("G4DNAPTBElasticModel::Initi <<
194 "not good modelData"); <<
195 } <<
196 else { <<
197 fpModelData = dataModel; <<
198 } <<
199 } <<
200 <<
201 if (verboseLevel > 2) { <<
202 G4cout << "Loaded cross section files for <<
203 } <<
204 <<
205 fParticleChangeForGamma = GetParticleChangeF <<
206 isInitialised = true; <<
207 } <<
208 <<
209 //....oooOO0OOooo........oooOO0OOooo........oo <<
210 <<
211 void G4DNAPTBElasticModel::ReadDiffCSFile(cons <<
212 cons <<
213 cons <<
214 { <<
215 // Method to read and save the information c <<
216 // This method is not yet standard. <<
217 <<
218 // get the path of the G4LEDATA data folder <<
219 const char* path = G4FindDataDir("G4LEDATA") <<
220 // if it is not found then quit and print er <<
221 if (path == nullptr) { <<
222 G4Exception("G4DNAPTBElasticModel::ReadAll <<
223 "G4LEDATA environment variable <<
224 return; <<
225 } <<
226 <<
227 // build the fullFileName path of the data f <<
228 std::ostringstream fullFileName; <<
229 fullFileName << path << "/" << file << ".dat <<
230 <<
231 // open the data file <<
232 std::ifstream diffCrossSection(fullFileName. <<
233 // error if file is not there <<
234 std::stringstream endPath; <<
235 if (!diffCrossSection) { <<
236 endPath << "Missing data file: " << file; <<
237 G4Exception("G4DNAPTBElasticModel::Initial <<
238 endPath.str().c_str()); <<
239 } <<
240 <<
241 tValuesVec[materialName][particleName].push_ <<
242 <<
243 G4String line; <<
244 <<
245 // read the file line by line until we reach <<
246 while (std::getline(diffCrossSection, line)) <<
247 // check if the line is comment or empty <<
248 // <<
249 std::istringstream testIss(line); <<
250 G4String test; <<
251 testIss >> test; <<
252 // check first caracter to determine if fo <<
253 if (test == "#") { <<
254 // skip the line by beginning a new whil <<
255 continue; <<
256 } <<
257 // check if line is empty <<
258 if (line.empty()) { <<
259 // skip the line by beginning a new whil <<
260 continue; <<
261 } <<
262 // <<
263 // end of the check <<
264 <<
265 // transform the line into a iss <<
266 std::istringstream iss(line); <<
267 <<
268 // Variables to be filled by the input fil <<
269 G4double tDummy; <<
270 G4double eDummy; <<
271 <<
272 // fill the variables with the content of <<
273 iss >> tDummy >> eDummy; <<
274 <<
275 // SI : mandatory Vecm initialization <<
276 <<
277 // Fill two vectors contained in maps of t <<
278 // [materialName][particleName]=vector <<
279 // [materialName][particleName][T]=vector <<
280 // to list all the incident energies (tVal <<
281 // file <<
282 // <<
283 // Check if we already have the current T <<
284 // If not then add it <<
285 if (tDummy != tValuesVec[materialName][par <<
286 // Add the current T value <<
287 tValuesVec[materialName][particleName].p <<
288 // Make it correspond to a default zero <<
289 eValuesVect[materialName][particleName][ <<
290 } <<
291 <<
292 // Put the differential cross section valu <<
293 // map <<
294 iss >> diffCrossSectionData[materialName][ <<
295 <<
296 // If the current E value (eDummy) is diff <<
297 // then add it to the vector <<
298 if (eDummy != eValuesVect[materialName][pa <<
299 eValuesVect[materialName][particleName][ <<
300 } <<
301 } <<
302 } <<
303 <<
304 //....oooOO0OOooo........oooOO0OOooo........oo <<
305 <<
306 G4double G4DNAPTBElasticModel::CrossSectionPer <<
307 <<
308 <<
309 { <<
310 if (verboseLevel > 3){ <<
311 G4cout << "Calling CrossSectionPerVolume() <<
312 } <<
313 <<
314 // Get the name of the current particle <<
315 const G4String& particleName = p->GetParticl <<
316 const std::size_t& materialID = pMaterial->G <<
317 <<
318 // set killBelowEnergy value for current mat <<
319 fKillBelowEnergy = fpModelData->GetLowELimit <<
320 // initialise the return value (cross sectio <<
321 G4double sigma = 0.; <<
322 <<
323 // check if we are below the high energy lim <<
324 if (ekin < fpModelData->GetHighELimit(materi <<
325 // This is used to kill the particle if it <<
326 // If the energy is lower then we return a <<
327 // method will be called for sure. SampleS <<
328 // simulation. <<
329 // <<
330 // SI : XS must not be zero otherwise samp <<
331 if (ekin < fKillBelowEnergy) { <<
332 return DBL_MAX; <<
333 } <<
334 <<
335 // Get the tables with the cross section d <<
336 auto tableData = fpModelData->GetData(); <<
337 if ((*tableData)[materialID][p] == nullptr <<
338 G4Exception("G4DNAPTBElasticModel::Cross <<
339 "No model is registered"); <<
340 } <<
341 // Retrieve the cross section value <<
342 sigma = (*tableData)[materialID][p]->FindV <<
343 } <<
344 <<
345 if (verboseLevel > 2) { <<
346 G4cout << "_______________________________ <<
347 G4cout << "°°° G4DNAPTBElasticModel - X <<
348 G4cout << "°°° Kinetic energy(eV)=" << <<
349 G4cout << "°°° Cross section per molecu <<
350 G4cout << "°°° G4DNAPTBElasticModel - X <<
351 } <<
352 <<
353 // Return the cross section <<
354 auto MolDensity = <<
355 (*G4DNAMolecularMaterial::Instance()->GetN <<
356 return sigma * MolDensity; <<
357 } 63 }
358 64
359 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
360 66
361 void G4DNAPTBElasticModel::SampleSecondaries(s << 67 void G4DNAPTBElasticModel::Initialise(const G4ParticleDefinition* particle,
362 c << 68 const G4DataVector& /*cuts*/, G4ParticleChangeForGamma*)
363 c <<
364 G <<
365 { 69 {
366 if (verboseLevel > 3) { << 70 if (verboseLevel > 3)
367 G4cout << "Calling SampleSecondaries() of << 71 G4cout << "Calling G4DNAPTBElasticModel::Initialise()" << G4endl;
368 } <<
369 72
370 G4double electronEnergy0 = aDynamicElectron- << 73 G4double scaleFactor = 1e-16*cm*cm;
371 const std::size_t& materialID = couple->GetM <<
372 auto p = aDynamicElectron->GetParticleDefini <<
373 74
374 // set killBelowEnergy value for material << 75 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
375 fKillBelowEnergy = fpModelData->GetLowELimit <<
376 76
377 // If the particle (electron here) energy is << 77 //*******************************************************
378 // simulation << 78 // Cross section data
379 if (electronEnergy0 < fKillBelowEnergy) { << 79 //*******************************************************
380 fParticleChangeForGamma->SetProposedKineti << 80
381 fParticleChangeForGamma->ProposeTrackStatu << 81 if(particle == electronDef)
382 fParticleChangeForGamma->ProposeLocalEnerg << 82 {
383 } << 83 G4String particleName = particle->GetParticleName();
384 // If we are above the kill limite and below << 84
385 else if (electronEnergy0 >= fKillBelowEnergy << 85 AddCrossSectionData("THF",
386 // Random sampling of the cosTheta << 86 particleName,
387 G4double cosTheta = fpModelData->Randomize << 87 "dna/sigma_elastic_e-_PTB_THF",
>> 88 "dna/sigmadiff_cumulated_elastic_e-_PTB_THF",
>> 89 scaleFactor);
>> 90 SetLowELimit("THF", particleName, 10*eV);
>> 91 SetHighELimit("THF", particleName, 1*keV);
>> 92
>> 93 AddCrossSectionData("PY",
>> 94 particleName,
>> 95 "dna/sigma_elastic_e-_PTB_PY",
>> 96 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
>> 97 scaleFactor);
>> 98 SetLowELimit("PY", particleName, 10*eV);
>> 99 SetHighELimit("PY", particleName, 1*keV);
>> 100
>> 101 AddCrossSectionData("PU",
>> 102 particleName,
>> 103 "dna/sigma_elastic_e-_PTB_PU",
>> 104 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
>> 105 scaleFactor);
>> 106 SetLowELimit("PU", particleName, 10*eV);
>> 107 SetHighELimit("PU", particleName, 1*keV);
>> 108
>> 109 AddCrossSectionData("TMP",
>> 110 particleName,
>> 111 "dna/sigma_elastic_e-_PTB_TMP",
>> 112 "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP",
>> 113 scaleFactor);
>> 114 SetLowELimit("TMP", particleName, 10*eV);
>> 115 SetHighELimit("TMP", particleName, 1*keV);
>> 116
>> 117 AddCrossSectionData("G4_WATER",
>> 118 particleName,
>> 119 "dna/sigma_elastic_e_champion",
>> 120 "dna/sigmadiff_cumulated_elastic_e_champion",
>> 121 scaleFactor);
>> 122 SetLowELimit("G4_WATER", particleName, 10*eV);
>> 123 SetHighELimit("G4_WATER", particleName, 1*keV);
>> 124
>> 125 // DNA materials
>> 126 //
>> 127 AddCrossSectionData("backbone_THF",
>> 128 particleName,
>> 129 "dna/sigma_elastic_e-_PTB_THF",
>> 130 "dna/sigmadiff_cumulated_elastic_e-_PTB_THF",
>> 131 scaleFactor*33./30);
>> 132 SetLowELimit("backbone_THF", particleName, 10*eV);
>> 133 SetHighELimit("backbone_THF", particleName, 1*keV);
>> 134
>> 135 AddCrossSectionData("cytosine_PY",
>> 136 particleName,
>> 137 "dna/sigma_elastic_e-_PTB_PY",
>> 138 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
>> 139 scaleFactor*42./30);
>> 140 SetLowELimit("cytosine_PY", particleName, 10*eV);
>> 141 SetHighELimit("cytosine_PY", particleName, 1*keV);
>> 142
>> 143 AddCrossSectionData("thymine_PY",
>> 144 particleName,
>> 145 "dna/sigma_elastic_e-_PTB_PY",
>> 146 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
>> 147 scaleFactor*48./30);
>> 148 SetLowELimit("thymine_PY", particleName, 10*eV);
>> 149 SetHighELimit("thymine_PY", particleName, 1*keV);
>> 150
>> 151 AddCrossSectionData("adenine_PU",
>> 152 particleName,
>> 153 "dna/sigma_elastic_e-_PTB_PU",
>> 154 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
>> 155 scaleFactor*50./44);
>> 156 SetLowELimit("adenine_PU", particleName, 10*eV);
>> 157 SetHighELimit("adenine_PU", particleName, 1*keV);
>> 158
>> 159 AddCrossSectionData("guanine_PU",
>> 160 particleName,
>> 161 "dna/sigma_elastic_e-_PTB_PU",
>> 162 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
>> 163 scaleFactor*56./44);
>> 164 SetLowELimit("guanine_PU", particleName, 10*eV);
>> 165 SetHighELimit("guanine_PU", particleName, 1*keV);
>> 166
>> 167 AddCrossSectionData("backbone_TMP",
>> 168 particleName,
>> 169 "dna/sigma_elastic_e-_PTB_TMP",
>> 170 "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP",
>> 171 scaleFactor*33./50);
>> 172 SetLowELimit("backbone_TMP", particleName, 10*eV);
>> 173 SetHighELimit("backbone_TMP", particleName, 1*keV);
>> 174 }
388 175
389 // Random sampling of phi << 176 //*******************************************************
390 G4double phi = 2. * CLHEP::pi * G4UniformR << 177 // Load the data
>> 178 //*******************************************************
391 179
392 auto zVers = aDynamicElectron->GetMomentum << 180 LoadCrossSectionData(particle->GetParticleName() );
393 auto xVers = zVers.orthogonal(); <<
394 auto yVers = zVers.cross(xVers); <<
395 181
396 G4double xDir = std::sqrt(1. - cosTheta * << 182 //*******************************************************
397 G4double yDir = xDir; << 183 // Verbose output
398 xDir *= std::cos(phi); << 184 //*******************************************************
399 yDir *= std::sin(phi); <<
400 185
401 // Particle direction after ModelInterface << 186 if (verboseLevel > 2)
402 G4ThreeVector zPrikeVers((xDir * xVers + y << 187 G4cout << "Loaded cross section files for PTB Elastic model" << G4endl;
403 188
404 // Give the new direction << 189 if( verboseLevel>0 )
405 fParticleChangeForGamma->ProposeMomentumDi << 190 {
>> 191 G4cout << "PTB Elastic model is initialized " << G4endl;
>> 192 }
>> 193 }
406 194
407 // Update the energy which does not change << 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408 fParticleChangeForGamma->SetProposedKineti << 196
409 } << 197 void G4DNAPTBElasticModel::ReadDiffCSFile(const G4String& materialName,
>> 198 const G4String& particleName,
>> 199 const G4String& file,
>> 200 const G4double)
>> 201 {
>> 202 // Method to read and save the information contained within the differential cross section files.
>> 203 // This method is not yet standard.
>> 204
>> 205 // get the path of the G4LEDATA data folder
>> 206 char *path = std::getenv("G4LEDATA");
>> 207 // if it is not found then quit and print error message
>> 208 if(!path)
>> 209 {
>> 210 G4Exception("G4DNAPTBElasticModel::ReadAllDiffCSFiles","em0006",
>> 211 FatalException,"G4LEDATA environment variable not set.");
>> 212 return;
>> 213 }
>> 214
>> 215 // build the fullFileName path of the data file
>> 216 std::ostringstream fullFileName;
>> 217 fullFileName << path <<"/"<< file<<".dat";
>> 218
>> 219 // open the data file
>> 220 std::ifstream diffCrossSection (fullFileName.str().c_str());
>> 221 // error if file is not there
>> 222 std::stringstream endPath;
>> 223 if (!diffCrossSection)
>> 224 {
>> 225 endPath << "Missing data file: "<<file;
>> 226 G4Exception("G4DNAPTBElasticModel::Initialise","em0003",
>> 227 FatalException, endPath.str().c_str());
>> 228 }
>> 229
>> 230 tValuesVec[materialName][particleName].push_back(0.);
>> 231
>> 232 G4String line;
>> 233
>> 234 // read the file line by line until we reach the end of file point
>> 235 while(std::getline(diffCrossSection, line))
>> 236 {
>> 237 // check if the line is comment or empty
>> 238 //
>> 239 std::istringstream testIss(line);
>> 240 G4String test;
>> 241 testIss >> test;
>> 242 // check first caracter to determine if following information is data or comments
>> 243 if(test=="#")
>> 244 {
>> 245 // skip the line by beginning a new while loop.
>> 246 continue;
>> 247 }
>> 248 // check if line is empty
>> 249 else if(line.empty())
>> 250 {
>> 251 // skip the line by beginning a new while loop.
>> 252 continue;
>> 253 }
>> 254 //
>> 255 // end of the check
>> 256
>> 257 // transform the line into a iss
>> 258 std::istringstream iss(line);
>> 259
>> 260 // Variables to be filled by the input file
>> 261 double tDummy;
>> 262 double eDummy;
>> 263
>> 264 // fill the variables with the content of the line
>> 265 iss>>tDummy>>eDummy;
>> 266
>> 267 // SI : mandatory Vecm initialization
>> 268
>> 269 // Fill two vectors contained in maps of types:
>> 270 // [materialName][particleName]=vector
>> 271 // [materialName][particleName][T]=vector
>> 272 // to list all the incident energies (tValues) and all the output energies (eValues) within the file
>> 273 //
>> 274 // Check if we already have the current T value in the vector.
>> 275 // If not then add it
>> 276 if (tDummy != tValuesVec[materialName][particleName].back())
>> 277 {
>> 278 // Add the current T value
>> 279 tValuesVec[materialName][particleName].push_back(tDummy);
>> 280
>> 281 // Make it correspond to a default zero E value
>> 282 eValuesVect[materialName][particleName][tDummy].push_back(0.);
>> 283 }
>> 284
>> 285 // Put the differential cross section value of the input file within the diffCrossSectionData map
>> 286 iss>>diffCrossSectionData[materialName][particleName][tDummy][eDummy];
>> 287
>> 288 // If the current E value (eDummy) is different from the one already registered in the eVector then add it to the vector
>> 289 if (eDummy != eValuesVect[materialName][particleName][tDummy].back()) eValuesVect[materialName][particleName][tDummy].push_back(eDummy);
>> 290 }
410 } 291 }
411 292
412 //....oooOO0OOooo........oooOO0OOooo........oo 293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413 294
414 G4double G4DNAPTBElasticModel::Theta(const G4P << 295 G4double G4DNAPTBElasticModel::CrossSectionPerVolume(const G4Material* /*material*/,
415 const std << 296 const G4String& materialName,
416 { << 297 const G4ParticleDefinition* p,
417 G4double theta = 0.; << 298 G4double ekin,
418 G4double valueT1 = 0; << 299 G4double /*emin*/,
419 G4double valueT2 = 0; << 300 G4double /*emax*/)
420 G4double valueE21 = 0; << 301 {
421 G4double valueE22 = 0; << 302 if (verboseLevel > 3)
422 G4double valueE12 = 0; << 303 G4cout << "Calling CrossSectionPerVolume() of G4DNAPTBElasticModel" << G4endl;
423 G4double valueE11 = 0; << 304
424 G4double xs11 = 0; << 305 // Get the name of the current particle
425 G4double xs12 = 0; << 306 const G4String& particleName = p->GetParticleName();
426 G4double xs21 = 0; << 307
427 G4double xs22 = 0; << 308 // set killBelowEnergy value for current material
428 if (p == G4Electron::ElectronDefinition()) { << 309 fKillBelowEnergy = GetLowELimit(materialName, particleName);
429 auto t2 = << 310
430 std::upper_bound(tValuesVec[materialID][ << 311 // initialise the return value (cross section) to zero
431 auto t1 = t2 - 1; << 312 G4double sigma(0);
>> 313
>> 314 // check if we are below the high energy limit
>> 315 if (ekin < GetHighELimit(materialName, particleName) )
>> 316 {
>> 317 // This is used to kill the particle if its kinetic energy is below fKillBelowEnergy.
>> 318 // If the energy is lower then we return a maximum cross section and thus the SampleSecondaries method will be called for sure.
>> 319 // SampleSecondaries will remove the particle from the simulation.
>> 320 //
>> 321 //SI : XS must not be zero otherwise sampling of secondaries method ignored
>> 322 if (ekin < fKillBelowEnergy) return DBL_MAX;
432 323
433 auto e12 = std::upper_bound(eValuesVect[ma << 324 // Get the tables with the cross section data
434 eValuesVect[ma << 325 TableMapData* tableData = GetTableData();
435 auto e11 = e12 - 1; <<
436 326
437 auto e22 = std::upper_bound(eValuesVect[ma << 327 // Retrieve the cross section value
438 eValuesVect[ma << 328 sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
439 auto e21 = e22 - 1; << 329 }
>> 330
>> 331 if (verboseLevel > 2)
>> 332 {
>> 333 G4cout << "__________________________________" << G4endl;
>> 334 G4cout << "°°° G4DNAPTBElasticModel - XS INFO START" << G4endl;
>> 335 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
>> 336 G4cout << "°°° Cross section per molecule (cm^2)=" << sigma/cm/cm << G4endl;
>> 337 G4cout << "°°° G4DNAPTBElasticModel - XS INFO END" << G4endl;
>> 338 }
440 339
441 valueT1 = *t1; << 340 // Return the cross section
442 valueT2 = *t2; << 341 return sigma;
443 valueE21 = *e21; << 342 }
444 valueE22 = *e22; << 343
445 valueE12 = *e12; << 344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
446 valueE11 = *e11; <<
447 345
448 xs11 = diffCrossSectionData[materialID][p] << 346 void G4DNAPTBElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
449 xs12 = diffCrossSectionData[materialID][p] << 347 const G4MaterialCutsCouple* /*couple*/,
450 xs21 = diffCrossSectionData[materialID][p] << 348 const G4String& materialName,
451 xs22 = diffCrossSectionData[materialID][p] << 349 const G4DynamicParticle* aDynamicElectron,
452 } << 350 G4ParticleChangeForGamma* particleChangeForGamma,
>> 351 G4double /*tmin*/,
>> 352 G4double /*tmax*/)
>> 353 {
>> 354 if (verboseLevel > 3)
>> 355 G4cout << "Calling SampleSecondaries() of G4DNAPTBElasticModel" << G4endl;
>> 356
>> 357 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
>> 358
>> 359 const G4String& particleName = aDynamicElectron->GetParticleDefinition()->GetParticleName();
>> 360
>> 361 // set killBelowEnergy value for material
>> 362 fKillBelowEnergy = GetLowELimit(materialName, particleName);
>> 363
>> 364 // If the particle (electron here) energy is below the kill limit then we remove it from the simulation
>> 365 if (electronEnergy0 < fKillBelowEnergy)
>> 366 {
>> 367 particleChangeForGamma->SetProposedKineticEnergy(0.);
>> 368 particleChangeForGamma->ProposeTrackStatus(fStopAndKill);
>> 369 particleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
>> 370 }
>> 371 // If we are above the kill limite and below the high limit then we proceed
>> 372 else if (electronEnergy0>= fKillBelowEnergy && electronEnergy0 < GetHighELimit(materialName, particleName) )
>> 373 {
>> 374 // Random sampling of the cosTheta
>> 375 G4double cosTheta = RandomizeCosTheta(electronEnergy0, materialName);
>> 376
>> 377 // Random sampling of phi
>> 378 G4double phi = 2. * pi * G4UniformRand();
>> 379
>> 380 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
>> 381 G4ThreeVector xVers = zVers.orthogonal();
>> 382 G4ThreeVector yVers = zVers.cross(xVers);
>> 383
>> 384 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
>> 385 G4double yDir = xDir;
>> 386 xDir *= std::cos(phi);
>> 387 yDir *= std::sin(phi);
453 388
454 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x << 389 // Particle direction after ModelInterface
455 return (0.); << 390 G4ThreeVector zPrikeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
456 } <<
457 391
458 theta = QuadInterpolator(valueE11, valueE12, << 392 // Give the new direction
459 valueT2, k, integrD << 393 particleChangeForGamma->ProposeMomentumDirection(zPrikeVers.unit()) ;
460 394
461 return theta; << 395 // Update the energy which does not change here
>> 396 particleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
>> 397 }
462 } 398 }
463 399
464 //....oooOO0OOooo........oooOO0OOooo........oo 400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 401
466 G4double G4DNAPTBElasticModel::LinLogInterpola << 402 G4double G4DNAPTBElasticModel::Theta
467 << 403 (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff, const G4String& materialName)
468 { << 404 {
469 G4double d1 = std::log(xs1); << 405 G4double theta = 0.;
470 G4double d2 = std::log(xs2); << 406 G4double valueT1 = 0;
471 G4double value = std::exp(d1 + (d2 - d1) * ( << 407 G4double valueT2 = 0;
472 return value; << 408 G4double valueE21 = 0;
>> 409 G4double valueE22 = 0;
>> 410 G4double valueE12 = 0;
>> 411 G4double valueE11 = 0;
>> 412 G4double xs11 = 0;
>> 413 G4double xs12 = 0;
>> 414 G4double xs21 = 0;
>> 415 G4double xs22 = 0;
>> 416 G4String particleName = particleDefinition->GetParticleName();
>> 417
>> 418 if (particleDefinition == G4Electron::ElectronDefinition())
>> 419 {
>> 420 std::vector<double>::iterator t2 = std::upper_bound(tValuesVec[materialName][particleName].begin(),tValuesVec[materialName][particleName].end(), k);
>> 421 std::vector<double>::iterator t1 = t2-1;
>> 422
>> 423 std::vector<double>::iterator e12 = std::upper_bound(eValuesVect[materialName][particleName][(*t1)].begin(),eValuesVect[materialName][particleName][(*t1)].end(), integrDiff);
>> 424 std::vector<double>::iterator e11 = e12-1;
>> 425
>> 426 std::vector<double>::iterator e22 = std::upper_bound(eValuesVect[materialName][particleName][(*t2)].begin(),eValuesVect[materialName][particleName][(*t2)].end(), integrDiff);
>> 427 std::vector<double>::iterator e21 = e22-1;
>> 428
>> 429 valueT1 =*t1;
>> 430 valueT2 =*t2;
>> 431 valueE21 =*e21;
>> 432 valueE22 =*e22;
>> 433 valueE12 =*e12;
>> 434 valueE11 =*e11;
>> 435
>> 436 xs11 = diffCrossSectionData[materialName][particleName][valueT1][valueE11];
>> 437 xs12 = diffCrossSectionData[materialName][particleName][valueT1][valueE12];
>> 438 xs21 = diffCrossSectionData[materialName][particleName][valueT2][valueE21];
>> 439 xs22 = diffCrossSectionData[materialName][particleName][valueT2][valueE22];
>> 440 }
>> 441
>> 442 if (xs11==0 && xs12==0 && xs21==0 && xs22==0) return (0.);
>> 443
>> 444 theta = QuadInterpolator ( valueE11, valueE12,
>> 445 valueE21, valueE22,
>> 446 xs11, xs12,
>> 447 xs21, xs22,
>> 448 valueT1, valueT2,
>> 449 k, integrDiff );
>> 450
>> 451 return theta;
473 } 452 }
474 453
475 //....oooOO0OOooo........oooOO0OOooo........oo 454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476 455
477 G4double G4DNAPTBElasticModel::LinLinInterpola << 456 G4double G4DNAPTBElasticModel::LinLogInterpolate(G4double e1,
>> 457 G4double e2,
>> 458 G4double e,
>> 459 G4double xs1,
478 460 G4double xs2)
479 { 461 {
480 G4double d1 = xs1; << 462 G4double d1 = std::log(xs1);
481 G4double d2 = xs2; << 463 G4double d2 = std::log(xs2);
482 G4double value = (d1 + (d2 - d1) * (e - e1) << 464 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
483 return value; << 465 return value;
484 } 466 }
485 467
486 //....oooOO0OOooo........oooOO0OOooo........oo 468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
487 469
488 G4double G4DNAPTBElasticModel::LogLogInterpola << 470 G4double G4DNAPTBElasticModel::LinLinInterpolate(G4double e1,
>> 471 G4double e2,
>> 472 G4double e,
>> 473 G4double xs1,
489 474 G4double xs2)
490 { 475 {
491 G4double a = (std::log10(xs2) - std::log10(x << 476 G4double d1 = xs1;
492 G4double b = std::log10(xs2) - a * std::log1 << 477 G4double d2 = xs2;
493 G4double sigma = a * std::log10(e) + b; << 478 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
494 G4double value = (std::pow(10., sigma)); << 479 return value;
495 return value; <<
496 } 480 }
497 481
498 //....oooOO0OOooo........oooOO0OOooo........oo 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
499 483
500 G4double G4DNAPTBElasticModel::QuadInterpolato << 484 G4double G4DNAPTBElasticModel::LogLogInterpolate(G4double e1,
501 << 485 G4double e2,
502 << 486 G4double e,
503 << 487 G4double xs1,
>> 488 G4double xs2)
504 { 489 {
505 // Log-Log << 490 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
506 /* << 491 G4double b = std::log10(xs2) - a*std::log10(e2);
507 G4double interpolatedvalue1 = LogLogInterpo << 492 G4double sigma = a*std::log10(e) + b;
508 G4double interpolatedvalue2 = LogLogInterpo << 493 G4double value = (std::pow(10.,sigma));
509 G4double value = LogLogInterpolate(t1, t2, << 494 return value;
510 << 495 }
511 496
512 // Lin-Log << 497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
513 G4double interpolatedvalue1 = LinLogInterpo <<
514 G4double interpolatedvalue2 = LinLogInterpo <<
515 G4double value = LinLogInterpolate(t1, t2, <<
516 */ <<
517 498
518 // Lin-Lin << 499 G4double G4DNAPTBElasticModel::QuadInterpolator(G4double e11, G4double e12,
519 G4double interpolatedvalue1 = LinLinInterpol << 500 G4double e21, G4double e22,
520 G4double interpolatedvalue2 = LinLinInterpol << 501 G4double xs11, G4double xs12,
521 G4double value = LinLinInterpolate(t1, t2, t << 502 G4double xs21, G4double xs22,
>> 503 G4double t1, G4double t2,
>> 504 G4double t, G4double e)
>> 505 {
>> 506 // Log-Log
>> 507 /*
>> 508 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
>> 509 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
>> 510 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
>> 511
>> 512
>> 513 // Lin-Log
>> 514 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
>> 515 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
>> 516 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
>> 517 */
>> 518
>> 519 // Lin-Lin
>> 520 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
>> 521 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
>> 522 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
522 523
523 return value; << 524 return value;
524 } 525 }
525 526
526 //....oooOO0OOooo........oooOO0OOooo........oo 527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
527 528
528 G4double G4DNAPTBElasticModel::RandomizeCosThe << 529 G4double G4DNAPTBElasticModel::RandomizeCosTheta(G4double k, const G4String& materialName)
529 { 530 {
530 G4double integrdiff = 0; << 531 G4double integrdiff=0;
531 G4double uniformRand = G4UniformRand(); << 532 G4double uniformRand=G4UniformRand();
532 integrdiff = uniformRand; << 533 integrdiff = uniformRand;
533 534
534 G4double theta = 0.; << 535 G4double theta=0.;
535 G4double cosTheta = 0.; << 536 G4double cosTheta=0.;
536 theta = Theta(G4Electron::ElectronDefinition << 537 theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff, materialName);
537 538
538 cosTheta = std::cos(theta * CLHEP::pi / 180) << 539 cosTheta= std::cos(theta*pi/180);
539 540
540 return cosTheta; << 541 return cosTheta;
541 } 542 }
>> 543
>> 544
>> 545
>> 546
542 547