00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "G4DNAMolecularMaterial.hh"
00029 #include "G4Material.hh"
00030 #include <utility>
00031 #include "G4StateManager.hh"
00032
00033 using namespace std;
00034
00035
00036 bool CompareMaterial::operator() (const G4Material* mat1, const G4Material* mat2) const
00037 {
00038 if(mat1==0 && mat2==0) return false;
00039 if(mat1==0) return true;
00040 if(mat2==0) return false;
00041
00042 const G4Material* baseMat1 = mat1->GetBaseMaterial();
00043 const G4Material* baseMat2 = mat2->GetBaseMaterial();
00044
00045 if((baseMat1 || baseMat2) == 0)
00046 {
00047 return mat1 < mat2;
00048 }
00049 else if(baseMat1 && baseMat2)
00050 {
00051 return baseMat1 < baseMat2;
00052 }
00053
00054 else if(baseMat1 && (baseMat2 == 0))
00055 {
00056 return baseMat1 < mat2;
00057 }
00058
00059 return mat1 < baseMat2;
00060 }
00061
00062 G4DNAMolecularMaterial* G4DNAMolecularMaterial::fInstance(0);
00063
00064 G4DNAMolecularMaterial* G4DNAMolecularMaterial::Instance()
00065 {
00066 if(! fInstance) new G4DNAMolecularMaterial();
00067 return fInstance;
00068 }
00069
00070 void G4DNAMolecularMaterial::DeleteInstance()
00071 {
00072 delete fInstance;
00073 fInstance = 0;
00074 }
00075
00076 void G4DNAMolecularMaterial::Create()
00077 {
00078 fpCompFractionTable = 0;
00079 fpCompDensityTable = 0;
00080 fpCompNumMolPerVolTable = 0;
00081 fIsInitialized = false;
00082 fInstance = this;
00083 }
00084
00085 G4DNAMolecularMaterial::G4DNAMolecularMaterial() :G4VStateDependent()
00086 {
00087 Create();
00088 fInstance = this;
00089 }
00090
00091 G4bool G4DNAMolecularMaterial::Notify(G4ApplicationState requestedState)
00092 {
00093 if(requestedState == G4State_Idle) Initialize();
00094 return true;
00095 }
00096
00097 G4DNAMolecularMaterial::G4DNAMolecularMaterial(const G4DNAMolecularMaterial& ) : G4VStateDependent()
00098 {
00099 Create();
00100 }
00101
00102 G4DNAMolecularMaterial& G4DNAMolecularMaterial::operator=(const G4DNAMolecularMaterial& rhs)
00103 {
00104 if(this == &rhs) return *this;
00105 Create();
00106 return *this;
00107 }
00108
00109 G4DNAMolecularMaterial::~G4DNAMolecularMaterial()
00110 {
00111 if(fpCompFractionTable)
00112 {
00113 fpCompFractionTable->clear();
00114 delete fpCompFractionTable;
00115 fpCompFractionTable = 0;
00116 }
00117 if(fpCompDensityTable)
00118 {
00119 fpCompDensityTable->clear();
00120 delete fpCompDensityTable;
00121 fpCompDensityTable = 0;
00122 }
00123 if(fpCompNumMolPerVolTable)
00124 {
00125 fpCompNumMolPerVolTable->clear();
00126 delete fpCompNumMolPerVolTable;
00127 fpCompNumMolPerVolTable = 0;
00128 }
00129
00130 std::map<const G4Material*,std::vector<double>*,CompareMaterial>::iterator it;
00131
00132 for(it= fAskedDensityTable.begin() ; it != fAskedDensityTable.end() ;it++)
00133 {
00134 if(it->second)
00135 {
00136 delete it->second;
00137 it->second = 0;
00138 }
00139 }
00140
00141 for(it= fAskedNumPerVolTable.begin() ; it != fAskedNumPerVolTable.end() ;it++)
00142 {
00143 if(it->second)
00144 {
00145 delete it->second;
00146 it->second = 0;
00147 }
00148 }
00149 }
00150
00151 void G4DNAMolecularMaterial::RecordMolecularMaterial(G4Material* parentMaterial, G4Material* molecularMaterial, G4double fraction)
00152 {
00153 ComponentMap& matComponent = (*fpCompFractionTable)[parentMaterial->GetIndex()];
00154
00155 if(matComponent.empty())
00156 {
00157 matComponent[molecularMaterial] = fraction;
00158 return;
00159 }
00160
00161 ComponentMap::iterator it = matComponent.find(molecularMaterial);
00162
00163 if(it == matComponent.end())
00164 {
00165 matComponent[molecularMaterial] = fraction;
00166 }
00167 else
00168 {
00169 matComponent[molecularMaterial] = it->second + fraction;
00170 }
00171 }
00172
00173 void G4DNAMolecularMaterial::SearchMolecularMaterial(G4Material* parentMaterial, G4Material* material, double currentFraction)
00174 {
00175 if(material->GetMassOfMolecule() != 0.0)
00176 {
00177 RecordMolecularMaterial(parentMaterial,material,currentFraction);
00178 return;
00179 }
00180
00181 G4Material* compMat(0);
00182 G4double fraction = -1;
00183 std::map<G4Material*,G4double> matComponent = material->GetMatComponents();
00184 std::map<G4Material*,G4double>::iterator it = matComponent.begin();
00185
00186 for( ; it!=matComponent.end() ; it++)
00187 {
00188 compMat = it->first;
00189 fraction = it->second;
00190 if(compMat->GetMassOfMolecule() == 0.0)
00191 {
00192 SearchMolecularMaterial(parentMaterial,compMat,currentFraction*fraction);
00193 }
00194 else
00195 {
00196 RecordMolecularMaterial(parentMaterial,compMat,currentFraction*fraction);
00197 }
00198
00199 compMat = 0;
00200 fraction = -1;
00201 }
00202 }
00203
00204 void G4DNAMolecularMaterial::InitializeDensity()
00205 {
00206 if(fpCompFractionTable)
00207 {
00208 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00209 fpCompDensityTable = new vector<ComponentMap>(G4Material::GetMaterialTable()->size());
00210
00211 G4Material* parentMat;
00212 const G4Material* compMat(0);
00213 double massFraction = -1;
00214 double parentDensity = -1;
00215
00216 for(int i = 0 ; i < int(materialTable->size()) ; i++)
00217 {
00218 parentMat = materialTable->at(i);
00219 ComponentMap& massFractionComp = (*fpCompFractionTable)[i];
00220 ComponentMap& densityComp = (*fpCompDensityTable)[i];
00221
00222 parentDensity = parentMat->GetDensity();
00223
00224 for(ComponentMap::iterator it = massFractionComp.begin() ; it!=massFractionComp.end() ; it++)
00225 {
00226 compMat = it->first;
00227 massFraction = it->second;
00228 densityComp[compMat] = massFraction*parentDensity;
00229 compMat = 0;
00230 massFraction = -1;
00231 }
00232 }
00233 }
00234 else
00235 {
00236 G4ExceptionDescription exceptionDescription;
00237 exceptionDescription << "The pointer fpCompFractionTable is not initialized" << G4endl;
00238 G4Exception("G4DNAMolecularMaterial::InitializeDensity","G4DNAMolecularMaterial001",
00239 FatalException,exceptionDescription);
00240 }
00241 }
00242
00243 void G4DNAMolecularMaterial::InitializeNumMolPerVol()
00244 {
00245 if(fpCompDensityTable)
00246 {
00247 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00248 fpCompNumMolPerVolTable = new vector<ComponentMap>(G4Material::GetMaterialTable()->size());
00249
00250 const G4Material* compMat(0);
00251
00252 for(int i = 0 ; i < int(materialTable->size()) ; i++)
00253 {
00254 ComponentMap& massFractionComp = (*fpCompFractionTable)[i];
00255 ComponentMap& densityComp = (*fpCompDensityTable)[i];
00256 ComponentMap& numMolPerVol = (*fpCompNumMolPerVolTable)[i];
00257
00258 for(ComponentMap::iterator it = massFractionComp.begin() ; it!=massFractionComp.end() ; it++)
00259 {
00260 compMat = it->first;
00261 numMolPerVol[compMat] = densityComp[compMat]/ compMat->GetMassOfMolecule();
00262 compMat = 0;
00263 }
00264 }
00265 }
00266 else
00267 {
00268 G4ExceptionDescription exceptionDescription;
00269 exceptionDescription << "The pointer fpCompDensityTable is not initialized" << G4endl;
00270 G4Exception("G4DNAMolecularMaterial::InitializeNumMolPerVol","G4DNAMolecularMaterial002",
00271 FatalException,exceptionDescription);
00272 }
00273 }
00274
00275 void G4DNAMolecularMaterial::Initialize()
00276 {
00277 if(fIsInitialized) return;
00278
00279 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00280
00281 if(fpCompFractionTable==0)
00282 {
00283 fpCompFractionTable = new vector<ComponentMap>(materialTable->size());
00284 }
00285
00286 G4Material* mat(0);
00287
00288 for(int i = 0 ; i < int(materialTable->size()) ; i++)
00289 {
00290 mat = materialTable->at(i);
00291 SearchMolecularMaterial(mat,mat,1);
00292
00293 mat = 0;
00294 }
00295
00296 InitializeDensity();
00297 InitializeNumMolPerVol();
00298
00299 fIsInitialized = true;
00300 }
00301
00302 const std::vector<double>* G4DNAMolecularMaterial::GetDensityTableFor(const G4Material* lookForMaterial) const
00303 {
00304 if(!fpCompDensityTable)
00305 {
00306 if(fIsInitialized)
00307 {
00308 G4ExceptionDescription exceptionDescription;
00309 exceptionDescription << "The pointer fpCompDensityTable is not initialized will the singleton of G4DNAMolecularMaterial "
00310 << "has already been initialized."<< G4endl;
00311 G4Exception("G4DNAMolecularMaterial::GetDensityTableFor","G4DNAMolecularMaterial003",
00312 FatalException,exceptionDescription);
00313 }
00314
00315 if(G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle)
00316 const_cast<G4DNAMolecularMaterial*>(this)->Initialize();
00317 else
00318 {
00319 G4ExceptionDescription exceptionDescription;
00320 exceptionDescription << "The geant4 application is at the wrong state. State must be: G4State_Idle."<< G4endl;
00321 G4Exception("G4DNAMolecularMaterial::GetDensityTableFor",
00322 "G4DNAMolecularMaterial_WRONG_STATE_APPLICATION",FatalException,exceptionDescription);
00323 }
00324 }
00325
00326 std::map<const G4Material*,std::vector<double>*,CompareMaterial>::const_iterator it_askedDensityTable = fAskedDensityTable.find(lookForMaterial);
00327 if(it_askedDensityTable != fAskedDensityTable.end())
00328 {
00329 return it_askedDensityTable->second;
00330 }
00331
00332 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00333
00334 std::vector<double>* output = new std::vector<double>(materialTable->size());
00335
00336 ComponentMap::const_iterator it;
00337
00338 G4bool materialWasNotFound = true;
00339
00340 for(int i = 0 ; i < int(materialTable->size()) ; i++)
00341 {
00342 ComponentMap& densityTable = (*fpCompDensityTable)[i];
00343
00344 it = densityTable.find(lookForMaterial);
00345
00346 if(it==densityTable.end())
00347 {
00348 (*output)[i] = 0.0;
00349 }
00350 else
00351 {
00352 materialWasNotFound = false;
00353 (*output)[i] = it->second;
00354 }
00355 }
00356
00357 if(materialWasNotFound)
00358 {
00359 PrintNotAMolecularMaterial("G4DNAMolecularMaterial::GetDensityTableFor",lookForMaterial);
00360 }
00361
00362 fAskedDensityTable.insert(make_pair(lookForMaterial, output));
00363
00364 return output;
00365 }
00366
00367 const std::vector<double>* G4DNAMolecularMaterial::GetNumMolPerVolTableFor(const G4Material* lookForMaterial) const
00368 {
00369 if(!fpCompNumMolPerVolTable)
00370 {
00371 if(fIsInitialized)
00372 {
00373 G4ExceptionDescription exceptionDescription;
00374 exceptionDescription << "The pointer fpCompNumMolPerVolTable is not initialized will the singleton of G4DNAMolecularMaterial "
00375 << "has already been initialized."<< G4endl;
00376 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolTableFor","G4DNAMolecularMaterial005",
00377 FatalException,exceptionDescription);
00378 }
00379
00380 if(G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle)
00381 {
00382 const_cast<G4DNAMolecularMaterial*>(this)->Initialize();
00383 }
00384 else
00385 {
00386 G4ExceptionDescription exceptionDescription;
00387 exceptionDescription << "The geant4 application is at the wrong state. State must be : G4State_Idle."<< G4endl;
00388 G4Exception("G4DNAMolecularMaterial::GetNumMolPerVolTableFor",
00389 "G4DNAMolecularMaterial_WRONG_STATE_APPLICATION",FatalException,exceptionDescription);
00390 }
00391 }
00392
00393 std::map<const G4Material*,std::vector<double>*,CompareMaterial>::const_iterator it_askedNumMolPerVolTable = fAskedNumPerVolTable.find(lookForMaterial);
00394 if(it_askedNumMolPerVolTable != fAskedNumPerVolTable.end())
00395 {
00396 return it_askedNumMolPerVolTable->second;
00397 }
00398
00399 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00400
00401 std::vector<double>* output = new std::vector<double>(materialTable->size());
00402
00403 ComponentMap::const_iterator it;
00404
00405 G4bool materialWasNotFound = true;
00406
00407 for(int i = 0 ; i < int(materialTable->size()) ; i++)
00408 {
00409 ComponentMap& densityTable = (*fpCompNumMolPerVolTable)[i];
00410
00411 it = densityTable.find(lookForMaterial);
00412
00413 if(it==densityTable.end())
00414 {
00415 (*output)[i] = 0.0;
00416 }
00417 else
00418 {
00419 materialWasNotFound = false;
00420 (*output)[i] = it->second;
00421 }
00422 }
00423
00424 if(materialWasNotFound)
00425 {
00426 PrintNotAMolecularMaterial("G4DNAMolecularMaterial::GetNumMolPerVolTableFor",lookForMaterial);
00427 }
00428
00429 fAskedNumPerVolTable.insert(make_pair(lookForMaterial, output));
00430
00431 return output;
00432 }
00433
00434 void G4DNAMolecularMaterial::PrintNotAMolecularMaterial(const char* methodName, const G4Material* lookForMaterial) const
00435 {
00436 std::map<const G4Material*,bool,CompareMaterial>::iterator it = fWarningPrinted.find(lookForMaterial);
00437
00438 if(it == fWarningPrinted.end())
00439 {
00440 G4ExceptionDescription exceptionDescription;
00441 exceptionDescription
00442 << "The material " << lookForMaterial->GetName()
00443 << " is not defined as a molecular material."<< G4endl
00444 << "Meaning: The elements should be added to the material using atom count rather than mass fraction (cf. G4Material)"
00445 << G4endl
00446 << "If you want to use DNA processes on liquid water, you should better use the NistManager to create the water material."
00447 << G4endl
00448 << "Since this message is displayed, it means that the DNA models will not be called."
00449 << "Please note that this message will only appear once even if you are using other methods of G4DNAMolecularMaterial."
00450 << G4endl;
00451
00452 G4Exception(methodName,"MATERIAL_NOT_DEFINE_USING_ATOM_COUNT",JustWarning,exceptionDescription);
00453 fWarningPrinted[lookForMaterial] = true;
00454 }
00455 }