| 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 //
28 // email: sylvain.meylan@symalgo-tech.com, car 29 // email: sylvain.meylan@symalgo-tech.com, carmen.villagrasa@irsn.fr
29 // updated : Hoang Tran : 6/1/2023 clean code <<
30 30
31 #include "G4DNAModelInterface.hh" 31 #include "G4DNAModelInterface.hh"
32 << 32 #include "G4PhysicalConstants.hh"
>> 33 #include "G4SystemOfUnits.hh"
33 #include "G4DNAMolecularMaterial.hh" 34 #include "G4DNAMolecularMaterial.hh"
34 #include "G4LossTableManager.hh" << 35
35 #include "G4ParticleChangeForGamma.hh" << 36
36 #include "G4VDNAModel.hh" << 37 G4DNAModelInterface::G4DNAModelInterface(const G4String &nam)
37 #include "G4VEmModel.hh" << 38 : G4VEmModel(nam), fName(nam), fpParticleChangeForGamma(0), fSampledMat("")
38 G4DNAModelInterface::G4DNAModelInterface(const << 39 {
>> 40
>> 41 }
39 42
40 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 44
42 void G4DNAModelInterface::Initialise(const G4P << 45 G4DNAModelInterface::~G4DNAModelInterface()
43 { 46 {
44 // Those two statements are necessary to ove << 47 // Loop on all the registered models to properly delete them (free the memory)
45 // (ionisation, elastic, etc...). Indeed, wi << 48 for(unsigned int i=0, ie = fRegisteredModels.size(); i<ie; ++i)
46 // themselves their energy limits per materi << 49 {
47 // in the G4DNAProcess classes. << 50 if(fRegisteredModels.at(i) != nullptr) delete fRegisteredModels.at(i);
48 // << 51 }
>> 52 }
49 53
50 fpG4_WATER = G4Material::GetMaterial("G4_WAT << 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 55
52 SetLowEnergyLimit(0.); << 56 void G4DNAModelInterface::Initialise(const G4ParticleDefinition* particle,
53 SetHighEnergyLimit(DBL_MAX); << 57 const G4DataVector& cuts)
>> 58 {
>> 59 // Those two statements are necessary to override the energy limits set in the G4DNAProcesses (ionisation, elastic, etc...).
>> 60 // Indeed, with the ModelInterface system, the model define themselves their energy limits per material and particle.
>> 61 // Therefore, such a limit should not be in the G4DNAProcess classes.
>> 62 //
>> 63 SetLowEnergyLimit(0.);
>> 64 SetHighEnergyLimit(DBL_MAX);
54 65
55 fpParticleChangeForGamma = GetParticleChange << 66 fpParticleChangeForGamma = GetParticleChangeForGamma();
56 67
57 // Loop on all the registered models to init << 68 // Loop on all the registered models to initialise them
58 for (auto & fRegisteredModel : fRegisteredMo << 69 for(unsigned int i=0, ie = fRegisteredModels.size(); i<ie; ++i)
59 fRegisteredModel->SetParticleChange(fpPart << 70 {
60 fRegisteredModel->Initialise(particle, cut << 71 fRegisteredModels.at(i)->Initialise(particle, cuts, fpParticleChangeForGamma);
61 } << 72 }
62 // used to retrieve the model corresponding <<
63 BuildMaterialParticleModelTable(particle); <<
64 73
65 BuildMaterialMolPerVolTable(); <<
66 74
67 StreamInfo(G4cout); << 75 // Build the [material][particle]=Models table
>> 76 // used to retrieve the model corresponding to the current material/particle couple
>> 77 BuildMaterialParticleModelTable(particle);
68 78
>> 79 BuildMaterialMolPerVolTable();
69 } 80 }
70 81
71 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 83
73 G4double G4DNAModelInterface::CrossSectionPerV 84 G4double G4DNAModelInterface::CrossSectionPerVolume(const G4Material* material,
74 << 85 const G4ParticleDefinition* p,
75 << 86 G4double ekin,
76 { << 87 G4double emin,
77 // Method to return the crossSection * nbMol << 88 G4double emax)
78 // Process class then calculates the path. << 89 {
79 // The cross section is calculated in the re << 90 // Method to return the crossSection * nbMoleculePerUnitVolume to the process class.
80 // Two cases are handled here: normal materi << 91 // Process class then calculates the path.
81 // << 92 // The cross section is calculated in the registered model(s) and this class just call the method
82 // Idea: << 93 // Two cases are handled here: normal material and composite material.
83 // *** Simple material *** << 94 //
84 // Ask for the cross section of the chosen m << 95 // Idea:
85 // Multiply it by the number of medium molec << 96 // *** Simple material ***
86 // Return the value. << 97 // Ask for the cross section of the chosen model.
87 // *** Composite material *** << 98 // Multiply it by the number of medium molecules per volume unit.
88 // Ask for the cross section of the chosen m << 99 // Return the value.
89 // Apply a factor to each cross section and << 100 // *** Composite material ***
90 // component per composite volume unit. The << 101 // Ask for the cross section of the chosen model for each component.
91 << 102 // Apply a factor to each cross section and sum the results. The factor is the molecule number of component per composite volume unit.
92 // To reset the sampledMat variable. << 103 // The total cross section is returned.
93 // Can be used by user to retrieve current c << 104
94 fSampledMat = 0; << 105 // To reset the sampledMat variable.
95 << 106 // Can be used by user to retrieve current component
96 // This is the value to be sum up and to be << 107 fSampledMat = "";
97 G4double crossSectionTimesNbMolPerVol(0.); << 108
98 << 109 // This is the value to be sum up and to be returned at then end
99 // Reset the map saving the material and the << 110 G4double crossSectionTimesNbMolPerVol (0);
100 // Used in SampleSecondaries if the interact << 111
101 // composite << 112 // Reset the map saving the material and the cumulated corresponding cross section
102 fMaterialCS.clear(); << 113 // Used in SampleSecondaries if the interaction is selected for the step and if the material is a composite
103 << 114 fMaterialCS.clear();
104 // This is the value to be used by SampleSec << 115
105 fCSsumTot = 0; << 116 // This is the value to be used by SampleSecondaries
106 << 117 fCSsumTot = 0;
107 // ***************************** << 118
108 // Material is not a composite << 119 // *****************************
109 // ***************************** << 120 // Material is not a composite
110 // << 121 // *****************************
111 if (material->GetMatComponents().empty()) { << 122 //
112 // Get the material name << 123 if(material->GetMatComponents().empty())
113 const size_t & materialID = material->GetI << 124 {
114 << 125 // Get the material name
115 // Use the table to get the model << 126 const G4String& materialName = material->GetName();
116 auto model = SelectModel(materialID, p, ek << 127
117 << 128 // Use the table to get the model
118 // Get the nunber of molecules per volume << 129 G4VDNAModel* model = GetDNAModel(materialName, p->GetParticleName(), ekin);
119 << 130
120 // Calculate the cross section times the n << 131 // Get the nunber of molecules per volume unit for that material
121 if (model != nullptr) { << 132 G4double nbOfMoleculePerVolumeUnit = GetNumMoleculePerVolumeUnitForMaterial(material);
122 if (dynamic_cast<G4VDNAModel*>(model) == << 133
123 // water material models only << 134 // Calculate the cross section times the number of molecules
124 crossSectionTimesNbMolPerVol = model-> << 135 if(model != 0)
125 } << 136 crossSectionTimesNbMolPerVol = nbOfMoleculePerVolumeUnit * model->CrossSectionPerVolume(material, materialName, p, ekin, emin, emax);
126 else { << 137 else // no model was selected, we are out of the energy ranges
127 crossSectionTimesNbMolPerVol = model-> << 138 crossSectionTimesNbMolPerVol = 0.;
128 } << 139 }
129 } << 140
130 else // no model was selected, we are out << 141 // ********************************
131 crossSectionTimesNbMolPerVol = 0.; << 142 // Material is a composite
132 } << 143 // ********************************
133 << 144 //
134 // ******************************** << 145 else
135 // Material is a composite << 146 {
136 // ******************************** << 147 // Copy the map in a local variable
137 // << 148 // Otherwise we get segmentation fault and iterator pointing to nowhere: do not know why...
138 else { << 149 // Maybe MatComponents map is overrided by something somewhere ?
139 // Copy the map in a local variable << 150 std::map<G4Material*, G4double> componentsMap = material->GetMatComponents();
140 // Otherwise we get segmentation fault and << 151
141 // Maybe MatComponents map is overrided by << 152 // Retrieve the iterator
142 auto componentsMap = material->GetMatCompo << 153 std::map<G4Material*, G4double>::const_iterator it = componentsMap.begin();
143 << 154
144 G4cout << material->GetName() << G4endl; << 155 // Get the size
145 << 156 unsigned int componentNumber = componentsMap.size();
146 // Loop on all the components << 157
147 for (const auto& it : componentsMap) { << 158 // Loop on all the components
148 // Get the current component << 159 //for(it = material->GetMatComponents().begin(); it!=material->GetMatComponents().end();++it)
149 auto component = it.first; << 160 for(unsigned int i=0; i<componentNumber; ++i)
150 // Get the current component mass fracti << 161 {
151 // G4double massFraction = it->second; << 162 // Get the current component
152 << 163 G4Material* component = it->first;
153 // Get the number of component molecules << 164
154 G4double nbMoleculeOfComponentInComposit << 165 // Get the current component mass fraction
155 GetNumMolPerVolUnitForComponentInCompo << 166 //G4double massFraction = it->second;
156 G4cout << " ==========>component : " << << 167
157 << " nbMoleculeOfComponentInCompo << 168 // Get the number of component molecules in a volume unit of composite material
158 << G4endl; << 169 G4double nbMoleculeOfComponentInCompositeMat = GetNumMolPerVolUnitForComponentInComposite(component, material);
159 << 170
160 // Get the current component name << 171 // Get the current component name
161 const std::size_t & componentID = compon << 172 const G4String componentName = component->GetName();
162 << 173
163 // Retrieve the model corresponding to t << 174 // Retrieve the model corresponding to the current component (ie material)
164 auto model = SelectModel(componentID, p, << 175 G4VDNAModel* model = GetDNAModel(componentName, p->GetParticleName(), ekin);
165 << 176
166 // Add the component part of the cross s << 177 // Add the component part of the cross section to the cross section variable.
167 // The component cross section is multip << 178 // The component cross section is multiplied by the total molecule number in the composite scaled by the mass fraction.
168 // scaled by the mass fraction. << 179 if(model != 0)
169 G4double crossSection; << 180 crossSectionTimesNbMolPerVol =
170 if (model != nullptr) { << 181 nbMoleculeOfComponentInCompositeMat * model->CrossSectionPerVolume(component, componentName, p, ekin, emin, emax);
171 if (dynamic_cast<G4VDNAModel*>(model) << 182 else // no model was selected, we are out of the energy ranges
172 // water models << 183 crossSectionTimesNbMolPerVol = 0.;
173 crossSection = << 184
174 model->CrossSectionPerVolume(compo << 185 // Save the component name and its calculated crossSectionTimesNbMolPerVol
175 / GetNumMoleculePerVolumeUnitForMa << 186 // To be used by sampling secondaries if the interaction is selected for the step
176 } << 187 fMaterialCS[componentName] = crossSectionTimesNbMolPerVol;
177 else { << 188
178 crossSection = model->CrossSectionPe << 189 // Save the component name and its calculated crossSectionTimesNbMolPerVol
179 / GetNumMoleculePerVo << 190 // To be used by sampling secondaries if the interaction is selected for the step
>> 191 fCSsumTot += crossSectionTimesNbMolPerVol;
>> 192
>> 193 // Move forward the iterator
>> 194 ++it;
180 } 195 }
181 crossSectionTimesNbMolPerVol = nbMolec << 196
182 } << 197 crossSectionTimesNbMolPerVol = fCSsumTot;
183 else // no model was selected, we are o << 198
184 { << 199 }
185 crossSectionTimesNbMolPerVol = 0.; << 200
186 } << 201 // return the cross section times the number of molecules
187 << 202 // the path of the interaction will be calculated using that value
188 // Save the component name and its calcu << 203 return crossSectionTimesNbMolPerVol;
189 // To be used by sampling secondaries if <<
190 fMaterialCS[componentID] = crossSectionT <<
191 <<
192 // Save the component name and its calcu <<
193 // To be used by sampling secondaries if <<
194 fCSsumTot += crossSectionTimesNbMolPerVo <<
195 } <<
196 crossSectionTimesNbMolPerVol = fCSsumTot; <<
197 } <<
198 <<
199 // return the cross section times the number <<
200 // the path of the interaction will be calcu <<
201 return crossSectionTimesNbMolPerVol; <<
202 } 204 }
203 205
204 //....oooOO0OOooo........oooOO0OOooo........oo 206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
205 207
206 void G4DNAModelInterface::SampleSecondaries(st 208 void G4DNAModelInterface::SampleSecondaries(std::vector<G4DynamicParticle*>* fVect,
207 co << 209 const G4MaterialCutsCouple* couple,
208 co << 210 const G4DynamicParticle* aDynamicParticle,
209 G4 << 211 G4double tmin,
210 { << 212 G4double tmax)
211 // To call the sampleSecondaries method of t << 213 {
212 // In the case of composite material, we nee << 214 // To call the sampleSecondaries method of the registered model(s)
213 // To do so we use a random sampling on the << 215 // In the case of composite material, we need to choose a component to call the method from.
214 // CrossSectionPerVolume method. If we enter << 216 // To do so we use a random sampling on the crossSectionTimesNbMolPerVol used in CrossSectionPerVolume method.
215 // (and process) has been chosen for the cur << 217 // If we enter that method it means the corresponding interaction (and process) has been chosen for the current step.
216 << 218
217 std::size_t materialID; << 219 G4String materialName;
218 << 220
219 // ******************************* << 221 // *******************************
220 // Material is not a composite << 222 // Material is not a composite
221 // ******************************* << 223 // *******************************
222 // << 224 //
223 if (couple->GetMaterial()->GetMatComponents( << 225 if(couple->GetMaterial()->GetMatComponents().empty())
224 materialID = couple->GetMaterial()->GetInd << 226 {
225 } << 227 materialName = couple->GetMaterial()->GetName();
226 << 228 }
227 // **************************** << 229
228 // Material is a composite << 230 // ****************************
229 // **************************** <<
230 // <<
231 else { <<
232 // Material is a composite 231 // Material is a composite
233 // We need to select a component << 232 // ****************************
>> 233 //
>> 234 else
>> 235 {
>> 236 // Material is a composite
>> 237 // We need to select a component
>> 238
>> 239 // We select a random number between 0 and fCSSumTot
>> 240 G4double rand = G4UniformRand()*fCSsumTot;
>> 241
>> 242 G4double cumulCS (0);
>> 243
>> 244 G4bool result = false;
>> 245
>> 246 // We loop on each component cumulated cross section
>> 247 //
>> 248 // Retrieve the iterators
>> 249 std::map<const G4String , G4double>::const_iterator it = fMaterialCS.begin();
>> 250 std::map<const G4String , G4double>::const_iterator ite = fMaterialCS.end();
>> 251 // While this is true we do not have found our component.
>> 252 while(rand>cumulCS)
>> 253 {
>> 254 // Check if the sampling is ok
>> 255 if(it==ite)
>> 256 {
>> 257 G4Exception("G4DNAModelManager::SampleSecondaries","em0006",
>> 258 FatalException,
>> 259 "The random component selection has failed: we ran into the end of the map without having a selected component");
>> 260 return; // to make some compilers happy
>> 261 }
>> 262
>> 263 // Set the cumulated value for the iteration
>> 264 cumulCS += it->second;
>> 265
>> 266 // Check if we have reach the material to be selected
>> 267 // The DBL_MAX is here to take into account a return DBL_MAX in CSPerVol for the elastic model
>> 268 // to force elastic sampleSecondaries where the particle can be killed.
>> 269 // Used when paticle energy is lower than limit.
>> 270 if(rand<cumulCS || cumulCS >= DBL_MAX)
>> 271 {
>> 272 // we have our selected material
>> 273 materialName = it->first;
>> 274 result = true;
>> 275 break;
>> 276 }
>> 277
>> 278 // make the iterator move forward
>> 279 ++it;
>> 280 }
234 281
235 // We select a random number between 0 and << 282 // Check that we get a result
236 G4double rand = G4UniformRand() * fCSsumTo << 283 if(!result)
>> 284 {
>> 285 // it is possible to end up here if the return DBL_MAX of CSPerVol in the elastic model is not taken into account
>> 286
>> 287 G4Exception("G4DNAModelManager::SampleSecondaries","em0006",
>> 288 FatalException,
>> 289 "The random component selection has failed: while loop ended without a selected component.");
>> 290 return; // to make some compilers happy
>> 291 }
237 292
238 G4double cumulCS(0); << 293 }
239 294
240 G4bool result = false; << 295 // **************************************
>> 296 // Call the SampleSecondaries method
>> 297 // **************************************
241 298
242 // We loop on each component cumulated cro << 299 // Rename material if modified NIST material
243 // << 300 // This is needed when material is obtained from G4MaterialCutsCouple
244 // Retrieve the iterators << 301 if(materialName.find("_MODIFIED")!=G4String::npos)
245 auto it = fMaterialCS.begin(); << 302 {
246 auto ite = fMaterialCS.end(); << 303 materialName = materialName.substr(0,materialName.size()-9);
247 // While this is true we do not have found << 304 }
248 while (rand > cumulCS) { <<
249 // Check if the sampling is ok <<
250 if (it == ite) { <<
251 G4Exception( <<
252 "G4DNAModelManager::SampleSecondarie <<
253 "The random component selection has <<
254 "having a selected component"); <<
255 return; // to make some compilers hap <<
256 } <<
257 <<
258 // Set the cumulated value for the itera <<
259 cumulCS += it->second; <<
260 <<
261 // Check if we have reach the material t <<
262 // The DBL_MAX is here to take into acco <<
263 // to force elastic sampleSecondaries wh <<
264 // Used when paticle energy is lower tha <<
265 if (rand < cumulCS || cumulCS >= DBL_MAX <<
266 // we have our selected material <<
267 materialID = it->first; <<
268 result = true; <<
269 break; <<
270 } <<
271 <<
272 // make the iterator move forward <<
273 ++it; <<
274 } <<
275 <<
276 // Check that we get a result <<
277 if (!result) { <<
278 // it is possible to end up here if the <<
279 // taken into account <<
280 <<
281 G4Exception("G4DNAModelManager::SampleSe <<
282 "The random component select <<
283 "component."); <<
284 return; // to make some compilers happy <<
285 } <<
286 } <<
287 <<
288 // ************************************** <<
289 // Call the SampleSecondaries method <<
290 // ************************************** <<
291 <<
292 // Rename material if modified NIST material <<
293 // This is needed when material is obtained <<
294 // if (materialName.find("_MODIFIED") != G4 <<
295 // materialName = materialName.substr(0, <<
296 // } <<
297 <<
298 fSampledMat = materialID; <<
299 <<
300 auto model = SelectModel(materialID, aDynami <<
301 aDynamicParticle->G <<
302 305
303 model->SampleSecondaries(fVect, couple, aDyn << 306 fSampledMat = materialName;
>> 307
>> 308 G4VDNAModel* model = GetDNAModel(materialName,
>> 309 aDynamicParticle->GetParticleDefinition()->GetParticleName(),
>> 310 aDynamicParticle->GetKineticEnergy() );
>> 311 //fMaterialParticleModelTable[materialName][aDynamicParticle->GetDefinition()->GetParticleName()][0];
>> 312
>> 313 model->SampleSecondaries(fVect, couple, materialName, aDynamicParticle, fpParticleChangeForGamma, tmin, tmax);
304 } 314 }
305 315
306 void G4DNAModelInterface::RegisterModel(G4VEmM << 316 void G4DNAModelInterface::RegisterModel(G4VDNAModel* model)
307 { 317 {
308 fRegisteredModels.push_back(model); << 318 fRegisteredModels.push_back(model);
>> 319 }
>> 320
>> 321 void G4DNAModelInterface::RegisterModel(G4VEmModel* model, const G4ParticleDefinition* particle)
>> 322 {
>> 323 G4DNADummyModel* dummyWrapper = new G4DNADummyModel("G4_WATER", particle, model->GetName(), model);
>> 324
>> 325 RegisterModel(dummyWrapper);
309 } 326 }
310 327
311 //....oooOO0OOooo........oooOO0OOooo........oo 328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
312 329
313 void G4DNAModelInterface::BuildMaterialParticl 330 void G4DNAModelInterface::BuildMaterialParticleModelTable(const G4ParticleDefinition* p)
314 { 331 {
315 // Method to build a map: [material][particl << 332 // Method to build a map: [material][particle] = Model*.
316 // The map is used to retrieve the correct m << 333 // The map is used to retrieve the correct model for the current particle/material couple.
>> 334
>> 335 // Get the current particle name
>> 336 const G4String& pName = p->GetParticleName();
>> 337
>> 338 // Retrieve the iterator
>> 339 G4MaterialTable::iterator it;
317 340
318 // Loop on all materials registered in the s << 341 // Loop on all materials registered in the simulation
319 for (auto it : *G4Material::GetMaterialTable << 342 for(it = G4Material::GetMaterialTable()->begin(); it!=G4Material::GetMaterialTable()->end(); ++it)
320 // Get the material pointer << 343 {
321 G4Material* mat = it; << 344 // Get the material pointer
322 // Get the map << 345 G4Material* mat = *it;
323 // Check that the material is not a compos << 346
324 auto componentMap = mat->GetMatComponents( << 347 // Get the map
325 if (componentMap.empty()) { << 348 std::map<G4Material*, G4double> componentMap = mat->GetMatComponents();
326 // Get the material name << 349
327 const std::size_t & matID = mat->GetInde << 350 // Get the number of component within the composite
328 InsertModelInTable(matID, p); << 351 unsigned int compositeSize = componentMap.size();
329 } << 352
330 // if the material is a composite material << 353 // Check that the material is not a composite material
331 // register them << 354 if(componentMap.empty())
332 else { << 355 {
333 // Loop on all the components of the mat << 356 // Get the material name
334 for (const auto& itComp : componentMap) << 357 const G4String& matName = mat->GetName();
335 G4Material* component = itComp.first; << 358
336 // Check that the component is not its << 359 // Insert the model in the table.
337 if (!component->GetMatComponents().emp << 360 InsertModelInTable(matName, pName);
338 std::ostringstream oss; << 361 }
339 oss << "Material " << mat->GetName() << 362 // if the material is a composite material then we need to loop on all its components to register them
340 oss << " " << component->GetName(); << 363 else
341 G4Exception("G4DNAModelManager::Buil << 364 {
342 FatalException, oss.str( << 365 // Retrieve the component map begin iterator
343 return; // to make some compilers h << 366 std::map<G4Material*, G4double>::const_iterator itComp = componentMap.begin();
>> 367
>> 368 // Loop on all the components of the material
>> 369 //for(itComp = mat->GetMatComponents().begin(); itComp != eitComp; ++itComp)
>> 370 for(unsigned int k=0; k<compositeSize; ++k)
>> 371 {
>> 372 G4Material* component = itComp->first;
>> 373
>> 374 // // Check that the component is not itself a composite
>> 375 // if(component->GetMatComponents().size()!=0)
>> 376 // {
>> 377 // std::ostringstream oss;
>> 378 // oss<<"Material "<<mat->GetName()<<" is a composite and its component ";
>> 379 // oss<<component->GetName()<<" is also a composite material. Building composite with other composites is not implemented yet";
>> 380 // oss<<G4endl;
>> 381 // G4Exception("G4DNAModelManager::BuildMaterialParticleModelTable","em0006",
>> 382 // FatalException, oss.str().c_str());
>> 383 // return; // to make some compilers happy
>> 384 // }
>> 385
>> 386 // Get the current component name
>> 387 const G4String compName = component->GetName();
>> 388
>> 389 // If there is a model then insert the model corresponding to the component in the table
>> 390 // contains a if statement to check we have not registered the material as a component or a normal material before.
>> 391 InsertModelInTable(compName, pName);
>> 392
>> 393 // move forward the iterator
>> 394 ++itComp;
>> 395 }
344 } 396 }
345 // Get the current component name <<
346 const std::size_t & compID = component <<
347 // If there is a model then insert the <<
348 // contains a if statement to check we <<
349 // normal material before. <<
350 InsertModelInTable(compID, p); <<
351 // move forward the iterator <<
352 } <<
353 } 397 }
354 } <<
355 } 398 }
356 399
357 //....oooOO0OOooo........oooOO0OOooo........oo 400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
358 401
359 void G4DNAModelInterface::BuildMaterialMolPerV 402 void G4DNAModelInterface::BuildMaterialMolPerVolTable()
360 { 403 {
361 // To be sure the G4DNAMolecularMaterial is << 404 // To be sure the G4DNAMolecularMaterial is initialized
362 G4DNAMolecularMaterial::Instance()->Initiali << 405 G4DNAMolecularMaterial::Instance()->Initialize();
363 406
364 G4MaterialTable* materialTable = G4Material: << 407 G4MaterialTable* materialTable = G4Material::GetMaterialTable();
365 408
366 // Loop on all the materials inside the "mat << 409 // Loop on all the materials inside the "materialTable"
367 for (auto currentMaterial : *materialTable) << 410 for(size_t i=0, ie=materialTable->size(); i<ie; i++)
368 // Current material << 411 {
369 // Current material name << 412 // Current material
370 const std::size_t & currentMatID = current << 413 G4Material* currentMaterial = materialTable->at(i);
371 << 414
372 // Will the material be used in this inter << 415 // Current material name
373 // Loop on all the materials that can be d << 416 const G4String& currentMatName = currentMaterial->GetName();
374 auto it = fMaterialParticleModelTable.begi << 417
375 auto ite = fMaterialParticleModelTable.end << 418 // Will the material be used in this interface instance ?
376 for (; it != ite; it++) { << 419 // Loop on all the materials that can be dealt with in this class
377 const std::size_t & materialID = it->fir << 420 MaterialParticleModelTable::iterator it = fMaterialParticleModelTable.begin();
378 << 421 MaterialParticleModelTable::iterator ite = fMaterialParticleModelTable.end();
379 if (materialID == currentMatID) { << 422 for(; it != ite; it++)
380 const std::vector<G4double>* numMolPer << 423 {
381 G4DNAMolecularMaterial::Instance()-> << 424 const G4String& materialName = it->first;
382 fMaterialMolPerVol[materialID] = numMo << 425
383 } << 426 if(materialName == currentMatName)
>> 427 {
>> 428 const std::vector<double>* numMolPerVolForMat = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(currentMaterial);
>> 429 fMaterialMolPerVol[materialName] = numMolPerVolForMat;
>> 430 }
>> 431 }
384 } 432 }
385 } <<
386 } 433 }
387 434
388 //....oooOO0OOooo........oooOO0OOooo........oo 435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389 436
390 void G4DNAModelInterface::InsertModelInTable(c << 437 void G4DNAModelInterface::InsertModelInTable(const G4String& matName, const G4String& pName)
391 { 438 {
392 // To insert the model(s) in the table Mater << 439 // To insert the model(s) in the table Material Particule -> Model(s)
393 440
394 // First, we need to check if the current ma << 441 // First, we need to check if the current material has already been inserted in the table.
395 // This is possible because of the composite << 442 // This is possible because of the composite material. We could add a component M1 and then try to add the independant M1 material.
396 // add the independant M1 material. This cas << 443 // This case must be avoided. Checking if M1 is already in the table is the way to avoid it.
397 // table is the way to avoid it. << 444 //
398 // << 445 // Chech if the current material and particle are already in the table.
399 // Check if the current material and particl << 446 // If they are: do nothing.
400 // If they are: do nothing. << 447 // If they are not: add the model(s)
401 // If they are not: add the model(s) << 448 //
402 // << 449 // Check for the material
403 // Check for the material << 450 if(fMaterialParticleModelTable.find(matName) == fMaterialParticleModelTable.end())
404 if (fMaterialParticleModelTable.find(matID) << 451 {
405 // Check for the particle << 452 // Check for the particle
406 if (fMaterialParticleModelTable[matID].fin << 453 if(fMaterialParticleModelTable[matName].find(pName) == fMaterialParticleModelTable[matName].end())
407 G4int modelNbForMaterial = 0; << 454 {
408 for (const auto& it : fRegisteredModels) << 455 G4int modelNbForMaterial (0);
409 auto model = dynamic_cast<G4VDNAModel* << 456
410 if (model != nullptr) { << 457 // Loop on all models registered in the simulation to check:
411 if (model->IsParticleExistingInModel << 458 // 1- if they can be applied to the current material
412 fMaterialParticleModelTable[matID] << 459 // 2- if they can be applied to the current particle
413 // and add one to the "there is a << 460 for(unsigned int i=0, ie=fRegisteredModels.size(); i<ie; ++i)
414 ++modelNbForMaterial; << 461 {
415 } << 462 // check if the model is correct for material and particle (previous 1 and 2)
>> 463 if(fRegisteredModels[i]->IsParticleExistingInModelForMaterial(pName, matName))
>> 464 {
>> 465 // if yes then add the model in the map
>> 466 fMaterialParticleModelTable[matName][pName].push_back(fRegisteredModels[i]);
>> 467
>> 468 // and add one to the "there is a model" material flag
>> 469 ++modelNbForMaterial;
>> 470 }
>> 471 }
>> 472
>> 473 // The model(s) applicable to the currently selected material should be in the map.
>> 474 // We check if there are several models for the material.
>> 475 if(modelNbForMaterial>1)
>> 476 {
>> 477 // If there are several models for a given material and particle couple it could be
>> 478 // because of the energy ranges. We will check if the energy ranges are coherent.
>> 479
>> 480 // Get the models (vector)
>> 481 std::vector<G4VDNAModel*>& models = fMaterialParticleModelTable[matName][pName];
>> 482
>> 483 // Declare a map to sort the limits (G4double) and a model "id" (G4int).
>> 484 // The model id is created on the fly here.
>> 485 // The idea is to fill a map with [limit] = counter. This map will be auto-sorted
>> 486 // and we will check by iterating on it that the counter order is maintained.
>> 487
>> 488 // Delcare the map
>> 489 std::map<G4double, G4int, std::less<G4double> > sortMap;
>> 490
>> 491 G4double smallDiff = 0.01 *eV;
>> 492
>> 493 // Loop on all the model for the current couple
>> 494 // and fill a map with [lim] = modelNumber
>> 495 for(unsigned int ii=0, em=models.size(); ii<em; ++ii)
>> 496 {
>> 497 G4double lowLim = models[ii]->GetLowELimit(matName, pName);
>> 498 G4double highLim = models[ii]->GetHighELimit(matName, pName);
>> 499
>> 500 if(sortMap.find(lowLim) != sortMap.end() )
>> 501 {
>> 502 lowLim += smallDiff;
>> 503 }
>> 504
>> 505 sortMap[lowLim] = ii;
>> 506
>> 507 if(sortMap.find(highLim) != sortMap.end() )
>> 508 {
>> 509 highLim -= smallDiff;
>> 510 }
>> 511
>> 512 sortMap[highLim] = ii;
>> 513 }
>> 514
>> 515 // The map has been created and ordered at this point.
>> 516 // We will check the map order.
>> 517
>> 518 // Loop on the sortMap with iterator and check the order is correct.
>> 519 std::map<G4double, G4int>::iterator it = sortMap.begin();
>> 520
>> 521 // First energy limit value
>> 522 G4double dummyLim = it->first - smallDiff;
>> 523
>> 524 // Loop on all the models again.
>> 525 // The goal is to check if for each limit pairs we have the same model number
>> 526 // and that the upper and lower limit are consistent.
>> 527 for(unsigned int ii=0, eii=models.size(); ii<eii; ++ii)
>> 528 {
>> 529 G4double lim1 = it->first - smallDiff;
>> 530 G4int count1 = it->second;
>> 531
>> 532 // Iterate
>> 533 ++it;
>> 534
>> 535 G4double lim2 = it->first + smallDiff;
>> 536 G4int count2 = it->second;
>> 537
>> 538 // Iterate
>> 539 ++it;
>> 540
>> 541 // Check model number and energy limit consistency
>> 542 // std::abs(dummyLim - lim1) > 1.*eV because we cannot do (dummyLim != lim1)
>> 543 // without experimenting precision loss. Therefore, the std::abs(...) > tolerance is the usual way of avoiding
>> 544 // the issue.
>> 545 if( (count1 != count2) || ( std::abs(dummyLim - lim1) > 1.*eV ) )
>> 546 {
>> 547 // Error
>> 548
>> 549 std::ostringstream oss;
>> 550 oss<<"The material "<<matName<<" and the particle "<<pName;
>> 551 oss<<" have several models registered for the "<<fName<<" interaction and their energy ranges ";
>> 552 oss<<"do not match. \nEnergy ranges: \n";
>> 553
>> 554 for(int iii=0, eiii=models.size(); iii<eiii; ++iii)
>> 555 {
>> 556 oss<<models[iii]->GetName()<<"\n";
>> 557 oss<<"low: "<<models[iii]->GetLowELimit(matName, pName)/eV<<" eV \n";
>> 558 oss<<"high: "<<models[iii]->GetHighELimit(matName, pName)/eV<<" eV \n";
>> 559 }
>> 560
>> 561 G4Exception("G4DNAModelManager::InsertModelInTable","em0006",
>> 562 FatalException, oss.str().c_str());
>> 563 return; // to make some compilers happy
>> 564 }
>> 565
>> 566 dummyLim = lim2;
>> 567 }
>> 568
>> 569 // If we are here then everything was ok.
>> 570 }
>> 571 // no model for the material case
>> 572 else if(modelNbForMaterial==0)
>> 573 {
>> 574 // std::ostringstream oss;
>> 575 // oss<<"The material "<<matName<<" and the particle "<<pName;
>> 576 // oss<<" does not have any model registered for the "<<fName<<" interaction. ";
>> 577
>> 578 // G4Exception("G4DNAModelManager::InsertModelInTable","em0006",
>> 579 // FatalException, oss.str().c_str());
>> 580 // return; // to make some compilers happy
>> 581 }
416 } 582 }
417 else { << 583 }
418 auto index = fpG4_WATER->GetIndex(); << 584 }
419 fMaterialParticleModelTable[index][p << 585
420 ++modelNbForMaterial; << 586 G4VDNAModel *G4DNAModelInterface::GetDNAModel(const G4String &material, const G4String &particle, G4double ekin)
>> 587 {
>> 588 // Output pointer
>> 589 G4VDNAModel* model = 0;
>> 590
>> 591 // Get a reference to all the models for the couple (material and particle)
>> 592 std::vector<G4VDNAModel*>& models = fMaterialParticleModelTable[material][particle];
>> 593
>> 594 // We must choose one of the model(s) accordingly to the particle energy and the model energy range(s)
>> 595
>> 596 //G4bool isOneModelSelected = false;
>> 597
>> 598 // Loop on all the models within the models vector and check if ekin is within the energy range.
>> 599 for(int i=0, ie=models.size(); i<ie; ++i)
>> 600 {
>> 601 // ekin is in the energy range: we select the model and stop the loop.
>> 602 if( ekin >= models[i]->GetLowELimit(material, particle)
>> 603 && ekin < models[i]->GetHighELimit(material, particle) )
>> 604 {
>> 605 // Select the model
>> 606 model = models[i];
>> 607
>> 608 // Boolean flag
>> 609 //isOneModelSelected = true;
>> 610
>> 611 // Quit the for loop
>> 612 break;
421 } 613 }
422 } << 614
423 if (modelNbForMaterial == 0) { << 615 // ekin is not in the energy range: we continue the loop.
424 std::ostringstream oss; << 616 }
425 oss << "The material " << (*G4Material << 617
426 << " and the particle " << p->GetP << 618 // // If no model was selected then fatal error
427 oss << " does not have any model regis << 619 // if(!isOneModelSelected)
428 G4Exception("G4DNAModelInterface::Inse << 620 // {
429 oss.str().c_str()); << 621 // G4String msg = "No model has ";
430 return; // to make some compilers hap << 622 // msg += ekin/eV;
431 } << 623 // msg += " eV in its energy range. Therefore nothing was selected.";
432 } << 624
433 } << 625 // G4Exception("G4DNAModelManager::GetDNAModel","em0006",
434 } << 626 // FatalException,
435 << 627 // msg);
436 G4VEmModel* G4DNAModelInterface::SelectModel(c << 628 // }
437 c << 629
438 c << 630 // Return a pointer to the selected model
439 { << 631 return model;
440 // Output pointer <<
441 G4VEmModel* model = nullptr; <<
442 <<
443 // Get a reference to all the models for the <<
444 auto modelData = fMaterialParticleModelTable <<
445 <<
446 // We must choose one of the model(s) accord <<
447 // range(s) <<
448 <<
449 // Loop on all the models within the models <<
450 auto DNAModel = dynamic_cast<G4VDNAModel*>(m <<
451 G4double lowL, highL; <<
452 if (DNAModel == nullptr) { <<
453 // ekin is in the energy range: we select <<
454 lowL = modelData->LowEnergyLimit(); <<
455 highL = modelData->HighEnergyLimit(); <<
456 if (ekin >= lowL && ekin < highL) { <<
457 // Select the model <<
458 <<
459 model = modelData; <<
460 // return model; <<
461 // Quit the for loop <<
462 // break; <<
463 } <<
464 // ekin is not in the energy range: we con <<
465 } <<
466 else { <<
467 // ekin is in the energy range: we select <<
468 lowL = DNAModel->GetLowELimit(materialID, <<
469 highL = DNAModel->GetHighELimit(materialID <<
470 if (ekin >= lowL && ekin < highL) { <<
471 // Select the model <<
472 model = modelData; <<
473 // return model; <<
474 // Quit the for loop <<
475 // break; <<
476 } <<
477 // ekin is not in the energy range: we con <<
478 } <<
479 //} <<
480 return model; <<
481 } 632 }
482 633
483 //....oooOO0OOooo........oooOO0OOooo........oo 634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
484 635
485 G4double G4DNAModelInterface::GetNumMoleculePe 636 G4double G4DNAModelInterface::GetNumMoleculePerVolumeUnitForMaterial(const G4Material* mat)
486 { 637 {
487 return fMaterialMolPerVol[mat->GetIndex()]-> << 638 return fMaterialMolPerVol[mat->GetName()]->at(mat->GetIndex() );
488 } 639 }
489 640
490 //....oooOO0OOooo........oooOO0OOooo........oo 641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
491 642
492 G4double << 643 G4double G4DNAModelInterface::GetNumMolPerVolUnitForComponentInComposite(const G4Material* component, const G4Material* composite)
493 G4DNAModelInterface::GetNumMolPerVolUnitForCom << 644 {
494 << 645 return fMaterialMolPerVol[component->GetName() ]->at(composite->GetIndex() );
495 { <<
496 return fMaterialMolPerVol[component->GetInde <<
497 } <<
498 <<
499 void G4DNAModelInterface::StreamInfo(std::ostr <<
500 { <<
501 G4long prec = os.precision(5); <<
502 os << "===================================== <<
503 << this->GetName() << " ============ <<
504 << "\n"; <<
505 os << std::setw(15) << "Material#" << std::s <<
506 << std::setw(17) << "LowLimit(MeV)" << st <<
507 << "Fast" << std::setw(13) << "Stationary <<
508 for (const auto& it1 : fMaterialParticleMode <<
509 os << std::setw(15) << (*G4Material::GetMa <<
510 for (const auto& it2 : it1.second) { <<
511 os << std::setw(13) << it2.first->GetPar <<
512 os << std::setw(35) << it2.second->GetNa <<
513 auto DNAModel = dynamic_cast<G4VDNAModel <<
514 if (DNAModel == nullptr) { <<
515 os << std::setw(17) << it2.second->Low <<
516 os << std::setw(17) << it2.second->Hig <<
517 } <<
518 else { <<
519 auto lowL = DNAModel->GetLowELimit(it1 <<
520 auto highL = DNAModel->GetHighELimit(i <<
521 os << std::setw(17) << lowL; <<
522 os << std::setw(17) << highL; <<
523 } <<
524 os << std::setw(13) << "no"; <<
525 os << std::setw(13) << "no"; <<
526 os << std::setw(13) << "no" << G4endl; <<
527 } <<
528 } <<
529 <<
530 os << "===================================== <<
531 "===================================== <<
532 << G4endl; <<
533 os.precision(prec); <<
534 } 646 }
535 647