| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Contact authors: S. Meylan, C. Villagrasa
28 // email: sylvain.meylan@symalgo-tech.com, car
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
39
40 //....oooOO0OOooo........oooOO0OOooo........oo
41
42 void G4DNAModelInterface::Initialise(const G4P
43 {
44 // Those two statements are necessary to ove
45 // (ionisation, elastic, etc...). Indeed, wi
46 // themselves their energy limits per materi
47 // in the G4DNAProcess classes.
48 //
49
50 fpG4_WATER = G4Material::GetMaterial("G4_WAT
51
52 SetLowEnergyLimit(0.);
53 SetHighEnergyLimit(DBL_MAX);
54
55 fpParticleChangeForGamma = GetParticleChange
56
57 // Loop on all the registered models to init
58 for (auto & fRegisteredModel : fRegisteredMo
59 fRegisteredModel->SetParticleChange(fpPart
60 fRegisteredModel->Initialise(particle, cut
61 }
62 // used to retrieve the model corresponding
63 BuildMaterialParticleModelTable(particle);
64
65 BuildMaterialMolPerVolTable();
66
67 StreamInfo(G4cout);
68
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oo
72
73 G4double G4DNAModelInterface::CrossSectionPerV
74
75
76 {
77 // Method to return the crossSection * nbMol
78 // Process class then calculates the path.
79 // The cross section is calculated in the re
80 // Two cases are handled here: normal materi
81 //
82 // Idea:
83 // *** Simple material ***
84 // Ask for the cross section of the chosen m
85 // Multiply it by the number of medium molec
86 // Return the value.
87 // *** Composite material ***
88 // Ask for the cross section of the chosen m
89 // Apply a factor to each cross section and
90 // component per composite volume unit. The
91
92 // To reset the sampledMat variable.
93 // Can be used by user to retrieve current c
94 fSampledMat = 0;
95
96 // This is the value to be sum up and to be
97 G4double crossSectionTimesNbMolPerVol(0.);
98
99 // Reset the map saving the material and the
100 // Used in SampleSecondaries if the interact
101 // composite
102 fMaterialCS.clear();
103
104 // This is the value to be used by SampleSec
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->GetI
114
115 // Use the table to get the model
116 auto model = SelectModel(materialID, p, ek
117
118 // Get the nunber of molecules per volume
119
120 // Calculate the cross section times the n
121 if (model != nullptr) {
122 if (dynamic_cast<G4VDNAModel*>(model) ==
123 // water material models only
124 crossSectionTimesNbMolPerVol = model->
125 }
126 else {
127 crossSectionTimesNbMolPerVol = model->
128 }
129 }
130 else // no model was selected, we are out
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
141 // Maybe MatComponents map is overrided by
142 auto componentsMap = material->GetMatCompo
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 fracti
151 // G4double massFraction = it->second;
152
153 // Get the number of component molecules
154 G4double nbMoleculeOfComponentInComposit
155 GetNumMolPerVolUnitForComponentInCompo
156 G4cout << " ==========>component : " <<
157 << " nbMoleculeOfComponentInCompo
158 << G4endl;
159
160 // Get the current component name
161 const std::size_t & componentID = compon
162
163 // Retrieve the model corresponding to t
164 auto model = SelectModel(componentID, p,
165
166 // Add the component part of the cross s
167 // The component cross section is multip
168 // scaled by the mass fraction.
169 G4double crossSection;
170 if (model != nullptr) {
171 if (dynamic_cast<G4VDNAModel*>(model)
172 // water models
173 crossSection =
174 model->CrossSectionPerVolume(compo
175 / GetNumMoleculePerVolumeUnitForMa
176 }
177 else {
178 crossSection = model->CrossSectionPe
179 / GetNumMoleculePerVo
180 }
181 crossSectionTimesNbMolPerVol = nbMolec
182 }
183 else // no model was selected, we are o
184 {
185 crossSectionTimesNbMolPerVol = 0.;
186 }
187
188 // Save the component name and its calcu
189 // To be used by sampling secondaries if
190 fMaterialCS[componentID] = crossSectionT
191
192 // Save the component name and its calcu
193 // To be used by sampling secondaries if
194 fCSsumTot += crossSectionTimesNbMolPerVo
195 }
196 crossSectionTimesNbMolPerVol = fCSsumTot;
197 }
198
199 // return the cross section times the number
200 // the path of the interaction will be calcu
201 return crossSectionTimesNbMolPerVol;
202 }
203
204 //....oooOO0OOooo........oooOO0OOooo........oo
205
206 void G4DNAModelInterface::SampleSecondaries(st
207 co
208 co
209 G4
210 {
211 // To call the sampleSecondaries method of t
212 // In the case of composite material, we nee
213 // To do so we use a random sampling on the
214 // CrossSectionPerVolume method. If we enter
215 // (and process) has been chosen for the cur
216
217 std::size_t materialID;
218
219 // *******************************
220 // Material is not a composite
221 // *******************************
222 //
223 if (couple->GetMaterial()->GetMatComponents(
224 materialID = couple->GetMaterial()->GetInd
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
236 G4double rand = G4UniformRand() * fCSsumTo
237
238 G4double cumulCS(0);
239
240 G4bool result = false;
241
242 // We loop on each component cumulated cro
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
248 while (rand > cumulCS) {
249 // Check if the sampling is ok
250 if (it == ite) {
251 G4Exception(
252 "G4DNAModelManager::SampleSecondarie
253 "The random component selection has
254 "having a selected component");
255 return; // to make some compilers hap
256 }
257
258 // Set the cumulated value for the itera
259 cumulCS += it->second;
260
261 // Check if we have reach the material t
262 // The DBL_MAX is here to take into acco
263 // to force elastic sampleSecondaries wh
264 // Used when paticle energy is lower tha
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
279 // taken into account
280
281 G4Exception("G4DNAModelManager::SampleSe
282 "The random component select
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
294 // if (materialName.find("_MODIFIED") != G4
295 // materialName = materialName.substr(0,
296 // }
297
298 fSampledMat = materialID;
299
300 auto model = SelectModel(materialID, aDynami
301 aDynamicParticle->G
302
303 model->SampleSecondaries(fVect, couple, aDyn
304 }
305
306 void G4DNAModelInterface::RegisterModel(G4VEmM
307 {
308 fRegisteredModels.push_back(model);
309 }
310
311 //....oooOO0OOooo........oooOO0OOooo........oo
312
313 void G4DNAModelInterface::BuildMaterialParticl
314 {
315 // Method to build a map: [material][particl
316 // The map is used to retrieve the correct m
317
318 // Loop on all materials registered in the s
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 compos
324 auto componentMap = mat->GetMatComponents(
325 if (componentMap.empty()) {
326 // Get the material name
327 const std::size_t & matID = mat->GetInde
328 InsertModelInTable(matID, p);
329 }
330 // if the material is a composite material
331 // register them
332 else {
333 // Loop on all the components of the mat
334 for (const auto& itComp : componentMap)
335 G4Material* component = itComp.first;
336 // Check that the component is not its
337 if (!component->GetMatComponents().emp
338 std::ostringstream oss;
339 oss << "Material " << mat->GetName()
340 oss << " " << component->GetName();
341 G4Exception("G4DNAModelManager::Buil
342 FatalException, oss.str(
343 return; // to make some compilers h
344 }
345 // Get the current component name
346 const std::size_t & compID = component
347 // If there is a model then insert the
348 // contains a if statement to check we
349 // normal material before.
350 InsertModelInTable(compID, p);
351 // move forward the iterator
352 }
353 }
354 }
355 }
356
357 //....oooOO0OOooo........oooOO0OOooo........oo
358
359 void G4DNAModelInterface::BuildMaterialMolPerV
360 {
361 // To be sure the G4DNAMolecularMaterial is
362 G4DNAMolecularMaterial::Instance()->Initiali
363
364 G4MaterialTable* materialTable = G4Material:
365
366 // Loop on all the materials inside the "mat
367 for (auto currentMaterial : *materialTable)
368 // Current material
369 // Current material name
370 const std::size_t & currentMatID = current
371
372 // Will the material be used in this inter
373 // Loop on all the materials that can be d
374 auto it = fMaterialParticleModelTable.begi
375 auto ite = fMaterialParticleModelTable.end
376 for (; it != ite; it++) {
377 const std::size_t & materialID = it->fir
378
379 if (materialID == currentMatID) {
380 const std::vector<G4double>* numMolPer
381 G4DNAMolecularMaterial::Instance()->
382 fMaterialMolPerVol[materialID] = numMo
383 }
384 }
385 }
386 }
387
388 //....oooOO0OOooo........oooOO0OOooo........oo
389
390 void G4DNAModelInterface::InsertModelInTable(c
391 {
392 // To insert the model(s) in the table Mater
393
394 // First, we need to check if the current ma
395 // This is possible because of the composite
396 // add the independant M1 material. This cas
397 // table is the way to avoid it.
398 //
399 // Check if the current material and particl
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)
405 // Check for the particle
406 if (fMaterialParticleModelTable[matID].fin
407 G4int modelNbForMaterial = 0;
408 for (const auto& it : fRegisteredModels)
409 auto model = dynamic_cast<G4VDNAModel*
410 if (model != nullptr) {
411 if (model->IsParticleExistingInModel
412 fMaterialParticleModelTable[matID]
413 // and add one to the "there is a
414 ++modelNbForMaterial;
415 }
416 }
417 else {
418 auto index = fpG4_WATER->GetIndex();
419 fMaterialParticleModelTable[index][p
420 ++modelNbForMaterial;
421 }
422 }
423 if (modelNbForMaterial == 0) {
424 std::ostringstream oss;
425 oss << "The material " << (*G4Material
426 << " and the particle " << p->GetP
427 oss << " does not have any model regis
428 G4Exception("G4DNAModelInterface::Inse
429 oss.str().c_str());
430 return; // to make some compilers hap
431 }
432 }
433 }
434 }
435
436 G4VEmModel* G4DNAModelInterface::SelectModel(c
437 c
438 c
439 {
440 // Output pointer
441 G4VEmModel* model = nullptr;
442
443 // Get a reference to all the models for the
444 auto modelData = fMaterialParticleModelTable
445
446 // We must choose one of the model(s) accord
447 // range(s)
448
449 // Loop on all the models within the models
450 auto DNAModel = dynamic_cast<G4VDNAModel*>(m
451 G4double lowL, highL;
452 if (DNAModel == nullptr) {
453 // ekin is in the energy range: we select
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 con
465 }
466 else {
467 // ekin is in the energy range: we select
468 lowL = DNAModel->GetLowELimit(materialID,
469 highL = DNAModel->GetHighELimit(materialID
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 con
478 }
479 //}
480 return model;
481 }
482
483 //....oooOO0OOooo........oooOO0OOooo........oo
484
485 G4double G4DNAModelInterface::GetNumMoleculePe
486 {
487 return fMaterialMolPerVol[mat->GetIndex()]->
488 }
489
490 //....oooOO0OOooo........oooOO0OOooo........oo
491
492 G4double
493 G4DNAModelInterface::GetNumMolPerVolUnitForCom
494
495 {
496 return fMaterialMolPerVol[component->GetInde
497 }
498
499 void G4DNAModelInterface::StreamInfo(std::ostr
500 {
501 G4long prec = os.precision(5);
502 os << "=====================================
503 << this->GetName() << " ============
504 << "\n";
505 os << std::setw(15) << "Material#" << std::s
506 << std::setw(17) << "LowLimit(MeV)" << st
507 << "Fast" << std::setw(13) << "Stationary
508 for (const auto& it1 : fMaterialParticleMode
509 os << std::setw(15) << (*G4Material::GetMa
510 for (const auto& it2 : it1.second) {
511 os << std::setw(13) << it2.first->GetPar
512 os << std::setw(35) << it2.second->GetNa
513 auto DNAModel = dynamic_cast<G4VDNAModel
514 if (DNAModel == nullptr) {
515 os << std::setw(17) << it2.second->Low
516 os << std::setw(17) << it2.second->Hig
517 }
518 else {
519 auto lowL = DNAModel->GetLowELimit(it1
520 auto highL = DNAModel->GetHighELimit(i
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