| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // This class is used to support PTB models that come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
29 //
30
31 #include "G4VDNAModel.hh"
32
33 #include "G4ParticleTable.hh"
34 #include "G4SystemOfUnits.hh"
35
36 G4VDNAModel::G4VDNAModel(const G4String& nam, const G4String& applyToMaterial)
37 : G4VEmModel(nam), fStringOfMaterials(applyToMaterial)
38 {}
39
40 G4VDNAModel::~G4VDNAModel() = default;
41
42 void G4VDNAModel::AddCrossSectionData(const std::size_t& materialID,
43 const G4ParticleDefinition* particleName,
44 const G4String& fileCS, const G4String& fileDiffCS,
45 const G4double& scaleFactor)
46 {
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
54 void G4VDNAModel::AddCrossSectionData(const std::size_t& materialID,
55 const G4ParticleDefinition* particleName,
56 const G4String& fileCS, const G4double& scaleFactor)
57 {
58 fModelMaterials.push_back(materialID);
59 fModelParticles.push_back(particleName);
60 fModelCSFiles.push_back(fileCS);
61 fModelScaleFactors.push_back(scaleFactor);
62 }
63
64 void G4VDNAModel::LoadCrossSectionData(const G4ParticleDefinition* particleName)
65 {
66 G4String fileElectron, fileDiffElectron = "";
67 G4String materialName, modelParticleName;
68 G4double scaleFactor;
69 std::size_t materialID;
70
71 const G4ParticleDefinition* pParticle;
72
73 // construct applyToMatVect with materials specified by the user
74 std::vector<G4String> applyToMatVect = BuildApplyToMatVect(fStringOfMaterials);
75
76 // iterate on each material contained into the fStringOfMaterials variable (through
77 // applyToMatVect)
78 for (const auto & i : applyToMatVect) {
79 auto pMat = G4Material::GetMaterial(i, false);
80 if (i != "all" && pMat == nullptr) {
81 continue;
82 }
83
84 // We have selected a material coming from applyToMatVect
85 // We try to find if this material correspond to a model registered material
86 // If it is, then isMatFound becomes true
87 G4bool isMatFound = false;
88
89 // We iterate on each model registered materials to load the CS data
90 // We have to do a for loop because of the "all" option
91 // applyToMatVect[i] == "all" implies applyToMatVect.size()=1 and we want to iterate on all
92 // registered materials
93 for (std::size_t j = 0; j < fModelMaterials.size(); ++j) {
94 if (i == "all" || pMat->GetIndex() == fModelMaterials[j]) {
95 isMatFound = true;
96 materialID = fModelMaterials[j];
97 pParticle = fModelParticles[j];
98 fileElectron = fModelCSFiles[j];
99 if (!fModelDiffCSFiles.empty()) fileDiffElectron = fModelDiffCSFiles[j];
100 scaleFactor = fModelScaleFactors[j];
101
102 ReadAndSaveCSFile(materialID, pParticle, fileElectron, scaleFactor);
103
104 if (!fileDiffElectron.empty())
105 ReadDiffCSFile(materialID, pParticle, fileDiffElectron, scaleFactor);
106 }
107 }
108
109 // check if we found a correspondance, if not: fatal error
110 if (!isMatFound) {
111 std::ostringstream oss;
112 oss << i
113 << " material was not found. It means the material specified in the UserPhysicsList is "
114 "not a model material for ";
115 oss << particleName;
116 G4Exception("G4VDNAModel::LoadCrossSectionData", "em0003", FatalException, oss.str().c_str());
117 return;
118 }
119 }
120 }
121
122 void G4VDNAModel::ReadDiffCSFile(const std::size_t&, const G4ParticleDefinition*, const G4String&,
123 const G4double&)
124 {
125 G4String text(
126 "ReadDiffCSFile must be implemented in the model class using a differential cross section data "
127 "file");
128
129 G4Exception("G4VDNAModel::ReadDiffCSFile", "em0003", FatalException, text);
130 }
131
132 void G4VDNAModel::EnableForMaterialAndParticle(const std::size_t& materialID,
133 const G4ParticleDefinition* p)
134 {
135 fData[materialID][p] = nullptr;
136 }
137
138 std::vector<G4String> G4VDNAModel::BuildApplyToMatVect(const G4String& materials)
139 {
140 // output material vector
141 std::vector<G4String> materialVect;
142
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) {
145 // we add the material to the output vector
146 materialVect.push_back(materials);
147 }
148 // if we have several materials listed in the string then we must retrieve them
149 else {
150 G4String materialsNonIdentified = materials;
151
152 while (materialsNonIdentified.find_first_of('/') != std::string::npos) {
153 // we select the first material and stop at the "/" caracter
154 const G4String& mat = materialsNonIdentified.substr(0, materialsNonIdentified.find_first_of('/'));
155 materialVect.push_back(mat);
156
157 // we remove the previous material from the materialsNonIdentified string
158 materialsNonIdentified = materialsNonIdentified.substr(
159 materialsNonIdentified.find_first_of('/') + 1,
160 materialsNonIdentified.size() - materialsNonIdentified.find_first_of('/'));
161 }
162
163 // we don't find "/" anymore, it means we only have one material string left
164 // we get it
165 materialVect.push_back(std::move(materialsNonIdentified));
166 }
167
168 return materialVect;
169 }
170
171 void G4VDNAModel::ReadAndSaveCSFile(const std::size_t& materialID, const G4ParticleDefinition* p,
172 const G4String& file, const G4double& scaleFactor)
173 {
174 fData[materialID][p] =
175 std::make_unique<G4DNACrossSectionDataSet>(new G4LogLogInterpolation, eV, scaleFactor);
176 fData[materialID][p]->LoadData(file);
177 }
178
179 G4int G4VDNAModel::RandomSelectShell(const G4double& k, const G4ParticleDefinition* particle,
180 const std::size_t& materialID)
181 {
182 G4int level = 0;
183
184 auto pos = fData[materialID].find(particle);
185
186 if (pos != fData[materialID].end()) {
187 G4DNACrossSectionDataSet* table = pos->second.get();
188
189 if (table != nullptr) {
190 auto valuesBuffer = new G4double[table->NumberOfComponents()];
191 auto n = (G4int)table->NumberOfComponents();
192 G4int i(n);
193 G4double value = 0.;
194
195 while (i > 0) {
196 --i;
197 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
198 value += valuesBuffer[i];
199 }
200
201 value *= G4UniformRand();
202
203 i = n;
204
205 while (i > 0) {
206 --i;
207
208 if (valuesBuffer[i] > value) {
209 delete[] valuesBuffer;
210 return i;
211 }
212 value -= valuesBuffer[i];
213 }
214
215 delete[] valuesBuffer;
216 }
217 }
218 else {
219 G4cout << "particle : " << particle->GetParticleName()
220 << " Materials : " << (*G4Material::GetMaterialTable())[materialID]->GetName() << " "
221 << this->GetName() << G4endl;
222 G4Exception("G4VDNAModel::RandomSelectShell", "em0002", FatalException,
223 "Model not applicable to particle type : ");
224 }
225 return level;
226 }
227
228 G4bool G4VDNAModel::IsMaterialDefine(const std::size_t& materialID)
229 {
230 // Check if the given material is defined in the simulation
231
232 G4bool exist(false);
233
234 G4double matTableSize = G4Material::GetMaterialTable()->size();
235
236 for (int i = 0; i < matTableSize; i++) {
237 if (materialID == G4Material::GetMaterialTable()->at(i)->GetIndex()) {
238 exist = true;
239 return exist;
240 }
241 }
242
243 G4Exception("G4VDNAModel::IsMaterialDefine", "em0003", FatalException,
244 "Materials are not defined!!");
245 return exist;
246 }
247
248 G4bool G4VDNAModel::IsMaterialExistingInModel(const std::size_t& materialID)
249 {
250 // Check if the given material is defined in the current model class
251
252 for (const auto& it : fModelMaterials) {
253 if (it == materialID) {
254 return true;
255 }
256 }
257 return false;
258 }
259
260 G4bool G4VDNAModel::IsParticleExistingInModelForMaterial(const G4ParticleDefinition* particleName,
261 const std::size_t& materialID)
262 {
263 // To check two things:
264 // 1- is the material existing in model ?
265 // 2- if yes, is the particle defined for that material ?
266
267 if (IsMaterialExistingInModel(materialID)) {
268 for (const auto& it : fModelParticles) {
269 if (it == particleName) {
270 return true;
271 }
272 }
273 }
274 return false;
275 }
276