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