Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // Contact authors: S. Meylan, C. Villagrasa 27 // Contact authors: S. Meylan, C. Villagrasa 28 // email: sylvain.meylan@symalgo-tech.com, car 28 // email: sylvain.meylan@symalgo-tech.com, carmen.villagrasa@irsn.fr 29 // updated : Hoang Tran : 6/1/2023 clean code 29 // updated : Hoang Tran : 6/1/2023 clean code 30 30 31 #include "G4DNAModelInterface.hh" 31 #include "G4DNAModelInterface.hh" 32 32 33 #include "G4DNAMolecularMaterial.hh" 33 #include "G4DNAMolecularMaterial.hh" 34 #include "G4LossTableManager.hh" 34 #include "G4LossTableManager.hh" 35 #include "G4ParticleChangeForGamma.hh" 35 #include "G4ParticleChangeForGamma.hh" 36 #include "G4VDNAModel.hh" 36 #include "G4VDNAModel.hh" 37 #include "G4VEmModel.hh" 37 #include "G4VEmModel.hh" 38 G4DNAModelInterface::G4DNAModelInterface(const 38 G4DNAModelInterface::G4DNAModelInterface(const G4String& nam) : G4VEmModel(nam), fName(nam) {} 39 39 40 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 41 41 42 void G4DNAModelInterface::Initialise(const G4P 42 void G4DNAModelInterface::Initialise(const G4ParticleDefinition* particle, const G4DataVector& cuts) 43 { 43 { 44 // Those two statements are necessary to ove 44 // Those two statements are necessary to override the energy limits set in the G4DNAProcesses 45 // (ionisation, elastic, etc...). Indeed, wi 45 // (ionisation, elastic, etc...). Indeed, with the ModelInterface system, the model define 46 // themselves their energy limits per materi 46 // themselves their energy limits per material and particle. Therefore, such a limit should not be 47 // in the G4DNAProcess classes. 47 // in the G4DNAProcess classes. 48 // 48 // 49 49 50 fpG4_WATER = G4Material::GetMaterial("G4_WAT 50 fpG4_WATER = G4Material::GetMaterial("G4_WATER", false); 51 51 52 SetLowEnergyLimit(0.); 52 SetLowEnergyLimit(0.); 53 SetHighEnergyLimit(DBL_MAX); 53 SetHighEnergyLimit(DBL_MAX); 54 54 55 fpParticleChangeForGamma = GetParticleChange 55 fpParticleChangeForGamma = GetParticleChangeForGamma(); 56 56 57 // Loop on all the registered models to init 57 // Loop on all the registered models to initialise them 58 for (auto & fRegisteredModel : fRegisteredMo 58 for (auto & fRegisteredModel : fRegisteredModels) { 59 fRegisteredModel->SetParticleChange(fpPart 59 fRegisteredModel->SetParticleChange(fpParticleChangeForGamma); 60 fRegisteredModel->Initialise(particle, cut 60 fRegisteredModel->Initialise(particle, cuts); 61 } 61 } 62 // used to retrieve the model corresponding 62 // used to retrieve the model corresponding to the current material/particle couple 63 BuildMaterialParticleModelTable(particle); 63 BuildMaterialParticleModelTable(particle); 64 64 65 BuildMaterialMolPerVolTable(); 65 BuildMaterialMolPerVolTable(); 66 66 67 StreamInfo(G4cout); 67 StreamInfo(G4cout); 68 68 69 } 69 } 70 70 71 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 72 72 73 G4double G4DNAModelInterface::CrossSectionPerV 73 G4double G4DNAModelInterface::CrossSectionPerVolume(const G4Material* material, 74 74 const G4ParticleDefinition* p, G4double ekin, 75 75 G4double emin, G4double emax) 76 { 76 { 77 // Method to return the crossSection * nbMol 77 // Method to return the crossSection * nbMoleculePerUnitVolume to the process class. 78 // Process class then calculates the path. 78 // Process class then calculates the path. 79 // The cross section is calculated in the re 79 // The cross section is calculated in the registered model(s) and this class just call the method 80 // Two cases are handled here: normal materi 80 // Two cases are handled here: normal material and composite material. 81 // 81 // 82 // Idea: 82 // Idea: 83 // *** Simple material *** 83 // *** Simple material *** 84 // Ask for the cross section of the chosen m 84 // Ask for the cross section of the chosen model. 85 // Multiply it by the number of medium molec 85 // Multiply it by the number of medium molecules per volume unit. 86 // Return the value. 86 // Return the value. 87 // *** Composite material *** 87 // *** Composite material *** 88 // Ask for the cross section of the chosen m 88 // Ask for the cross section of the chosen model for each component. 89 // Apply a factor to each cross section and 89 // Apply a factor to each cross section and sum the results. The factor is the molecule number of 90 // component per composite volume unit. The 90 // component per composite volume unit. The total cross section is returned. 91 91 92 // To reset the sampledMat variable. 92 // To reset the sampledMat variable. 93 // Can be used by user to retrieve current c 93 // Can be used by user to retrieve current component 94 fSampledMat = 0; 94 fSampledMat = 0; 95 95 96 // This is the value to be sum up and to be 96 // This is the value to be sum up and to be returned at then end 97 G4double crossSectionTimesNbMolPerVol(0.); 97 G4double crossSectionTimesNbMolPerVol(0.); 98 98 99 // Reset the map saving the material and the 99 // Reset the map saving the material and the cumulated corresponding cross section 100 // Used in SampleSecondaries if the interact 100 // Used in SampleSecondaries if the interaction is selected for the step and if the material is a 101 // composite 101 // composite 102 fMaterialCS.clear(); 102 fMaterialCS.clear(); 103 103 104 // This is the value to be used by SampleSec 104 // This is the value to be used by SampleSecondaries 105 fCSsumTot = 0; 105 fCSsumTot = 0; 106 106 107 // ***************************** 107 // ***************************** 108 // Material is not a composite 108 // Material is not a composite 109 // ***************************** 109 // ***************************** 110 // 110 // 111 if (material->GetMatComponents().empty()) { 111 if (material->GetMatComponents().empty()) { 112 // Get the material name 112 // Get the material name 113 const size_t & materialID = material->GetI 113 const size_t & materialID = material->GetIndex(); 114 114 115 // Use the table to get the model 115 // Use the table to get the model 116 auto model = SelectModel(materialID, p, ek 116 auto model = SelectModel(materialID, p, ekin); 117 117 118 // Get the nunber of molecules per volume 118 // Get the nunber of molecules per volume unit for that material 119 119 120 // Calculate the cross section times the n 120 // Calculate the cross section times the number of molecules 121 if (model != nullptr) { 121 if (model != nullptr) { 122 if (dynamic_cast<G4VDNAModel*>(model) == 122 if (dynamic_cast<G4VDNAModel*>(model) == nullptr) { 123 // water material models only 123 // water material models only 124 crossSectionTimesNbMolPerVol = model-> 124 crossSectionTimesNbMolPerVol = model->CrossSectionPerVolume(material, p, ekin, emin, emax); 125 } 125 } 126 else { 126 else { 127 crossSectionTimesNbMolPerVol = model-> 127 crossSectionTimesNbMolPerVol = model->CrossSectionPerVolume(material, p, ekin, emin, emax); 128 } 128 } 129 } 129 } 130 else // no model was selected, we are out 130 else // no model was selected, we are out of the energy ranges 131 crossSectionTimesNbMolPerVol = 0.; 131 crossSectionTimesNbMolPerVol = 0.; 132 } 132 } 133 133 134 // ******************************** 134 // ******************************** 135 // Material is a composite 135 // Material is a composite 136 // ******************************** 136 // ******************************** 137 // 137 // 138 else { 138 else { 139 // Copy the map in a local variable 139 // Copy the map in a local variable 140 // Otherwise we get segmentation fault and 140 // Otherwise we get segmentation fault and iterator pointing to nowhere: do not know why... 141 // Maybe MatComponents map is overrided by 141 // Maybe MatComponents map is overrided by something somewhere ? 142 auto componentsMap = material->GetMatCompo 142 auto componentsMap = material->GetMatComponents(); 143 143 144 G4cout << material->GetName() << G4endl; 144 G4cout << material->GetName() << G4endl; 145 145 146 // Loop on all the components 146 // Loop on all the components 147 for (const auto& it : componentsMap) { 147 for (const auto& it : componentsMap) { 148 // Get the current component 148 // Get the current component 149 auto component = it.first; 149 auto component = it.first; 150 // Get the current component mass fracti 150 // Get the current component mass fraction 151 // G4double massFraction = it->second; 151 // G4double massFraction = it->second; 152 152 153 // Get the number of component molecules 153 // Get the number of component molecules in a volume unit of composite material 154 G4double nbMoleculeOfComponentInComposit 154 G4double nbMoleculeOfComponentInCompositeMat = 155 GetNumMolPerVolUnitForComponentInCompo 155 GetNumMolPerVolUnitForComponentInComposite(component, material); 156 G4cout << " ==========>component : " << 156 G4cout << " ==========>component : " << component->GetName() 157 << " nbMoleculeOfComponentInCompo 157 << " nbMoleculeOfComponentInCompositeMat: " << nbMoleculeOfComponentInCompositeMat 158 << G4endl; 158 << G4endl; 159 159 160 // Get the current component name 160 // Get the current component name 161 const std::size_t & componentID = compon 161 const std::size_t & componentID = component->GetIndex(); 162 162 163 // Retrieve the model corresponding to t 163 // Retrieve the model corresponding to the current component (ie material) 164 auto model = SelectModel(componentID, p, 164 auto model = SelectModel(componentID, p, ekin); 165 165 166 // Add the component part of the cross s 166 // Add the component part of the cross section to the cross section variable. 167 // The component cross section is multip 167 // The component cross section is multiplied by the total molecule number in the composite 168 // scaled by the mass fraction. 168 // scaled by the mass fraction. 169 G4double crossSection; 169 G4double crossSection; 170 if (model != nullptr) { 170 if (model != nullptr) { 171 if (dynamic_cast<G4VDNAModel*>(model) 171 if (dynamic_cast<G4VDNAModel*>(model) == nullptr) { 172 // water models 172 // water models 173 crossSection = 173 crossSection = 174 model->CrossSectionPerVolume(compo 174 model->CrossSectionPerVolume(component, p, ekin, emin, emax) 175 / GetNumMoleculePerVolumeUnitForMa 175 / GetNumMoleculePerVolumeUnitForMaterial(fpG4_WATER); 176 } 176 } 177 else { 177 else { 178 crossSection = model->CrossSectionPe 178 crossSection = model->CrossSectionPerVolume(component, p, ekin, emin, emax) 179 / GetNumMoleculePerVo 179 / GetNumMoleculePerVolumeUnitForMaterial(component); 180 } 180 } 181 crossSectionTimesNbMolPerVol = nbMolec 181 crossSectionTimesNbMolPerVol = nbMoleculeOfComponentInCompositeMat * crossSection; 182 } 182 } 183 else // no model was selected, we are o 183 else // no model was selected, we are out of the energy ranges 184 { 184 { 185 crossSectionTimesNbMolPerVol = 0.; 185 crossSectionTimesNbMolPerVol = 0.; 186 } 186 } 187 187 188 // Save the component name and its calcu 188 // Save the component name and its calculated crossSectionTimesNbMolPerVol 189 // To be used by sampling secondaries if 189 // To be used by sampling secondaries if the interaction is selected for the step 190 fMaterialCS[componentID] = crossSectionT 190 fMaterialCS[componentID] = crossSectionTimesNbMolPerVol; 191 191 192 // Save the component name and its calcu 192 // Save the component name and its calculated crossSectionTimesNbMolPerVol 193 // To be used by sampling secondaries if 193 // To be used by sampling secondaries if the interaction is selected for the step 194 fCSsumTot += crossSectionTimesNbMolPerVo 194 fCSsumTot += crossSectionTimesNbMolPerVol; 195 } 195 } 196 crossSectionTimesNbMolPerVol = fCSsumTot; 196 crossSectionTimesNbMolPerVol = fCSsumTot; 197 } 197 } 198 198 199 // return the cross section times the number 199 // return the cross section times the number of molecules 200 // the path of the interaction will be calcu 200 // the path of the interaction will be calculated using that value 201 return crossSectionTimesNbMolPerVol; 201 return crossSectionTimesNbMolPerVol; 202 } 202 } 203 203 204 //....oooOO0OOooo........oooOO0OOooo........oo 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 205 205 206 void G4DNAModelInterface::SampleSecondaries(st 206 void G4DNAModelInterface::SampleSecondaries(std::vector<G4DynamicParticle*>* fVect, 207 co 207 const G4MaterialCutsCouple* couple, 208 co 208 const G4DynamicParticle* aDynamicParticle, 209 G4 209 G4double tmin, G4double tmax) 210 { 210 { 211 // To call the sampleSecondaries method of t 211 // To call the sampleSecondaries method of the registered model(s) 212 // In the case of composite material, we nee 212 // In the case of composite material, we need to choose a component to call the method from. 213 // To do so we use a random sampling on the 213 // To do so we use a random sampling on the crossSectionTimesNbMolPerVol used in 214 // CrossSectionPerVolume method. If we enter 214 // CrossSectionPerVolume method. If we enter that method it means the corresponding interaction 215 // (and process) has been chosen for the cur 215 // (and process) has been chosen for the current step. 216 216 217 std::size_t materialID; 217 std::size_t materialID; 218 218 219 // ******************************* 219 // ******************************* 220 // Material is not a composite 220 // Material is not a composite 221 // ******************************* 221 // ******************************* 222 // 222 // 223 if (couple->GetMaterial()->GetMatComponents( 223 if (couple->GetMaterial()->GetMatComponents().empty()) { 224 materialID = couple->GetMaterial()->GetInd 224 materialID = couple->GetMaterial()->GetIndex(); 225 } 225 } 226 226 227 // **************************** 227 // **************************** 228 // Material is a composite 228 // Material is a composite 229 // **************************** 229 // **************************** 230 // 230 // 231 else { 231 else { 232 // Material is a composite 232 // Material is a composite 233 // We need to select a component 233 // We need to select a component 234 234 235 // We select a random number between 0 and 235 // We select a random number between 0 and fCSSumTot 236 G4double rand = G4UniformRand() * fCSsumTo 236 G4double rand = G4UniformRand() * fCSsumTot; 237 237 238 G4double cumulCS(0); 238 G4double cumulCS(0); 239 239 240 G4bool result = false; 240 G4bool result = false; 241 241 242 // We loop on each component cumulated cro 242 // We loop on each component cumulated cross section 243 // 243 // 244 // Retrieve the iterators 244 // Retrieve the iterators 245 auto it = fMaterialCS.begin(); 245 auto it = fMaterialCS.begin(); 246 auto ite = fMaterialCS.end(); 246 auto ite = fMaterialCS.end(); 247 // While this is true we do not have found 247 // While this is true we do not have found our component. 248 while (rand > cumulCS) { 248 while (rand > cumulCS) { 249 // Check if the sampling is ok 249 // Check if the sampling is ok 250 if (it == ite) { 250 if (it == ite) { 251 G4Exception( 251 G4Exception( 252 "G4DNAModelManager::SampleSecondarie 252 "G4DNAModelManager::SampleSecondaries", "em0003", FatalException, 253 "The random component selection has 253 "The random component selection has failed: we ran into the end of the map without " 254 "having a selected component"); 254 "having a selected component"); 255 return; // to make some compilers hap 255 return; // to make some compilers happy 256 } 256 } 257 257 258 // Set the cumulated value for the itera 258 // Set the cumulated value for the iteration 259 cumulCS += it->second; 259 cumulCS += it->second; 260 260 261 // Check if we have reach the material t 261 // Check if we have reach the material to be selected 262 // The DBL_MAX is here to take into acco 262 // The DBL_MAX is here to take into account a return DBL_MAX in CSPerVol for the elastic model 263 // to force elastic sampleSecondaries wh 263 // to force elastic sampleSecondaries where the particle can be killed. 264 // Used when paticle energy is lower tha 264 // Used when paticle energy is lower than limit. 265 if (rand < cumulCS || cumulCS >= DBL_MAX 265 if (rand < cumulCS || cumulCS >= DBL_MAX) { 266 // we have our selected material 266 // we have our selected material 267 materialID = it->first; 267 materialID = it->first; 268 result = true; 268 result = true; 269 break; 269 break; 270 } 270 } 271 271 272 // make the iterator move forward 272 // make the iterator move forward 273 ++it; 273 ++it; 274 } 274 } 275 275 276 // Check that we get a result 276 // Check that we get a result 277 if (!result) { 277 if (!result) { 278 // it is possible to end up here if the 278 // it is possible to end up here if the return DBL_MAX of CSPerVol in the elastic model is not 279 // taken into account 279 // taken into account 280 280 281 G4Exception("G4DNAModelManager::SampleSe 281 G4Exception("G4DNAModelManager::SampleSecondaries", "em0005", FatalException, 282 "The random component select 282 "The random component selection has failed: while loop ended without a selected " 283 "component."); 283 "component."); 284 return; // to make some compilers happy 284 return; // to make some compilers happy 285 } 285 } 286 } 286 } 287 287 288 // ************************************** 288 // ************************************** 289 // Call the SampleSecondaries method 289 // Call the SampleSecondaries method 290 // ************************************** 290 // ************************************** 291 291 292 // Rename material if modified NIST material 292 // Rename material if modified NIST material 293 // This is needed when material is obtained 293 // This is needed when material is obtained from G4MaterialCutsCouple 294 // if (materialName.find("_MODIFIED") != G4 294 // if (materialName.find("_MODIFIED") != G4String::npos) { 295 // materialName = materialName.substr(0, 295 // materialName = materialName.substr(0, materialName.size() - 9); 296 // } 296 // } 297 297 298 fSampledMat = materialID; 298 fSampledMat = materialID; 299 299 300 auto model = SelectModel(materialID, aDynami 300 auto model = SelectModel(materialID, aDynamicParticle->GetParticleDefinition(), 301 aDynamicParticle->G 301 aDynamicParticle->GetKineticEnergy()); 302 302 303 model->SampleSecondaries(fVect, couple, aDyn 303 model->SampleSecondaries(fVect, couple, aDynamicParticle, tmin, tmax); 304 } 304 } 305 305 306 void G4DNAModelInterface::RegisterModel(G4VEmM 306 void G4DNAModelInterface::RegisterModel(G4VEmModel* model) 307 { 307 { 308 fRegisteredModels.push_back(model); 308 fRegisteredModels.push_back(model); 309 } 309 } 310 310 311 //....oooOO0OOooo........oooOO0OOooo........oo 311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 312 312 313 void G4DNAModelInterface::BuildMaterialParticl 313 void G4DNAModelInterface::BuildMaterialParticleModelTable(const G4ParticleDefinition* p) 314 { 314 { 315 // Method to build a map: [material][particl 315 // Method to build a map: [material][particle] = Model*. 316 // The map is used to retrieve the correct m 316 // The map is used to retrieve the correct model for the current particle/material couple. 317 317 318 // Loop on all materials registered in the s 318 // Loop on all materials registered in the simulation 319 for (auto it : *G4Material::GetMaterialTable 319 for (auto it : *G4Material::GetMaterialTable()) { 320 // Get the material pointer 320 // Get the material pointer 321 G4Material* mat = it; 321 G4Material* mat = it; 322 // Get the map 322 // Get the map 323 // Check that the material is not a compos 323 // Check that the material is not a composite material 324 auto componentMap = mat->GetMatComponents( 324 auto componentMap = mat->GetMatComponents(); 325 if (componentMap.empty()) { 325 if (componentMap.empty()) { 326 // Get the material name 326 // Get the material name 327 const std::size_t & matID = mat->GetInde 327 const std::size_t & matID = mat->GetIndex(); 328 InsertModelInTable(matID, p); 328 InsertModelInTable(matID, p); 329 } 329 } 330 // if the material is a composite material 330 // if the material is a composite material then we need to loop on all its components to 331 // register them 331 // register them 332 else { 332 else { 333 // Loop on all the components of the mat 333 // Loop on all the components of the material 334 for (const auto& itComp : componentMap) 334 for (const auto& itComp : componentMap) { 335 G4Material* component = itComp.first; 335 G4Material* component = itComp.first; 336 // Check that the component is not its 336 // Check that the component is not itself a composite 337 if (!component->GetMatComponents().emp 337 if (!component->GetMatComponents().empty()) { 338 std::ostringstream oss; 338 std::ostringstream oss; 339 oss << "Material " << mat->GetName() 339 oss << "Material " << mat->GetName() << " is a composite and its component"; 340 oss << " " << component->GetName(); 340 oss << " " << component->GetName(); 341 G4Exception("G4DNAModelManager::Buil 341 G4Exception("G4DNAModelManager::BuildMaterialParticleModelTable", "em0007", 342 FatalException, oss.str( 342 FatalException, oss.str().c_str()); 343 return; // to make some compilers h 343 return; // to make some compilers happy 344 } 344 } 345 // Get the current component name 345 // Get the current component name 346 const std::size_t & compID = component 346 const std::size_t & compID = component->GetIndex(); 347 // If there is a model then insert the 347 // If there is a model then insert the model corresponding to the component in the table 348 // contains a if statement to check we 348 // contains a if statement to check we have not registered the material as a component or a 349 // normal material before. 349 // normal material before. 350 InsertModelInTable(compID, p); 350 InsertModelInTable(compID, p); 351 // move forward the iterator 351 // move forward the iterator 352 } 352 } 353 } 353 } 354 } 354 } 355 } 355 } 356 356 357 //....oooOO0OOooo........oooOO0OOooo........oo 357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 358 358 359 void G4DNAModelInterface::BuildMaterialMolPerV 359 void G4DNAModelInterface::BuildMaterialMolPerVolTable() 360 { 360 { 361 // To be sure the G4DNAMolecularMaterial is 361 // To be sure the G4DNAMolecularMaterial is initialized 362 G4DNAMolecularMaterial::Instance()->Initiali 362 G4DNAMolecularMaterial::Instance()->Initialize(); 363 363 364 G4MaterialTable* materialTable = G4Material: 364 G4MaterialTable* materialTable = G4Material::GetMaterialTable(); 365 365 366 // Loop on all the materials inside the "mat 366 // Loop on all the materials inside the "materialTable" 367 for (auto currentMaterial : *materialTable) 367 for (auto currentMaterial : *materialTable) { 368 // Current material 368 // Current material 369 // Current material name 369 // Current material name 370 const std::size_t & currentMatID = current 370 const std::size_t & currentMatID = currentMaterial->GetIndex(); 371 371 372 // Will the material be used in this inter 372 // Will the material be used in this interface instance ? 373 // Loop on all the materials that can be d 373 // Loop on all the materials that can be dealt with in this class 374 auto it = fMaterialParticleModelTable.begi 374 auto it = fMaterialParticleModelTable.begin(); 375 auto ite = fMaterialParticleModelTable.end 375 auto ite = fMaterialParticleModelTable.end(); 376 for (; it != ite; it++) { 376 for (; it != ite; it++) { 377 const std::size_t & materialID = it->fir 377 const std::size_t & materialID = it->first; 378 378 379 if (materialID == currentMatID) { 379 if (materialID == currentMatID) { 380 const std::vector<G4double>* numMolPer 380 const std::vector<G4double>* numMolPerVolForMat = 381 G4DNAMolecularMaterial::Instance()-> 381 G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(currentMaterial); 382 fMaterialMolPerVol[materialID] = numMo 382 fMaterialMolPerVol[materialID] = numMolPerVolForMat; 383 } 383 } 384 } 384 } 385 } 385 } 386 } 386 } 387 387 388 //....oooOO0OOooo........oooOO0OOooo........oo 388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 389 389 390 void G4DNAModelInterface::InsertModelInTable(c 390 void G4DNAModelInterface::InsertModelInTable(const std::size_t& matID, const G4ParticleDefinition* p) 391 { 391 { 392 // To insert the model(s) in the table Mater 392 // To insert the model(s) in the table Material Particule -> Model(s) 393 393 394 // First, we need to check if the current ma 394 // First, we need to check if the current material has already been inserted in the table. 395 // This is possible because of the composite 395 // This is possible because of the composite material. We could add a component M1 and then try to 396 // add the independant M1 material. This cas 396 // add the independant M1 material. This case must be avoided. Checking if M1 is already in the 397 // table is the way to avoid it. 397 // table is the way to avoid it. 398 // 398 // 399 // Check if the current material and particl 399 // Check if the current material and particle are already in the table. 400 // If they are: do nothing. 400 // If they are: do nothing. 401 // If they are not: add the model(s) 401 // If they are not: add the model(s) 402 // 402 // 403 // Check for the material 403 // Check for the material 404 if (fMaterialParticleModelTable.find(matID) 404 if (fMaterialParticleModelTable.find(matID) == fMaterialParticleModelTable.end()) { 405 // Check for the particle 405 // Check for the particle 406 if (fMaterialParticleModelTable[matID].fin 406 if (fMaterialParticleModelTable[matID].find(p) == fMaterialParticleModelTable[matID].end()) { 407 G4int modelNbForMaterial = 0; 407 G4int modelNbForMaterial = 0; 408 for (const auto& it : fRegisteredModels) 408 for (const auto& it : fRegisteredModels) { 409 auto model = dynamic_cast<G4VDNAModel* 409 auto model = dynamic_cast<G4VDNAModel*>(it); 410 if (model != nullptr) { 410 if (model != nullptr) { 411 if (model->IsParticleExistingInModel 411 if (model->IsParticleExistingInModelForMaterial(p, matID)) { 412 fMaterialParticleModelTable[matID] 412 fMaterialParticleModelTable[matID][p] = it; 413 // and add one to the "there is a 413 // and add one to the "there is a model" material flag 414 ++modelNbForMaterial; 414 ++modelNbForMaterial; 415 } 415 } 416 } 416 } 417 else { 417 else { 418 auto index = fpG4_WATER->GetIndex(); 418 auto index = fpG4_WATER->GetIndex(); 419 fMaterialParticleModelTable[index][p 419 fMaterialParticleModelTable[index][p] = it; 420 ++modelNbForMaterial; 420 ++modelNbForMaterial; 421 } 421 } 422 } 422 } 423 if (modelNbForMaterial == 0) { 423 if (modelNbForMaterial == 0) { 424 std::ostringstream oss; 424 std::ostringstream oss; 425 oss << "The material " << (*G4Material 425 oss << "The material " << (*G4Material::GetMaterialTable())[matID]->GetName() 426 << " and the particle " << p->GetP 426 << " and the particle " << p->GetParticleName(); 427 oss << " does not have any model regis 427 oss << " does not have any model registered for the " << fName << " interaction."; 428 G4Exception("G4DNAModelInterface::Inse 428 G4Exception("G4DNAModelInterface::InsertModelInTable", "em0006", FatalException, 429 oss.str().c_str()); 429 oss.str().c_str()); 430 return; // to make some compilers hap 430 return; // to make some compilers happy 431 } 431 } 432 } 432 } 433 } 433 } 434 } 434 } 435 435 436 G4VEmModel* G4DNAModelInterface::SelectModel(c 436 G4VEmModel* G4DNAModelInterface::SelectModel(const std::size_t& materialID, 437 c 437 const G4ParticleDefinition* particle, 438 c 438 const G4double& ekin) 439 { 439 { 440 // Output pointer 440 // Output pointer 441 G4VEmModel* model = nullptr; 441 G4VEmModel* model = nullptr; 442 442 443 // Get a reference to all the models for the 443 // Get a reference to all the models for the couple (material and particle) 444 auto modelData = fMaterialParticleModelTable 444 auto modelData = fMaterialParticleModelTable[materialID][particle]; 445 445 446 // We must choose one of the model(s) accord 446 // We must choose one of the model(s) accordingly to the particle energy and the model energy 447 // range(s) 447 // range(s) 448 448 449 // Loop on all the models within the models 449 // Loop on all the models within the models vector and check if ekin is within the energy range. 450 auto DNAModel = dynamic_cast<G4VDNAModel*>(m 450 auto DNAModel = dynamic_cast<G4VDNAModel*>(modelData); 451 G4double lowL, highL; 451 G4double lowL, highL; 452 if (DNAModel == nullptr) { 452 if (DNAModel == nullptr) { 453 // ekin is in the energy range: we select 453 // ekin is in the energy range: we select the model and stop the loop. 454 lowL = modelData->LowEnergyLimit(); 454 lowL = modelData->LowEnergyLimit(); 455 highL = modelData->HighEnergyLimit(); 455 highL = modelData->HighEnergyLimit(); 456 if (ekin >= lowL && ekin < highL) { 456 if (ekin >= lowL && ekin < highL) { 457 // Select the model 457 // Select the model 458 458 459 model = modelData; 459 model = modelData; 460 // return model; 460 // return model; 461 // Quit the for loop 461 // Quit the for loop 462 // break; 462 // break; 463 } 463 } 464 // ekin is not in the energy range: we con 464 // ekin is not in the energy range: we continue the loop. 465 } 465 } 466 else { 466 else { 467 // ekin is in the energy range: we select 467 // ekin is in the energy range: we select the model and stop the loop. 468 lowL = DNAModel->GetLowELimit(materialID, 468 lowL = DNAModel->GetLowELimit(materialID, particle); 469 highL = DNAModel->GetHighELimit(materialID 469 highL = DNAModel->GetHighELimit(materialID, particle); 470 if (ekin >= lowL && ekin < highL) { 470 if (ekin >= lowL && ekin < highL) { 471 // Select the model 471 // Select the model 472 model = modelData; 472 model = modelData; 473 // return model; 473 // return model; 474 // Quit the for loop 474 // Quit the for loop 475 // break; 475 // break; 476 } 476 } 477 // ekin is not in the energy range: we con 477 // ekin is not in the energy range: we continue the loop. 478 } 478 } 479 //} 479 //} 480 return model; 480 return model; 481 } 481 } 482 482 483 //....oooOO0OOooo........oooOO0OOooo........oo 483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 484 484 485 G4double G4DNAModelInterface::GetNumMoleculePe 485 G4double G4DNAModelInterface::GetNumMoleculePerVolumeUnitForMaterial(const G4Material* mat) 486 { 486 { 487 return fMaterialMolPerVol[mat->GetIndex()]-> 487 return fMaterialMolPerVol[mat->GetIndex()]->at(mat->GetIndex()); 488 } 488 } 489 489 490 //....oooOO0OOooo........oooOO0OOooo........oo 490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 491 491 492 G4double 492 G4double 493 G4DNAModelInterface::GetNumMolPerVolUnitForCom 493 G4DNAModelInterface::GetNumMolPerVolUnitForComponentInComposite(const G4Material* component, 494 494 const G4Material* composite) 495 { 495 { 496 return fMaterialMolPerVol[component->GetInde 496 return fMaterialMolPerVol[component->GetIndex()]->at(composite->GetIndex()); 497 } 497 } 498 498 499 void G4DNAModelInterface::StreamInfo(std::ostr 499 void G4DNAModelInterface::StreamInfo(std::ostream& os) const 500 { 500 { 501 G4long prec = os.precision(5); 501 G4long prec = os.precision(5); 502 os << "===================================== 502 os << "======================================= Materials of " << std::setw(17) 503 << this->GetName() << " ============ 503 << this->GetName() << " ================================================" 504 << "\n"; 504 << "\n"; 505 os << std::setw(15) << "Material#" << std::s 505 os << std::setw(15) << "Material#" << std::setw(13) << "Particle" << std::setw(35) << "Model" 506 << std::setw(17) << "LowLimit(MeV)" << st 506 << std::setw(17) << "LowLimit(MeV)" << std::setw(17) << "HighLimit(MeV)" << std::setw(13) 507 << "Fast" << std::setw(13) << "Stationary 507 << "Fast" << std::setw(13) << "Stationary" << std::setw(13) << "Chemistry" << G4endl; 508 for (const auto& it1 : fMaterialParticleMode 508 for (const auto& it1 : fMaterialParticleModelTable) { 509 os << std::setw(15) << (*G4Material::GetMa 509 os << std::setw(15) << (*G4Material::GetMaterialTable())[it1.first]->GetName(); 510 for (const auto& it2 : it1.second) { 510 for (const auto& it2 : it1.second) { 511 os << std::setw(13) << it2.first->GetPar 511 os << std::setw(13) << it2.first->GetParticleName(); 512 os << std::setw(35) << it2.second->GetNa 512 os << std::setw(35) << it2.second->GetName(); 513 auto DNAModel = dynamic_cast<G4VDNAModel 513 auto DNAModel = dynamic_cast<G4VDNAModel*>(it2.second); 514 if (DNAModel == nullptr) { 514 if (DNAModel == nullptr) { 515 os << std::setw(17) << it2.second->Low 515 os << std::setw(17) << it2.second->LowEnergyLimit(); 516 os << std::setw(17) << it2.second->Hig 516 os << std::setw(17) << it2.second->HighEnergyLimit(); 517 } 517 } 518 else { 518 else { 519 auto lowL = DNAModel->GetLowELimit(it1 519 auto lowL = DNAModel->GetLowELimit(it1.first, it2.first); 520 auto highL = DNAModel->GetHighELimit(i 520 auto highL = DNAModel->GetHighELimit(it1.first, it2.first); 521 os << std::setw(17) << lowL; 521 os << std::setw(17) << lowL; 522 os << std::setw(17) << highL; 522 os << std::setw(17) << highL; 523 } 523 } 524 os << std::setw(13) << "no"; 524 os << std::setw(13) << "no"; 525 os << std::setw(13) << "no"; 525 os << std::setw(13) << "no"; 526 os << std::setw(13) << "no" << G4endl; 526 os << std::setw(13) << "no" << G4endl; 527 } 527 } 528 } 528 } 529 529 530 os << "===================================== 530 os << "========================================================================================" 531 "===================================== 531 "==================================================" 532 << G4endl; 532 << G4endl; 533 os.precision(prec); 533 os.precision(prec); 534 } 534 } 535 535