| 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 // Author: Mathieu Karamitros
28 //
29
30 #include "G4DNAMolecularMaterial.hh"
31
32 #include "G4AutoLock.hh"
33 #include "G4Material.hh"
34 #include "G4MoleculeTable.hh"
35 #include "G4StateManager.hh"
36 #include "G4Threading.hh"
37
38 #include <utility>
39
40 using namespace std;
41
42 G4DNAMolecularMaterial* G4DNAMolecularMaterial::fInstance(nullptr);
43
44 namespace
45 {
46 G4Mutex aMutex = G4MUTEX_INITIALIZER;
47 }
48
49 //------------------------------------------------------------------------------
50
51 bool CompareMaterial::operator()(const G4Material* mat1,
52 const G4Material* mat2) const
53 {
54 if (mat1 == nullptr && mat2 == nullptr) return false; //(mat1 == mat2)
55 if (mat1 == nullptr) return true; // mat1 < mat2
56 if (mat2 == nullptr) return false; //mat2 < mat1
57
58 const G4Material* baseMat1 = mat1->GetBaseMaterial();
59 const G4Material* baseMat2 = mat2->GetBaseMaterial();
60
61 if ((baseMat1 == nullptr) && (baseMat2 == nullptr)){
62 // None of the materials derives from a base material
63 return mat1 < mat2;
64 }
65 if ((baseMat1 != nullptr) && (baseMat2 != nullptr)){
66 // Both materials derive from a base material
67 return baseMat1 < baseMat2;
68 }
69
70 if ((baseMat1 != nullptr) && (baseMat2 == nullptr)){
71 // Only the material 1 derives from a base material
72 return baseMat1 < mat2;
73 }
74 // only case baseMat1==nullptr && baseMat2 remains
75 return mat1 < baseMat2;
76 }
77
78 //------------------------------------------------------------------------------
79
80 G4DNAMolecularMaterial* G4DNAMolecularMaterial::Instance()
81 {
82 if (fInstance == nullptr) fInstance = new G4DNAMolecularMaterial();
83 return fInstance;
84 }
85
86 //------------------------------------------------------------------------------
87
88 void G4DNAMolecularMaterial::Create()
89 {
90 fpCompFractionTable = nullptr;
91 fpCompDensityTable = nullptr;
92 fpCompNumMolPerVolTable = nullptr;
93 fIsInitialized = false;
94 fNMaterials = 0;
95 }
96
97 //------------------------------------------------------------------------------
98
99 void G4DNAMolecularMaterial::Clear()
100 {
101 G4AutoLock l2(&aMutex);
102 if (fpCompFractionTable != nullptr){
103 fpCompFractionTable->clear();
104 delete fpCompFractionTable;
105 fpCompFractionTable = nullptr;
106 }
107 if (fpCompDensityTable != nullptr){
108 fpCompDensityTable->clear();
109 delete fpCompDensityTable;
110 fpCompDensityTable = nullptr;
111 }
112 if (fpCompNumMolPerVolTable != nullptr){
113 fpCompNumMolPerVolTable->clear();
114 delete fpCompNumMolPerVolTable;
115 fpCompNumMolPerVolTable = nullptr;
116 }
117
118 std::map<const G4Material*, std::vector<G4double>*, CompareMaterial>::iterator it;
119
120 for (it = fAskedDensityTable.begin(); it != fAskedDensityTable.end(); ++it){
121 if (it->second != nullptr){
122 delete it->second;
123 it->second = nullptr;
124 }
125 }
126
127 for (it = fAskedNumPerVolTable.begin(); it != fAskedNumPerVolTable.end(); ++it){
128 if (it->second != nullptr){
129 delete it->second;
130 it->second = nullptr;
131 }
132 }
133 l2.unlock();
134 }
135
136
137 //------------------------------------------------------------------------------
138
139 G4DNAMolecularMaterial::G4DNAMolecularMaterial()
140 {
141 Create();
142 }
143
144 //------------------------------------------------------------------------------
145
146 G4bool G4DNAMolecularMaterial::Notify(G4ApplicationState requestedState)
147 {
148 if (requestedState == G4State_Idle && G4StateManager::GetStateManager()
149 ->GetPreviousState() == G4State_PreInit){
150 Initialize();
151 }
152 return true;
153 }
154
155 //------------------------------------------------------------------------------
156
157 G4DNAMolecularMaterial::~G4DNAMolecularMaterial()
158 {
159 Clear();
160 }
161
162 //------------------------------------------------------------------------------
163
164 void G4DNAMolecularMaterial::Initialize()
165 {
166 if (fIsInitialized){
167 return;
168 }
169
170 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
171
172 fNMaterials = materialTable->size();
173 // This is to prevent segment fault if materials are created later on
174 // Actually this creation should not be done
175
176 G4AutoLock l1(&aMutex);
177 if (fpCompFractionTable == nullptr){
178 fpCompFractionTable = new vector<ComponentMap>(materialTable->size());
179 }
180
181 G4Material* mat(nullptr);
182
183 for (std::size_t i = 0; i < fNMaterials; ++i){
184 mat = materialTable->at(i);
185 SearchMolecularMaterial(mat, mat, 1);
186 }
187
188 InitializeDensity();
189 InitializeNumMolPerVol();
190 l1.unlock();
191
192 fIsInitialized = true;
193 }
194
195 //------------------------------------------------------------------------------
196
197 void G4DNAMolecularMaterial::InitializeDensity()
198 {
199 if (fpCompFractionTable != nullptr){
200 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
201 fpCompDensityTable = new vector<ComponentMap>(
202 G4Material::GetMaterialTable()->size());
203
204 G4Material* parentMat;
205 const G4Material* compMat(nullptr);
206 G4double massFraction = -1;
207 G4double parentDensity = -1;
208
209 for (std::size_t i = 0; i < fNMaterials; ++i){
210 parentMat = materialTable->at(i);
211 ComponentMap& massFractionComp = (*fpCompFractionTable)[i];
212 ComponentMap& densityComp = (*fpCompDensityTable)[i];
213
214 parentDensity = parentMat->GetDensity();
215
216 for (const auto& it : massFractionComp){
217 compMat = it.first;
218 massFraction = it.second;
219 densityComp[compMat] = massFraction * parentDensity;
220 compMat = nullptr;
221 massFraction = -1;
222 }
223 }
224 }
225 else{
226 G4ExceptionDescription exceptionDescription;
227 exceptionDescription << "The pointer fpCompFractionTable is not initialized"
228 << G4endl;
229 G4Exception("G4DNAMolecularMaterial::InitializeDensity",
230 "G4DNAMolecularMaterial001", FatalException,
231 exceptionDescription);
232 }
233 }
234
235 //------------------------------------------------------------------------------
236
237 void G4DNAMolecularMaterial::InitializeNumMolPerVol()
238 {
239 if (fpCompDensityTable != nullptr){
240 fpCompNumMolPerVolTable = new vector<ComponentMap>(fNMaterials);
241
242 const G4Material* compMat(nullptr);
243
244 for (std::size_t i = 0; i < fNMaterials; ++i){
245 ComponentMap& massFractionComp = (*fpCompFractionTable)[i];
246 ComponentMap& densityComp = (*fpCompDensityTable)[i];
247 ComponentMap& numMolPerVol = (*fpCompNumMolPerVolTable)[i];
248
249 for (auto& it : massFractionComp){
250 compMat = it.first;
251 numMolPerVol[compMat] = densityComp[compMat]
252 / compMat->GetMassOfMolecule();
253 compMat = nullptr;
254 }
255 }
256 }
257 else{
258 G4ExceptionDescription exceptionDescription;
259 exceptionDescription << "The pointer fpCompDensityTable is not initialized"
260 << G4endl;
261 G4Exception("G4DNAMolecularMaterial::InitializeNumMolPerVol",
262 "G4DNAMolecularMaterial002", FatalException,
263 exceptionDescription);
264 }
265 }
266
267 //------------------------------------------------------------------------------
268
269 void
270 G4DNAMolecularMaterial::RecordMolecularMaterial(G4Material* parentMaterial,
271 G4Material* molecularMaterial,
272 G4double fraction)
273 {
274 ComponentMap& matComponent =
275 (*fpCompFractionTable)[parentMaterial->GetIndex()];
276
277 if (matComponent.empty()){
278 matComponent[molecularMaterial] = fraction;
279 return;
280 }
281
282 auto it = matComponent.find(molecularMaterial);
283
284 if (it == matComponent.cend()){
285 matComponent[molecularMaterial] = fraction;
286 }
287 else{
288 matComponent[molecularMaterial] = it->second + fraction;
289 // handle "base material"
290 }
291 }
292
293 //------------------------------------------------------------------------------
294
295 void G4DNAMolecularMaterial::SearchMolecularMaterial(G4Material* parentMaterial,
296 G4Material* material,
297 G4double currentFraction)
298 {
299 if (material->GetMassOfMolecule() != 0.0){ // is a molecular material
300 RecordMolecularMaterial(parentMaterial, material, currentFraction);
301 return;
302 }
303
304 G4Material* compMat(nullptr);
305 G4double fraction = -1.;
306 std::map<G4Material*, G4double> matComponent = material->GetMatComponents();
307 auto it = matComponent.cbegin();
308
309 for (; it != matComponent.cend(); ++it){
310 compMat = it->first;
311 fraction = it->second;
312 if (compMat->GetMassOfMolecule() == 0.0){ // is not a molecular material
313 SearchMolecularMaterial(parentMaterial, compMat,
314 currentFraction * fraction);
315 }
316 else{ // is a molecular material
317 RecordMolecularMaterial(parentMaterial, compMat,
318 currentFraction * fraction);
319 }
320 }
321 }
322
323 //------------------------------------------------------------------------------
324
325 const std::vector<G4double>*
326 G4DNAMolecularMaterial::
327 GetDensityTableFor(const G4Material* lookForMaterial) const
328 {
329 if (fpCompDensityTable == nullptr){
330 if (fIsInitialized){
331 G4ExceptionDescription exceptionDescription;
332 exceptionDescription
333 << "The pointer fpCompDensityTable is not initialized will the "
334 "singleton of G4DNAMolecularMaterial "
335 << "has already been initialized." << G4endl;
336 G4Exception("G4DNAMolecularMaterial::GetDensityTableFor",
337 "G4DNAMolecularMaterial003", FatalException,
338 exceptionDescription);
339 }
340
341 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Init){
342 const_cast<G4DNAMolecularMaterial*>(this)->Initialize();
343 }
344 else{
345 G4ExceptionDescription exceptionDescription;
346 exceptionDescription
347 << "The geant4 application is at the wrong state. State must be: "
348 "G4State_Init."
349 << G4endl;
350 G4Exception("G4DNAMolecularMaterial::GetDensityTableFor",
351 "G4DNAMolecularMaterial_WRONG_STATE_APPLICATION",
352 FatalException, exceptionDescription);
353 }
354 }
355
356 auto it_askedDensityTable = fAskedDensityTable.find(lookForMaterial);
357
358 if (it_askedDensityTable != fAskedDensityTable.cend()){
359 return it_askedDensityTable->second;
360 }
361
362 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
363
364 auto output = new std::vector<G4double>(materialTable->size());
365
366 ComponentMap::const_iterator it;
367
368 G4bool materialWasNotFound = true;
369
370 for (std::size_t i = 0; i < fNMaterials; ++i){
371 ComponentMap& densityTable = (*fpCompDensityTable)[i];
372
373 it = densityTable.find(lookForMaterial);
374
375 if (it == densityTable.cend()){
376 (*output)[i] = 0.0;
377 }
378 else{
379 materialWasNotFound = false;
380 (*output)[i] = it->second;
381 }
382 }
383
384 if (materialWasNotFound){
385 PrintNotAMolecularMaterial("G4DNAMolecularMaterial::GetDensityTableFor",
386 lookForMaterial);
387 }
388
389 fAskedDensityTable.insert(make_pair(lookForMaterial, output));
390
391 return output;
392 }
393
394 //------------------------------------------------------------------------------
395
396 const std::vector<G4double>* G4DNAMolecularMaterial::GetNumMolPerVolTableFor(
397 const G4Material* lookForMaterial) const
398 {
399 if(lookForMaterial==nullptr) return nullptr;
400
401 if (fpCompNumMolPerVolTable == nullptr){
402 if (fIsInitialized){
403 G4ExceptionDescription exceptionDescription;
404 exceptionDescription
405 << "The pointer fpCompNumMolPerVolTable is not initialized whereas "
406 "the singleton of G4DNAMolecularMaterial "
407 << "has already been initialized." << G4endl;
408 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolTableFor",
409 "G4DNAMolecularMaterial005", FatalException,
410 exceptionDescription);
411 }
412
413 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Init){
414 const_cast<G4DNAMolecularMaterial*>(this)->Initialize();
415 }
416 else{
417 G4ExceptionDescription exceptionDescription;
418 exceptionDescription
419 << "The geant4 application is at the wrong state. State must be : "
420 "G4State_Init."
421 << G4endl;
422 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolTableFor",
423 "G4DNAMolecularMaterial_WRONG_STATE_APPLICATION",
424 FatalException, exceptionDescription);
425 }
426 }
427
428 auto it_askedNumMolPerVolTable = fAskedNumPerVolTable.find(lookForMaterial);
429 if (it_askedNumMolPerVolTable != fAskedNumPerVolTable.cend()){
430 return it_askedNumMolPerVolTable->second;
431 }
432
433 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
434
435 auto output = new std::vector<G4double>(materialTable->size());
436
437 ComponentMap::const_iterator it;
438
439 G4bool materialWasNotFound = true;
440
441 for (std::size_t i = 0; i < fNMaterials; ++i){
442 ComponentMap& densityTable = (*fpCompNumMolPerVolTable)[i];
443
444 it = densityTable.find(lookForMaterial);
445
446 if (it == densityTable.cend()){
447 (*output)[i] = 0.0;
448 }
449 else{
450 materialWasNotFound = false;
451 (*output)[i] = it->second;
452 }
453 }
454
455 if (materialWasNotFound){
456 PrintNotAMolecularMaterial(
457 "G4DNAMolecularMaterial::GetNumMolPerVolTableFor", lookForMaterial);
458 }
459
460 fAskedNumPerVolTable.insert(make_pair(lookForMaterial, output));
461
462 return output;
463 }
464
465 //------------------------------------------------------------------------------
466
467 void G4DNAMolecularMaterial::
468 PrintNotAMolecularMaterial(const char* methodName,
469 const G4Material* lookForMaterial) const
470 {
471 auto it = fWarningPrinted.find(lookForMaterial);
472
473 if (it == fWarningPrinted.cend()){
474 G4ExceptionDescription exceptionDescription;
475 exceptionDescription << "The material " << lookForMaterial->GetName()
476 << " is not defined as a molecular material."
477 << G4endl
478 << "Meaning: The elements should be added to the "
479 "material using atom count rather than mass fraction "
480 "(cf. G4Material)"
481 << G4endl
482 << "If you want to use DNA processes on liquid water, you should better use "
483 "the NistManager to create the water material."
484 << G4endl
485 << "Since this message is displayed, it means that the DNA models will not "
486 "be called."
487 << "Please note that this message will only appear once even if you are "
488 "using other methods of G4DNAMolecularMaterial."
489 << G4endl;
490
491 G4Exception(methodName, "MATERIAL_NOT_DEFINE_USING_ATOM_COUNT", JustWarning,
492 exceptionDescription);
493 fWarningPrinted[lookForMaterial] = true;
494 }
495 }
496
497 //------------------------------------------------------------------------------
498
499 G4MolecularConfiguration*
500 G4DNAMolecularMaterial::
501 GetMolecularConfiguration(const G4Material* material) const
502 {
503 auto material_id = (G4int)material->GetIndex();
504 auto it = fMaterialToMolecularConf.find(material_id);
505 if(it == fMaterialToMolecularConf.cend()) return nullptr;
506 return it->second;
507 }
508
509 //------------------------------------------------------------------------------
510
511 void
512 G4DNAMolecularMaterial::
513 SetMolecularConfiguration(const G4Material* material,
514 G4MolecularConfiguration* molConf)
515 {
516 assert(material != nullptr);
517 auto material_id = (G4int)material->GetIndex();
518 fMaterialToMolecularConf[material_id] = molConf;
519 }
520
521 //------------------------------------------------------------------------------
522
523 void
524 G4DNAMolecularMaterial::SetMolecularConfiguration(const G4Material* material,
525 const G4String& molUserID)
526 {
527 assert(material != nullptr);
528 auto material_id = (G4int)material->GetIndex();
529 fMaterialToMolecularConf[material_id] =
530 G4MoleculeTable::Instance()->GetConfiguration(molUserID, true);
531 }
532
533 //------------------------------------------------------------------------------
534
535 void
536 G4DNAMolecularMaterial::SetMolecularConfiguration(const G4String& materialName,
537 const G4String& molUserID)
538 {
539 G4Material* material = G4Material::GetMaterial(materialName);
540
541 if(material == nullptr){
542 G4cout<< "Material " << materialName
543 << " was not found and therefore won't be linked to "
544 << molUserID << G4endl;
545 return;
546 }
547 SetMolecularConfiguration(material, molUserID);
548 }
549
550 //------------------------------------------------------------------------------
551
552 G4double
553 G4DNAMolecularMaterial::
554 GetNumMoleculePerVolumeUnitForMaterial(const G4Material*)
555 {
556 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolForComponentInComposite",
557 "DEPRECATED",
558 FatalException,"Use standard method: GetNumMolPerVolTableFor"
559 " at the run initialization to retrieve a read-only table used"
560 " during stepping. The method is thread-safe.");
561 return 0.;
562 }
563
564 //------------------------------------------------------------------------------
565
566 G4double
567 G4DNAMolecularMaterial::
568 GetNumMolPerVolForComponentInComposite(const G4Material*,
569 const G4Material*,
570 G4double)
571 {
572 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolForComponentInComposite",
573 "DEPRECATED",
574 FatalException,"Use standard method: GetNumMolPerVolTableFor"
575 " at the run initialization to retrieve a read-only table used"
576 " during stepping. The method is thread-safe.");
577 return 0.;
578 }
579