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