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