Geant4-11
G4DNAMolecularMaterial.cc
Go to the documentation of this file.
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 <utility>
31
33#include "G4Material.hh"
34#include "G4StateManager.hh"
35#include "G4Threading.hh"
36#include "G4AutoLock.hh"
37#include "G4StateManager.hh"
38#include "G4MoleculeTable.hh"
39
40using namespace std;
41
43
44namespace
45{
47}
48
49//------------------------------------------------------------------------------
50
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)) == false){
62 // None of the materials derives from a base material
63 return mat1 < mat2;
64 }
65 else if ((baseMat1 != nullptr) && (baseMat2 != nullptr)){
66 // Both materials derive from a base material
67 return baseMat1 < baseMat2;
68 }
69
70 else 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
81{
83 return fInstance;
84}
85
86//------------------------------------------------------------------------------
87
89{
90 fpCompFractionTable = nullptr;
91 fpCompDensityTable = nullptr;
93 fIsInitialized = false;
94 fNMaterials = 0;
95}
96
97//------------------------------------------------------------------------------
98
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){
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
141{
142 Create();
143}
144
145//------------------------------------------------------------------------------
146
148{
149 if (requestedState == G4State_Idle && G4StateManager::GetStateManager()
150 ->GetPreviousState() == G4State_PreInit){
151 Initialize();
152 }
153 return true;
154}
155
156//------------------------------------------------------------------------------
157
159 const G4DNAMolecularMaterial& /*rhs*/) :
161{
162 Create();
163}
164
165//------------------------------------------------------------------------------
166
169{
170 if (this == &rhs) return *this;
171 Create();
172 return *this;
173}
174
175//------------------------------------------------------------------------------
176
178{
179 Clear();
180}
181
182//------------------------------------------------------------------------------
183
185{
186 if (fIsInitialized){
187 return;
188 }
189
190 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
191
192 fNMaterials = materialTable->size();
193 // This is to prevent segment fault if materials are created later on
194 // Actually this creation should not be done
195
196 G4AutoLock l1(&aMutex);
197 if (fpCompFractionTable == nullptr){
198 fpCompFractionTable = new vector<ComponentMap>(materialTable->size());
199 }
200
201 G4Material* mat(nullptr);
202
203 for (std::size_t i = 0; i < fNMaterials; ++i){
204 mat = materialTable->at(i);
205 SearchMolecularMaterial(mat, mat, 1);
206 }
207
210 l1.unlock();
211
212 fIsInitialized = true;
213}
214
215//------------------------------------------------------------------------------
216
218{
220 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
221 fpCompDensityTable = new vector<ComponentMap>(
223
224 G4Material* parentMat;
225 const G4Material* compMat(nullptr);
226 G4double massFraction = -1;
227 G4double parentDensity = -1;
228
229 for (std::size_t i = 0; i < fNMaterials; ++i){
230 parentMat = materialTable->at(i);
231 ComponentMap& massFractionComp = (*fpCompFractionTable)[i];
232 ComponentMap& densityComp = (*fpCompDensityTable)[i];
233
234 parentDensity = parentMat->GetDensity();
235
236 for (auto it = massFractionComp.cbegin();
237 it != massFractionComp.cend(); ++it){
238 compMat = it->first;
239 massFraction = it->second;
240 densityComp[compMat] = massFraction * parentDensity;
241 compMat = nullptr;
242 massFraction = -1;
243 }
244 }
245 }
246 else{
247 G4ExceptionDescription exceptionDescription;
248 exceptionDescription << "The pointer fpCompFractionTable is not initialized"
249 << G4endl;
250 G4Exception("G4DNAMolecularMaterial::InitializeDensity",
251 "G4DNAMolecularMaterial001", FatalException,
252 exceptionDescription);
253 }
254}
255
256//------------------------------------------------------------------------------
257
259{
261 fpCompNumMolPerVolTable = new vector<ComponentMap>(fNMaterials);
262
263 const G4Material* compMat(nullptr);
264
265 for (std::size_t i = 0; i < fNMaterials; ++i){
266 ComponentMap& massFractionComp = (*fpCompFractionTable)[i];
267 ComponentMap& densityComp = (*fpCompDensityTable)[i];
268 ComponentMap& numMolPerVol = (*fpCompNumMolPerVolTable)[i];
269
270 for (auto it = massFractionComp.cbegin();
271 it != massFractionComp.cend(); ++it){
272 compMat = it->first;
273 numMolPerVol[compMat] = densityComp[compMat]
274 / compMat->GetMassOfMolecule();
275 compMat = nullptr;
276 }
277 }
278 }
279 else{
280 G4ExceptionDescription exceptionDescription;
281 exceptionDescription << "The pointer fpCompDensityTable is not initialized"
282 << G4endl;
283 G4Exception("G4DNAMolecularMaterial::InitializeNumMolPerVol",
284 "G4DNAMolecularMaterial002", FatalException,
285 exceptionDescription);
286 }
287}
288
289//------------------------------------------------------------------------------
290
291void
293 G4Material* molecularMaterial,
294 G4double fraction)
295{
296 ComponentMap& matComponent =
297 (*fpCompFractionTable)[parentMaterial->GetIndex()];
298
299 if (matComponent.empty()){
300 matComponent[molecularMaterial] = fraction;
301 return;
302 }
303
304 auto it = matComponent.find(molecularMaterial);
305
306 if (it == matComponent.cend()){
307 matComponent[molecularMaterial] = fraction;
308 }
309 else{
310 matComponent[molecularMaterial] = it->second + fraction;
311 // handle "base material"
312 }
313}
314
315//------------------------------------------------------------------------------
316
319 G4double currentFraction)
320{
321 if (material->GetMassOfMolecule() != 0.0){ // is a molecular material
322 RecordMolecularMaterial(parentMaterial, material, currentFraction);
323 return;
324 }
325
326 G4Material* compMat(nullptr);
327 G4double fraction = -1.;
328 std::map<G4Material*, G4double> matComponent = material->GetMatComponents();
329 auto it = matComponent.cbegin();
330
331 for (; it != matComponent.cend(); ++it){
332 compMat = it->first;
333 fraction = it->second;
334 if (compMat->GetMassOfMolecule() == 0.0){ // is not a molecular material
335 SearchMolecularMaterial(parentMaterial, compMat,
336 currentFraction * fraction);
337 }
338 else{ // is a molecular material
339 RecordMolecularMaterial(parentMaterial, compMat,
340 currentFraction * fraction);
341 }
342 }
343}
344
345//------------------------------------------------------------------------------
346
347const std::vector<G4double>*
349GetDensityTableFor(const G4Material* lookForMaterial) const
350{
351 if (!fpCompDensityTable){
352 if (fIsInitialized){
353 G4ExceptionDescription exceptionDescription;
354 exceptionDescription
355 << "The pointer fpCompDensityTable is not initialized will the "
356 "singleton of G4DNAMolecularMaterial "
357 << "has already been initialized." << G4endl;
358 G4Exception("G4DNAMolecularMaterial::GetDensityTableFor",
359 "G4DNAMolecularMaterial003", FatalException,
360 exceptionDescription);
361 }
362
363 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Init){
364 const_cast<G4DNAMolecularMaterial*>(this)->Initialize();
365 }
366 else{
367 G4ExceptionDescription exceptionDescription;
368 exceptionDescription
369 << "The geant4 application is at the wrong state. State must be: "
370 "G4State_Init."
371 << G4endl;
372 G4Exception("G4DNAMolecularMaterial::GetDensityTableFor",
373 "G4DNAMolecularMaterial_WRONG_STATE_APPLICATION",
374 FatalException, exceptionDescription);
375 }
376 }
377
378 auto it_askedDensityTable = fAskedDensityTable.find(lookForMaterial);
379
380 if (it_askedDensityTable != fAskedDensityTable.cend()){
381 return it_askedDensityTable->second;
382 }
383
384 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
385
386 std::vector<G4double>* output = new std::vector<G4double>(materialTable->size());
387
388 ComponentMap::const_iterator it;
389
390 G4bool materialWasNotFound = true;
391
392 for (std::size_t i = 0; i < fNMaterials; ++i){
393 ComponentMap& densityTable = (*fpCompDensityTable)[i];
394
395 it = densityTable.find(lookForMaterial);
396
397 if (it == densityTable.cend()){
398 (*output)[i] = 0.0;
399 }
400 else{
401 materialWasNotFound = false;
402 (*output)[i] = it->second;
403 }
404 }
405
406 if (materialWasNotFound){
407 PrintNotAMolecularMaterial("G4DNAMolecularMaterial::GetDensityTableFor",
408 lookForMaterial);
409 }
410
411 fAskedDensityTable.insert(make_pair(lookForMaterial, output));
412
413 return output;
414}
415
416//------------------------------------------------------------------------------
417
419 const G4Material* lookForMaterial) const
420{
421 if(lookForMaterial==nullptr) return nullptr;
422
424 if (fIsInitialized){
425 G4ExceptionDescription exceptionDescription;
426 exceptionDescription
427 << "The pointer fpCompNumMolPerVolTable is not initialized whereas "
428 "the singleton of G4DNAMolecularMaterial "
429 << "has already been initialized." << G4endl;
430 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolTableFor",
431 "G4DNAMolecularMaterial005", FatalException,
432 exceptionDescription);
433 }
434
435 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Init){
436 const_cast<G4DNAMolecularMaterial*>(this)->Initialize();
437 }
438 else{
439 G4ExceptionDescription exceptionDescription;
440 exceptionDescription
441 << "The geant4 application is at the wrong state. State must be : "
442 "G4State_Init."
443 << G4endl;
444 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolTableFor",
445 "G4DNAMolecularMaterial_WRONG_STATE_APPLICATION",
446 FatalException, exceptionDescription);
447 }
448 }
449
450 auto it_askedNumMolPerVolTable = fAskedNumPerVolTable.find(lookForMaterial);
451 if (it_askedNumMolPerVolTable != fAskedNumPerVolTable.cend()){
452 return it_askedNumMolPerVolTable->second;
453 }
454
455 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
456
457 std::vector<G4double>* output = new std::vector<G4double>(materialTable->size());
458
459 ComponentMap::const_iterator it;
460
461 G4bool materialWasNotFound = true;
462
463 for (std::size_t i = 0; i < fNMaterials; ++i){
464 ComponentMap& densityTable = (*fpCompNumMolPerVolTable)[i];
465
466 it = densityTable.find(lookForMaterial);
467
468 if (it == densityTable.cend()){
469 (*output)[i] = 0.0;
470 }
471 else{
472 materialWasNotFound = false;
473 (*output)[i] = it->second;
474 }
475 }
476
477 if (materialWasNotFound){
479 "G4DNAMolecularMaterial::GetNumMolPerVolTableFor", lookForMaterial);
480 }
481
482 fAskedNumPerVolTable.insert(make_pair(lookForMaterial, output));
483
484 return output;
485}
486
487//------------------------------------------------------------------------------
488
490PrintNotAMolecularMaterial(const char* methodName,
491 const G4Material* lookForMaterial) const
492{
493 auto it = fWarningPrinted.find(lookForMaterial);
494
495 if (it == fWarningPrinted.cend()){
496 G4ExceptionDescription exceptionDescription;
497 exceptionDescription << "The material " << lookForMaterial->GetName()
498 << " is not defined as a molecular material."
499 << G4endl
500 << "Meaning: The elements should be added to the "
501 "material using atom count rather than mass fraction "
502 "(cf. G4Material)"
503 << G4endl
504 << "If you want to use DNA processes on liquid water, you should better use "
505 "the NistManager to create the water material."
506 << G4endl
507 << "Since this message is displayed, it means that the DNA models will not "
508 "be called."
509 << "Please note that this message will only appear once even if you are "
510 "using other methods of G4DNAMolecularMaterial."
511 << G4endl;
512
513 G4Exception(methodName, "MATERIAL_NOT_DEFINE_USING_ATOM_COUNT", JustWarning,
514 exceptionDescription);
515 fWarningPrinted[lookForMaterial] = true;
516 }
517}
518
519//------------------------------------------------------------------------------
520
524{
525 int material_id = material->GetIndex();
526 auto it = fMaterialToMolecularConf.find(material_id);
527 if(it == fMaterialToMolecularConf.cend()) return nullptr;
528 return it->second;
529}
530
531//------------------------------------------------------------------------------
532
533void
537{
538 assert(material != nullptr);
539 int material_id = material->GetIndex();
540 fMaterialToMolecularConf[material_id] = molConf;
541}
542
543//------------------------------------------------------------------------------
544
545void
547 const G4String& molUserID)
548{
549 assert(material != nullptr);
550 int material_id = material->GetIndex();
551 fMaterialToMolecularConf[material_id] =
553}
554
555//------------------------------------------------------------------------------
556
557void
559 const G4String& molUserID)
560{
562
563 if(material == nullptr){
564 G4cout<< "Material " << materialName
565 << " was not found and therefore won't be linked to "
566 << molUserID << G4endl;
567 return;
568 }
570}
571
572//------------------------------------------------------------------------------
573
577{
578 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolForComponentInComposite",
579 "DEPRECATED",
580 FatalException,"Use standard method: GetNumMolPerVolTableFor"
581 " at the run initialization to retrieve a read-only table used"
582 " during stepping. The method is thread-safe.");
583 return 0.;
584}
585
586//------------------------------------------------------------------------------
587
591 const G4Material*,
592 G4double)
593{
594 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolForComponentInComposite",
595 "DEPRECATED",
596 FatalException,"Use standard method: GetNumMolPerVolTableFor"
597 " at the run initialization to retrieve a read-only table used"
598 " during stepping. The method is thread-safe.");
599 return 0.;
600}
G4ApplicationState
@ G4State_Init
@ G4State_Idle
@ G4State_PreInit
std::map< const G4Material *, G4double, CompareMaterial > ComponentMap
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
std::vector< G4Material * > G4MaterialTable
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:85
std::mutex G4Mutex
Definition: G4Threading.hh:81
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
G4DNAMolecularMaterial builds tables of molecular densities for chosen molecular materials....
std::map< int, G4MolecularConfiguration * > fMaterialToMolecularConf
void SetMolecularConfiguration(const G4Material *, G4MolecularConfiguration *)
Associate a molecular configuration to a G4material.
std::map< const G4Material *, std::vector< G4double > *, CompareMaterial > fAskedDensityTable
void RecordMolecularMaterial(G4Material *parentMaterial, G4Material *molecularMaterial, G4double fraction)
G4MolecularConfiguration * GetMolecularConfiguration(const G4Material *) const
void PrintNotAMolecularMaterial(const char *methodName, const G4Material *lookForMaterial) const
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
virtual G4bool Notify(G4ApplicationState requestedState)
G4double GetNumMoleculePerVolumeUnitForMaterial(const G4Material *mat)
Deprecated.
std::vector< ComponentMap > * fpCompFractionTable
std::vector< ComponentMap > * fpCompNumMolPerVolTable
static G4DNAMolecularMaterial * Instance()
std::map< const G4Material *, bool, CompareMaterial > fWarningPrinted
G4double GetNumMolPerVolForComponentInComposite(const G4Material *composite, const G4Material *component, G4double massFraction)
Deprecated.
const std::vector< G4double > * GetDensityTableFor(const G4Material *) const
Retrieve a table of volumetric mass densities (mass per unit volume) in the G4 unit system for chosen...
static G4DNAMolecularMaterial * fInstance
void SearchMolecularMaterial(G4Material *parentMaterial, G4Material *material, G4double currentFraction)
std::map< const G4Material *, std::vector< G4double > *, CompareMaterial > fAskedNumPerVolTable
G4DNAMolecularMaterial & operator=(const G4DNAMolecularMaterial &)
std::vector< ComponentMap > * fpCompDensityTable
G4double GetDensity() const
Definition: G4Material.hh:176
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:229
G4double GetMassOfMolecule() const
Definition: G4Material.hh:237
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:672
const G4String & GetName() const
Definition: G4Material.hh:173
size_t GetIndex() const
Definition: G4Material.hh:256
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:686
G4MolecularConfiguration * GetConfiguration(const G4String &, bool mustExist=true)
static G4MoleculeTable * Instance()
static G4StateManager * GetStateManager()
string material
Definition: eplot.py:19
Materials can be described as a derivation of existing "parent" materials in order to alter few of th...
bool operator()(const G4Material *mat1, const G4Material *mat2) const