| 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 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 070523 Try to limit sum of secondary photon << 30 // 070523 Try to limit sum of secondary photon energy while keeping distribution shape
31 // in the of nDiscrete = 1 an nPartial << 31 // in the of nDiscrete = 1 an nPartial = 1. Most case are satisfied.
32 // T. Koi 32 // T. Koi
33 // 070606 Add Partial case by T. Koi << 33 // 070606 Add Partial case by T. Koi
34 // 070618 fix memory leaking by T. Koi 34 // 070618 fix memory leaking by T. Koi
35 // 080801 fix memory leaking by T. Koi 35 // 080801 fix memory leaking by T. Koi
36 // 080801 Correcting data disorder which happe << 36 // 080801 Correcting data disorder which happened when both InitPartial
37 // and InitAnglurar methods was called 37 // and InitAnglurar methods was called in a same instance by T. Koi
38 // 090514 Fix bug in IC electron emission case << 38 // 090514 Fix bug in IC electron emission case
39 // Contribution from Chao Zhang (Chao.Z 39 // Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
40 // But it looks like never cause real e 40 // But it looks like never cause real effect in G4NDL3.13 (at least Natural elements) TK
41 // 101111 Change warning message for "repFlag << 41 // 101111 Change warning message for "repFlag == 2 && isoFlag != 1" case
42 // 42 //
43 // there is a lot of unused (and undebugged) c 43 // there is a lot of unused (and undebugged) code in this file. Kept for the moment just in case. @@
44 // P. Arce, June-2014 Conversion neutron_hp to 44 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
45 // 45 //
46 #include "G4ParticleHPPhotonDist.hh" << 46 #include <numeric>
47 47
48 #include "G4Electron.hh" << 48 #include "G4ParticleHPPhotonDist.hh"
49 #include "G4ParticleHPLegendreStore.hh" <<
50 #include "G4PhysicalConstants.hh" 49 #include "G4PhysicalConstants.hh"
51 #include "G4Poisson.hh" <<
52 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
>> 51 #include "G4ParticleHPLegendreStore.hh"
>> 52 #include "G4Electron.hh"
>> 53 #include "G4Poisson.hh"
53 54
54 #include <numeric> << 55 G4bool G4ParticleHPPhotonDist::InitMean(std::istream & aDataFile)
55 <<
56 G4bool G4ParticleHPPhotonDist::InitMean(std::i <<
57 { 56 {
58 G4bool result = true; 57 G4bool result = true;
59 if (aDataFile >> repFlag) { << 58 if(aDataFile >> repFlag)
>> 59 {
>> 60
60 aDataFile >> targetMass; 61 aDataFile >> targetMass;
61 if (repFlag == 1) { << 62 if(repFlag==1)
>> 63 {
62 // multiplicities 64 // multiplicities
63 aDataFile >> nDiscrete; 65 aDataFile >> nDiscrete;
64 const std::size_t msize = nDiscrete > 0 << 66 disType = new G4int[nDiscrete];
65 disType = new G4int[msize]; << 67 energy = new G4double[nDiscrete];
66 energy = new G4double[msize]; << 68 //actualMult = new G4int[nDiscrete];
67 // actualMult = new G4int[msize]; << 69 theYield = new G4ParticleHPVector[nDiscrete];
68 theYield = new G4ParticleHPVector[msize] << 70 for (G4int i=0; i<nDiscrete; ++i)
69 for (std::size_t i = 0; i < msize; ++i) << 71 {
70 aDataFile >> disType[i] >> energy[i]; << 72 aDataFile >> disType[i]>>energy[i];
71 energy[i] *= eV; << 73 energy[i]*=eV;
72 theYield[i].Init(aDataFile, eV); << 74 theYield[i].Init(aDataFile, eV);
73 } <<
74 } <<
75 else if (repFlag == 2) { <<
76 aDataFile >> theInternalConversionFlag; <<
77 aDataFile >> theBaseEnergy; <<
78 theBaseEnergy *= eV; <<
79 aDataFile >> theInternalConversionFlag; <<
80 aDataFile >> nGammaEnergies; <<
81 const std::size_t esize = nGammaEnergies <<
82 theLevelEnergies = new G4double[esize]; <<
83 theTransitionProbabilities = new G4doubl <<
84 if (theInternalConversionFlag == 2) <<
85 thePhotonTransitionFraction = new G4do <<
86 for (std::size_t ii = 0; ii < esize; ++i <<
87 if (theInternalConversionFlag == 1) { <<
88 aDataFile >> theLevelEnergies[ii] >> <<
89 theLevelEnergies[ii] *= eV; <<
90 } <<
91 else if (theInternalConversionFlag == <<
92 aDataFile >> theLevelEnergies[ii] >> <<
93 >> thePhotonTransitionFraction[ii] <<
94 theLevelEnergies[ii] *= eV; <<
95 } <<
96 else { <<
97 throw G4HadronicException(__FILE__, <<
98 "G4Particl <<
99 } <<
100 } 75 }
101 } 76 }
102 else { << 77 else if(repFlag == 2)
103 G4cout << "Data representation in G4Part << 78 {
104 throw G4HadronicException( << 79 aDataFile >> theInternalConversionFlag;
105 __FILE__, __LINE__, "G4ParticleHPPhoto << 80 aDataFile >> theBaseEnergy;
>> 81 theBaseEnergy*=eV;
>> 82 aDataFile >> theInternalConversionFlag;
>> 83 aDataFile >> nGammaEnergies;
>> 84 theLevelEnergies = new G4double[nGammaEnergies];
>> 85 theTransitionProbabilities = new G4double[nGammaEnergies];
>> 86 if(theInternalConversionFlag == 2) thePhotonTransitionFraction = new G4double[nGammaEnergies];
>> 87 for(G4int ii=0; ii<nGammaEnergies; ++ii)
>> 88 {
>> 89 if(theInternalConversionFlag == 1)
>> 90 {
>> 91 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii];
>> 92 theLevelEnergies[ii]*=eV;
>> 93 }
>> 94 else if(theInternalConversionFlag == 2)
>> 95 {
>> 96 aDataFile >> theLevelEnergies[ii] >> theTransitionProbabilities[ii] >> thePhotonTransitionFraction[ii];
>> 97 theLevelEnergies[ii]*=eV;
>> 98 }
>> 99 else
>> 100 {
>> 101 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPPhotonDist: Unknown conversion flag");
>> 102 }
>> 103 }
>> 104 }
>> 105 else
>> 106 {
>> 107 G4cout << "Data representation in G4ParticleHPPhotonDist: "<<repFlag<<G4endl;
>> 108 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPPhotonDist: This data representation is not implemented.");
106 } 109 }
107 } 110 }
108 else { << 111 else
>> 112 {
109 result = false; 113 result = false;
110 } 114 }
111 return result; 115 return result;
112 } 116 }
113 117
114 void G4ParticleHPPhotonDist::InitAngular(std:: << 118 void G4ParticleHPPhotonDist::InitAngular(std::istream & aDataFile)
115 { 119 {
116 G4int i, ii; 120 G4int i, ii;
117 // angular distributions 121 // angular distributions
118 aDataFile >> isoFlag; 122 aDataFile >> isoFlag;
119 if (isoFlag != 1) { << 123 if (isoFlag != 1)
120 if (repFlag == 2) << 124 {
121 G4cout << "G4ParticleHPPhotonDist: repFl << 125 if (repFlag == 2) G4cout << "G4ParticleHPPhotonDist: repFlag == 2 && isoFlag != 1 is unexpected! If you use G4ND3.x, then please report to Geant4 HyperNews. " << G4endl;
122 "G4ND3.x, then please report t <<
123 << G4endl; <<
124 aDataFile >> tabulationType >> nDiscrete2 126 aDataFile >> tabulationType >> nDiscrete2 >> nIso;
125 if (theGammas != nullptr && nDiscrete2 != << 127 if (theGammas != NULL && nDiscrete2 != nDiscrete)
126 G4cout << "080731c G4ParticleHPPhotonDis << 128 G4cout << "080731c G4ParticleHPPhotonDist nDiscrete2 != nDiscrete, It looks like something wrong in your NDL files. Please update the latest. If you still have this messages after the update, then please report to Geant4 Hyper News." << G4endl;
127 "wrong in your NDL files. Plea <<
128 "messages after the update, th <<
129 << G4endl; <<
130 129
131 // The order of cross section (InitPartial 130 // The order of cross section (InitPartials) and distribution
132 // (InitAngular here) data are different, 131 // (InitAngular here) data are different, we have to re-coordinate
133 // consistent data order. 132 // consistent data order.
134 std::vector<G4double> vct_gammas_par; << 133 std::vector < G4double > vct_gammas_par;
135 std::vector<G4double> vct_shells_par; << 134 std::vector < G4double > vct_shells_par;
136 std::vector<G4int> vct_primary_par; << 135 std::vector < G4int > vct_primary_par;
137 std::vector<G4int> vct_distype_par; << 136 std::vector < G4int > vct_distype_par;
138 std::vector<G4ParticleHPVector*> vct_pXS_p << 137 std::vector < G4ParticleHPVector* > vct_pXS_par;
139 if (theGammas != nullptr && theShells != n << 138 if ( theGammas != nullptr && theShells != nullptr )
140 // copy the cross section data << 139 {
141 for (i = 0; i < nDiscrete; ++i) { << 140 //copy the cross section data
142 vct_gammas_par.push_back(theGammas[i]) << 141 for ( i = 0 ; i < nDiscrete ; ++i )
143 vct_shells_par.push_back(theShells[i]) << 142 {
144 vct_primary_par.push_back(isPrimary[i] << 143 vct_gammas_par.push_back( theGammas[ i ] );
145 vct_distype_par.push_back(disType[i]); << 144 vct_shells_par.push_back( theShells[ i ] );
146 auto hpv = new G4ParticleHPVector; << 145 vct_primary_par.push_back( isPrimary[ i ] );
147 *hpv = thePartialXsec[i]; << 146 vct_distype_par.push_back( disType[ i ] );
148 vct_pXS_par.push_back(hpv); << 147 G4ParticleHPVector* hpv = new G4ParticleHPVector;
149 } << 148 *hpv = thePartialXsec[ i ];
150 } << 149 vct_pXS_par.push_back( hpv );
151 const std::size_t psize = nDiscrete2 > 0 ? << 150 }
152 if (theGammas == nullptr) theGammas = new << 151 }
153 if (theShells == nullptr) theShells = new << 152 if ( theGammas == nullptr ) theGammas = new G4double[nDiscrete2];
>> 153 if ( theShells == nullptr ) theShells = new G4double[nDiscrete2];
154 154
155 for (i = 0; i < nIso; ++i) // isotropic p << 155 for (i=0; i< nIso; ++i) // isotropic photons
156 { 156 {
157 aDataFile >> theGammas[i] >> theShells[i 157 aDataFile >> theGammas[i] >> theShells[i];
158 theGammas[i] *= eV; << 158 theGammas[i]*=eV;
159 theShells[i] *= eV; << 159 theShells[i]*=eV;
160 } 160 }
161 const std::size_t tsize = nDiscrete2 - nIs << 161 nNeu = new G4int [nDiscrete2-nIso];
162 nNeu = new G4int[tsize]; << 162 if(tabulationType==1)theLegendre=new G4ParticleHPLegendreTable *[nDiscrete2-nIso];
163 if (tabulationType == 1) theLegendre = new << 163 if(tabulationType==2)theAngular =new G4ParticleHPAngularP *[nDiscrete2-nIso];
164 if (tabulationType == 2) theAngular = new << 164 for(i=nIso; i< nDiscrete2; ++i)
165 for (i = nIso; i < nDiscrete2; ++i) { << 165 {
166 if (tabulationType == 1) { << 166 if(tabulationType==1)
167 aDataFile >> theGammas[i] >> theShells << 167 {
168 theGammas[i] *= eV; << 168 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
169 theShells[i] *= eV; << 169 theGammas[i]*=eV;
170 const std::size_t lsize = nNeu[i - nIs << 170 theShells[i]*=eV;
171 theLegendre[i - nIso] = new G4Particle << 171 theLegendre[i-nIso]=new G4ParticleHPLegendreTable[nNeu[i-nIso]];
172 theLegendreManager.Init(aDataFile); << 172 theLegendreManager.Init(aDataFile);
173 for (ii = 0; ii < nNeu[i - nIso]; ++ii << 173 for (ii=0; ii<nNeu[i-nIso]; ++ii)
174 theLegendre[i - nIso][ii].Init(aData << 174 {
175 } << 175 theLegendre[i-nIso][ii].Init(aDataFile);
176 } << 176 }
177 else if (tabulationType == 2) { << 177 }
178 aDataFile >> theGammas[i] >> theShells << 178 else if(tabulationType==2)
179 theGammas[i] *= eV; << 179 {
180 theShells[i] *= eV; << 180 aDataFile >> theGammas[i] >> theShells[i] >> nNeu[i-nIso];
181 const std::size_t asize = nNeu[i - nIs << 181 theGammas[i]*=eV;
182 theAngular[i - nIso] = new G4ParticleH << 182 theShells[i]*=eV;
183 for (ii = 0; ii < nNeu[i - nIso]; ++ii << 183 theAngular[i-nIso]=new G4ParticleHPAngularP[nNeu[i-nIso]];
184 theAngular[i - nIso][ii].Init(aDataF << 184 for (ii=0; ii<nNeu[i-nIso]; ++ii)
185 } << 185 {
186 } << 186 theAngular[i-nIso][ii].Init(aDataFile);
187 else { << 187 }
188 G4cout << "tabulation type: tabulation << 188 }
189 throw G4HadronicException( << 189 else
190 __FILE__, __LINE__, "cannot deal wit << 190 {
>> 191 G4cout << "tabulation type: tabulationType"<<G4endl;
>> 192 throw G4HadronicException(__FILE__, __LINE__, "cannot deal with this tabulation type for angular distributions.");
191 } 193 }
192 } 194 }
193 195
194 if (!vct_gammas_par.empty()) { << 196 if ( vct_gammas_par.size() > 0 )
195 // Reordering cross section data to corr << 197 {
196 for (i = 0; i < nDiscrete; ++i) { << 198 // Reordering cross section data to corrsponding distribution data
197 for (G4int j = 0; j < nDiscrete; ++j) << 199 for ( i = 0 ; i < nDiscrete ; ++i )
198 // Checking gamma and shell to ident << 200 {
199 if (theGammas[i] == vct_gammas_par[j << 201 for ( G4int j = 0 ; j < nDiscrete ; ++j )
200 isPrimary[i] = vct_primary_par[j]; << 202 {
201 disType[i] = vct_distype_par[j]; << 203 // Checking gamma and shell to identification
202 thePartialXsec[i] = (*(vct_pXS_par << 204 if ( theGammas[ i ] == vct_gammas_par [ j ] && theShells [ i ] == vct_shells_par[ j ] )
>> 205 {
>> 206 isPrimary [ i ] = vct_primary_par [ j ];
>> 207 disType [ i ] = vct_distype_par [ j ];
>> 208 thePartialXsec[ i ] = ( *( vct_pXS_par[ j ] ) );
203 } 209 }
204 } 210 }
205 } 211 }
206 // Garbage collection << 212 // Garbage collection
207 for (auto it = vct_pXS_par.cbegin(); it << 213 for ( auto it = vct_pXS_par.cbegin() ; it != vct_pXS_par.cend() ; ++it )
>> 214 {
208 delete *it; 215 delete *it;
209 } 216 }
210 } 217 }
211 } 218 }
212 } 219 }
213 220
214 void G4ParticleHPPhotonDist::InitEnergies(std: << 221 void G4ParticleHPPhotonDist::InitEnergies(std::istream & aDataFile)
215 { 222 {
216 G4int i, energyDistributionsNeeded = 0; 223 G4int i, energyDistributionsNeeded = 0;
217 for (i = 0; i < nDiscrete; ++i) { << 224 for (i=0; i<nDiscrete; ++i)
218 if (disType[i] == 1) energyDistributionsNe << 225 {
219 } << 226 if( disType[i]==1) energyDistributionsNeeded =1;
220 if (energyDistributionsNeeded == 0) return; << 227 }
221 aDataFile >> nPartials; << 228 if(!energyDistributionsNeeded) return;
222 const std::size_t dsize = nPartials > 0 ? nP << 229 aDataFile >> nPartials;
223 distribution = new G4int[dsize]; << 230 distribution = new G4int[nPartials];
224 probs = new G4ParticleHPVector[dsize]; << 231 probs = new G4ParticleHPVector[nPartials];
225 partials = new G4ParticleHPPartial*[dsize]; << 232 partials = new G4ParticleHPPartial * [nPartials];
226 G4int nen; 233 G4int nen;
227 G4int dummy; 234 G4int dummy;
228 for (i = 0; i < nPartials; ++i) { << 235 for (i=0; i<nPartials; ++i)
>> 236 {
229 aDataFile >> dummy; 237 aDataFile >> dummy;
230 probs[i].Init(aDataFile, eV); 238 probs[i].Init(aDataFile, eV);
231 aDataFile >> nen; 239 aDataFile >> nen;
232 partials[i] = new G4ParticleHPPartial(nen) 240 partials[i] = new G4ParticleHPPartial(nen);
233 partials[i]->InitInterpolation(aDataFile); 241 partials[i]->InitInterpolation(aDataFile);
234 partials[i]->Init(aDataFile); 242 partials[i]->Init(aDataFile);
235 } 243 }
236 } 244 }
237 245
238 void G4ParticleHPPhotonDist::InitPartials(std: << 246 void G4ParticleHPPhotonDist::InitPartials(std::istream& aDataFile,
>> 247 G4ParticleHPVector* theXsec)
239 { 248 {
240 if (theXsec != nullptr) theReactionXsec = th << 249 if (theXsec) theReactionXsec = theXsec;
241 250
242 aDataFile >> nDiscrete >> targetMass; 251 aDataFile >> nDiscrete >> targetMass;
243 if (nDiscrete != 1) { << 252 if(nDiscrete != 1)
>> 253 {
244 theTotalXsec.Init(aDataFile, eV); 254 theTotalXsec.Init(aDataFile, eV);
245 } 255 }
246 const std::size_t dsize = nDiscrete > 0 ? nD << 256 G4int i;
247 theGammas = new G4double[dsize]; << 257 theGammas = new G4double[nDiscrete];
248 theShells = new G4double[dsize]; << 258 theShells = new G4double[nDiscrete];
249 isPrimary = new G4int[dsize]; << 259 isPrimary = new G4int[nDiscrete];
250 disType = new G4int[dsize]; << 260 disType = new G4int[nDiscrete];
251 thePartialXsec = new G4ParticleHPVector[dsiz << 261 thePartialXsec = new G4ParticleHPVector[nDiscrete];
252 for (std::size_t i = 0; i < dsize; ++i) { << 262 for(i=0; i<nDiscrete; ++i)
253 aDataFile >> theGammas[i] >> theShells[i] << 263 {
254 theGammas[i] *= eV; << 264 aDataFile>>theGammas[i]>>theShells[i]>>isPrimary[i]>>disType[i];
255 theShells[i] *= eV; << 265 theGammas[i]*=eV;
>> 266 theShells[i]*=eV;
256 thePartialXsec[i].Init(aDataFile, eV); 267 thePartialXsec[i].Init(aDataFile, eV);
257 } << 268 }
258 } 269 }
259 270
260 G4ReactionProductVector* G4ParticleHPPhotonDis << 271 G4ReactionProductVector * G4ParticleHPPhotonDist::GetPhotons(G4double anEnergy)
261 { 272 {
262 // the partial cross-section case is not all 273 // the partial cross-section case is not all in this yet.
263 if (actualMult.Get() == nullptr) { << 274 if ( actualMult.Get() == nullptr ) {
264 actualMult.Get() = new std::vector<G4int>( << 275 actualMult.Get() = new std::vector<G4int>( nDiscrete );
265 } 276 }
266 G4int i, ii, iii; 277 G4int i, ii, iii;
267 G4int nSecondaries = 0; 278 G4int nSecondaries = 0;
268 auto thePhotons = new G4ReactionProductVecto << 279 G4ReactionProductVector* thePhotons = new G4ReactionProductVector;
269 280
270 if (repFlag == 1) { << 281 if (repFlag==1) {
271 G4double current = 0; << 282 G4double current=0;
272 for (i = 0; i < nDiscrete; ++i) { 283 for (i = 0; i < nDiscrete; ++i) {
273 current = theYield[i].GetY(anEnergy); 284 current = theYield[i].GetY(anEnergy);
274 actualMult.Get()->at(i) = (G4int)G4Poiss << 285 actualMult.Get()->at(i) = (G4int)G4Poisson(current); // max cut-off still missing @@@
275 if (nDiscrete == 1 && current < 1.0001) 286 if (nDiscrete == 1 && current < 1.0001) {
276 actualMult.Get()->at(i) = static_cast< 287 actualMult.Get()->at(i) = static_cast<G4int>(current);
277 if (current < 1) { << 288 if(current<1)
>> 289 {
278 actualMult.Get()->at(i) = 0; 290 actualMult.Get()->at(i) = 0;
279 if (G4UniformRand() < current) actua << 291 if(G4UniformRand()<current) actualMult.Get()->at(i) = 1;
280 } 292 }
281 } 293 }
282 nSecondaries += actualMult.Get()->at(i); 294 nSecondaries += actualMult.Get()->at(i);
283 } 295 }
284 for (i = 0; i < nSecondaries; ++i) { 296 for (i = 0; i < nSecondaries; ++i) {
285 auto theOne = new G4ReactionProduct; << 297 G4ReactionProduct * theOne = new G4ReactionProduct;
286 theOne->SetDefinition(G4Gamma::Gamma()); 298 theOne->SetDefinition(G4Gamma::Gamma());
287 thePhotons->push_back(theOne); 299 thePhotons->push_back(theOne);
288 } 300 }
289 301
290 G4int count = 0; 302 G4int count = 0;
291 303
292 if (nDiscrete == 1 && nPartials == 1) { 304 if (nDiscrete == 1 && nPartials == 1) {
293 if (actualMult.Get()->at(0) > 0) { 305 if (actualMult.Get()->at(0) > 0) {
294 if (disType[0] == 1) { 306 if (disType[0] == 1) {
295 // continuum 307 // continuum
296 G4ParticleHPVector* temp; 308 G4ParticleHPVector* temp;
297 temp = partials[0]->GetY(anEnergy); << 309 temp = partials[ 0 ]->GetY(anEnergy); //@@@ look at, seems fishy
298 G4double maximumE = temp->GetX(temp- << 310 G4double maximumE = temp->GetX( temp->GetVectorLength()-1 ); // This is an assumption.
299 311
300 std::vector<G4double> photons_e_best << 312 std::vector< G4double > photons_e_best( actualMult.Get()->at(0) , 0.0 );
301 G4double best = DBL_MAX; 313 G4double best = DBL_MAX;
302 G4int maxTry = 1000; << 314 G4int maxTry = 1000;
303 for (G4int j = 0; j < maxTry; ++j) { 315 for (G4int j = 0; j < maxTry; ++j) {
304 std::vector<G4double> photons_e(ac 316 std::vector<G4double> photons_e(actualMult.Get()->at(0), 0.0);
305 for (auto it = photons_e.begin(); 317 for (auto it = photons_e.begin(); it < photons_e.end(); ++it) {
306 *it = temp->Sample(); 318 *it = temp->Sample();
307 } 319 }
308 320
309 if (std::accumulate(photons_e.cbeg 321 if (std::accumulate(photons_e.cbegin(), photons_e.cend(), 0.0) > maximumE) {
310 if (std::accumulate(photons_e.cb 322 if (std::accumulate(photons_e.cbegin(), photons_e.cend(), 0.0) < best)
311 photons_e_best = std::move(pho << 323 photons_e_best = photons_e;
312 continue; 324 continue;
313 } <<
314 G4int iphot = 0; <<
315 for (auto it = photons_e.cbegin(); <<
316 thePhotons->operator[](iphot)->S <<
317 *it); // Replace index count, <<
318 // with iphot, which is <<
319 // bug report 2167 <<
320 ++iphot; <<
321 } <<
322 325
323 break; << 326 } else {
>> 327 G4int iphot = 0;
>> 328 for (auto it = photons_e.cbegin(); it < photons_e.cend(); ++it) {
>> 329 thePhotons->operator[](iphot)->SetKineticEnergy(*it); // Replace index count, which was not incremented,
>> 330 // with iphot, which is, as per Artem Zontikov,
>> 331 // bug report 2167
>> 332 ++iphot;
>> 333 }
>> 334
>> 335 break;
>> 336 }
324 } 337 }
325 delete temp; 338 delete temp;
326 } << 339
327 else { << 340 } else {
328 // discrete 341 // discrete
329 thePhotons->operator[](count)->SetKi 342 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
330 } 343 }
331 ++count; 344 ++count;
332 if (count > nSecondaries) << 345 if (count > nSecondaries)
333 throw G4HadronicException(__FILE__, << 346 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPPhotonDist::GetPhotons inconsistency");
334 "G4Particl <<
335 } 347 }
336 } << 348
337 else { // nDiscrete != 1 or nPartials != << 349 } else { // nDiscrete != 1 or nPartials != 1
338 for (i = 0; i < nDiscrete; ++i) { << 350 for (i=0; i<nDiscrete; ++i) {
339 for (ii = 0; ii < actualMult.Get()->at << 351 for (ii=0; ii< actualMult.Get()->at(i); ++ii) {
340 if (disType[i] == 1) { 352 if (disType[i] == 1) {
341 // continuum 353 // continuum
342 G4double sum = 0, run = 0; << 354 G4double sum=0, run=0;
343 for (iii = 0; iii < nPartials; ++i 355 for (iii = 0; iii < nPartials; ++iii)
344 sum += probs[iii].GetY(anEnergy) << 356 sum+=probs[iii].GetY(anEnergy);
345 G4double random = G4UniformRand(); 357 G4double random = G4UniformRand();
346 G4int theP = 0; 358 G4int theP = 0;
347 for (iii = 0; iii < nPartials; ++i 359 for (iii = 0; iii < nPartials; ++iii) {
348 run += probs[iii].GetY(anEnergy) << 360 run+=probs[iii].GetY(anEnergy);
349 theP = iii; 361 theP = iii;
350 if (random < run / sum) break; << 362 if(random<run/sum) break;
351 } 363 }
352 364
353 if (theP == nPartials) theP = nPar << 365 if (theP == nPartials) theP=nPartials-1; // das sortiert J aus.
354 sum = 0; << 366 sum = 0;
355 G4ParticleHPVector* temp; << 367 G4ParticleHPVector * temp;
356 temp = partials[theP]->GetY(anEner << 368 temp = partials[theP]->GetY(anEnergy); //@@@ look at, seems fishy
357 G4double eGamm = temp->Sample(); 369 G4double eGamm = temp->Sample();
358 thePhotons->operator[](count)->Set 370 thePhotons->operator[](count)->SetKineticEnergy(eGamm);
359 delete temp; 371 delete temp;
360 } << 372
361 else { << 373 } else {
362 // discrete 374 // discrete
363 thePhotons->operator[](count)->Set 375 thePhotons->operator[](count)->SetKineticEnergy(energy[i]);
364 } << 376 }
365 ++count; 377 ++count;
366 if (count > nSecondaries) 378 if (count > nSecondaries)
367 throw G4HadronicException(__FILE__ << 379 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPPhotonDist::GetPhotons inconsistency");
368 "G4Parti <<
369 } 380 }
370 } 381 }
371 } 382 }
372 383
373 // now do the angular distributions... 384 // now do the angular distributions...
374 if (isoFlag == 1) { 385 if (isoFlag == 1) {
375 for (i = 0; i < nSecondaries; ++i) { << 386 for (i=0; i< nSecondaries; ++i)
376 G4double costheta = 2. * G4UniformRand << 387 {
377 G4double theta = std::acos(costheta); << 388 G4double costheta = 2.*G4UniformRand()-1;
378 G4double phi = twopi * G4UniformRand() << 389 G4double theta = std::acos(costheta);
379 G4double sinth = std::sin(theta); << 390 G4double phi = twopi*G4UniformRand();
380 G4double en = thePhotons->operator[](i << 391 G4double sinth = std::sin(theta);
381 G4ThreeVector temp(en * sinth * std::c << 392 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
382 en * std::cos(theta << 393 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
383 thePhotons->operator[](i)->SetMomentum << 394 thePhotons->operator[](i)->SetMomentum( temp ) ;
384 } << 395 }
385 } << 396 }
386 else { << 397 else
387 for (i = 0; i < nSecondaries; ++i) { << 398 {
388 G4double currentEnergy = thePhotons->o << 399 for(i=0; i<nSecondaries; ++i)
389 for (ii = 0; ii < nDiscrete2; ++ii) { << 400 {
390 if (std::abs(currentEnergy - theGamm << 401 G4double currentEnergy = thePhotons->operator[](i)->GetTotalEnergy();
391 } << 402 for(ii=0; ii<nDiscrete2; ++ii)
392 if (ii == nDiscrete2) << 403 {
393 --ii; // fix for what seems an (fil << 404 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
394 // data. @@ << 405 }
395 if (ii < nIso) { << 406 if(ii==nDiscrete2) --ii; // fix for what seems an (file12 vs file 14) inconsistency found in the ENDF 7N14 data. @@
>> 407 if(ii<nIso)
>> 408 {
396 // isotropic distribution 409 // isotropic distribution
397 // 410 //
398 // Fix Bugzilla report #1745 << 411 //Fix Bugzilla report #1745
399 // G4double theta = pi*G4UniformRand << 412 //G4double theta = pi*G4UniformRand();
400 G4double costheta = 2. * G4UniformRa << 413 G4double costheta = 2.*G4UniformRand()-1;
401 G4double theta = std::acos(costheta) << 414 G4double theta = std::acos(costheta);
402 G4double phi = twopi * G4UniformRand << 415 G4double phi = twopi*G4UniformRand();
403 G4double sinth = std::sin(theta); 416 G4double sinth = std::sin(theta);
404 G4double en = thePhotons->operator[] 417 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
405 // DHW G4ThreeVector tempVector(en*s << 418 // DHW G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
406 // en*std::cos(theta) ); << 419 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costheta );
407 G4ThreeVector tempVector(en * sinth << 420 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
408 en * costhe << 421 }
409 thePhotons->operator[](i)->SetMoment << 422 else if(tabulationType==1)
410 } << 423 {
411 else if (tabulationType == 1) { <<
412 // legendre polynomials 424 // legendre polynomials
413 G4int it(0); 425 G4int it(0);
414 for (iii = 0; iii < nNeu[ii - nIso]; << 426 for (iii=0; iii<nNeu[ii-nIso]; ++iii) // find the neutron energy
415 { 427 {
416 it = iii; 428 it = iii;
417 if (theLegendre[ii - nIso][iii].Ge << 429 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
>> 430 break;
418 } 431 }
419 G4ParticleHPLegendreStore aStore(2); 432 G4ParticleHPLegendreStore aStore(2);
420 aStore.SetCoeff(1, &(theLegendre[ii << 433 aStore.SetCoeff(1, &(theLegendre[ii-nIso][it]));
421 if (it > 0) { << 434 if ( it > 0 )
422 aStore.SetCoeff(0, &(theLegendre[i << 435 {
>> 436 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
423 } 437 }
424 else { << 438 else
425 aStore.SetCoeff(0, &(theLegendre[i << 439 {
>> 440 aStore.SetCoeff(0, &(theLegendre[ii-nIso][it]));
426 } 441 }
427 G4double cosTh = aStore.SampleMax(an 442 G4double cosTh = aStore.SampleMax(anEnergy);
428 G4double theta = std::acos(cosTh); 443 G4double theta = std::acos(cosTh);
429 G4double phi = twopi * G4UniformRand << 444 G4double phi = twopi*G4UniformRand();
430 G4double sinth = std::sin(theta); 445 G4double sinth = std::sin(theta);
431 G4double en = thePhotons->operator[] 446 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
432 G4ThreeVector tempVector(en * sinth << 447 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
433 en * std::c << 448 thePhotons->operator[](i)->SetMomentum( tempVector ) ;
434 thePhotons->operator[](i)->SetMoment << 449 }
435 } << 450 else
436 else { << 451 {
437 // tabulation of probabilities. 452 // tabulation of probabilities.
438 G4int it(0); 453 G4int it(0);
439 for (iii = 0; iii < nNeu[ii - nIso]; << 454 for (iii=0; iii<nNeu[ii-nIso]; ++iii) // find the neutron energy
440 { 455 {
441 it = iii; 456 it = iii;
442 if (theAngular[ii - nIso][iii].Get << 457 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
>> 458 break;
443 } 459 }
444 G4double costh = theAngular[ii - nIs << 460 G4double costh = theAngular[ii-nIso][it].GetCosTh(); // no interpolation yet @@
445 G4double theta = std::acos(costh); 461 G4double theta = std::acos(costh);
446 G4double phi = twopi * G4UniformRand << 462 G4double phi = twopi*G4UniformRand();
447 G4double sinth = std::sin(theta); 463 G4double sinth = std::sin(theta);
448 G4double en = thePhotons->operator[] 464 G4double en = thePhotons->operator[](i)->GetTotalEnergy();
449 G4ThreeVector tmpVector(en * sinth * << 465 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
450 en * costh); << 466 thePhotons->operator[](i)->SetMomentum( tmpVector ) ;
451 thePhotons->operator[](i)->SetMoment << 467 }
452 } << 468 }
453 } << 469 }
454 } << 470
455 } << 471 } else if (repFlag == 2) {
456 else if (repFlag == 2) { << 472 G4double * running = new G4double[nGammaEnergies];
457 auto running = new G4double[nGammaEnergies << 473 running[0]=theTransitionProbabilities[0];
458 running[0] = theTransitionProbabilities[0] << 474 for(i=1; i<nGammaEnergies; ++i)
459 for (i = 1; i < nGammaEnergies; ++i) { << 475 {
460 running[i] = running[i - 1] + theTransit << 476 running[i]=running[i-1]+theTransitionProbabilities[i];
461 } 477 }
462 G4double random = G4UniformRand(); 478 G4double random = G4UniformRand();
463 G4int it = 0; << 479 G4int it=0;
464 for (i = 0; i < nGammaEnergies; ++i) { << 480 for(i=0; i<nGammaEnergies; ++i)
>> 481 {
465 it = i; 482 it = i;
466 if (random < running[i] / running[nGamma << 483 if(random < running[i]/running[nGammaEnergies-1]) break;
467 } 484 }
468 delete[] running; << 485 delete [] running;
469 G4double totalEnergy = theBaseEnergy - the 486 G4double totalEnergy = theBaseEnergy - theLevelEnergies[it];
470 auto theOne = new G4ReactionProduct; << 487 G4ReactionProduct * theOne = new G4ReactionProduct;
471 theOne->SetDefinition(G4Gamma::Gamma()); 488 theOne->SetDefinition(G4Gamma::Gamma());
472 random = G4UniformRand(); 489 random = G4UniformRand();
473 if (theInternalConversionFlag == 2 && rand << 490 if(theInternalConversionFlag==2 && random>thePhotonTransitionFraction[it])
>> 491 {
474 theOne->SetDefinition(G4Electron::Electr 492 theOne->SetDefinition(G4Electron::Electron());
475 // Bug reported Chao Zhang (Chao.Zhang@u << 493 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
476 // 2009 But never enter at least with G4 << 494 //But never enter at least with G4NDL3.13
477 totalEnergy += << 495 totalEnergy += G4Electron::Electron()->GetPDGMass(); // proposed correction: add this line for electron
478 G4Electron::Electron()->GetPDGMass(); <<
479 } 496 }
480 theOne->SetTotalEnergy(totalEnergy); 497 theOne->SetTotalEnergy(totalEnergy);
481 if (isoFlag == 1) { << 498 if( isoFlag == 1 )
482 G4double costheta = 2. * G4UniformRand() << 499 {
>> 500 G4double costheta = 2.*G4UniformRand()-1;
483 G4double theta = std::acos(costheta); 501 G4double theta = std::acos(costheta);
484 G4double phi = twopi * G4UniformRand(); << 502 G4double phi = twopi*G4UniformRand();
485 G4double sinth = std::sin(theta); 503 G4double sinth = std::sin(theta);
486 // Bug reported Chao Zhang (Chao.Zhang@u << 504 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
487 // 2009 G4double en = theOne->GetTotalEn << 505 //G4double en = theOne->GetTotalEnergy();
488 G4double en = theOne->GetTotalMomentum() 506 G4double en = theOne->GetTotalMomentum();
489 // But never cause real effect at least << 507 //But never cause real effect at least with G4NDL3.13 TK
490 G4ThreeVector temp(en * sinth * std::cos << 508 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
491 en * std::cos(theta)) << 509 theOne->SetMomentum( temp ) ;
492 theOne->SetMomentum(temp); <<
493 } 510 }
494 else { << 511 else
>> 512 {
495 G4double currentEnergy = theOne->GetTota 513 G4double currentEnergy = theOne->GetTotalEnergy();
496 for (ii = 0; ii < nDiscrete2; ++ii) { << 514 for(ii=0; ii<nDiscrete2; ++ii)
497 if (std::abs(currentEnergy - theGammas << 515 {
>> 516 if (std::abs(currentEnergy-theGammas[ii])<0.1*keV) break;
498 } 517 }
499 if (ii == nDiscrete2) << 518 if(ii==nDiscrete2) --ii; // fix for what seems an (file12 vs file 14) inconsistency found in the ENDF 7N14 data. @@
500 --ii; // fix for what seems an (file1 << 519 if(ii<nIso)
501 // data. @@ << 520 {
502 if (ii < nIso) { << 521 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
503 // Bug reported Chao Zhang (Chao.Zhang << 522 // isotropic distribution
504 // 2009 << 523 //G4double theta = pi*G4UniformRand();
505 // isotropic distribution << 524 G4double theta = std::acos(2.*G4UniformRand()-1.);
506 // G4double theta = pi*G4UniformRand() << 525 //But this is alos never cause real effect at least with G4NDL3.13 TK not repFlag == 2 AND isoFlag != 1
507 G4double theta = std::acos(2. * G4Unif << 526 G4double phi = twopi*G4UniformRand();
508 // But this is alos never cause real e <<
509 // isoFlag != 1 <<
510 G4double phi = twopi * G4UniformRand() <<
511 G4double sinth = std::sin(theta); 527 G4double sinth = std::sin(theta);
512 // Bug reported Chao Zhang (Chao.Zhang << 528 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
513 // 2009 G4double en = theOne->GetTotal << 529 //G4double en = theOne->GetTotalEnergy();
514 G4double en = theOne->GetTotalMomentum 530 G4double en = theOne->GetTotalMomentum();
515 // But never cause real effect at leas << 531 //But never cause real effect at least with G4NDL3.13 TK
516 G4ThreeVector tempVector(en * sinth * << 532 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
517 en * std::cos << 533 theOne->SetMomentum( tempVector ) ;
518 theOne->SetMomentum(tempVector); <<
519 } 534 }
520 else if (tabulationType == 1) { << 535 else if(tabulationType==1)
>> 536 {
521 // legendre polynomials 537 // legendre polynomials
522 G4int itt(0); 538 G4int itt(0);
523 for (iii = 0; iii < nNeu[ii - nIso]; + << 539 for (iii=0; iii<nNeu[ii-nIso]; ++iii) // find the neutron energy
524 { 540 {
525 itt = iii; 541 itt = iii;
526 if (theLegendre[ii - nIso][iii].GetE << 542 if(theLegendre[ii-nIso][iii].GetEnergy()>anEnergy)
>> 543 break;
527 } 544 }
528 G4ParticleHPLegendreStore aStore(2); 545 G4ParticleHPLegendreStore aStore(2);
529 aStore.SetCoeff(1, &(theLegendre[ii - << 546 aStore.SetCoeff(1, &(theLegendre[ii-nIso][itt]));
530 // aStore.SetCoeff(0, &(theLegendre[ii << 547 //aStore.SetCoeff(0, &(theLegendre[ii-nIso][it-1]));
531 // TKDB 110512 << 548 //TKDB 110512
532 if (itt > 0) { << 549 if ( itt > 0 )
533 aStore.SetCoeff(0, &(theLegendre[ii << 550 {
>> 551 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt-1]));
534 } 552 }
535 else { << 553 else
536 aStore.SetCoeff(0, &(theLegendre[ii << 554 {
>> 555 aStore.SetCoeff(0, &(theLegendre[ii-nIso][itt]));
537 } 556 }
538 G4double cosTh = aStore.SampleMax(anEn 557 G4double cosTh = aStore.SampleMax(anEnergy);
539 G4double theta = std::acos(cosTh); 558 G4double theta = std::acos(cosTh);
540 G4double phi = twopi * G4UniformRand() << 559 G4double phi = twopi*G4UniformRand();
541 G4double sinth = std::sin(theta); 560 G4double sinth = std::sin(theta);
542 // Bug reported Chao Zhang (Chao.Zhang << 561 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
543 // 2009 G4double en = theOne->GetTotal << 562 //G4double en = theOne->GetTotalEnergy();
544 G4double en = theOne->GetTotalMomentum 563 G4double en = theOne->GetTotalMomentum();
545 // But never cause real effect at leas << 564 //But never cause real effect at least with G4NDL3.13 TK
546 G4ThreeVector tempVector(en * sinth * << 565 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
547 en * std::cos << 566 theOne->SetMomentum( tempVector ) ;
548 theOne->SetMomentum(tempVector); <<
549 } 567 }
550 else { << 568 else
>> 569 {
551 // tabulation of probabilities. 570 // tabulation of probabilities.
552 G4int itt(0); 571 G4int itt(0);
553 for (iii = 0; iii < nNeu[ii - nIso]; + << 572 for (iii=0; iii<nNeu[ii-nIso]; ++iii) // find the neutron energy
554 { 573 {
555 itt = iii; 574 itt = iii;
556 if (theAngular[ii - nIso][iii].GetEn << 575 if(theAngular[ii-nIso][iii].GetEnergy()>anEnergy)
>> 576 break;
557 } 577 }
558 G4double costh = theAngular[ii - nIso] << 578 G4double costh = theAngular[ii-nIso][itt].GetCosTh(); // no interpolation yet @@
559 G4double theta = std::acos(costh); 579 G4double theta = std::acos(costh);
560 G4double phi = twopi * G4UniformRand() << 580 G4double phi = twopi*G4UniformRand();
561 G4double sinth = std::sin(theta); 581 G4double sinth = std::sin(theta);
562 // Bug reported Chao Zhang (Chao.Zhang << 582 //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009
563 // 2009 G4double en = theOne->GetTotal << 583 //G4double en = theOne->GetTotalEnergy();
564 G4double en = theOne->GetTotalMomentum 584 G4double en = theOne->GetTotalMomentum();
565 // But never cause real effect at leas << 585 //But never cause real effect at least with G4NDL3.13 TK
566 G4ThreeVector tmpVector(en * sinth * s << 586 G4ThreeVector tmpVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costh );
567 theOne->SetMomentum(tmpVector); << 587 theOne->SetMomentum( tmpVector ) ;
568 } 588 }
569 } 589 }
570 thePhotons->push_back(theOne); 590 thePhotons->push_back(theOne);
571 } 591 }
572 else if (repFlag == 0) { << 592 else if( repFlag==0 )
573 if (thePartialXsec == nullptr) { << 593 {
574 return thePhotons; << 594 if ( thePartialXsec == 0 )
575 } << 595 {
576 << 596 return thePhotons;
577 // Partial Case << 597 }
578 << 598
579 auto theOne = new G4ReactionProduct; << 599 // Partial Case
580 theOne->SetDefinition(G4Gamma::Gamma()); << 600
581 thePhotons->push_back(theOne); << 601 G4ReactionProduct * theOne = new G4ReactionProduct;
582 << 602 theOne->SetDefinition( G4Gamma::Gamma() );
583 // Energy << 603 thePhotons->push_back( theOne );
584 << 604
585 G4double sum = 0.0; << 605 // Energy
586 std::vector<G4double> dif(nDiscrete, 0.0); << 606
587 for (G4int j = 0; j < nDiscrete; ++j) { << 607 G4double sum = 0.0;
588 G4double x = thePartialXsec[j].GetXsec(a << 608 std::vector < G4double > dif( nDiscrete , 0.0 );
589 if (x > 0) { << 609 for ( G4int j = 0 ; j < nDiscrete ; ++j )
590 sum += x; << 610 {
591 } << 611 G4double x = thePartialXsec[ j ].GetXsec( anEnergy ); // x in barn
592 dif[j] = sum; << 612 if ( x > 0 )
593 } << 613 {
594 << 614 sum += x;
595 G4double rand = G4UniformRand(); << 615 }
596 << 616 dif [ j ] = sum;
597 G4int iphoton = 0; <<
598 for (G4int j = 0; j < nDiscrete; ++j) { <<
599 G4double y = rand * sum; <<
600 if (dif[j] > y) { <<
601 iphoton = j; <<
602 break; <<
603 } 617 }
604 } << 618
>> 619 G4double rand = G4UniformRand();
605 620
606 // Statistically suppress the photon accor << 621 G4int iphoton = 0;
607 // Fix proposed by Artem Zontikov, Bug rep << 622 for ( G4int j = 0 ; j < nDiscrete ; ++j )
608 if (theReactionXsec != nullptr) { <<
609 if (thePartialXsec[iphoton].GetXsec(anEn <<
610 < G4UniformRand()) <<
611 { 623 {
612 delete thePhotons; << 624 G4double y = rand*sum;
613 thePhotons = nullptr; << 625 if ( dif [ j ] > y )
614 return thePhotons; << 626 {
>> 627 iphoton = j;
>> 628 break;
>> 629 }
615 } 630 }
616 } <<
617 631
618 // Angle << 632 // Statistically suppress the photon according to reaction cross section
619 G4double cosTheta = 0.0; // mu << 633 // Fix proposed by Artem Zontikov, Bug report #1824
620 << 634 if (theReactionXsec) {
621 if (isoFlag == 1) { << 635 if (thePartialXsec[iphoton].GetXsec(anEnergy)/theReactionXsec->GetXsec(anEnergy) < G4UniformRand() ) {
622 // Isotropic Case << 636 delete thePhotons;
623 << 637 thePhotons = nullptr;
624 cosTheta = 2. * G4UniformRand() - 1; << 638 return thePhotons;
625 } << 639 }
626 else { <<
627 if (iphoton < nIso) { <<
628 // still Isotropic <<
629 <<
630 cosTheta = 2. * G4UniformRand() - 1; <<
631 } 640 }
632 else { <<
633 if (tabulationType == 1) { <<
634 // Legendre polynomials <<
635 <<
636 G4int iangle = 0; <<
637 for (G4int j = 0; j < nNeu[iphoton - <<
638 iangle = j; <<
639 if (theLegendre[iphoton - nIso][j] <<
640 } <<
641 641
642 G4ParticleHPLegendreStore aStore(2); << 642 // Angle
643 aStore.SetCoeff(1, &(theLegendre[iph << 643 G4double cosTheta = 0.0; // mu
644 aStore.SetCoeff(0, &(theLegendre[iph <<
645 644
646 cosTheta = aStore.SampleMax(anEnergy << 645 if ( isoFlag == 1 )
647 } << 646 {
648 else if (tabulationType == 2) { << 647 // Isotropic Case
649 // tabulation of probabilities. <<
650 648
651 G4int iangle = 0; << 649 cosTheta = 2.*G4UniformRand()-1;
652 for (G4int j = 0; j < nNeu[iphoton - <<
653 iangle = j; <<
654 if (theAngular[iphoton - nIso][j]. <<
655 } <<
656 cosTheta = theAngular[iphoton - nIso <<
657 // no interpolation yet @@ <<
658 } <<
659 } 650 }
660 } << 651 else
661 << 652 {
662 // Set << 653 if ( iphoton < nIso )
663 G4double phi = twopi * G4UniformRand(); << 654 {
664 G4double theta = std::acos(cosTheta); << 655 // still Isotropic
665 G4double sinTheta = std::sin(theta); << 656
>> 657 cosTheta = 2.*G4UniformRand()-1;
>> 658 }
>> 659 else
>> 660 {
>> 661 if ( tabulationType == 1 )
>> 662 {
>> 663 // Legendre polynomials
>> 664
>> 665 G4int iangle = 0;
>> 666 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; ++j )
>> 667 {
>> 668 iangle = j;
>> 669 if ( theLegendre[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
>> 670 }
>> 671
>> 672 G4ParticleHPLegendreStore aStore( 2 );
>> 673 aStore.SetCoeff( 1 , &( theLegendre[ iphoton - nIso ][ iangle ] ) );
>> 674 aStore.SetCoeff( 0 , &( theLegendre[ iphoton - nIso ][ iangle - 1 ] ) );
666 675
667 G4double photonE = theGammas[iphoton]; << 676 cosTheta = aStore.SampleMax( anEnergy );
668 G4ThreeVector direction(sinTheta * std::co << 677 }
669 G4ThreeVector photonP = photonE * directio << 678 else if ( tabulationType == 2 )
670 thePhotons->operator[](0)->SetMomentum(pho << 679 {
>> 680 // tabulation of probabilities.
>> 681
>> 682 G4int iangle = 0;
>> 683 for ( G4int j = 0 ; j < nNeu [ iphoton - nIso ] ; ++j )
>> 684 {
>> 685 iangle = j;
>> 686 if ( theAngular[ iphoton - nIso ][ j ].GetEnergy() > anEnergy ) break;
>> 687 }
>> 688 cosTheta = theAngular[iphoton-nIso][ iangle ].GetCosTh();
>> 689 // no interpolation yet @@
>> 690 }
>> 691 }
>> 692 }
>> 693
>> 694 // Set
>> 695 G4double phi = twopi*G4UniformRand();
>> 696 G4double theta = std::acos( cosTheta );
>> 697 G4double sinTheta = std::sin( theta );
>> 698
>> 699 G4double photonE = theGammas[ iphoton ];
>> 700 G4ThreeVector direction ( sinTheta*std::cos( phi ) , sinTheta * std::sin( phi ) , cosTheta );
>> 701 G4ThreeVector photonP = photonE * direction;
>> 702 thePhotons->operator[]( 0 )->SetMomentum( photonP ) ;
671 } 703 }
672 else { << 704 else
>> 705 {
673 delete thePhotons; 706 delete thePhotons;
674 thePhotons = nullptr; // no gamma data av << 707 thePhotons = nullptr; // no gamma data available; some work needed @@@@@@@
675 } << 708 }
676 return thePhotons; 709 return thePhotons;
677 } 710 }
678 711