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