| 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 // This class is used to support PTB models th 27 // This class is used to support PTB models that 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 "G4VDNAModel.hh" 31 #include "G4VDNAModel.hh"
32 32
33 #include "G4ParticleTable.hh" 33 #include "G4ParticleTable.hh"
34 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh"
35 35
36 G4VDNAModel::G4VDNAModel(const G4String& nam, 36 G4VDNAModel::G4VDNAModel(const G4String& nam, const G4String& applyToMaterial)
37 : G4VEmModel(nam), fStringOfMaterials(applyT 37 : G4VEmModel(nam), fStringOfMaterials(applyToMaterial)
38 {} 38 {}
39 39
40 G4VDNAModel::~G4VDNAModel() = default; 40 G4VDNAModel::~G4VDNAModel() = default;
41 41
42 void G4VDNAModel::AddCrossSectionData(const st 42 void G4VDNAModel::AddCrossSectionData(const std::size_t& materialID,
43 const G4 43 const G4ParticleDefinition* particleName,
44 const G4 44 const G4String& fileCS, const G4String& fileDiffCS,
45 const G4 45 const G4double& scaleFactor)
46 { 46 {
47 fModelMaterials.push_back(materialID); 47 fModelMaterials.push_back(materialID);
48 fModelParticles.push_back(particleName); 48 fModelParticles.push_back(particleName);
49 fModelCSFiles.push_back(fileCS); 49 fModelCSFiles.push_back(fileCS);
50 fModelDiffCSFiles.push_back(fileDiffCS); 50 fModelDiffCSFiles.push_back(fileDiffCS);
51 fModelScaleFactors.push_back(scaleFactor); 51 fModelScaleFactors.push_back(scaleFactor);
52 } 52 }
53 53
54 void G4VDNAModel::AddCrossSectionData(const st 54 void G4VDNAModel::AddCrossSectionData(const std::size_t& materialID,
55 const G4 55 const G4ParticleDefinition* particleName,
56 const G4 56 const G4String& fileCS, const G4double& scaleFactor)
57 { 57 {
58 fModelMaterials.push_back(materialID); 58 fModelMaterials.push_back(materialID);
59 fModelParticles.push_back(particleName); 59 fModelParticles.push_back(particleName);
60 fModelCSFiles.push_back(fileCS); 60 fModelCSFiles.push_back(fileCS);
61 fModelScaleFactors.push_back(scaleFactor); 61 fModelScaleFactors.push_back(scaleFactor);
62 } 62 }
63 63
64 void G4VDNAModel::LoadCrossSectionData(const G 64 void G4VDNAModel::LoadCrossSectionData(const G4ParticleDefinition* particleName)
65 { 65 {
66 G4String fileElectron, fileDiffElectron = "" 66 G4String fileElectron, fileDiffElectron = "";
67 G4String materialName, modelParticleName; 67 G4String materialName, modelParticleName;
68 G4double scaleFactor; 68 G4double scaleFactor;
69 std::size_t materialID; 69 std::size_t materialID;
70 70
71 const G4ParticleDefinition* pParticle; 71 const G4ParticleDefinition* pParticle;
72 72
73 // construct applyToMatVect with materials s 73 // construct applyToMatVect with materials specified by the user
74 std::vector<G4String> applyToMatVect = Build 74 std::vector<G4String> applyToMatVect = BuildApplyToMatVect(fStringOfMaterials);
75 75
76 // iterate on each material contained into t 76 // iterate on each material contained into the fStringOfMaterials variable (through
77 // applyToMatVect) 77 // applyToMatVect)
78 for (const auto & i : applyToMatVect) { 78 for (const auto & i : applyToMatVect) {
79 auto pMat = G4Material::GetMaterial(i, fal 79 auto pMat = G4Material::GetMaterial(i, false);
80 if (i != "all" && pMat == nullptr) { 80 if (i != "all" && pMat == nullptr) {
81 continue; 81 continue;
82 } 82 }
83 83
84 // We have selected a material coming from 84 // We have selected a material coming from applyToMatVect
85 // We try to find if this material corresp 85 // We try to find if this material correspond to a model registered material
86 // If it is, then isMatFound becomes true 86 // If it is, then isMatFound becomes true
87 G4bool isMatFound = false; 87 G4bool isMatFound = false;
88 88
89 // We iterate on each model registered mat 89 // We iterate on each model registered materials to load the CS data
90 // We have to do a for loop because of the 90 // We have to do a for loop because of the "all" option
91 // applyToMatVect[i] == "all" implies appl 91 // applyToMatVect[i] == "all" implies applyToMatVect.size()=1 and we want to iterate on all
92 // registered materials 92 // registered materials
93 for (std::size_t j = 0; j < fModelMaterial 93 for (std::size_t j = 0; j < fModelMaterials.size(); ++j) {
94 if (i == "all" || pMat->GetIndex() == fM 94 if (i == "all" || pMat->GetIndex() == fModelMaterials[j]) {
95 isMatFound = true; 95 isMatFound = true;
96 materialID = fModelMaterials[j]; 96 materialID = fModelMaterials[j];
97 pParticle = fModelParticles[j]; 97 pParticle = fModelParticles[j];
98 fileElectron = fModelCSFiles[j]; 98 fileElectron = fModelCSFiles[j];
99 if (!fModelDiffCSFiles.empty()) fileDi 99 if (!fModelDiffCSFiles.empty()) fileDiffElectron = fModelDiffCSFiles[j];
100 scaleFactor = fModelScaleFactors[j]; 100 scaleFactor = fModelScaleFactors[j];
101 101
102 ReadAndSaveCSFile(materialID, pParticl 102 ReadAndSaveCSFile(materialID, pParticle, fileElectron, scaleFactor);
103 103
104 if (!fileDiffElectron.empty()) 104 if (!fileDiffElectron.empty())
105 ReadDiffCSFile(materialID, pParticle 105 ReadDiffCSFile(materialID, pParticle, fileDiffElectron, scaleFactor);
106 } 106 }
107 } 107 }
108 108
109 // check if we found a correspondance, if 109 // check if we found a correspondance, if not: fatal error
110 if (!isMatFound) { 110 if (!isMatFound) {
111 std::ostringstream oss; 111 std::ostringstream oss;
112 oss << i 112 oss << i
113 << " material was not found. It mean 113 << " material was not found. It means the material specified in the UserPhysicsList is "
114 "not a model material for "; 114 "not a model material for ";
115 oss << particleName; 115 oss << particleName;
116 G4Exception("G4VDNAModel::LoadCrossSecti 116 G4Exception("G4VDNAModel::LoadCrossSectionData", "em0003", FatalException, oss.str().c_str());
117 return; 117 return;
118 } 118 }
119 } 119 }
120 } 120 }
121 121
122 void G4VDNAModel::ReadDiffCSFile(const std::si 122 void G4VDNAModel::ReadDiffCSFile(const std::size_t&, const G4ParticleDefinition*, const G4String&,
123 const G4doubl 123 const G4double&)
124 { 124 {
125 G4String text( 125 G4String text(
126 "ReadDiffCSFile must be implemented in the 126 "ReadDiffCSFile must be implemented in the model class using a differential cross section data "
127 "file"); 127 "file");
128 128
129 G4Exception("G4VDNAModel::ReadDiffCSFile", " 129 G4Exception("G4VDNAModel::ReadDiffCSFile", "em0003", FatalException, text);
130 } 130 }
131 131
132 void G4VDNAModel::EnableForMaterialAndParticle 132 void G4VDNAModel::EnableForMaterialAndParticle(const std::size_t& materialID,
133 133 const G4ParticleDefinition* p)
134 { 134 {
135 fData[materialID][p] = nullptr; 135 fData[materialID][p] = nullptr;
136 } 136 }
137 137
138 std::vector<G4String> G4VDNAModel::BuildApplyT 138 std::vector<G4String> G4VDNAModel::BuildApplyToMatVect(const G4String& materials)
139 { 139 {
140 // output material vector 140 // output material vector
141 std::vector<G4String> materialVect; 141 std::vector<G4String> materialVect;
142 142
143 // if we don't find any "/" then it means we 143 // if we don't find any "/" then it means we only have one "material" (could be the "all" option)
144 if (materials.find('/') == std::string::npos 144 if (materials.find('/') == std::string::npos) {
145 // we add the material to the output vecto 145 // we add the material to the output vector
146 materialVect.push_back(materials); 146 materialVect.push_back(materials);
147 } 147 }
148 // if we have several materials listed in th 148 // if we have several materials listed in the string then we must retrieve them
149 else { 149 else {
150 G4String materialsNonIdentified = material 150 G4String materialsNonIdentified = materials;
151 151
152 while (materialsNonIdentified.find_first_o 152 while (materialsNonIdentified.find_first_of('/') != std::string::npos) {
153 // we select the first material and stop 153 // we select the first material and stop at the "/" caracter
154 const G4String& mat = materialsNonIdenti << 154 G4String mat = materialsNonIdentified.substr(0, materialsNonIdentified.find_first_of('/'));
155 materialVect.push_back(mat); 155 materialVect.push_back(mat);
156 156
157 // we remove the previous material from 157 // we remove the previous material from the materialsNonIdentified string
158 materialsNonIdentified = materialsNonIde 158 materialsNonIdentified = materialsNonIdentified.substr(
159 materialsNonIdentified.find_first_of(' 159 materialsNonIdentified.find_first_of('/') + 1,
160 materialsNonIdentified.size() - materi 160 materialsNonIdentified.size() - materialsNonIdentified.find_first_of('/'));
161 } 161 }
162 162
163 // we don't find "/" anymore, it means we 163 // we don't find "/" anymore, it means we only have one material string left
164 // we get it 164 // we get it
165 materialVect.push_back(std::move(materials << 165 materialVect.push_back(materialsNonIdentified);
166 } 166 }
167 167
168 return materialVect; 168 return materialVect;
169 } 169 }
170 170
171 void G4VDNAModel::ReadAndSaveCSFile(const std: 171 void G4VDNAModel::ReadAndSaveCSFile(const std::size_t& materialID, const G4ParticleDefinition* p,
172 const G4St 172 const G4String& file, const G4double& scaleFactor)
173 { 173 {
174 fData[materialID][p] = 174 fData[materialID][p] =
175 std::make_unique<G4DNACrossSectionDataSet> 175 std::make_unique<G4DNACrossSectionDataSet>(new G4LogLogInterpolation, eV, scaleFactor);
176 fData[materialID][p]->LoadData(file); 176 fData[materialID][p]->LoadData(file);
177 } 177 }
178 178
179 G4int G4VDNAModel::RandomSelectShell(const G4d 179 G4int G4VDNAModel::RandomSelectShell(const G4double& k, const G4ParticleDefinition* particle,
180 const std 180 const std::size_t& materialID)
181 { 181 {
182 G4int level = 0; 182 G4int level = 0;
183 183
184 auto pos = fData[materialID].find(particle); 184 auto pos = fData[materialID].find(particle);
185 185
186 if (pos != fData[materialID].end()) { 186 if (pos != fData[materialID].end()) {
187 G4DNACrossSectionDataSet* table = pos->sec 187 G4DNACrossSectionDataSet* table = pos->second.get();
188 188
189 if (table != nullptr) { 189 if (table != nullptr) {
190 auto valuesBuffer = new G4double[table-> 190 auto valuesBuffer = new G4double[table->NumberOfComponents()];
191 auto n = (G4int)table->NumberOfComponent 191 auto n = (G4int)table->NumberOfComponents();
192 G4int i(n); 192 G4int i(n);
193 G4double value = 0.; 193 G4double value = 0.;
194 194
195 while (i > 0) { 195 while (i > 0) {
196 --i; 196 --i;
197 valuesBuffer[i] = table->GetComponent( 197 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
198 value += valuesBuffer[i]; 198 value += valuesBuffer[i];
199 } 199 }
200 200
201 value *= G4UniformRand(); 201 value *= G4UniformRand();
202 202
203 i = n; 203 i = n;
204 204
205 while (i > 0) { 205 while (i > 0) {
206 --i; 206 --i;
207 207
208 if (valuesBuffer[i] > value) { 208 if (valuesBuffer[i] > value) {
209 delete[] valuesBuffer; 209 delete[] valuesBuffer;
210 return i; 210 return i;
211 } 211 }
212 value -= valuesBuffer[i]; 212 value -= valuesBuffer[i];
213 } 213 }
214 214
215 delete[] valuesBuffer; 215 delete[] valuesBuffer;
216 } 216 }
217 } 217 }
218 else { 218 else {
219 G4cout << "particle : " << particle->GetPa 219 G4cout << "particle : " << particle->GetParticleName()
220 << " Materials : " << (*G4Material: 220 << " Materials : " << (*G4Material::GetMaterialTable())[materialID]->GetName() << " "
221 << this->GetName() << G4endl; 221 << this->GetName() << G4endl;
222 G4Exception("G4VDNAModel::RandomSelectShel 222 G4Exception("G4VDNAModel::RandomSelectShell", "em0002", FatalException,
223 "Model not applicable to parti 223 "Model not applicable to particle type : ");
224 } 224 }
225 return level; 225 return level;
226 } 226 }
227 227
228 G4bool G4VDNAModel::IsMaterialDefine(const std 228 G4bool G4VDNAModel::IsMaterialDefine(const std::size_t& materialID)
229 { 229 {
230 // Check if the given material is defined in 230 // Check if the given material is defined in the simulation
231 231
232 G4bool exist(false); 232 G4bool exist(false);
233 233
234 G4double matTableSize = G4Material::GetMater 234 G4double matTableSize = G4Material::GetMaterialTable()->size();
235 235
236 for (int i = 0; i < matTableSize; i++) { 236 for (int i = 0; i < matTableSize; i++) {
237 if (materialID == G4Material::GetMaterialT 237 if (materialID == G4Material::GetMaterialTable()->at(i)->GetIndex()) {
238 exist = true; 238 exist = true;
239 return exist; 239 return exist;
240 } 240 }
241 } 241 }
242 242
243 G4Exception("G4VDNAModel::IsMaterialDefine", 243 G4Exception("G4VDNAModel::IsMaterialDefine", "em0003", FatalException,
244 "Materials are not defined!!"); 244 "Materials are not defined!!");
245 return exist; 245 return exist;
246 } 246 }
247 247
248 G4bool G4VDNAModel::IsMaterialExistingInModel( 248 G4bool G4VDNAModel::IsMaterialExistingInModel(const std::size_t& materialID)
249 { 249 {
250 // Check if the given material is defined in 250 // Check if the given material is defined in the current model class
251 251
252 for (const auto& it : fModelMaterials) { 252 for (const auto& it : fModelMaterials) {
253 if (it == materialID) { 253 if (it == materialID) {
254 return true; 254 return true;
255 } 255 }
256 } 256 }
257 return false; 257 return false;
258 } 258 }
259 259
260 G4bool G4VDNAModel::IsParticleExistingInModelF 260 G4bool G4VDNAModel::IsParticleExistingInModelForMaterial(const G4ParticleDefinition* particleName,
261 261 const std::size_t& materialID)
262 { 262 {
263 // To check two things: 263 // To check two things:
264 // 1- is the material existing in model ? 264 // 1- is the material existing in model ?
265 // 2- if yes, is the particle defined for th 265 // 2- if yes, is the particle defined for that material ?
266 266
267 if (IsMaterialExistingInModel(materialID)) { 267 if (IsMaterialExistingInModel(materialID)) {
268 for (const auto& it : fModelParticles) { 268 for (const auto& it : fModelParticles) {
269 if (it == particleName) { 269 if (it == particleName) {
270 return true; 270 return true;
271 } 271 }
272 } 272 }
273 } 273 }
274 return false; 274 return false;
275 } 275 }
276 276