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