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