| 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 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4CrossSectionDataStore 32 // File name: G4CrossSectionDataStore
33 // 33 //
34 // Modifications: 34 // Modifications:
35 // 23.01.2009 V.Ivanchenko add destruction of 35 // 23.01.2009 V.Ivanchenko add destruction of data sets
36 // 29.04.2010 G.Folger modifictaions for i 36 // 29.04.2010 G.Folger modifictaions for integer A & Z
37 // 14.03.2011 V.Ivanchenko fixed DumpPhysicsTa 37 // 14.03.2011 V.Ivanchenko fixed DumpPhysicsTable
38 // 15.08.2011 G.Folger, V.Ivanchenko, T.Koi, D 38 // 15.08.2011 G.Folger, V.Ivanchenko, T.Koi, D.Wright redesign the class
39 // 07.03.2013 M.Maire cosmetic in DumpPhysicsT 39 // 07.03.2013 M.Maire cosmetic in DumpPhysicsTable
40 // 40 //
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
43 43
44 #include "G4CrossSectionDataStore.hh" 44 #include "G4CrossSectionDataStore.hh"
45 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
46 #include "G4UnitsTable.hh" 46 #include "G4UnitsTable.hh"
47 #include "Randomize.hh" 47 #include "Randomize.hh"
48 #include "G4Nucleus.hh" 48 #include "G4Nucleus.hh"
49 49
50 #include "G4DynamicParticle.hh" 50 #include "G4DynamicParticle.hh"
51 #include "G4Isotope.hh" 51 #include "G4Isotope.hh"
52 #include "G4Element.hh" 52 #include "G4Element.hh"
53 #include "G4Material.hh" 53 #include "G4Material.hh"
54 #include "G4MaterialTable.hh" <<
55 #include "G4NistManager.hh" 54 #include "G4NistManager.hh"
56 #include "G4HadronicParameters.hh" <<
57 #include <algorithm> 55 #include <algorithm>
58 #include <typeinfo> << 56
59 57
60 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
61 59
62 G4CrossSectionDataStore::G4CrossSectionDataSto 60 G4CrossSectionDataStore::G4CrossSectionDataStore()
63 : nist(G4NistManager::Instance()) 61 : nist(G4NistManager::Instance())
>> 62 , currentMaterial(nullptr)
>> 63 , matParticle(nullptr)
>> 64 , matKinEnergy(0.0)
>> 65 , matCrossSection(0.0)
>> 66 , nDataSetList(0)
>> 67 , verboseLevel(0)
>> 68 , fastPathFlags()
>> 69 , fastPathParams()
>> 70 , counters()
>> 71 , fastPathCache()
64 {} 72 {}
65 73
66 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
67 75
>> 76 G4CrossSectionDataStore::~G4CrossSectionDataStore()
>> 77 {}
>> 78
>> 79
>> 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
>> 81 G4double
>> 82 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
>> 83 const G4Material* mat , G4bool requiresSlowPath)
>> 84 {
>> 85 //The fast-path algorithm
>> 86 // requiresSlowPath == true => Use slow path independently of other conditions
>> 87 //A. Dotti: modifications to this algorithm following the studies of P. Ruth and R. Fowler
>> 88 // on speeding up the cross-sections calculations. Their algorithm is called in the
>> 89 // following "fast-path" while the normal approach to calcualte cross-section is
>> 90 // referred to as "slow-path".
>> 91 //Some details on the fast-path algorithm:
>> 92 //The idea is to use a cached approximation of the material cross-section.
>> 93 //Starting points:
>> 94 //1- We need the material cross-section for navigation purposes: e.g. to calculate the PIL.
>> 95 //2- If this interaction occurs at the end of a step we need to select the G4Element on which
>> 96 // nucleus the interaction is actually happening. This is done calculating the element cross-section
>> 97 // throwing a random number and selecting the appropriate nucleus (see SampleZandA function)
>> 98 //3- To calculate the material cross section needed for 1- we use the G4Element cross sections.
>> 99 // Since the material cross-section is simply the weighted sum of the element cross-sections.
>> 100 //4- The slow path algorithm here accomplishes two use cases (not very good design IMHO): it
>> 101 // calculates the material cross-section and it updates the xsecelem array that contains the
>> 102 // relative cross-sections used by SampleZandA to select the element on which the interaction
>> 103 // occurs.
>> 104 //The idea of the fast-path algorithm is to replace the request for 1- with a faster calculation
>> 105 //of the material cross-section from an approximation contained in cache.
>> 106 //There two complications to take into account:
>> 107 //A- If I use a fast parametrization for material cross-section, I still need to do the full
>> 108 // calculations if an interaction occurs, because I need to calculate the element cross-sections.
>> 109 // Since the function that updates xsecelem is the same (this one) I need to be sure that
>> 110 // if I call this method for SampleAandZ the xsecelem is updated.
>> 111 //B- It exists the possibility to be even fast the the fast-path algorithm: this happens when
>> 112 // to select the element of the interaction via SampleZandI I call again this method exactly
>> 113 // with the same conditions as when I called this method to calculate the material cross-section.
>> 114 // In such a case xsecelem is updated and we do not need to do much more. This happens when
>> 115 // for example a neutron undergoes an interaction at the end of the step.
>> 116 //Dealing with B- complicates a bit the algorithm.
>> 117 //In summary:
>> 118 // If no fast-path algo is available (or user does not want that), go with the old plain algorithm
>> 119 // If a fast-path algo is avilable, use it whenever possible.
>> 120 //
>> 121 // In general we expect user to selectively decide for which processes, materials and particle combinations
>> 122 // we want to use the fast-path. If this is activated we also expect that the cross-section fast-path
>> 123 // cache is created during the run initialization via calls to this method.
>> 124 //
>> 125 //fastPathFlags contains control flags for the fast-path algorithm:
>> 126 // .prevCalcUsedFastPath == true => Previous call to GetCrossSection used the fast-path
>> 127 // it is used in the decision to assess if xsecelem is correctly set-up
>> 128 // .useFastPathIfAvailable == true => User requested the use of fast-path algorithm
>> 129 // .initializationPhase == true => If true we are in Geant4 Init phase before the event-loop
>> 130
>> 131 //Check user-request, does he want fast-path? if not
>> 132 // OR
>> 133 // we want fast-path and we are in initialization phase?
>> 134 if ( !fastPathFlags.useFastPathIfAvailable
>> 135 || (fastPathFlags.useFastPathIfAvailable&&fastPathFlags.initializationPhase) ) {
>> 136 //Traditional algorithm is requested
>> 137 requiresSlowPath=true;
>> 138 }
>> 139
>> 140 //Logging for performance calculations and counter, active only in FPDEBUG mode
>> 141 counters.MethodCalled();
>> 142 //Measure number of cycles
>> 143 G4FastPathHadronicCrossSection::logStartCountCycles(timing);
>> 144
>> 145 //This is the cache entry of the fast-path cross-section parametrization
>> 146 G4FastPathHadronicCrossSection::cycleCountEntry* entry = nullptr;
>> 147 //Did user request fast-path in first place and are we not in the initialization phase
>> 148 if ( fastPathFlags.useFastPathIfAvailable && !fastPathFlags.initializationPhase ) {
>> 149 //Important: if it is in initialization phase we should NOT use fast path: we are going to build it
>> 150 //G4FastPathHadronicCrossSection::G4CrossSectionDataStore_Key searchkey = {part->GetParticleDefinition(),mat};
>> 151 entry = fastPathCache[{part->GetParticleDefinition(),mat}];
>> 152 }
>> 153
>> 154 //Super-fast-path: are we calling again this method for exactly the same conditions
>> 155 //of the triplet {particle,material,energy}?
>> 156 if(mat == currentMaterial && part->GetDefinition() == matParticle
>> 157 && part->GetKineticEnergy() == matKinEnergy)
>> 158 {
>> 159 G4FastPathHadronicCrossSection::logInvocationTriedOneLine(entry);
>> 160 //If there is no user-request for the fast-path in first place?
>> 161 //It means we built the xsecelem for sure, let's return immediately
>> 162 if ( !fastPathFlags.useFastPathIfAvailable ) {
>> 163 return matCrossSection;
>> 164 } else {
>> 165 //Check that the last time we called this method we used the slow
>> 166 //path: we need the data-member xsecelem to be setup correctly for the current
>> 167 //interaction. This is ensured only if: we will do the slow path right now or we
>> 168 //did it exactly for the same conditions of the last call.
>> 169 if ( !fastPathFlags.prevCalcUsedFastPath && ! requiresSlowPath ) {
>> 170 counters.HitOneLine();
>> 171 G4FastPathHadronicCrossSection::logInvocationOneLine(entry);
>> 172 //Good everything is setup correctly, exit!
>> 173 return matCrossSection;
>> 174 } else {
>> 175 //We need to follow the slow-path because
>> 176 //xsecelem is not calculated correctly
>> 177 requiresSlowPath = true;
>> 178 }
>> 179 }
>> 180 }
>> 181
>> 182 //Ok, now check if we have cached for this {particle,material,energy} the cross-section
>> 183 //in this case let's return immediately, if we are not forced to take the slow path
>> 184 //(e.g. as before if the xsecelem is not up-to-date we need to take the slow-path).
>> 185 //Note that this is not equivalent to the previous ultra-fast check: we now have a map here
>> 186 //So we can have for example a different particle.
>> 187 if ( entry != nullptr && entry->energy == part->GetKineticEnergy() ) {
>> 188 G4FastPathHadronicCrossSection::logHit(entry);
>> 189 if ( !requiresSlowPath ) {
>> 190 return entry->crossSection;
>> 191 }
>> 192 }
>> 193
>> 194 currentMaterial = mat;
>> 195 matParticle = part->GetDefinition();
>> 196 matKinEnergy = part->GetKineticEnergy();
>> 197 matCrossSection = 0;
>> 198
>> 199 //Now check if the cache entry has a fast-path cross-section calculation available
>> 200 G4FastPathHadronicCrossSection::fastPathEntry* fast_entry = nullptr;
>> 201 if ( entry != nullptr && ! requiresSlowPath ) {
>> 202 fast_entry = entry->fastPath;
>> 203 assert(fast_entry!=nullptr && !fastPathFlags.initializationPhase);
>> 204 }
>> 205
>> 206 //Each fast-path cross-section has a minimum value of validity, if energy is below
>> 207 //that skip fast-path algorithm
>> 208 if ( fast_entry != nullptr && part->GetKineticEnergy() < fast_entry->min_cutoff )
>> 209 {
>> 210 assert(requiresSlowPath==false);
>> 211 requiresSlowPath = true;
>> 212 }
>> 213
>> 214 //Ready to use the fast-path calculation
>> 215 if ( !requiresSlowPath && fast_entry != nullptr ) {
>> 216 counters.FastPath();
>> 217 //Retrieve cross-section from fast-path cache
>> 218 matCrossSection = fast_entry->GetCrossSection(part->GetKineticEnergy());
>> 219 fastPathFlags.prevCalcUsedFastPath=true;
>> 220 } else {
>> 221 counters.SlowPath();
>> 222 //Remember that we are now doing the full calculation: xsecelem will
>> 223 //be made valid
>> 224 fastPathFlags.prevCalcUsedFastPath=false;
>> 225
>> 226 G4int nElements = mat->GetNumberOfElements();
>> 227 const G4double* nAtomsPerVolume = mat->GetVecNbOfAtomsPerVolume();
>> 228
>> 229 if(G4int(xsecelm.size()) < nElements) { xsecelm.resize(nElements); }
>> 230
>> 231 for(G4int i=0; i<nElements; ++i) {
>> 232 matCrossSection += nAtomsPerVolume[i] *
>> 233 GetCrossSection(part, (*mat->GetElementVector())[i], mat);
>> 234 xsecelm[i] = matCrossSection;
>> 235 }
>> 236 }
>> 237 //Stop measurement of cpu cycles
>> 238 G4FastPathHadronicCrossSection::logStopCountCycles(timing);
>> 239
>> 240 if ( entry != nullptr ) {
>> 241 entry->energy = part->GetKineticEnergy();
>> 242 entry->crossSection = matCrossSection;
>> 243 }
>> 244 //Some logging of timing
>> 245 G4FastPathHadronicCrossSection::logTiming(entry,fast_entry,timing);
>> 246 return matCrossSection;
>> 247 }
>> 248
>> 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
>> 250
>> 251 void
>> 252 G4CrossSectionDataStore::DumpFastPath(const G4ParticleDefinition* pd, const G4Material* mat,std::ostream& os)
>> 253 {
>> 254 const G4FastPathHadronicCrossSection::cycleCountEntry* entry = fastPathCache[{pd,mat}];
>> 255 if ( entry != nullptr ) {
>> 256 if ( entry->fastPath != nullptr ) {
>> 257 os<<*entry->fastPath;
>> 258 } else {
>> 259 os<<"#Cache entry for {"<<(pd!=nullptr?pd->GetParticleName():"UNDEFINED")<<",";
>> 260 os<<(mat!=nullptr?mat->GetName():"UNDEFINED")<<"} found, but no fast path defined";
>> 261 }
>> 262 } else {
>> 263 os<<"#Cache entry for {"<<(pd!=nullptr?pd->GetParticleName():"UNDEFINED")<<",";
>> 264 os<<(mat!=nullptr?mat->GetName():"UNDEFINED")<<"} not found.";
>> 265 }
>> 266 }
>> 267
>> 268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
>> 269
68 G4double 270 G4double
69 G4CrossSectionDataStore::ComputeCrossSection(c << 271 G4CrossSectionDataStore::ComputeCrossSection(const G4DynamicParticle* part,
70 const G4Material* mat) 272 const G4Material* mat)
71 { 273 {
>> 274 if(part->GetKineticEnergy() == matKinEnergy && mat == currentMaterial &&
>> 275 part->GetDefinition() == matParticle)
>> 276 return matCrossSection;
>> 277
72 currentMaterial = mat; 278 currentMaterial = mat;
73 matParticle = dp->GetDefinition(); << 279 matParticle = part->GetDefinition();
74 matKinEnergy = dp->GetKineticEnergy(); << 280 matKinEnergy = part->GetKineticEnergy();
75 matCrossSection = 0.0; 281 matCrossSection = 0.0;
76 282
77 std::size_t nElements = mat->GetNumberOfElem << 283 size_t nElements = mat->GetNumberOfElements();
78 const G4double* nAtomsPerVolume = mat->GetVe 284 const G4double* nAtomsPerVolume = mat->GetVecNbOfAtomsPerVolume();
79 285
80 if(xsecelm.size() < nElements) { xsecelm.res 286 if(xsecelm.size() < nElements) { xsecelm.resize(nElements); }
81 287
82 for(G4int i=0; i<(G4int)nElements; ++i) { << 288 for(size_t i=0; i<nElements; ++i) {
83 G4double xs = << 289 matCrossSection += nAtomsPerVolume[i] *
84 nAtomsPerVolume[i]*GetCrossSection(dp, m << 290 GetCrossSection(part, mat->GetElement(i), mat);
85 matCrossSection += std::max(xs, 0.0); <<
86 xsecelm[i] = matCrossSection; 291 xsecelm[i] = matCrossSection;
87 } 292 }
88 return matCrossSection; 293 return matCrossSection;
89 } 294 }
90 295
91 //....oooOO0OOooo........oooOO0OOooo........oo 296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
92 297
93 G4double G4CrossSectionDataStore::GetCrossSect << 298 G4double G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
94 299 const G4Element* elm,
95 300 const G4Material* mat)
96 { 301 {
97 // first check the most last cross section << 302 G4int i = nDataSetList - 1;
98 G4int i = nDataSetList-1; <<
99 G4int Z = elm->GetZasInt(); 303 G4int Z = elm->GetZasInt();
100 304
101 if(elm->GetNaturalAbundanceFlag() && 305 if(elm->GetNaturalAbundanceFlag() &&
102 dataSetList[i]->IsElementApplicable(dp, Z << 306 dataSetList[i]->IsElementApplicable(part, Z, mat))
103 { 307 {
104 // element wise cross section 308 // element wise cross section
105 return dataSetList[i]->GetElementCrossSect << 309 return dataSetList[i]->GetElementCrossSection(part, Z, mat);
106 } 310 }
107 311
108 // isotope wise cross section 312 // isotope wise cross section
109 G4int nIso = (G4int)elm->GetNumberOfIsotopes << 313 size_t nIso = elm->GetNumberOfIsotopes();
110 314
111 // user-defined isotope abundances 315 // user-defined isotope abundances
112 const G4double* abundVector = elm->GetRelati 316 const G4double* abundVector = elm->GetRelativeAbundanceVector();
113 317
114 G4double sigma = 0.0; 318 G4double sigma = 0.0;
115 319
116 // isotope and element wise cross sections << 320 for(size_t j = 0; j < nIso; ++j)
117 for(G4int j = 0; j < nIso; ++j) <<
118 { 321 {
119 const G4Isotope* iso = elm->GetIsotope(j); 322 const G4Isotope* iso = elm->GetIsotope(j);
120 sigma += abundVector[j] * 323 sigma += abundVector[j] *
121 GetIsoCrossSection(dp, Z, iso->GetN(), i << 324 GetIsoCrossSection(part, Z, iso->GetN(), iso, elm, mat, i);
122 } 325 }
>> 326
123 return sigma; 327 return sigma;
124 } 328 }
125 329
126 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
127 331
128 G4double 332 G4double
129 G4CrossSectionDataStore::GetIsoCrossSection(co << 333 G4CrossSectionDataStore::GetIsoCrossSection(const G4DynamicParticle* part,
130 G4int Z, G4int A, 334 G4int Z, G4int A,
131 const G4Isotope* iso, 335 const G4Isotope* iso,
132 const G4Element* elm, 336 const G4Element* elm,
133 const G4Material* mat, 337 const G4Material* mat,
134 G4int idx) 338 G4int idx)
135 { 339 {
136 // this methods is called after the check th 340 // this methods is called after the check that dataSetList[idx]
137 // depend on isotopes, so first isotopes are << 341 // depend on isotopes, so for this DataSet only isotopes are checked
138 if(dataSetList[idx]->IsIsoApplicable(dp, Z, <<
139 return dataSetList[idx]->GetIsoCrossSectio <<
140 } <<
141 342
142 // no isotope wise cross section - check oth << 343 // isotope-wise cross section does exist
143 for (G4int j = nDataSetList-1; j >= 0; --j) << 344 if(dataSetList[idx]->IsIsoApplicable(part, Z, A, elm, mat) ) {
144 if(dataSetList[j]->IsElementApplicable(dp, << 345 return dataSetList[idx]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
145 return dataSetList[j]->GetElementCrossSe << 346
146 } else if (dataSetList[j]->IsIsoApplicable << 347 } else {
147 return dataSetList[j]->GetIsoCrossSectio << 348 // seach for other dataSet
>> 349 for (G4int j = nDataSetList-1; j >= 0; --j) {
>> 350 if (dataSetList[j]->IsElementApplicable(part, Z, mat)) {
>> 351 return dataSetList[j]->GetElementCrossSection(part, Z, mat);
>> 352 } else if (dataSetList[j]->IsIsoApplicable(part, Z, A, elm, mat)) {
>> 353 return dataSetList[j]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
>> 354 }
148 } 355 }
149 } 356 }
150 G4ExceptionDescription ed; 357 G4ExceptionDescription ed;
151 ed << "No isotope cross section found for " 358 ed << "No isotope cross section found for "
152 << dp->GetDefinition()->GetParticleName() << 359 << part->GetDefinition()->GetParticleName()
153 << " off target Element " << elm->GetName << 360 << " off Element " << elm->GetName()
154 << " Z= " << Z << " A= " << A; << 361 << " in " << mat->GetName() << " Z= " << Z << " A= " << A
155 if(nullptr != mat) ed << " from " << mat->Ge << 362 << " E(MeV)= " << part->GetKineticEnergy()/MeV << G4endl;
156 ed << " E(MeV)=" << dp->GetKineticEnergy()/M <<
157 G4Exception("G4CrossSectionDataStore::GetIso 363 G4Exception("G4CrossSectionDataStore::GetIsoCrossSection", "had001",
158 FatalException, ed); 364 FatalException, ed);
159 return 0.0; 365 return 0.0;
160 } 366 }
161 367
162 //....oooOO0OOooo........oooOO0OOooo........oo 368 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
163 369
164 G4double 370 G4double
165 G4CrossSectionDataStore::GetCrossSection(const << 371 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
166 G4int 372 G4int Z, G4int A,
167 const G4Isotope* iso, 373 const G4Isotope* iso,
168 const 374 const G4Element* elm,
169 const G4Material* mat) 375 const G4Material* mat)
170 { 376 {
171 for (G4int i = nDataSetList-1; i >= 0; --i) 377 for (G4int i = nDataSetList-1; i >= 0; --i) {
172 if (dataSetList[i]->IsIsoApplicable(dp, Z, << 378 if (dataSetList[i]->IsIsoApplicable(part, Z, A, elm, mat) ) {
173 return dataSetList[i]->GetIsoCrossSectio << 379 return dataSetList[i]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
174 } else if(dataSetList[i]->IsElementApplica <<
175 return dataSetList[i]->GetElementCrossSe <<
176 } 380 }
177 } 381 }
178 G4ExceptionDescription ed; 382 G4ExceptionDescription ed;
179 ed << "No isotope cross section found for " 383 ed << "No isotope cross section found for "
180 << dp->GetDefinition()->GetParticleName() << 384 << part->GetDefinition()->GetParticleName()
181 << " off target Element " << elm->GetName << 385 << " off Element " << elm->GetName()
182 << " Z= " << Z << " A= " << A; << 386 << " in " << mat->GetName() << " Z= " << Z << " A= " << A
183 if(nullptr != mat) ed << " from " << mat->Ge << 387 << " E(MeV)= " << part->GetKineticEnergy()/MeV << G4endl;
184 ed << " E(MeV)=" << dp->GetKineticEnergy()/M <<
185 G4Exception("G4CrossSectionDataStore::GetCro 388 G4Exception("G4CrossSectionDataStore::GetCrossSection", "had001",
186 FatalException, ed); 389 FatalException, ed);
187 return 0.0; 390 return 0.0;
188 } 391 }
189 392
190 //....oooOO0OOooo........oooOO0OOooo........oo 393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
191 394
192 const G4Element* 395 const G4Element*
193 G4CrossSectionDataStore::SampleZandA(const G4D << 396 G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* part,
194 const G4M 397 const G4Material* mat,
195 G4Nucleus& target) 398 G4Nucleus& target)
196 { 399 {
197 if(nullptr != forcedElement) { return forced << 400 size_t nElements = mat->GetNumberOfElements();
198 std::size_t nElements = mat->GetNumberOfElem <<
199 const G4Element* anElement = mat->GetElement 401 const G4Element* anElement = mat->GetElement(0);
200 402
201 // select element from a compound 403 // select element from a compound
202 if(1 < nElements) { 404 if(1 < nElements) {
203 G4double cross = matCrossSection*G4Uniform 405 G4double cross = matCrossSection*G4UniformRand();
204 for(G4int i=0; i<(G4int)nElements; ++i) { << 406 for(size_t i=0; i<nElements; ++i) {
205 if(cross <= xsecelm[i]) { 407 if(cross <= xsecelm[i]) {
206 anElement = mat->GetElement(i); << 408 anElement = mat->GetElement(i);
207 break; 409 break;
208 } 410 }
209 } 411 }
210 } 412 }
211 413
212 G4int Z = anElement->GetZasInt(); 414 G4int Z = anElement->GetZasInt();
213 const G4Isotope* iso = nullptr; 415 const G4Isotope* iso = nullptr;
214 416
215 G4int i = nDataSetList-1; 417 G4int i = nDataSetList-1;
216 if (dataSetList[i]->IsElementApplicable(dp, << 418 if (dataSetList[i]->IsElementApplicable(part, Z, mat)) {
217 419
218 //---------------------------------------- 420 //----------------------------------------------------------------
219 // element-wise cross section 421 // element-wise cross section
220 // isotope cross section is not computed 422 // isotope cross section is not computed
221 //---------------------------------------- 423 //----------------------------------------------------------------
222 std::size_t nIso = anElement->GetNumberOfI << 424 size_t nIso = anElement->GetNumberOfIsotopes();
223 iso = anElement->GetIsotope(0); 425 iso = anElement->GetIsotope(0);
224 426
225 // more than 1 isotope 427 // more than 1 isotope
226 if(1 < nIso) { 428 if(1 < nIso) {
227 iso = dataSetList[i]->SelectIsotope(anEl 429 iso = dataSetList[i]->SelectIsotope(anElement,
228 dp-> << 430 part->GetKineticEnergy(),
229 dp->GetLogKineticEnergy()); << 431 part->GetLogKineticEnergy());
230 } 432 }
231 } else { 433 } else {
232 434
233 //---------------------------------------- 435 //----------------------------------------------------------------
234 // isotope-wise cross section 436 // isotope-wise cross section
235 // isotope cross section is computed 437 // isotope cross section is computed
236 //---------------------------------------- 438 //----------------------------------------------------------------
237 std::size_t nIso = anElement->GetNumberOfI << 439 size_t nIso = anElement->GetNumberOfIsotopes();
238 iso = anElement->GetIsotope(0); 440 iso = anElement->GetIsotope(0);
239 441
240 // more than 1 isotope 442 // more than 1 isotope
241 if(1 < nIso) { 443 if(1 < nIso) {
242 const G4double* abundVector = anElement- 444 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
243 if(xseciso.size() < nIso) { xseciso.resi 445 if(xseciso.size() < nIso) { xseciso.resize(nIso); }
244 446
245 G4double cross = 0.0; 447 G4double cross = 0.0;
246 G4int j; << 448 size_t j;
247 for (j = 0; j<(G4int)nIso; ++j) { << 449 for (j = 0; j<nIso; ++j) {
248 G4double xsec = 0.0; 450 G4double xsec = 0.0;
249 if(abundVector[j] > 0.0) { 451 if(abundVector[j] > 0.0) {
250 iso = anElement->GetIsotope(j); 452 iso = anElement->GetIsotope(j);
251 xsec = abundVector[j]* 453 xsec = abundVector[j]*
252 GetIsoCrossSection(dp, Z, iso->GetN(), i << 454 GetIsoCrossSection(part, Z, iso->GetN(), iso, anElement, mat, i);
253 } 455 }
254 cross += xsec; 456 cross += xsec;
255 xseciso[j] = cross; 457 xseciso[j] = cross;
256 } 458 }
257 cross *= G4UniformRand(); 459 cross *= G4UniformRand();
258 for (j = 0; j<(G4int)nIso; ++j) { << 460 for (j = 0; j<nIso; ++j) {
259 if(cross <= xseciso[j]) { 461 if(cross <= xseciso[j]) {
260 iso = anElement->GetIsotope(j); 462 iso = anElement->GetIsotope(j);
261 break; 463 break;
262 } 464 }
263 } 465 }
264 } 466 }
265 } 467 }
266 target.SetIsotope(iso); 468 target.SetIsotope(iso);
267 return anElement; 469 return anElement;
268 } 470 }
269 471
270 //....oooOO0OOooo........oooOO0OOooo........oo 472 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
271 473
272 void 474 void
273 G4CrossSectionDataStore::BuildPhysicsTable(con << 475 G4CrossSectionDataStore::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
274 { 476 {
275 if (nDataSetList == 0) { 477 if (nDataSetList == 0) {
276 G4ExceptionDescription ed; 478 G4ExceptionDescription ed;
277 ed << "No cross section is registered for 479 ed << "No cross section is registered for "
278 << part.GetParticleName() << G4endl; << 480 << aParticleType.GetParticleName() << G4endl;
279 G4Exception("G4CrossSectionDataStore::Buil 481 G4Exception("G4CrossSectionDataStore::BuildPhysicsTable", "had001",
280 FatalException, ed); 482 FatalException, ed);
281 return; 483 return;
282 } 484 }
283 matParticle = ∂ <<
284 for (G4int i=0; i<nDataSetList; ++i) { 485 for (G4int i=0; i<nDataSetList; ++i) {
285 dataSetList[i]->BuildPhysicsTable(part); << 486 dataSetList[i]->BuildPhysicsTable(aParticleType);
>> 487 }
>> 488 //A.Dotti: if fast-path has been requested we can now create the surrogate
>> 489 // model for fast path.
>> 490 if ( fastPathFlags.useFastPathIfAvailable ) {
>> 491 fastPathFlags.initializationPhase = true;
>> 492 using my_value_type=G4FastPathHadronicCrossSection::G4CrossSectionDataStore_Requests::value_type;
>> 493 //Loop on all requests, if particle matches create the corresponding fsat-path
>> 494 std::for_each( requests.begin() , requests.end() ,
>> 495 [&aParticleType,this](const my_value_type& req) {
>> 496 if ( aParticleType == *req.part_mat.first ) {
>> 497 G4FastPathHadronicCrossSection::cycleCountEntry* entry =
>> 498 new G4FastPathHadronicCrossSection::cycleCountEntry(aParticleType.GetParticleName(),req.part_mat.second);
>> 499 entry->fastPath =
>> 500 new G4FastPathHadronicCrossSection::fastPathEntry(&aParticleType,req.part_mat.second,req.min_cutoff);
>> 501 entry->fastPath->Initialize(this);
>> 502 fastPathCache[req.part_mat] = entry;
>> 503 }
>> 504 }
>> 505 );
>> 506 fastPathFlags.initializationPhase = false;
286 } 507 }
287 const G4MaterialTable* theMatTable = G4Mater << 508 }
288 std::size_t nelm = 0; << 509
289 std::size_t niso = 0; << 510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
290 for(auto mat : *theMatTable) { << 511
291 std::size_t nElements = mat->GetNumberOfEl << 512 void G4CrossSectionDataStore::ActivateFastPath( const G4ParticleDefinition* pdef,
292 nelm = std::max(nelm, nElements); << 513 const G4Material* mat, G4double min_cutoff)
293 for(G4int j=0; j<(G4int)nElements; ++j) { << 514 {
294 niso = std::max(niso, mat->GetElement(j) << 515 assert(pdef!=nullptr&&mat!=nullptr);
295 } << 516 G4FastPathHadronicCrossSection::G4CrossSectionDataStore_Key key={pdef,mat};
>> 517 if ( requests.insert( { key , min_cutoff } ).second ) {
>> 518 G4ExceptionDescription ed;
>> 519 ed << "Attempting to request FastPath for couple: <"
>> 520 << pdef->GetParticleName() << ", " <<mat->GetName()
>> 521 << "> but combination already exists" << G4endl;
>> 522 G4Exception("G4CrossSectionDataStore::ActivateFastPath", "had001",
>> 523 FatalException, ed);
296 } 524 }
297 // define vectors for a run <<
298 xsecelm.resize(nelm, 0.0); <<
299 xseciso.resize(niso, 0.0); <<
300 } 525 }
301 526
302 //....oooOO0OOooo........oooOO0OOooo........oo 527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
303 528
304 void 529 void
305 G4CrossSectionDataStore::DumpPhysicsTable(cons << 530 G4CrossSectionDataStore::DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
306 { 531 {
307 // Print out all cross section data sets use 532 // Print out all cross section data sets used and the energies at
308 // which they apply 533 // which they apply
309 534
310 if (nDataSetList == 0) { 535 if (nDataSetList == 0) {
311 G4cout << "WARNING - G4CrossSectionDataSto 536 G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: "
312 << " no data sets registered" << G4endl; 537 << " no data sets registered" << G4endl;
313 return; 538 return;
314 } 539 }
315 540
316 for (G4int i = nDataSetList-1; i >= 0; --i) 541 for (G4int i = nDataSetList-1; i >= 0; --i) {
317 G4double e1 = dataSetList[i]->GetMinKinEne 542 G4double e1 = dataSetList[i]->GetMinKinEnergy();
318 G4double e2 = dataSetList[i]->GetMaxKinEne 543 G4double e2 = dataSetList[i]->GetMaxKinEnergy();
319 G4cout << 544 G4cout
320 << " Cr_sctns: " << std::setw(25) << 545 << " Cr_sctns: " << std::setw(25) << dataSetList[i]->GetName() << ": "
321 << G4BestUnit(e1, "Energy") << " ---> " << 546 << G4BestUnit(e1, "Energy")
322 << G4BestUnit(e2, "Energy") << "\n"; << 547 << " ---> "
323 if (dataSetList[i]->GetName() == "G4CrossS << 548 << G4BestUnit(e2, "Energy") << "\n";
324 dataSetList[i]->DumpPhysicsTable(part); << 549 if (dataSetList[i]->GetName() == "G4CrossSectionPairGG") {
325 G4cout << G4endl; << 550 dataSetList[i]->DumpPhysicsTable(aParticleType);
326 } << 551 }
327 } 552 }
328 } 553 }
329 554
330 //....oooOO0OOooo........oooOO0OOooo........oo 555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
331 << 556 #include <typeinfo>
332 void G4CrossSectionDataStore::DumpHtml(const G 557 void G4CrossSectionDataStore::DumpHtml(const G4ParticleDefinition& /* pD */,
333 std::of 558 std::ofstream& outFile) const
334 { 559 {
335 // Write cross section data set info to html 560 // Write cross section data set info to html physics list
336 // documentation page 561 // documentation page
337 562
338 G4double ehi = 0; 563 G4double ehi = 0;
339 G4double elo = 0; 564 G4double elo = 0;
340 auto param = G4HadronicParameters::Instance( << 565 G4String physListName(std::getenv("G4PhysListName"));
341 G4String physListName = param->GetPhysListNa <<
342 G4String dirName = param->GetPhysListDocDir( <<
343 <<
344 for (G4int i = nDataSetList-1; i > 0; i--) { 566 for (G4int i = nDataSetList-1; i > 0; i--) {
345 elo = dataSetList[i]->GetMinKinEnergy()/Ge 567 elo = dataSetList[i]->GetMinKinEnergy()/GeV;
346 ehi = dataSetList[i]->GetMaxKinEnergy()/Ge 568 ehi = dataSetList[i]->GetMaxKinEnergy()/GeV;
347 outFile << " <li><b><a href=\"" << ph 569 outFile << " <li><b><a href=\"" << physListName << "_"
348 << dataSetList[i]->GetName() << ".html\" << 570 << dataSetList[i]->GetName() << ".html\"> "
349 << dataSetList[i]->GetName() << "< 571 << dataSetList[i]->GetName() << "</a> from "
350 << elo << " GeV to " << ehi << " G 572 << elo << " GeV to " << ehi << " GeV </b></li>\n";
351 PrintCrossSectionHtml(dataSetList[i], phys << 573 //G4cerr << i << ": XS for " << pD.GetParticleName() << " : " << dataSetList[i]->GetName()
>> 574 // << " typeid : " << typeid(dataSetList[i]).name()<< G4endl;
>> 575 PrintCrossSectionHtml(dataSetList[i]);
352 } 576 }
353 577
354 G4double defaultHi = dataSetList[0]->GetMaxK 578 G4double defaultHi = dataSetList[0]->GetMaxKinEnergy()/GeV;
355 if (ehi < defaultHi) { 579 if (ehi < defaultHi) {
356 outFile << " <li><b><a href=\"" << da << 580 outFile << " <li><b><a href=\"" << dataSetList[0]->GetName() << ".html\"> "
357 << ".html\"> " <<
358 << dataSetList[0]->GetName() << "< 581 << dataSetList[0]->GetName() << "</a> from "
359 << ehi << " GeV to " << defaultHi 582 << ehi << " GeV to " << defaultHi << " GeV </b></li>\n";
360 PrintCrossSectionHtml(dataSetList[0], phys << 583 PrintCrossSectionHtml(dataSetList[0]);
361 } 584 }
362 } 585 }
363 586
364 //....oooOO0OOooo........oooOO0OOooo........oo 587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
365 588
366 void G4CrossSectionDataStore::PrintCrossSectio << 589 void G4CrossSectionDataStore::PrintCrossSectionHtml(const G4VCrossSectionDataSet *cs) const
367 <<
368 <<
369 { 590 {
>> 591 G4String dirName(std::getenv("G4PhysListDocDir"));
>> 592 G4String physListName(std::getenv("G4PhysListName"));
370 593
371 G4String pathName = dirName + "/" + physList << 594 G4String pathName = dirName + "/" + physListName + "_" + HtmlFileName(cs->GetName());
372 std::ofstream outCS; << 595 std::ofstream outCS;
373 outCS.open(pathName); << 596 outCS.open(pathName);
374 outCS << "<html>\n"; << 597 outCS << "<html>\n";
375 outCS << "<head>\n"; << 598 outCS << "<head>\n";
376 outCS << "<title>Description of " << cs->Get << 599 outCS << "<title>Description of " << cs->GetName()
377 << "</title>\n"; << 600 << "</title>\n";
378 outCS << "</head>\n"; << 601 outCS << "</head>\n";
379 outCS << "<body>\n"; << 602 outCS << "<body>\n";
380 <<
381 cs->CrossSectionDescription(outCS); <<
382 <<
383 outCS << "</body>\n"; <<
384 outCS << "</html>\n"; <<
385 } <<
386 603
>> 604 cs->CrossSectionDescription(outCS);
>> 605
>> 606 outCS << "</body>\n";
>> 607 outCS << "</html>\n";
>> 608
>> 609 }
387 //....oooOO0OOooo........oooOO0OOooo........oo 610 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
388 611
389 G4String G4CrossSectionDataStore::HtmlFileName 612 G4String G4CrossSectionDataStore::HtmlFileName(const G4String & in) const
390 { 613 {
391 G4String str(in); 614 G4String str(in);
392 // replace blanks by _ C++11 version: 615 // replace blanks by _ C++11 version:
393 std::transform(str.begin(), str.end(), str. 616 std::transform(str.begin(), str.end(), str.begin(), [](char ch) {
394 return ch == ' ' ? '_' : ch; 617 return ch == ' ' ? '_' : ch;
395 }); 618 });
396 str=str + ".html"; 619 str=str + ".html";
397 return str; 620 return str;
398 } 621 }
399 622
400 //....oooOO0OOooo........oooOO0OOooo........oo 623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
401 624
402 void G4CrossSectionDataStore::AddDataSet(G4VCr 625 void G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* p)
403 { 626 {
404 if(p->ForAllAtomsAndEnergies()) { 627 if(p->ForAllAtomsAndEnergies()) {
405 dataSetList.clear(); 628 dataSetList.clear();
406 nDataSetList = 0; 629 nDataSetList = 0;
407 } 630 }
408 dataSetList.push_back(p); 631 dataSetList.push_back(p);
409 ++nDataSetList; 632 ++nDataSetList;
410 } 633 }
411 634
>> 635
>> 636
412 //....oooOO0OOooo........oooOO0OOooo........oo 637 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
413 638
414 void G4CrossSectionDataStore::AddDataSet(G4VCr << 639 void G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* p, size_t i )
415 { 640 {
416 if(p->ForAllAtomsAndEnergies()) { 641 if(p->ForAllAtomsAndEnergies()) {
417 dataSetList.clear(); 642 dataSetList.clear();
418 dataSetList.push_back(p); 643 dataSetList.push_back(p);
419 nDataSetList = 1; 644 nDataSetList = 1;
420 } else if ( i >= dataSetList.size() ) { <<
421 dataSetList.push_back(p); <<
422 ++nDataSetList; <<
423 } else { 645 } else {
>> 646 if ( i > dataSetList.size() ) i = dataSetList.size();
424 std::vector< G4VCrossSectionDataSet* >::it 647 std::vector< G4VCrossSectionDataSet* >::iterator it = dataSetList.end() - i;
425 dataSetList.insert(it , p); 648 dataSetList.insert(it , p);
426 ++nDataSetList; 649 ++nDataSetList;
427 } 650 }
428 } 651 }
429 <<
430 //....oooOO0OOooo........oooOO0OOooo........oo 652 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
431 653