| 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 // 26 //
27 // Contact authors: S. Meylan, C. Villagrasa 27 // Contact authors: S. Meylan, C. Villagrasa
28 // email: sylvain.meylan@symalgo-tech.com, car 28 // email: sylvain.meylan@symalgo-tech.com, carmen.villagrasa@irsn.fr
29 // updated : Hoang Tran : 6/1/2023 clean code 29 // updated : Hoang Tran : 6/1/2023 clean code
30 30
31 #include "G4DNAModelInterface.hh" 31 #include "G4DNAModelInterface.hh"
32 32
33 #include "G4DNAMolecularMaterial.hh" 33 #include "G4DNAMolecularMaterial.hh"
34 #include "G4LossTableManager.hh" 34 #include "G4LossTableManager.hh"
35 #include "G4ParticleChangeForGamma.hh" 35 #include "G4ParticleChangeForGamma.hh"
36 #include "G4VDNAModel.hh" 36 #include "G4VDNAModel.hh"
37 #include "G4VEmModel.hh" 37 #include "G4VEmModel.hh"
38 G4DNAModelInterface::G4DNAModelInterface(const 38 G4DNAModelInterface::G4DNAModelInterface(const G4String& nam) : G4VEmModel(nam), fName(nam) {}
39 39
40 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 41
42 void G4DNAModelInterface::Initialise(const G4P 42 void G4DNAModelInterface::Initialise(const G4ParticleDefinition* particle, const G4DataVector& cuts)
43 { 43 {
44 // Those two statements are necessary to ove 44 // Those two statements are necessary to override the energy limits set in the G4DNAProcesses
45 // (ionisation, elastic, etc...). Indeed, wi 45 // (ionisation, elastic, etc...). Indeed, with the ModelInterface system, the model define
46 // themselves their energy limits per materi 46 // themselves their energy limits per material and particle. Therefore, such a limit should not be
47 // in the G4DNAProcess classes. 47 // in the G4DNAProcess classes.
48 // 48 //
49 49
50 fpG4_WATER = G4Material::GetMaterial("G4_WAT 50 fpG4_WATER = G4Material::GetMaterial("G4_WATER", false);
51 51
52 SetLowEnergyLimit(0.); 52 SetLowEnergyLimit(0.);
53 SetHighEnergyLimit(DBL_MAX); 53 SetHighEnergyLimit(DBL_MAX);
54 54
55 fpParticleChangeForGamma = GetParticleChange 55 fpParticleChangeForGamma = GetParticleChangeForGamma();
56 56
57 // Loop on all the registered models to init 57 // Loop on all the registered models to initialise them
58 for (auto & fRegisteredModel : fRegisteredMo 58 for (auto & fRegisteredModel : fRegisteredModels) {
59 fRegisteredModel->SetParticleChange(fpPart 59 fRegisteredModel->SetParticleChange(fpParticleChangeForGamma);
60 fRegisteredModel->Initialise(particle, cut 60 fRegisteredModel->Initialise(particle, cuts);
61 } 61 }
62 // used to retrieve the model corresponding 62 // used to retrieve the model corresponding to the current material/particle couple
63 BuildMaterialParticleModelTable(particle); 63 BuildMaterialParticleModelTable(particle);
64 64
65 BuildMaterialMolPerVolTable(); 65 BuildMaterialMolPerVolTable();
66 66
67 StreamInfo(G4cout); 67 StreamInfo(G4cout);
68 68
69 } 69 }
70 70
71 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 72
73 G4double G4DNAModelInterface::CrossSectionPerV 73 G4double G4DNAModelInterface::CrossSectionPerVolume(const G4Material* material,
74 74 const G4ParticleDefinition* p, G4double ekin,
75 75 G4double emin, G4double emax)
76 { 76 {
77 // Method to return the crossSection * nbMol 77 // Method to return the crossSection * nbMoleculePerUnitVolume to the process class.
78 // Process class then calculates the path. 78 // Process class then calculates the path.
79 // The cross section is calculated in the re 79 // The cross section is calculated in the registered model(s) and this class just call the method
80 // Two cases are handled here: normal materi 80 // Two cases are handled here: normal material and composite material.
81 // 81 //
82 // Idea: 82 // Idea:
83 // *** Simple material *** 83 // *** Simple material ***
84 // Ask for the cross section of the chosen m 84 // Ask for the cross section of the chosen model.
85 // Multiply it by the number of medium molec 85 // Multiply it by the number of medium molecules per volume unit.
86 // Return the value. 86 // Return the value.
87 // *** Composite material *** 87 // *** Composite material ***
88 // Ask for the cross section of the chosen m 88 // Ask for the cross section of the chosen model for each component.
89 // Apply a factor to each cross section and 89 // Apply a factor to each cross section and sum the results. The factor is the molecule number of
90 // component per composite volume unit. The 90 // component per composite volume unit. The total cross section is returned.
91 91
92 // To reset the sampledMat variable. 92 // To reset the sampledMat variable.
93 // Can be used by user to retrieve current c 93 // Can be used by user to retrieve current component
94 fSampledMat = 0; 94 fSampledMat = 0;
95 95
96 // This is the value to be sum up and to be 96 // This is the value to be sum up and to be returned at then end
97 G4double crossSectionTimesNbMolPerVol(0.); 97 G4double crossSectionTimesNbMolPerVol(0.);
98 98
99 // Reset the map saving the material and the 99 // Reset the map saving the material and the cumulated corresponding cross section
100 // Used in SampleSecondaries if the interact 100 // Used in SampleSecondaries if the interaction is selected for the step and if the material is a
101 // composite 101 // composite
102 fMaterialCS.clear(); 102 fMaterialCS.clear();
103 103
104 // This is the value to be used by SampleSec 104 // This is the value to be used by SampleSecondaries
105 fCSsumTot = 0; 105 fCSsumTot = 0;
106 106
107 // ***************************** 107 // *****************************
108 // Material is not a composite 108 // Material is not a composite
109 // ***************************** 109 // *****************************
110 // 110 //
111 if (material->GetMatComponents().empty()) { 111 if (material->GetMatComponents().empty()) {
112 // Get the material name 112 // Get the material name
113 const size_t & materialID = material->GetI 113 const size_t & materialID = material->GetIndex();
114 114
115 // Use the table to get the model 115 // Use the table to get the model
116 auto model = SelectModel(materialID, p, ek 116 auto model = SelectModel(materialID, p, ekin);
117 117
118 // Get the nunber of molecules per volume 118 // Get the nunber of molecules per volume unit for that material
119 119
120 // Calculate the cross section times the n 120 // Calculate the cross section times the number of molecules
121 if (model != nullptr) { 121 if (model != nullptr) {
122 if (dynamic_cast<G4VDNAModel*>(model) == 122 if (dynamic_cast<G4VDNAModel*>(model) == nullptr) {
123 // water material models only 123 // water material models only
124 crossSectionTimesNbMolPerVol = model-> 124 crossSectionTimesNbMolPerVol = model->CrossSectionPerVolume(material, p, ekin, emin, emax);
125 } 125 }
126 else { 126 else {
127 crossSectionTimesNbMolPerVol = model-> 127 crossSectionTimesNbMolPerVol = model->CrossSectionPerVolume(material, p, ekin, emin, emax);
128 } 128 }
129 } 129 }
130 else // no model was selected, we are out 130 else // no model was selected, we are out of the energy ranges
131 crossSectionTimesNbMolPerVol = 0.; 131 crossSectionTimesNbMolPerVol = 0.;
132 } 132 }
133 133
134 // ******************************** 134 // ********************************
135 // Material is a composite 135 // Material is a composite
136 // ******************************** 136 // ********************************
137 // 137 //
138 else { 138 else {
139 // Copy the map in a local variable 139 // Copy the map in a local variable
140 // Otherwise we get segmentation fault and 140 // Otherwise we get segmentation fault and iterator pointing to nowhere: do not know why...
141 // Maybe MatComponents map is overrided by 141 // Maybe MatComponents map is overrided by something somewhere ?
142 auto componentsMap = material->GetMatCompo 142 auto componentsMap = material->GetMatComponents();
143 143
144 G4cout << material->GetName() << G4endl; 144 G4cout << material->GetName() << G4endl;
145 145
146 // Loop on all the components 146 // Loop on all the components
147 for (const auto& it : componentsMap) { 147 for (const auto& it : componentsMap) {
148 // Get the current component 148 // Get the current component
149 auto component = it.first; 149 auto component = it.first;
150 // Get the current component mass fracti 150 // Get the current component mass fraction
151 // G4double massFraction = it->second; 151 // G4double massFraction = it->second;
152 152
153 // Get the number of component molecules 153 // Get the number of component molecules in a volume unit of composite material
154 G4double nbMoleculeOfComponentInComposit 154 G4double nbMoleculeOfComponentInCompositeMat =
155 GetNumMolPerVolUnitForComponentInCompo 155 GetNumMolPerVolUnitForComponentInComposite(component, material);
156 G4cout << " ==========>component : " << 156 G4cout << " ==========>component : " << component->GetName()
157 << " nbMoleculeOfComponentInCompo 157 << " nbMoleculeOfComponentInCompositeMat: " << nbMoleculeOfComponentInCompositeMat
158 << G4endl; 158 << G4endl;
159 159
160 // Get the current component name 160 // Get the current component name
161 const std::size_t & componentID = compon 161 const std::size_t & componentID = component->GetIndex();
162 162
163 // Retrieve the model corresponding to t 163 // Retrieve the model corresponding to the current component (ie material)
164 auto model = SelectModel(componentID, p, 164 auto model = SelectModel(componentID, p, ekin);
165 165
166 // Add the component part of the cross s 166 // Add the component part of the cross section to the cross section variable.
167 // The component cross section is multip 167 // The component cross section is multiplied by the total molecule number in the composite
168 // scaled by the mass fraction. 168 // scaled by the mass fraction.
169 G4double crossSection; 169 G4double crossSection;
170 if (model != nullptr) { 170 if (model != nullptr) {
171 if (dynamic_cast<G4VDNAModel*>(model) 171 if (dynamic_cast<G4VDNAModel*>(model) == nullptr) {
172 // water models 172 // water models
173 crossSection = 173 crossSection =
174 model->CrossSectionPerVolume(compo 174 model->CrossSectionPerVolume(component, p, ekin, emin, emax)
175 / GetNumMoleculePerVolumeUnitForMa 175 / GetNumMoleculePerVolumeUnitForMaterial(fpG4_WATER);
176 } 176 }
177 else { 177 else {
178 crossSection = model->CrossSectionPe 178 crossSection = model->CrossSectionPerVolume(component, p, ekin, emin, emax)
179 / GetNumMoleculePerVo 179 / GetNumMoleculePerVolumeUnitForMaterial(component);
180 } 180 }
181 crossSectionTimesNbMolPerVol = nbMolec 181 crossSectionTimesNbMolPerVol = nbMoleculeOfComponentInCompositeMat * crossSection;
182 } 182 }
183 else // no model was selected, we are o 183 else // no model was selected, we are out of the energy ranges
184 { 184 {
185 crossSectionTimesNbMolPerVol = 0.; 185 crossSectionTimesNbMolPerVol = 0.;
186 } 186 }
187 187
188 // Save the component name and its calcu 188 // Save the component name and its calculated crossSectionTimesNbMolPerVol
189 // To be used by sampling secondaries if 189 // To be used by sampling secondaries if the interaction is selected for the step
190 fMaterialCS[componentID] = crossSectionT 190 fMaterialCS[componentID] = crossSectionTimesNbMolPerVol;
191 191
192 // Save the component name and its calcu 192 // Save the component name and its calculated crossSectionTimesNbMolPerVol
193 // To be used by sampling secondaries if 193 // To be used by sampling secondaries if the interaction is selected for the step
194 fCSsumTot += crossSectionTimesNbMolPerVo 194 fCSsumTot += crossSectionTimesNbMolPerVol;
195 } 195 }
196 crossSectionTimesNbMolPerVol = fCSsumTot; 196 crossSectionTimesNbMolPerVol = fCSsumTot;
197 } 197 }
198 198
199 // return the cross section times the number 199 // return the cross section times the number of molecules
200 // the path of the interaction will be calcu 200 // the path of the interaction will be calculated using that value
201 return crossSectionTimesNbMolPerVol; 201 return crossSectionTimesNbMolPerVol;
202 } 202 }
203 203
204 //....oooOO0OOooo........oooOO0OOooo........oo 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
205 205
206 void G4DNAModelInterface::SampleSecondaries(st 206 void G4DNAModelInterface::SampleSecondaries(std::vector<G4DynamicParticle*>* fVect,
207 co 207 const G4MaterialCutsCouple* couple,
208 co 208 const G4DynamicParticle* aDynamicParticle,
209 G4 209 G4double tmin, G4double tmax)
210 { 210 {
211 // To call the sampleSecondaries method of t 211 // To call the sampleSecondaries method of the registered model(s)
212 // In the case of composite material, we nee 212 // In the case of composite material, we need to choose a component to call the method from.
213 // To do so we use a random sampling on the 213 // To do so we use a random sampling on the crossSectionTimesNbMolPerVol used in
214 // CrossSectionPerVolume method. If we enter 214 // CrossSectionPerVolume method. If we enter that method it means the corresponding interaction
215 // (and process) has been chosen for the cur 215 // (and process) has been chosen for the current step.
216 216
217 std::size_t materialID; 217 std::size_t materialID;
218 218
219 // ******************************* 219 // *******************************
220 // Material is not a composite 220 // Material is not a composite
221 // ******************************* 221 // *******************************
222 // 222 //
223 if (couple->GetMaterial()->GetMatComponents( 223 if (couple->GetMaterial()->GetMatComponents().empty()) {
224 materialID = couple->GetMaterial()->GetInd 224 materialID = couple->GetMaterial()->GetIndex();
225 } 225 }
226 226
227 // **************************** 227 // ****************************
228 // Material is a composite 228 // Material is a composite
229 // **************************** 229 // ****************************
230 // 230 //
231 else { 231 else {
232 // Material is a composite 232 // Material is a composite
233 // We need to select a component 233 // We need to select a component
234 234
235 // We select a random number between 0 and 235 // We select a random number between 0 and fCSSumTot
236 G4double rand = G4UniformRand() * fCSsumTo 236 G4double rand = G4UniformRand() * fCSsumTot;
237 237
238 G4double cumulCS(0); 238 G4double cumulCS(0);
239 239
240 G4bool result = false; 240 G4bool result = false;
241 241
242 // We loop on each component cumulated cro 242 // We loop on each component cumulated cross section
243 // 243 //
244 // Retrieve the iterators 244 // Retrieve the iterators
245 auto it = fMaterialCS.begin(); 245 auto it = fMaterialCS.begin();
246 auto ite = fMaterialCS.end(); 246 auto ite = fMaterialCS.end();
247 // While this is true we do not have found 247 // While this is true we do not have found our component.
248 while (rand > cumulCS) { 248 while (rand > cumulCS) {
249 // Check if the sampling is ok 249 // Check if the sampling is ok
250 if (it == ite) { 250 if (it == ite) {
251 G4Exception( 251 G4Exception(
252 "G4DNAModelManager::SampleSecondarie 252 "G4DNAModelManager::SampleSecondaries", "em0003", FatalException,
253 "The random component selection has 253 "The random component selection has failed: we ran into the end of the map without "
254 "having a selected component"); 254 "having a selected component");
255 return; // to make some compilers hap 255 return; // to make some compilers happy
256 } 256 }
257 257
258 // Set the cumulated value for the itera 258 // Set the cumulated value for the iteration
259 cumulCS += it->second; 259 cumulCS += it->second;
260 260
261 // Check if we have reach the material t 261 // Check if we have reach the material to be selected
262 // The DBL_MAX is here to take into acco 262 // The DBL_MAX is here to take into account a return DBL_MAX in CSPerVol for the elastic model
263 // to force elastic sampleSecondaries wh 263 // to force elastic sampleSecondaries where the particle can be killed.
264 // Used when paticle energy is lower tha 264 // Used when paticle energy is lower than limit.
265 if (rand < cumulCS || cumulCS >= DBL_MAX 265 if (rand < cumulCS || cumulCS >= DBL_MAX) {
266 // we have our selected material 266 // we have our selected material
267 materialID = it->first; 267 materialID = it->first;
268 result = true; 268 result = true;
269 break; 269 break;
270 } 270 }
271 271
272 // make the iterator move forward 272 // make the iterator move forward
273 ++it; 273 ++it;
274 } 274 }
275 275
276 // Check that we get a result 276 // Check that we get a result
277 if (!result) { 277 if (!result) {
278 // it is possible to end up here if the 278 // it is possible to end up here if the return DBL_MAX of CSPerVol in the elastic model is not
279 // taken into account 279 // taken into account
280 280
281 G4Exception("G4DNAModelManager::SampleSe 281 G4Exception("G4DNAModelManager::SampleSecondaries", "em0005", FatalException,
282 "The random component select 282 "The random component selection has failed: while loop ended without a selected "
283 "component."); 283 "component.");
284 return; // to make some compilers happy 284 return; // to make some compilers happy
285 } 285 }
286 } 286 }
287 287
288 // ************************************** 288 // **************************************
289 // Call the SampleSecondaries method 289 // Call the SampleSecondaries method
290 // ************************************** 290 // **************************************
291 291
292 // Rename material if modified NIST material 292 // Rename material if modified NIST material
293 // This is needed when material is obtained 293 // This is needed when material is obtained from G4MaterialCutsCouple
294 // if (materialName.find("_MODIFIED") != G4 294 // if (materialName.find("_MODIFIED") != G4String::npos) {
295 // materialName = materialName.substr(0, 295 // materialName = materialName.substr(0, materialName.size() - 9);
296 // } 296 // }
297 297
298 fSampledMat = materialID; 298 fSampledMat = materialID;
299 299
300 auto model = SelectModel(materialID, aDynami 300 auto model = SelectModel(materialID, aDynamicParticle->GetParticleDefinition(),
301 aDynamicParticle->G 301 aDynamicParticle->GetKineticEnergy());
302 302
303 model->SampleSecondaries(fVect, couple, aDyn 303 model->SampleSecondaries(fVect, couple, aDynamicParticle, tmin, tmax);
304 } 304 }
305 305
306 void G4DNAModelInterface::RegisterModel(G4VEmM 306 void G4DNAModelInterface::RegisterModel(G4VEmModel* model)
307 { 307 {
308 fRegisteredModels.push_back(model); 308 fRegisteredModels.push_back(model);
309 } 309 }
310 310
311 //....oooOO0OOooo........oooOO0OOooo........oo 311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
312 312
313 void G4DNAModelInterface::BuildMaterialParticl 313 void G4DNAModelInterface::BuildMaterialParticleModelTable(const G4ParticleDefinition* p)
314 { 314 {
315 // Method to build a map: [material][particl 315 // Method to build a map: [material][particle] = Model*.
316 // The map is used to retrieve the correct m 316 // The map is used to retrieve the correct model for the current particle/material couple.
317 317
318 // Loop on all materials registered in the s 318 // Loop on all materials registered in the simulation
319 for (auto it : *G4Material::GetMaterialTable 319 for (auto it : *G4Material::GetMaterialTable()) {
320 // Get the material pointer 320 // Get the material pointer
321 G4Material* mat = it; 321 G4Material* mat = it;
322 // Get the map 322 // Get the map
323 // Check that the material is not a compos 323 // Check that the material is not a composite material
324 auto componentMap = mat->GetMatComponents( 324 auto componentMap = mat->GetMatComponents();
325 if (componentMap.empty()) { 325 if (componentMap.empty()) {
326 // Get the material name 326 // Get the material name
327 const std::size_t & matID = mat->GetInde 327 const std::size_t & matID = mat->GetIndex();
328 InsertModelInTable(matID, p); 328 InsertModelInTable(matID, p);
329 } 329 }
330 // if the material is a composite material 330 // if the material is a composite material then we need to loop on all its components to
331 // register them 331 // register them
332 else { 332 else {
333 // Loop on all the components of the mat 333 // Loop on all the components of the material
334 for (const auto& itComp : componentMap) 334 for (const auto& itComp : componentMap) {
335 G4Material* component = itComp.first; 335 G4Material* component = itComp.first;
336 // Check that the component is not its 336 // Check that the component is not itself a composite
337 if (!component->GetMatComponents().emp 337 if (!component->GetMatComponents().empty()) {
338 std::ostringstream oss; 338 std::ostringstream oss;
339 oss << "Material " << mat->GetName() 339 oss << "Material " << mat->GetName() << " is a composite and its component";
340 oss << " " << component->GetName(); 340 oss << " " << component->GetName();
341 G4Exception("G4DNAModelManager::Buil 341 G4Exception("G4DNAModelManager::BuildMaterialParticleModelTable", "em0007",
342 FatalException, oss.str( 342 FatalException, oss.str().c_str());
343 return; // to make some compilers h 343 return; // to make some compilers happy
344 } 344 }
345 // Get the current component name 345 // Get the current component name
346 const std::size_t & compID = component 346 const std::size_t & compID = component->GetIndex();
347 // If there is a model then insert the 347 // If there is a model then insert the model corresponding to the component in the table
348 // contains a if statement to check we 348 // contains a if statement to check we have not registered the material as a component or a
349 // normal material before. 349 // normal material before.
350 InsertModelInTable(compID, p); 350 InsertModelInTable(compID, p);
351 // move forward the iterator 351 // move forward the iterator
352 } 352 }
353 } 353 }
354 } 354 }
355 } 355 }
356 356
357 //....oooOO0OOooo........oooOO0OOooo........oo 357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
358 358
359 void G4DNAModelInterface::BuildMaterialMolPerV 359 void G4DNAModelInterface::BuildMaterialMolPerVolTable()
360 { 360 {
361 // To be sure the G4DNAMolecularMaterial is 361 // To be sure the G4DNAMolecularMaterial is initialized
362 G4DNAMolecularMaterial::Instance()->Initiali 362 G4DNAMolecularMaterial::Instance()->Initialize();
363 363
364 G4MaterialTable* materialTable = G4Material: 364 G4MaterialTable* materialTable = G4Material::GetMaterialTable();
365 365
366 // Loop on all the materials inside the "mat 366 // Loop on all the materials inside the "materialTable"
367 for (auto currentMaterial : *materialTable) 367 for (auto currentMaterial : *materialTable) {
368 // Current material 368 // Current material
369 // Current material name 369 // Current material name
370 const std::size_t & currentMatID = current 370 const std::size_t & currentMatID = currentMaterial->GetIndex();
371 371
372 // Will the material be used in this inter 372 // Will the material be used in this interface instance ?
373 // Loop on all the materials that can be d 373 // Loop on all the materials that can be dealt with in this class
374 auto it = fMaterialParticleModelTable.begi 374 auto it = fMaterialParticleModelTable.begin();
375 auto ite = fMaterialParticleModelTable.end 375 auto ite = fMaterialParticleModelTable.end();
376 for (; it != ite; it++) { 376 for (; it != ite; it++) {
377 const std::size_t & materialID = it->fir 377 const std::size_t & materialID = it->first;
378 378
379 if (materialID == currentMatID) { 379 if (materialID == currentMatID) {
380 const std::vector<G4double>* numMolPer 380 const std::vector<G4double>* numMolPerVolForMat =
381 G4DNAMolecularMaterial::Instance()-> 381 G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(currentMaterial);
382 fMaterialMolPerVol[materialID] = numMo 382 fMaterialMolPerVol[materialID] = numMolPerVolForMat;
383 } 383 }
384 } 384 }
385 } 385 }
386 } 386 }
387 387
388 //....oooOO0OOooo........oooOO0OOooo........oo 388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389 389
390 void G4DNAModelInterface::InsertModelInTable(c 390 void G4DNAModelInterface::InsertModelInTable(const std::size_t& matID, const G4ParticleDefinition* p)
391 { 391 {
392 // To insert the model(s) in the table Mater 392 // To insert the model(s) in the table Material Particule -> Model(s)
393 393
394 // First, we need to check if the current ma 394 // First, we need to check if the current material has already been inserted in the table.
395 // This is possible because of the composite 395 // This is possible because of the composite material. We could add a component M1 and then try to
396 // add the independant M1 material. This cas 396 // add the independant M1 material. This case must be avoided. Checking if M1 is already in the
397 // table is the way to avoid it. 397 // table is the way to avoid it.
398 // 398 //
399 // Check if the current material and particl 399 // Check if the current material and particle are already in the table.
400 // If they are: do nothing. 400 // If they are: do nothing.
401 // If they are not: add the model(s) 401 // If they are not: add the model(s)
402 // 402 //
403 // Check for the material 403 // Check for the material
404 if (fMaterialParticleModelTable.find(matID) 404 if (fMaterialParticleModelTable.find(matID) == fMaterialParticleModelTable.end()) {
405 // Check for the particle 405 // Check for the particle
406 if (fMaterialParticleModelTable[matID].fin 406 if (fMaterialParticleModelTable[matID].find(p) == fMaterialParticleModelTable[matID].end()) {
407 G4int modelNbForMaterial = 0; 407 G4int modelNbForMaterial = 0;
408 for (const auto& it : fRegisteredModels) 408 for (const auto& it : fRegisteredModels) {
409 auto model = dynamic_cast<G4VDNAModel* 409 auto model = dynamic_cast<G4VDNAModel*>(it);
410 if (model != nullptr) { 410 if (model != nullptr) {
411 if (model->IsParticleExistingInModel 411 if (model->IsParticleExistingInModelForMaterial(p, matID)) {
412 fMaterialParticleModelTable[matID] 412 fMaterialParticleModelTable[matID][p] = it;
413 // and add one to the "there is a 413 // and add one to the "there is a model" material flag
414 ++modelNbForMaterial; 414 ++modelNbForMaterial;
415 } 415 }
416 } 416 }
417 else { 417 else {
418 auto index = fpG4_WATER->GetIndex(); 418 auto index = fpG4_WATER->GetIndex();
419 fMaterialParticleModelTable[index][p 419 fMaterialParticleModelTable[index][p] = it;
420 ++modelNbForMaterial; 420 ++modelNbForMaterial;
421 } 421 }
422 } 422 }
423 if (modelNbForMaterial == 0) { 423 if (modelNbForMaterial == 0) {
424 std::ostringstream oss; 424 std::ostringstream oss;
425 oss << "The material " << (*G4Material 425 oss << "The material " << (*G4Material::GetMaterialTable())[matID]->GetName()
426 << " and the particle " << p->GetP 426 << " and the particle " << p->GetParticleName();
427 oss << " does not have any model regis 427 oss << " does not have any model registered for the " << fName << " interaction.";
428 G4Exception("G4DNAModelInterface::Inse 428 G4Exception("G4DNAModelInterface::InsertModelInTable", "em0006", FatalException,
429 oss.str().c_str()); 429 oss.str().c_str());
430 return; // to make some compilers hap 430 return; // to make some compilers happy
431 } 431 }
432 } 432 }
433 } 433 }
434 } 434 }
435 435
436 G4VEmModel* G4DNAModelInterface::SelectModel(c 436 G4VEmModel* G4DNAModelInterface::SelectModel(const std::size_t& materialID,
437 c 437 const G4ParticleDefinition* particle,
438 c 438 const G4double& ekin)
439 { 439 {
440 // Output pointer 440 // Output pointer
441 G4VEmModel* model = nullptr; 441 G4VEmModel* model = nullptr;
442 442
443 // Get a reference to all the models for the 443 // Get a reference to all the models for the couple (material and particle)
444 auto modelData = fMaterialParticleModelTable 444 auto modelData = fMaterialParticleModelTable[materialID][particle];
445 445
446 // We must choose one of the model(s) accord 446 // We must choose one of the model(s) accordingly to the particle energy and the model energy
447 // range(s) 447 // range(s)
448 448
449 // Loop on all the models within the models 449 // Loop on all the models within the models vector and check if ekin is within the energy range.
450 auto DNAModel = dynamic_cast<G4VDNAModel*>(m 450 auto DNAModel = dynamic_cast<G4VDNAModel*>(modelData);
451 G4double lowL, highL; 451 G4double lowL, highL;
452 if (DNAModel == nullptr) { 452 if (DNAModel == nullptr) {
453 // ekin is in the energy range: we select 453 // ekin is in the energy range: we select the model and stop the loop.
454 lowL = modelData->LowEnergyLimit(); 454 lowL = modelData->LowEnergyLimit();
455 highL = modelData->HighEnergyLimit(); 455 highL = modelData->HighEnergyLimit();
456 if (ekin >= lowL && ekin < highL) { 456 if (ekin >= lowL && ekin < highL) {
457 // Select the model 457 // Select the model
458 458
459 model = modelData; 459 model = modelData;
460 // return model; 460 // return model;
461 // Quit the for loop 461 // Quit the for loop
462 // break; 462 // break;
463 } 463 }
464 // ekin is not in the energy range: we con 464 // ekin is not in the energy range: we continue the loop.
465 } 465 }
466 else { 466 else {
467 // ekin is in the energy range: we select 467 // ekin is in the energy range: we select the model and stop the loop.
468 lowL = DNAModel->GetLowELimit(materialID, 468 lowL = DNAModel->GetLowELimit(materialID, particle);
469 highL = DNAModel->GetHighELimit(materialID 469 highL = DNAModel->GetHighELimit(materialID, particle);
470 if (ekin >= lowL && ekin < highL) { 470 if (ekin >= lowL && ekin < highL) {
471 // Select the model 471 // Select the model
472 model = modelData; 472 model = modelData;
473 // return model; 473 // return model;
474 // Quit the for loop 474 // Quit the for loop
475 // break; 475 // break;
476 } 476 }
477 // ekin is not in the energy range: we con 477 // ekin is not in the energy range: we continue the loop.
478 } 478 }
479 //} 479 //}
480 return model; 480 return model;
481 } 481 }
482 482
483 //....oooOO0OOooo........oooOO0OOooo........oo 483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
484 484
485 G4double G4DNAModelInterface::GetNumMoleculePe 485 G4double G4DNAModelInterface::GetNumMoleculePerVolumeUnitForMaterial(const G4Material* mat)
486 { 486 {
487 return fMaterialMolPerVol[mat->GetIndex()]-> 487 return fMaterialMolPerVol[mat->GetIndex()]->at(mat->GetIndex());
488 } 488 }
489 489
490 //....oooOO0OOooo........oooOO0OOooo........oo 490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
491 491
492 G4double 492 G4double
493 G4DNAModelInterface::GetNumMolPerVolUnitForCom 493 G4DNAModelInterface::GetNumMolPerVolUnitForComponentInComposite(const G4Material* component,
494 494 const G4Material* composite)
495 { 495 {
496 return fMaterialMolPerVol[component->GetInde 496 return fMaterialMolPerVol[component->GetIndex()]->at(composite->GetIndex());
497 } 497 }
498 498
499 void G4DNAModelInterface::StreamInfo(std::ostr 499 void G4DNAModelInterface::StreamInfo(std::ostream& os) const
500 { 500 {
501 G4long prec = os.precision(5); 501 G4long prec = os.precision(5);
502 os << "===================================== 502 os << "======================================= Materials of " << std::setw(17)
503 << this->GetName() << " ============ 503 << this->GetName() << " ================================================"
504 << "\n"; 504 << "\n";
505 os << std::setw(15) << "Material#" << std::s 505 os << std::setw(15) << "Material#" << std::setw(13) << "Particle" << std::setw(35) << "Model"
506 << std::setw(17) << "LowLimit(MeV)" << st 506 << std::setw(17) << "LowLimit(MeV)" << std::setw(17) << "HighLimit(MeV)" << std::setw(13)
507 << "Fast" << std::setw(13) << "Stationary 507 << "Fast" << std::setw(13) << "Stationary" << std::setw(13) << "Chemistry" << G4endl;
508 for (const auto& it1 : fMaterialParticleMode 508 for (const auto& it1 : fMaterialParticleModelTable) {
509 os << std::setw(15) << (*G4Material::GetMa 509 os << std::setw(15) << (*G4Material::GetMaterialTable())[it1.first]->GetName();
510 for (const auto& it2 : it1.second) { 510 for (const auto& it2 : it1.second) {
511 os << std::setw(13) << it2.first->GetPar 511 os << std::setw(13) << it2.first->GetParticleName();
512 os << std::setw(35) << it2.second->GetNa 512 os << std::setw(35) << it2.second->GetName();
513 auto DNAModel = dynamic_cast<G4VDNAModel 513 auto DNAModel = dynamic_cast<G4VDNAModel*>(it2.second);
514 if (DNAModel == nullptr) { 514 if (DNAModel == nullptr) {
515 os << std::setw(17) << it2.second->Low 515 os << std::setw(17) << it2.second->LowEnergyLimit();
516 os << std::setw(17) << it2.second->Hig 516 os << std::setw(17) << it2.second->HighEnergyLimit();
517 } 517 }
518 else { 518 else {
519 auto lowL = DNAModel->GetLowELimit(it1 519 auto lowL = DNAModel->GetLowELimit(it1.first, it2.first);
520 auto highL = DNAModel->GetHighELimit(i 520 auto highL = DNAModel->GetHighELimit(it1.first, it2.first);
521 os << std::setw(17) << lowL; 521 os << std::setw(17) << lowL;
522 os << std::setw(17) << highL; 522 os << std::setw(17) << highL;
523 } 523 }
524 os << std::setw(13) << "no"; 524 os << std::setw(13) << "no";
525 os << std::setw(13) << "no"; 525 os << std::setw(13) << "no";
526 os << std::setw(13) << "no" << G4endl; 526 os << std::setw(13) << "no" << G4endl;
527 } 527 }
528 } 528 }
529 529
530 os << "===================================== 530 os << "========================================================================================"
531 "===================================== 531 "=================================================="
532 << G4endl; 532 << G4endl;
533 os.precision(prec); 533 os.precision(prec);
534 } 534 }
535 535