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