| 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 // $Id: G4CrossSectionDataStore.cc,v 1.15.2.1 2009/03/03 11:48:00 gcosmo Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-01 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class file 31 // GEANT4 Class file
30 // 32 //
31 // 33 //
32 // File name: G4CrossSectionDataStore 34 // File name: G4CrossSectionDataStore
33 // 35 //
>> 36 //
34 // Modifications: 37 // Modifications:
35 // 23.01.2009 V.Ivanchenko add destruction of 38 // 23.01.2009 V.Ivanchenko add destruction of data sets
36 // 29.04.2010 G.Folger modifictaions for i <<
37 // 14.03.2011 V.Ivanchenko fixed DumpPhysicsTa <<
38 // 15.08.2011 G.Folger, V.Ivanchenko, T.Koi, D <<
39 // 07.03.2013 M.Maire cosmetic in DumpPhysicsT <<
40 // 39 //
41 //....oooOO0OOooo........oooOO0OOooo........oo << 40 //
42 //....oooOO0OOooo........oooOO0OOooo........oo <<
43 41
44 #include "G4CrossSectionDataStore.hh" 42 #include "G4CrossSectionDataStore.hh"
45 #include "G4SystemOfUnits.hh" << 43 #include "G4HadronicException.hh"
46 #include "G4UnitsTable.hh" << 44 #include "G4StableIsotopes.hh"
>> 45 #include "G4HadTmpUtil.hh"
47 #include "Randomize.hh" 46 #include "Randomize.hh"
48 #include "G4Nucleus.hh" << 47 #include "G4Nucleus.hh"
49 <<
50 #include "G4DynamicParticle.hh" <<
51 #include "G4Isotope.hh" <<
52 #include "G4Element.hh" <<
53 #include "G4Material.hh" <<
54 #include "G4MaterialTable.hh" <<
55 #include "G4NistManager.hh" <<
56 #include "G4HadronicParameters.hh" <<
57 #include <algorithm> <<
58 #include <typeinfo> <<
59 48
60 //....oooOO0OOooo........oooOO0OOooo........oo << 49 G4CrossSectionDataStore::G4CrossSectionDataStore() :
61 << 50 NDataSetList(0), verboseLevel(0)
62 G4CrossSectionDataStore::G4CrossSectionDataSto <<
63 : nist(G4NistManager::Instance()) <<
64 {} 51 {}
65 52
66 //....oooOO0OOooo........oooOO0OOooo........oo << 53 G4CrossSectionDataStore::~G4CrossSectionDataStore()
>> 54 {}
67 55
68 G4double << 56 G4double
69 G4CrossSectionDataStore::ComputeCrossSection(c << 57 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
70 const G4Material* mat) << 58 const G4Element* anElement,
>> 59 G4double aTemperature)
71 { 60 {
72 currentMaterial = mat; << 61 if (NDataSetList == 0)
73 matParticle = dp->GetDefinition(); << 62 {
74 matKinEnergy = dp->GetKineticEnergy(); << 63 throw G4HadronicException(__FILE__, __LINE__,
75 matCrossSection = 0.0; << 64 "G4CrossSectionDataStore: no data sets registered");
76 << 65 return DBL_MIN;
77 std::size_t nElements = mat->GetNumberOfElem << 66 }
78 const G4double* nAtomsPerVolume = mat->GetVe << 67 for (G4int i = NDataSetList-1; i >= 0; i--) {
79 << 68 if (DataSetList[i]->IsApplicable(aParticle, anElement))
80 if(xsecelm.size() < nElements) { xsecelm.res << 69 return DataSetList[i]->GetCrossSection(aParticle,anElement,aTemperature);
81 <<
82 for(G4int i=0; i<(G4int)nElements; ++i) { <<
83 G4double xs = <<
84 nAtomsPerVolume[i]*GetCrossSection(dp, m <<
85 matCrossSection += std::max(xs, 0.0); <<
86 xsecelm[i] = matCrossSection; <<
87 } 70 }
88 return matCrossSection; << 71 throw G4HadronicException(__FILE__, __LINE__,
>> 72 "G4CrossSectionDataStore: no applicable data set found "
>> 73 "for particle/element");
>> 74 return DBL_MIN;
89 } 75 }
90 76
91 //....oooOO0OOooo........oooOO0OOooo........oo << 77 G4VCrossSectionDataSet*
92 << 78 G4CrossSectionDataStore::whichDataSetInCharge(const G4DynamicParticle* aParticle,
93 G4double G4CrossSectionDataStore::GetCrossSect << 79 const G4Element* anElement)
94 <<
95 <<
96 { 80 {
97 // first check the most last cross section << 81 if (NDataSetList == 0)
98 G4int i = nDataSetList-1; <<
99 G4int Z = elm->GetZasInt(); <<
100 <<
101 if(elm->GetNaturalAbundanceFlag() && <<
102 dataSetList[i]->IsElementApplicable(dp, Z <<
103 { 82 {
104 // element wise cross section << 83 throw G4HadronicException(__FILE__, __LINE__,
105 return dataSetList[i]->GetElementCrossSect << 84 "G4CrossSectionDataStore: no data sets registered");
>> 85 return 0;
>> 86 }
>> 87 for (G4int i = NDataSetList-1; i >= 0; i--) {
>> 88 if (DataSetList[i]->IsApplicable(aParticle, anElement) )
>> 89 {
>> 90 return DataSetList[i];
>> 91 }
106 } 92 }
>> 93 throw G4HadronicException(__FILE__, __LINE__,
>> 94 "G4CrossSectionDataStore: no applicable data set found "
>> 95 "for particle/element");
>> 96 return 0;
>> 97 }
107 98
108 // isotope wise cross section <<
109 G4int nIso = (G4int)elm->GetNumberOfIsotopes <<
110 99
111 // user-defined isotope abundances << 100 G4double
112 const G4double* abundVector = elm->GetRelati << 101 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
>> 102 const G4Isotope* anIsotope,
>> 103 G4double aTemperature)
>> 104 {
>> 105 if (NDataSetList == 0)
>> 106 {
>> 107 throw G4HadronicException(__FILE__, __LINE__,
>> 108 "G4CrossSectionDataStore: no data sets registered");
>> 109 return DBL_MIN;
>> 110 }
>> 111 for (G4int i = NDataSetList-1; i >= 0; i--) {
>> 112 if (DataSetList[i]->IsZAApplicable(aParticle, anIsotope->GetZ(), anIsotope->GetN()))
>> 113 return DataSetList[i]->GetIsoCrossSection(aParticle,anIsotope,aTemperature);
>> 114 }
>> 115 throw G4HadronicException(__FILE__, __LINE__,
>> 116 "G4CrossSectionDataStore: no applicable data set found "
>> 117 "for particle/element");
>> 118 return DBL_MIN;
>> 119 }
113 120
114 G4double sigma = 0.0; <<
115 121
116 // isotope and element wise cross sections << 122 G4double
117 for(G4int j = 0; j < nIso; ++j) << 123 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
>> 124 G4double Z, G4double A,
>> 125 G4double aTemperature)
>> 126 {
>> 127 if (NDataSetList == 0)
118 { 128 {
119 const G4Isotope* iso = elm->GetIsotope(j); << 129 throw G4HadronicException(__FILE__, __LINE__,
120 sigma += abundVector[j] * << 130 "G4CrossSectionDataStore: no data sets registered");
121 GetIsoCrossSection(dp, Z, iso->GetN(), i << 131 return DBL_MIN;
122 } << 132 }
123 return sigma; << 133 for (G4int i = NDataSetList-1; i >= 0; i--) {
>> 134 if (DataSetList[i]->IsZAApplicable(aParticle, Z, A))
>> 135 return DataSetList[i]->GetIsoZACrossSection(aParticle,Z,A,aTemperature);
>> 136 }
>> 137 throw G4HadronicException(__FILE__, __LINE__,
>> 138 "G4CrossSectionDataStore: no applicable data set found "
>> 139 "for particle/element");
>> 140 return DBL_MIN;
124 } 141 }
125 142
126 //....oooOO0OOooo........oooOO0OOooo........oo <<
127 143
128 G4double 144 G4double
129 G4CrossSectionDataStore::GetIsoCrossSection(co << 145 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
130 G4int Z, G4int A, << 146 const G4Material* aMaterial)
131 const G4Isotope* iso, << 147 {
132 const G4Element* elm, << 148 G4double sigma(0);
133 const G4Material* mat, << 149 G4double aTemp = aMaterial->GetTemperature();
134 G4int idx) << 150 G4int nElements = aMaterial->GetNumberOfElements();
135 { << 151 const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
136 // this methods is called after the check th << 152 G4double xSection(0);
137 // depend on isotopes, so first isotopes are << 153
138 if(dataSetList[idx]->IsIsoApplicable(dp, Z, << 154 for(G4int i=0; i<nElements; i++) {
139 return dataSetList[idx]->GetIsoCrossSectio << 155 xSection = GetCrossSection( aParticle, (*aMaterial->GetElementVector())[i], aTemp);
140 } << 156 sigma += theAtomsPerVolumeVector[i] * xSection;
141 <<
142 // no isotope wise cross section - check oth <<
143 for (G4int j = nDataSetList-1; j >= 0; --j) <<
144 if(dataSetList[j]->IsElementApplicable(dp, <<
145 return dataSetList[j]->GetElementCrossSe <<
146 } else if (dataSetList[j]->IsIsoApplicable <<
147 return dataSetList[j]->GetIsoCrossSectio <<
148 } <<
149 } 157 }
150 G4ExceptionDescription ed; << 158
151 ed << "No isotope cross section found for " << 159 return sigma;
152 << dp->GetDefinition()->GetParticleName() <<
153 << " off target Element " << elm->GetName <<
154 << " Z= " << Z << " A= " << A; <<
155 if(nullptr != mat) ed << " from " << mat->Ge <<
156 ed << " E(MeV)=" << dp->GetKineticEnergy()/M <<
157 G4Exception("G4CrossSectionDataStore::GetIso <<
158 FatalException, ed); <<
159 return 0.0; <<
160 } 160 }
161 161
162 //....oooOO0OOooo........oooOO0OOooo........oo <<
163 162
164 G4double << 163 G4Element* G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* particle,
165 G4CrossSectionDataStore::GetCrossSection(const << 164 const G4Material* aMaterial,
166 G4int << 165 G4Nucleus& target)
167 const G4Isotope* iso, << 166 {
168 const << 167 G4double aTemp = aMaterial->GetTemperature();
169 const G4Material* mat) << 168 const G4int nElements = aMaterial->GetNumberOfElements();
170 { << 169 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
171 for (G4int i = nDataSetList-1; i >= 0; --i) << 170 G4Element* anElement = (*theElementVector)[0];
172 if (dataSetList[i]->IsIsoApplicable(dp, Z, << 171 G4VCrossSectionDataSet* inCharge;
173 return dataSetList[i]->GetIsoCrossSectio << 172 G4int i;
174 } else if(dataSetList[i]->IsElementApplica <<
175 return dataSetList[i]->GetElementCrossSe <<
176 } <<
177 } <<
178 G4ExceptionDescription ed; <<
179 ed << "No isotope cross section found for " <<
180 << dp->GetDefinition()->GetParticleName() <<
181 << " off target Element " << elm->GetName <<
182 << " Z= " << Z << " A= " << A; <<
183 if(nullptr != mat) ed << " from " << mat->Ge <<
184 ed << " E(MeV)=" << dp->GetKineticEnergy()/M <<
185 G4Exception("G4CrossSectionDataStore::GetCro <<
186 FatalException, ed); <<
187 return 0.0; <<
188 } <<
189 <<
190 //....oooOO0OOooo........oooOO0OOooo........oo <<
191 <<
192 const G4Element* <<
193 G4CrossSectionDataStore::SampleZandA(const G4D <<
194 const G4M <<
195 G4Nucleus& target) <<
196 { <<
197 if(nullptr != forcedElement) { return forced <<
198 std::size_t nElements = mat->GetNumberOfElem <<
199 const G4Element* anElement = mat->GetElement <<
200 173
201 // select element from a compound << 174 // compounds
202 if(1 < nElements) { 175 if(1 < nElements) {
203 G4double cross = matCrossSection*G4Uniform << 176 G4double* xsec = new G4double [nElements];
204 for(G4int i=0; i<(G4int)nElements; ++i) { << 177 const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
205 if(cross <= xsecelm[i]) { << 178 G4double cross = 0.0;
206 anElement = mat->GetElement(i); << 179 for(i=0; i<nElements; i++) {
207 break; << 180 anElement= (*theElementVector)[i];
>> 181 inCharge = whichDataSetInCharge(particle, anElement);
>> 182 cross += theAtomsPerVolumeVector[i]*
>> 183 inCharge->GetCrossSection(particle, anElement, aTemp);
>> 184 xsec[i] = cross;
>> 185 }
>> 186 cross *= G4UniformRand();
>> 187
>> 188 for(i=0; i<nElements; i++) {
>> 189 if( cross <= xsec[i] ) {
>> 190 anElement = (*theElementVector)[i];
>> 191 break;
208 } 192 }
209 } 193 }
>> 194 delete [] xsec;
210 } 195 }
211 196
212 G4int Z = anElement->GetZasInt(); << 197 // element have been selected
213 const G4Isotope* iso = nullptr; << 198 inCharge = whichDataSetInCharge(particle, anElement);
>> 199 G4double ZZ = anElement->GetZ();
>> 200 G4double AA;
>> 201
>> 202 // Collect abundance weighted cross sections and A values for each isotope
>> 203 // in each element
>> 204
>> 205 const G4int nIsoPerElement = anElement->GetNumberOfIsotopes();
214 206
215 G4int i = nDataSetList-1; << 207 // user-defined isotope abundances
216 if (dataSetList[i]->IsElementApplicable(dp, << 208 if (0 < nIsoPerElement) {
>> 209 G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
>> 210 AA = G4double((*isoVector)[0]->GetN());
>> 211 if(1 < nIsoPerElement) {
>> 212
>> 213 G4double* xsec = new G4double [nIsoPerElement];
>> 214 G4double iso_xs = 0.0;
>> 215 G4double cross = 0.0;
>> 216
>> 217 G4double* abundVector = anElement->GetRelativeAbundanceVector();
>> 218 G4bool elementXS = false;
>> 219 for (i = 0; i<nIsoPerElement; i++) {
>> 220 if (inCharge->IsZAApplicable(particle, ZZ, G4double((*isoVector)[i]->GetN()))) {
>> 221 iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[i], aTemp);
>> 222 } else if (elementXS == false) {
>> 223 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
>> 224 elementXS = true;
>> 225 }
217 226
218 //---------------------------------------- << 227 cross += abundVector[i]*iso_xs;
219 // element-wise cross section << 228 xsec[i] = cross;
220 // isotope cross section is not computed << 229 }
221 //---------------------------------------- << 230 cross *= G4UniformRand();
222 std::size_t nIso = anElement->GetNumberOfI << 231 for (i = 0; i<nIsoPerElement; i++) {
223 iso = anElement->GetIsotope(0); << 232 if(cross <= xsec[i]) {
224 << 233 AA = G4double((*isoVector)[i]->GetN());
225 // more than 1 isotope << 234 break;
226 if(1 < nIso) { << 235 }
227 iso = dataSetList[i]->SelectIsotope(anEl << 236 }
228 dp-> << 237 delete [] xsec;
229 dp->GetLogKineticEnergy()); <<
230 } 238 }
231 } else { << 239 // natural abundances
>> 240 } else {
232 241
233 //---------------------------------------- << 242 G4StableIsotopes theDefaultIsotopes;
234 // isotope-wise cross section << 243 G4int Z = G4int(ZZ + 0.5);
235 // isotope cross section is computed << 244 const G4int nIso = theDefaultIsotopes.GetNumberOfIsotopes(Z);
236 //---------------------------------------- << 245 G4int index = theDefaultIsotopes.GetFirstIsotope(Z);
237 std::size_t nIso = anElement->GetNumberOfI << 246 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index));
238 iso = anElement->GetIsotope(0); << 247
239 <<
240 // more than 1 isotope <<
241 if(1 < nIso) { 248 if(1 < nIso) {
242 const G4double* abundVector = anElement- <<
243 if(xseciso.size() < nIso) { xseciso.resi <<
244 249
245 G4double cross = 0.0; << 250 G4double* xsec = new G4double [nIso];
246 G4int j; << 251 G4double iso_xs = 0.0;
247 for (j = 0; j<(G4int)nIso; ++j) { << 252 G4double cross = 0.0;
248 G4double xsec = 0.0; << 253 G4bool elementXS= false;
249 if(abundVector[j] > 0.0) { << 254
250 iso = anElement->GetIsotope(j); << 255 for (i = 0; i<nIso; i++) {
251 xsec = abundVector[j]* << 256 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+i));
252 GetIsoCrossSection(dp, Z, iso->GetN(), i << 257 if (inCharge->IsZAApplicable(particle, ZZ, AA )) {
>> 258 iso_xs = inCharge->GetIsoZACrossSection(particle, ZZ, AA, aTemp);
>> 259 } else if (elementXS == false) {
>> 260 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
>> 261 elementXS = true;
253 } 262 }
254 cross += xsec; << 263 cross += theDefaultIsotopes.GetAbundance(index+i)*iso_xs;
255 xseciso[j] = cross; << 264 xsec[i] = cross;
256 } 265 }
257 cross *= G4UniformRand(); 266 cross *= G4UniformRand();
258 for (j = 0; j<(G4int)nIso; ++j) { << 267 for (i = 0; i<nIso; i++) {
259 if(cross <= xseciso[j]) { << 268 if(cross <= xsec[i]) {
260 iso = anElement->GetIsotope(j); << 269 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+i));
261 break; 270 break;
262 } 271 }
263 } 272 }
>> 273 delete [] xsec;
264 } 274 }
265 } 275 }
266 target.SetIsotope(iso); << 276 //G4cout << "XS: " << particle->GetDefinition()->GetParticleName()
>> 277 // << " e(GeV)= " << particle->GetKineticEnergy()/GeV
>> 278 // << " in " << aMaterial->GetName()
>> 279 // << " ZZ= " << ZZ << " AA= " << AA << " " << anElement->GetName() << G4endl;
>> 280
>> 281 target.SetParameters(AA, ZZ);
267 return anElement; 282 return anElement;
268 } 283 }
269 284
270 //....oooOO0OOooo........oooOO0OOooo........oo << 285 /*
>> 286 G4Element* G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* particle,
>> 287 const G4Material* aMaterial,
>> 288 G4Nucleus& target)
>> 289 {
>> 290 static G4StableIsotopes theDefaultIsotopes; // natural abundances and
>> 291 // stable isotopes
>> 292 G4double aTemp = aMaterial->GetTemperature();
>> 293 G4int nElements = aMaterial->GetNumberOfElements();
>> 294 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 295 const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
>> 296 std::vector<std::vector<G4double> > awicsPerElement;
>> 297 std::vector<std::vector<G4double> > AvaluesPerElement;
>> 298 G4Element* anElement;
>> 299
>> 300 // Collect abundance weighted cross sections and A values for each isotope
>> 301 // in each element
>> 302
>> 303 for (G4int i = 0; i < nElements; i++) {
>> 304 anElement = (*theElementVector)[i];
>> 305 G4int nIsoPerElement = anElement->GetNumberOfIsotopes();
>> 306 std::vector<G4double> isoholder;
>> 307 std::vector<G4double> aholder;
>> 308 G4double iso_xs = DBL_MIN;
>> 309
>> 310 if (nIsoPerElement) { // user-defined isotope abundances
>> 311 G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
>> 312 G4double* abundVector = anElement->GetRelativeAbundanceVector();
>> 313 G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
>> 314 G4bool elementXS = false;
>> 315 for (G4int j = 0; j < nIsoPerElement; j++) {
>> 316 if (inCharge->IsZAApplicable(particle, (*isoVector)[j]->GetZ(),
>> 317 (*isoVector)[j]->GetN() ) ) {
>> 318 iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[j], aTemp);
>> 319 } else if (elementXS == false) {
>> 320 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
>> 321 elementXS = true;
>> 322 }
271 323
272 void << 324 isoholder.push_back(abundVector[j]*iso_xs);
273 G4CrossSectionDataStore::BuildPhysicsTable(con << 325 aholder.push_back(G4double((*isoVector)[j]->GetN()));
274 { << 326 }
275 if (nDataSetList == 0) { <<
276 G4ExceptionDescription ed; <<
277 ed << "No cross section is registered for <<
278 << part.GetParticleName() << G4endl; <<
279 G4Exception("G4CrossSectionDataStore::Buil <<
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 } <<
296 } <<
297 // define vectors for a run <<
298 xsecelm.resize(nelm, 0.0); <<
299 xseciso.resize(niso, 0.0); <<
300 } <<
301 <<
302 //....oooOO0OOooo........oooOO0OOooo........oo <<
303 327
304 void << 328 } else { // natural abundances
305 G4CrossSectionDataStore::DumpPhysicsTable(cons << 329 G4int ZZ = G4lrint(anElement->GetZ());
306 { << 330 nIsoPerElement = theDefaultIsotopes.GetNumberOfIsotopes(ZZ);
307 // Print out all cross section data sets use << 331 G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ);
308 // which they apply << 332 G4double AA;
>> 333 G4double abundance;
>> 334
>> 335 G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
>> 336 G4bool elementXS = false;
>> 337
>> 338 for (G4int j = 0; j < nIsoPerElement; j++) {
>> 339 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+j));
>> 340 aholder.push_back(AA);
>> 341 if (inCharge->IsZAApplicable(particle, G4double(ZZ), AA) ) {
>> 342 iso_xs = inCharge->GetIsoZACrossSection(particle, G4double(ZZ), AA, aTemp);
>> 343 } else if (elementXS == false) {
>> 344 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
>> 345 elementXS = true;
>> 346 }
>> 347 abundance = theDefaultIsotopes.GetAbundance(index+j)/100.0;
>> 348 isoholder.push_back(abundance*iso_xs);
>> 349 }
>> 350 }
309 351
310 if (nDataSetList == 0) { << 352 awicsPerElement.push_back(isoholder);
311 G4cout << "WARNING - G4CrossSectionDataSto << 353 AvaluesPerElement.push_back(aholder);
312 << " no data sets registered" << G4endl; <<
313 return; <<
314 } 354 }
315 355
316 for (G4int i = nDataSetList-1; i >= 0; --i) << 356 // Calculate running sums for isotope selection
317 G4double e1 = dataSetList[i]->GetMinKinEne << 357
318 G4double e2 = dataSetList[i]->GetMaxKinEne << 358 G4double crossSectionTotal = 0;
319 G4cout << 359 G4double xSectionPerElement;
320 << " Cr_sctns: " << std::setw(25) << << 360 std::vector<G4double> runningSum;
321 << G4BestUnit(e1, "Energy") << " ---> " << 361
322 << G4BestUnit(e2, "Energy") << "\n"; << 362 for (G4int i=0; i < nElements; i++) {
323 if (dataSetList[i]->GetName() == "G4CrossS << 363 xSectionPerElement = 0;
324 dataSetList[i]->DumpPhysicsTable(part); << 364 for (G4int j=0; j < G4int(awicsPerElement[i].size()); j++)
325 G4cout << G4endl; << 365 xSectionPerElement += awicsPerElement[i][j];
>> 366 runningSum.push_back(theAtomsPerVolumeVector[i]*xSectionPerElement);
>> 367 crossSectionTotal += runningSum[i];
>> 368 }
>> 369
>> 370 // Compare random number to running sum over element xc to choose Z
>> 371
>> 372 // Initialize Z and A to first element and first isotope in case
>> 373 // cross section is zero
>> 374
>> 375 anElement = (*theElementVector)[0];
>> 376 G4double ZZ = anElement->GetZ();
>> 377 G4double AA = AvaluesPerElement[0][0];
>> 378 if (crossSectionTotal != 0.) {
>> 379 G4double random = G4UniformRand();
>> 380 for(G4int i=0; i < nElements; i++) {
>> 381 if(i!=0) runningSum[i] += runningSum[i-1];
>> 382 if(random <= runningSum[i]/crossSectionTotal) {
>> 383 anElement = (*theElementVector)[i];
>> 384 ZZ = anElement->GetZ();
>> 385
>> 386 // Compare random number to running sum over isotope xc to choose A
>> 387
>> 388 G4int nIso = awicsPerElement[i].size();
>> 389 G4double* running = new G4double[nIso];
>> 390 for (G4int j=0; j < nIso; j++) {
>> 391 running[j] = awicsPerElement[i][j];
>> 392 if(j!=0) running[j] += running[j-1];
>> 393 }
>> 394
>> 395 G4double trial = G4UniformRand();
>> 396 for (G4int j=0; j < nIso; j++) {
>> 397 AA = AvaluesPerElement[i][j];
>> 398 if (trial <= running[j]/running[nIso-1]) break;
>> 399 }
>> 400 delete [] running;
>> 401 break;
>> 402 }
326 } 403 }
327 } 404 }
328 } <<
329 405
330 //....oooOO0OOooo........oooOO0OOooo........oo << 406 //G4cout << "XS: " << particle->GetDefinition()->GetParticleName()
>> 407 // << " e(GeV)= " << particle->GetKineticEnergy()/GeV
>> 408 // << " in " << aMaterial->GetName()
>> 409 // << " ZZ= " << ZZ << " AA= " << AA << " " << anElement->GetName() << G4endl;
331 410
332 void G4CrossSectionDataStore::DumpHtml(const G << 411 target.SetParameters(AA, ZZ);
333 std::of << 412 return anElement;
334 { << 413 }
335 // Write cross section data set info to html <<
336 // documentation page <<
337 414
338 G4double ehi = 0; << 415 std::pair<G4double, G4double>
339 G4double elo = 0; << 416 G4CrossSectionDataStore::SelectRandomIsotope(const G4DynamicParticle* particle,
340 auto param = G4HadronicParameters::Instance( << 417 const G4Material* aMaterial)
341 G4String physListName = param->GetPhysListNa << 418 {
342 G4String dirName = param->GetPhysListDocDir( << 419 static G4StableIsotopes theDefaultIsotopes; // natural abundances and
>> 420 // stable isotopes
>> 421 G4double aTemp = aMaterial->GetTemperature();
>> 422 G4int nElements = aMaterial->GetNumberOfElements();
>> 423 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 424 const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume();
>> 425 std::vector<std::vector<G4double> > awicsPerElement;
>> 426 std::vector<std::vector<G4double> > AvaluesPerElement;
>> 427 G4Element* anElement;
>> 428
>> 429 // Collect abundance weighted cross sections and A values for each isotope
>> 430 // in each element
>> 431
>> 432 for (G4int i = 0; i < nElements; i++) {
>> 433 anElement = (*theElementVector)[i];
>> 434 G4int nIsoPerElement = anElement->GetNumberOfIsotopes();
>> 435 std::vector<G4double> isoholder;
>> 436 std::vector<G4double> aholder;
>> 437 G4double iso_xs = DBL_MIN;
>> 438
>> 439 if (nIsoPerElement) { // user-defined isotope abundances
>> 440 G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
>> 441 G4double* abundVector = anElement->GetRelativeAbundanceVector();
>> 442 G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
>> 443 G4bool elementXS = false;
>> 444 for (G4int j = 0; j < nIsoPerElement; j++) {
>> 445 if (inCharge->IsZAApplicable(particle, (*isoVector)[j]->GetZ(),
>> 446 (*isoVector)[j]->GetN() ) ) {
>> 447 iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[j], aTemp);
>> 448 } else if (elementXS == false) {
>> 449 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
>> 450 elementXS = true;
>> 451 }
343 452
344 for (G4int i = nDataSetList-1; i > 0; i--) { << 453 isoholder.push_back(abundVector[j]*iso_xs);
345 elo = dataSetList[i]->GetMinKinEnergy()/Ge << 454 aholder.push_back(G4double((*isoVector)[j]->GetN()));
346 ehi = dataSetList[i]->GetMaxKinEnergy()/Ge << 455 }
347 outFile << " <li><b><a href=\"" << ph <<
348 << dataSetList[i]->GetName() << ".html\" <<
349 << dataSetList[i]->GetName() << "< <<
350 << elo << " GeV to " << ehi << " G <<
351 PrintCrossSectionHtml(dataSetList[i], phys <<
352 } <<
353 456
354 G4double defaultHi = dataSetList[0]->GetMaxK << 457 } else { // natural abundances
355 if (ehi < defaultHi) { << 458 G4int ZZ = G4lrint(anElement->GetZ());
356 outFile << " <li><b><a href=\"" << da << 459 nIsoPerElement = theDefaultIsotopes.GetNumberOfIsotopes(ZZ);
357 << ".html\"> " << 460 G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ);
358 << dataSetList[0]->GetName() << "< << 461 G4double AA;
359 << ehi << " GeV to " << defaultHi << 462 G4double abundance;
360 PrintCrossSectionHtml(dataSetList[0], phys << 463
361 } << 464 G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement);
362 } << 465 G4bool elementXS = false;
>> 466
>> 467 for (G4int j = 0; j < nIsoPerElement; j++) {
>> 468 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+j));
>> 469 aholder.push_back(AA);
>> 470 if (inCharge->IsZAApplicable(particle, G4double(ZZ), AA) ) {
>> 471 iso_xs = inCharge->GetIsoZACrossSection(particle, G4double(ZZ), AA, aTemp);
>> 472 } else if (elementXS == false) {
>> 473 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp);
>> 474 elementXS = true;
>> 475 }
>> 476 abundance = theDefaultIsotopes.GetAbundance(index+j)/100.0;
>> 477 isoholder.push_back(abundance*iso_xs);
>> 478 }
>> 479 }
363 480
364 //....oooOO0OOooo........oooOO0OOooo........oo << 481 awicsPerElement.push_back(isoholder);
>> 482 AvaluesPerElement.push_back(aholder);
>> 483 }
365 484
366 void G4CrossSectionDataStore::PrintCrossSectio << 485 // Calculate running sums for isotope selection
367 <<
368 <<
369 { <<
370 486
371 G4String pathName = dirName + "/" + physList << 487 G4double crossSectionTotal = 0;
372 std::ofstream outCS; << 488 G4double xSectionPerElement;
373 outCS.open(pathName); << 489 std::vector<G4double> runningSum;
374 outCS << "<html>\n"; << 490
375 outCS << "<head>\n"; << 491 for (G4int i=0; i < nElements; i++) {
376 outCS << "<title>Description of " << cs->Get << 492 xSectionPerElement = 0;
377 << "</title>\n"; << 493 for (G4int j=0; j < G4int(awicsPerElement[i].size()); j++)
378 outCS << "</head>\n"; << 494 xSectionPerElement += awicsPerElement[i][j];
379 outCS << "<body>\n"; << 495 runningSum.push_back(theAtomsPerVolumeVector[i]*xSectionPerElement);
380 << 496 crossSectionTotal += runningSum[i];
381 cs->CrossSectionDescription(outCS); << 497 }
382 << 498
383 outCS << "</body>\n"; << 499 // Compare random number to running sum over element xc to choose Z
384 outCS << "</html>\n"; << 500
>> 501 // Initialize Z and A to first element and first isotope in case
>> 502 // cross section is zero
>> 503
>> 504 G4double ZZ = (*theElementVector)[0]->GetZ();
>> 505 G4double AA = AvaluesPerElement[0][0];
>> 506 if (crossSectionTotal != 0.) {
>> 507 G4double random = G4UniformRand();
>> 508 for(G4int i=0; i < nElements; i++) {
>> 509 if(i!=0) runningSum[i] += runningSum[i-1];
>> 510 if(random <= runningSum[i]/crossSectionTotal) {
>> 511 ZZ = ((*theElementVector)[i])->GetZ();
>> 512
>> 513 // Compare random number to running sum over isotope xc to choose A
>> 514
>> 515 G4int nIso = awicsPerElement[i].size();
>> 516 G4double* running = new G4double[nIso];
>> 517 for (G4int j=0; j < nIso; j++) {
>> 518 running[j] = awicsPerElement[i][j];
>> 519 if(j!=0) running[j] += running[j-1];
>> 520 }
>> 521
>> 522 G4double trial = G4UniformRand();
>> 523 for (G4int j=0; j < nIso; j++) {
>> 524 AA = AvaluesPerElement[i][j];
>> 525 if (trial <= running[j]/running[nIso-1]) break;
>> 526 }
>> 527 delete [] running;
>> 528 break;
>> 529 }
>> 530 }
>> 531 }
>> 532 return std::pair<G4double, G4double>(ZZ, AA);
385 } 533 }
>> 534 */
386 535
387 //....oooOO0OOooo........oooOO0OOooo........oo << 536 void
388 << 537 G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* aDataSet)
389 G4String G4CrossSectionDataStore::HtmlFileName <<
390 { 538 {
391 G4String str(in); << 539 DataSetList.push_back(aDataSet);
392 // replace blanks by _ C++11 version: << 540 NDataSetList++;
393 std::transform(str.begin(), str.end(), str. <<
394 return ch == ' ' ? '_' : ch; <<
395 }); <<
396 str=str + ".html"; <<
397 return str; <<
398 } 541 }
399 542
400 //....oooOO0OOooo........oooOO0OOooo........oo << 543 void
401 << 544 G4CrossSectionDataStore::
402 void G4CrossSectionDataStore::AddDataSet(G4VCr << 545 BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
403 { 546 {
404 if(p->ForAllAtomsAndEnergies()) { << 547 if(NDataSetList > 0) {
405 dataSetList.clear(); << 548 for (G4int i=0; i<NDataSetList; i++) {
406 nDataSetList = 0; << 549 DataSetList[i]->BuildPhysicsTable(aParticleType);
>> 550 }
407 } 551 }
408 dataSetList.push_back(p); <<
409 ++nDataSetList; <<
410 } 552 }
411 553
412 //....oooOO0OOooo........oooOO0OOooo........oo <<
413 554
414 void G4CrossSectionDataStore::AddDataSet(G4VCr << 555 void
>> 556 G4CrossSectionDataStore::
>> 557 DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
415 { 558 {
416 if(p->ForAllAtomsAndEnergies()) { << 559 if (NDataSetList == 0) {
417 dataSetList.clear(); << 560 G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: "
418 dataSetList.push_back(p); << 561 << " no data sets registered"<<G4endl;
419 nDataSetList = 1; << 562 return;
420 } else if ( i >= dataSetList.size() ) { << 563 }
421 dataSetList.push_back(p); << 564 for (G4int i = NDataSetList-1; i >= 0; i--) {
422 ++nDataSetList; << 565 DataSetList[i]->DumpPhysicsTable(aParticleType);
423 } else { <<
424 std::vector< G4VCrossSectionDataSet* >::it <<
425 dataSetList.insert(it , p); <<
426 ++nDataSetList; <<
427 } 566 }
428 } 567 }
429 <<
430 //....oooOO0OOooo........oooOO0OOooo........oo <<
431 568