| 1 | //
|
|---|
| 2 | // ********************************************************************
|
|---|
| 3 | // * License and Disclaimer *
|
|---|
| 4 | // * *
|
|---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of *
|
|---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and *
|
|---|
| 7 | // * conditions of the Geant4 Software License, included in the file *
|
|---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These *
|
|---|
| 9 | // * include a list of copyright holders. *
|
|---|
| 10 | // * *
|
|---|
| 11 | // * Neither the authors of this software system, nor their employing *
|
|---|
| 12 | // * institutes,nor the agencies providing financial support for this *
|
|---|
| 13 | // * work make any representation or warranty, express or implied, *
|
|---|
| 14 | // * regarding this software system or assume any liability for its *
|
|---|
| 15 | // * use. Please see the license in the file LICENSE and URL above *
|
|---|
| 16 | // * for the full disclaimer and the limitation of liability. *
|
|---|
| 17 | // * *
|
|---|
| 18 | // * This code implementation is the result of the scientific and *
|
|---|
| 19 | // * technical work of the GEANT4 collaboration. *
|
|---|
| 20 | // * By using, copying, modifying or distributing the software (or *
|
|---|
| 21 | // * any work based on the software) you agree to acknowledge its *
|
|---|
| 22 | // * use in resulting scientific publications, and indicate your *
|
|---|
| 23 | // * acceptance of all terms of the Geant4 Software license. *
|
|---|
| 24 | // ********************************************************************
|
|---|
| 25 | //
|
|---|
| 26 | // neutron_hp -- source file
|
|---|
| 27 | // J.P. Wellisch, Nov-1996
|
|---|
| 28 | // A prototype of the low energy neutron transport model.
|
|---|
| 29 | //
|
|---|
| 30 | // 09-May-06 fix in Sample by T. Koi
|
|---|
| 31 | // 080318 Fix Compilation warnings - gcc-4.3.0 by T. Koi
|
|---|
| 32 | // (This fix has a real effect to the code.)
|
|---|
| 33 | // 080409 Fix div0 error with G4FPE by T. Koi
|
|---|
| 34 | //
|
|---|
| 35 |
|
|---|
| 36 | #include "G4NeutronHPContAngularPar.hh"
|
|---|
| 37 | #include "G4NeutronHPLegendreStore.hh"
|
|---|
| 38 | #include "G4Gamma.hh"
|
|---|
| 39 | #include "G4Electron.hh"
|
|---|
| 40 | #include "G4Positron.hh"
|
|---|
| 41 | #include "G4Neutron.hh"
|
|---|
| 42 | #include "G4Proton.hh"
|
|---|
| 43 | #include "G4Deuteron.hh"
|
|---|
| 44 | #include "G4Triton.hh"
|
|---|
| 45 | #include "G4He3.hh"
|
|---|
| 46 | #include "G4Alpha.hh"
|
|---|
| 47 | #include "G4NeutronHPVector.hh"
|
|---|
| 48 | #include "G4NucleiPropertiesTable.hh"
|
|---|
| 49 | #include "G4NeutronHPKallbachMannSyst.hh"
|
|---|
| 50 | #include "G4ParticleTable.hh"
|
|---|
| 51 |
|
|---|
| 52 | void G4NeutronHPContAngularPar::Init(std::ifstream & aDataFile)
|
|---|
| 53 | {
|
|---|
| 54 | aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
|
|---|
| 55 | theEnergy *= eV;
|
|---|
| 56 | theAngular = new G4NeutronHPList [nEnergies];
|
|---|
| 57 | for(G4int i=0; i<nEnergies; i++)
|
|---|
| 58 | {
|
|---|
| 59 | G4double sEnergy;
|
|---|
| 60 | aDataFile >> sEnergy;
|
|---|
| 61 | sEnergy*=eV;
|
|---|
| 62 | theAngular[i].SetLabel(sEnergy);
|
|---|
| 63 | theAngular[i].Init(aDataFile, nAngularParameters, 1.);
|
|---|
| 64 | }
|
|---|
| 65 | }
|
|---|
| 66 |
|
|---|
| 67 | G4ReactionProduct *
|
|---|
| 68 | G4NeutronHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /*targetMass*/,
|
|---|
| 69 | G4int angularRep, G4int /*interpolE*/)
|
|---|
| 70 | {
|
|---|
| 71 | G4ReactionProduct * result = new G4ReactionProduct;
|
|---|
| 72 | G4int Z = static_cast<G4int>(massCode/1000);
|
|---|
| 73 | G4int A = static_cast<G4int>(massCode-1000*Z);
|
|---|
| 74 | if(massCode==0)
|
|---|
| 75 | {
|
|---|
| 76 | result->SetDefinition(G4Gamma::Gamma());
|
|---|
| 77 | }
|
|---|
| 78 | else if(A==0)
|
|---|
| 79 | {
|
|---|
| 80 | result->SetDefinition(G4Electron::Electron());
|
|---|
| 81 | if(Z==1) result->SetDefinition(G4Positron::Positron());
|
|---|
| 82 | }
|
|---|
| 83 | else if(A==1)
|
|---|
| 84 | {
|
|---|
| 85 | result->SetDefinition(G4Neutron::Neutron());
|
|---|
| 86 | if(Z==1) result->SetDefinition(G4Proton::Proton());
|
|---|
| 87 | }
|
|---|
| 88 | else if(A==2)
|
|---|
| 89 | {
|
|---|
| 90 | result->SetDefinition(G4Deuteron::Deuteron());
|
|---|
| 91 | }
|
|---|
| 92 | else if(A==3)
|
|---|
| 93 | {
|
|---|
| 94 | result->SetDefinition(G4Triton::Triton());
|
|---|
| 95 | if(Z==2) result->SetDefinition(G4He3::He3());
|
|---|
| 96 | }
|
|---|
| 97 | else if(A==4)
|
|---|
| 98 | {
|
|---|
| 99 | result->SetDefinition(G4Alpha::Alpha());
|
|---|
| 100 | if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unknown ion case 1");
|
|---|
| 101 | }
|
|---|
| 102 | else
|
|---|
| 103 | {
|
|---|
| 104 | result->SetDefinition(G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z));
|
|---|
| 105 | }
|
|---|
| 106 | G4int i(0);
|
|---|
| 107 | G4int it(0);
|
|---|
| 108 | G4double fsEnergy(0);
|
|---|
| 109 | G4double cosTh(0);
|
|---|
| 110 | if(angularRep==1)
|
|---|
| 111 | {
|
|---|
| 112 | G4double random = G4UniformRand();
|
|---|
| 113 | G4double * running = new G4double[nEnergies];
|
|---|
| 114 | running[0]=0;
|
|---|
| 115 | G4double weighted = 0;
|
|---|
| 116 | for(i=1; i<nEnergies; i++)
|
|---|
| 117 | {
|
|---|
| 118 | if(i!=0)
|
|---|
| 119 | {
|
|---|
| 120 | running[i]=running[i-1];
|
|---|
| 121 | }
|
|---|
| 122 | running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 123 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 124 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 125 | weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 126 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 127 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 128 | }
|
|---|
| 129 | // cash the mean energy in this distribution
|
|---|
| 130 | //080409 TKDB
|
|---|
| 131 | if ( nEnergies == 1 )
|
|---|
| 132 | currentMeanEnergy = 0.0;
|
|---|
| 133 | else
|
|---|
| 134 | currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 135 |
|
|---|
| 136 | //080409 TKDB
|
|---|
| 137 | if ( nEnergies == 1 ) it = 0;
|
|---|
| 138 | //for(i=1; i<nEnergies; i++)
|
|---|
| 139 | for(i=1; i<nEnergies; i++)
|
|---|
| 140 | {
|
|---|
| 141 | it = i;
|
|---|
| 142 | if(random<running[i]/running[nEnergies-1]) break;
|
|---|
| 143 | }
|
|---|
| 144 | if(it<nDiscreteEnergies||it==0)
|
|---|
| 145 | {
|
|---|
| 146 | if(it == 0)
|
|---|
| 147 | {
|
|---|
| 148 | fsEnergy = theAngular[it].GetLabel();
|
|---|
| 149 | G4NeutronHPLegendreStore theStore(1);
|
|---|
| 150 | theStore.Init(0,fsEnergy,nAngularParameters);
|
|---|
| 151 | for(i=0;i<nAngularParameters;i++)
|
|---|
| 152 | {
|
|---|
| 153 | theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
|
|---|
| 154 | }
|
|---|
| 155 | // use it to sample.
|
|---|
| 156 | cosTh = theStore.SampleMax(fsEnergy);
|
|---|
| 157 | }
|
|---|
| 158 | else
|
|---|
| 159 | {
|
|---|
| 160 | G4double e1, e2;
|
|---|
| 161 | e1 = theAngular[it-1].GetLabel();
|
|---|
| 162 | e2 = theAngular[it].GetLabel();
|
|---|
| 163 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
|
|---|
| 164 | random,
|
|---|
| 165 | running[it-1]/running[nEnergies-1],
|
|---|
| 166 | running[it]/running[nEnergies-1],
|
|---|
| 167 | e1, e2);
|
|---|
| 168 | // fill a Legendrestore
|
|---|
| 169 | G4NeutronHPLegendreStore theStore(2);
|
|---|
| 170 | theStore.Init(0,e1,nAngularParameters);
|
|---|
| 171 | theStore.Init(1,e2,nAngularParameters);
|
|---|
| 172 | for(i=0;i<nAngularParameters;i++)
|
|---|
| 173 | {
|
|---|
| 174 | theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
|
|---|
| 175 | theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
|
|---|
| 176 | }
|
|---|
| 177 | // use it to sample.
|
|---|
| 178 | theStore.SetManager(theManager);
|
|---|
| 179 | cosTh = theStore.SampleMax(fsEnergy);
|
|---|
| 180 | }
|
|---|
| 181 | }
|
|---|
| 182 | else // continuum contribution
|
|---|
| 183 | {
|
|---|
| 184 | G4double x1 = running[it-1]/running[nEnergies-1];
|
|---|
| 185 | G4double x2 = running[it]/running[nEnergies-1];
|
|---|
| 186 | G4double y1 = theAngular[it-1].GetLabel();
|
|---|
| 187 | G4double y2 = theAngular[it].GetLabel();
|
|---|
| 188 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
|
|---|
| 189 | random,x1,x2,y1,y2);
|
|---|
| 190 | G4NeutronHPLegendreStore theStore(2);
|
|---|
| 191 | theStore.Init(0,y1,nAngularParameters);
|
|---|
| 192 | theStore.Init(1,y2,nAngularParameters);
|
|---|
| 193 | theStore.SetManager(theManager);
|
|---|
| 194 | for(i=0;i<nAngularParameters;i++)
|
|---|
| 195 | {
|
|---|
| 196 | theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
|
|---|
| 197 | theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
|
|---|
| 198 | }
|
|---|
| 199 | // use it to sample.
|
|---|
| 200 | cosTh = theStore.SampleMax(fsEnergy);
|
|---|
| 201 | }
|
|---|
| 202 | delete [] running;
|
|---|
| 203 | }
|
|---|
| 204 | else if(angularRep==2)
|
|---|
| 205 | {
|
|---|
| 206 | // first get the energy (already the right for this incoming energy)
|
|---|
| 207 | G4int i;
|
|---|
| 208 | G4double * running = new G4double[nEnergies];
|
|---|
| 209 | running[0]=0;
|
|---|
| 210 | G4double weighted = 0;
|
|---|
| 211 | for(i=1; i<nEnergies; i++)
|
|---|
| 212 | {
|
|---|
| 213 | if(i!=0) running[i]=running[i-1];
|
|---|
| 214 | running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 215 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 216 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 217 | weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 218 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 219 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 220 | }
|
|---|
| 221 | // cash the mean energy in this distribution
|
|---|
| 222 | //080409 TKDB
|
|---|
| 223 | //currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 224 | if ( nEnergies == 1 )
|
|---|
| 225 | currentMeanEnergy = 0.0;
|
|---|
| 226 | else
|
|---|
| 227 | currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 228 |
|
|---|
| 229 | G4int it(0);
|
|---|
| 230 | G4double randkal = G4UniformRand();
|
|---|
| 231 | //080409 TKDB
|
|---|
| 232 | //for(i=0; i<nEnergies; i++)
|
|---|
| 233 | for(i=1; i<nEnergies; i++)
|
|---|
| 234 | {
|
|---|
| 235 | it = i;
|
|---|
| 236 | if(randkal<running[i]/running[nEnergies-1]) break;
|
|---|
| 237 | }
|
|---|
| 238 |
|
|---|
| 239 | // interpolate the secondary energy.
|
|---|
| 240 | G4double x, x1,x2,y1,y2;
|
|---|
| 241 | if(it==0) it=1;
|
|---|
| 242 | x = randkal*running[nEnergies-1];
|
|---|
| 243 | x1 = running[it-1];
|
|---|
| 244 | x2 = running[it];
|
|---|
| 245 | G4double compoundFraction;
|
|---|
| 246 | // interpolate energy
|
|---|
| 247 | y1 = theAngular[it-1].GetLabel();
|
|---|
| 248 | y2 = theAngular[it].GetLabel();
|
|---|
| 249 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it-1),
|
|---|
| 250 | x, x1,x2,y1,y2);
|
|---|
| 251 | // for theta interpolate the compoundFractions
|
|---|
| 252 | G4double cLow = theAngular[it-1].GetValue(1);
|
|---|
| 253 | G4double cHigh = theAngular[it].GetValue(1);
|
|---|
| 254 | compoundFraction = theInt.Interpolate(theManager.GetScheme(it),
|
|---|
| 255 | fsEnergy, y1, y2, cLow,cHigh);
|
|---|
| 256 | delete [] running;
|
|---|
| 257 |
|
|---|
| 258 | // get cosTh
|
|---|
| 259 | G4double incidentEnergy = anEnergy;
|
|---|
| 260 | G4double incidentMass = G4Neutron::Neutron()->GetPDGMass();
|
|---|
| 261 | G4double productEnergy = fsEnergy;
|
|---|
| 262 | G4double productMass = result->GetMass();
|
|---|
| 263 | G4int targetZ = G4int(theTargetCode/1000);
|
|---|
| 264 | G4int targetA = G4int(theTargetCode-1000*targetZ);
|
|---|
| 265 | // To correspond to natural composition (-nat-) data files.
|
|---|
| 266 | if ( targetA == 0 )
|
|---|
| 267 | targetA = int ( theTarget->GetMass()/amu_c2 + 0.5 );
|
|---|
| 268 | G4double targetMass = theTarget->GetMass();
|
|---|
| 269 | G4int residualA = targetA+1-A;
|
|---|
| 270 | G4int residualZ = targetZ-Z;
|
|---|
| 271 | G4double residualMass = residualZ*G4Proton::Proton()->GetPDGMass();
|
|---|
| 272 | residualMass +=(residualA-residualZ)*G4Neutron::Neutron()->GetPDGMass();
|
|---|
| 273 | residualMass -= G4NucleiPropertiesTable::GetBindingEnergy(residualZ, residualA);
|
|---|
| 274 | G4NeutronHPKallbachMannSyst theKallbach(compoundFraction,
|
|---|
| 275 | incidentEnergy, incidentMass,
|
|---|
| 276 | productEnergy, productMass,
|
|---|
| 277 | residualMass, residualA, residualZ,
|
|---|
| 278 | targetMass, targetA, targetZ);
|
|---|
| 279 | cosTh = theKallbach.Sample(anEnergy);
|
|---|
| 280 | }
|
|---|
| 281 | else if(angularRep>10&&angularRep<16)
|
|---|
| 282 | {
|
|---|
| 283 | G4double random = G4UniformRand();
|
|---|
| 284 | G4double * running = new G4double[nEnergies];
|
|---|
| 285 | running[0]=0;
|
|---|
| 286 | G4double weighted = 0;
|
|---|
| 287 | for(i=1; i<nEnergies; i++)
|
|---|
| 288 | {
|
|---|
| 289 | if(i!=0) running[i]=running[i-1];
|
|---|
| 290 | running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 291 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 292 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 293 | weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 294 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 295 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 296 | }
|
|---|
| 297 | // cash the mean energy in this distribution
|
|---|
| 298 | //currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 299 | if ( nEnergies == 1 )
|
|---|
| 300 | currentMeanEnergy = 0.0;
|
|---|
| 301 | else
|
|---|
| 302 | currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 303 |
|
|---|
| 304 | //080409 TKDB
|
|---|
| 305 | if ( nEnergies == 1 ) it = 0;
|
|---|
| 306 | //for(i=0; i<nEnergies; i++)
|
|---|
| 307 | for(i=1; i<nEnergies; i++)
|
|---|
| 308 | {
|
|---|
| 309 | it = i;
|
|---|
| 310 | if(random<running[i]/running[nEnergies-1]) break;
|
|---|
| 311 | }
|
|---|
| 312 | if(it<nDiscreteEnergies||it==0)
|
|---|
| 313 | {
|
|---|
| 314 | if(it==0)
|
|---|
| 315 | {
|
|---|
| 316 | fsEnergy = theAngular[0].GetLabel();
|
|---|
| 317 | G4NeutronHPVector theStore;
|
|---|
| 318 | G4int aCounter = 0;
|
|---|
| 319 | for(G4int i=1; i<nAngularParameters; i+=2)
|
|---|
| 320 | {
|
|---|
| 321 | theStore.SetX(aCounter, theAngular[0].GetValue(i));
|
|---|
| 322 | theStore.SetY(aCounter, theAngular[0].GetValue(i+1));
|
|---|
| 323 | aCounter++;
|
|---|
| 324 | }
|
|---|
| 325 | G4InterpolationManager aMan;
|
|---|
| 326 | aMan.Init(angularRep-10, nAngularParameters-1);
|
|---|
| 327 | theStore.SetInterpolationManager(aMan);
|
|---|
| 328 | cosTh = theStore.Sample();
|
|---|
| 329 | }
|
|---|
| 330 | else
|
|---|
| 331 | {
|
|---|
| 332 | fsEnergy = theAngular[it].GetLabel();
|
|---|
| 333 | G4NeutronHPVector theStore;
|
|---|
| 334 | G4InterpolationManager aMan;
|
|---|
| 335 | aMan.Init(angularRep-10, nAngularParameters-1);
|
|---|
| 336 | theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 337 | G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
|
|---|
| 338 | G4int aCounter = 0;
|
|---|
| 339 | for(G4int i=1; i<nAngularParameters; i+=2)
|
|---|
| 340 | {
|
|---|
| 341 | theStore.SetX(aCounter, theAngular[it].GetValue(i));
|
|---|
| 342 | theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
|
|---|
| 343 | random,
|
|---|
| 344 | running[it-1]/running[nEnergies-1],
|
|---|
| 345 | running[it]/running[nEnergies-1],
|
|---|
| 346 | theAngular[it-1].GetValue(i+1),
|
|---|
| 347 | theAngular[it].GetValue(i+1)));
|
|---|
| 348 | aCounter++;
|
|---|
| 349 | }
|
|---|
| 350 | cosTh = theStore.Sample();
|
|---|
| 351 | }
|
|---|
| 352 | }
|
|---|
| 353 | else
|
|---|
| 354 | {
|
|---|
| 355 | G4double x1 = running[it-1]/running[nEnergies-1];
|
|---|
| 356 | G4double x2 = running[it]/running[nEnergies-1];
|
|---|
| 357 | G4double y1 = theAngular[it-1].GetLabel();
|
|---|
| 358 | G4double y2 = theAngular[it].GetLabel();
|
|---|
| 359 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
|
|---|
| 360 | random,x1,x2,y1,y2);
|
|---|
| 361 | G4NeutronHPVector theBuff1;
|
|---|
| 362 | G4NeutronHPVector theBuff2;
|
|---|
| 363 | G4InterpolationManager aMan;
|
|---|
| 364 | aMan.Init(angularRep-10, nAngularParameters-1);
|
|---|
| 365 | // theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 366 | // theBuff2.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 367 | for(i=0; i<nAngularParameters; i++) // i=1 ist wichtig!
|
|---|
| 368 | {
|
|---|
| 369 | theBuff1.SetX(i, theAngular[it-1].GetValue(i));
|
|---|
| 370 | theBuff1.SetY(i, theAngular[it-1].GetValue(i+1));
|
|---|
| 371 | theBuff2.SetX(i, theAngular[it].GetValue(i));
|
|---|
| 372 | theBuff2.SetY(i, theAngular[it].GetValue(i+1));
|
|---|
| 373 | i++;
|
|---|
| 374 | }
|
|---|
| 375 | G4NeutronHPVector theStore;
|
|---|
| 376 | theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 377 | x1 = y1;
|
|---|
| 378 | x2 = y2;
|
|---|
| 379 | G4double x, y;
|
|---|
| 380 | //for(i=0;i<theBuff1.GetVectorLength(); i++);
|
|---|
| 381 | for(i=0;i<theBuff1.GetVectorLength(); i++)
|
|---|
| 382 | {
|
|---|
| 383 | x = theBuff1.GetX(i); // costh binning identical
|
|---|
| 384 | y1 = theBuff1.GetY(i);
|
|---|
| 385 | y2 = theBuff2.GetY(i);
|
|---|
| 386 | y = theInt.Interpolate(theManager.GetScheme(it),
|
|---|
| 387 | fsEnergy, theAngular[it-1].GetLabel(),
|
|---|
| 388 | theAngular[it].GetLabel(), y1, y2);
|
|---|
| 389 | theStore.SetX(i, x);
|
|---|
| 390 | theStore.SetY(i, y);
|
|---|
| 391 | }
|
|---|
| 392 | cosTh = theStore.Sample();
|
|---|
| 393 | }
|
|---|
| 394 | delete [] running;
|
|---|
| 395 | }
|
|---|
| 396 | else
|
|---|
| 397 | {
|
|---|
| 398 | throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar::Sample: Unknown angular representation");
|
|---|
| 399 | }
|
|---|
| 400 | result->SetKineticEnergy(fsEnergy);
|
|---|
| 401 | G4double phi = twopi*G4UniformRand();
|
|---|
| 402 | G4double theta = std::acos(cosTh);
|
|---|
| 403 | G4double sinth = std::sin(theta);
|
|---|
| 404 | G4double mtot = result->GetTotalMomentum();
|
|---|
| 405 | G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
|
|---|
| 406 | result->SetMomentum(tempVector);
|
|---|
| 407 | // return the result.
|
|---|
| 408 | return result;
|
|---|
| 409 | }
|
|---|
