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