| 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 | // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
|
|---|
| 35 | // 080714 Limiting the sum of energy of secondary particles by T. Koi
|
|---|
| 36 | // 080801 Fix div0 error wiht G4FPE and memory leak by T. Koi
|
|---|
| 37 | // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
|
|---|
| 38 | //
|
|---|
| 39 |
|
|---|
| 40 | #include "G4NeutronHPContAngularPar.hh"
|
|---|
| 41 | #include "G4NeutronHPLegendreStore.hh"
|
|---|
| 42 | #include "G4Gamma.hh"
|
|---|
| 43 | #include "G4Electron.hh"
|
|---|
| 44 | #include "G4Positron.hh"
|
|---|
| 45 | #include "G4Neutron.hh"
|
|---|
| 46 | #include "G4Proton.hh"
|
|---|
| 47 | #include "G4Deuteron.hh"
|
|---|
| 48 | #include "G4Triton.hh"
|
|---|
| 49 | #include "G4He3.hh"
|
|---|
| 50 | #include "G4Alpha.hh"
|
|---|
| 51 | #include "G4NeutronHPVector.hh"
|
|---|
| 52 | #include "G4NucleiProperties.hh"
|
|---|
| 53 | #include "G4NeutronHPKallbachMannSyst.hh"
|
|---|
| 54 | #include "G4ParticleTable.hh"
|
|---|
| 55 |
|
|---|
| 56 | void G4NeutronHPContAngularPar::Init(std::ifstream & aDataFile)
|
|---|
| 57 | {
|
|---|
| 58 | aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
|
|---|
| 59 | theEnergy *= eV;
|
|---|
| 60 | theAngular = new G4NeutronHPList [nEnergies];
|
|---|
| 61 | for(G4int i=0; i<nEnergies; i++)
|
|---|
| 62 | {
|
|---|
| 63 | G4double sEnergy;
|
|---|
| 64 | aDataFile >> sEnergy;
|
|---|
| 65 | sEnergy*=eV;
|
|---|
| 66 | theAngular[i].SetLabel(sEnergy);
|
|---|
| 67 | theAngular[i].Init(aDataFile, nAngularParameters, 1.);
|
|---|
| 68 | }
|
|---|
| 69 | }
|
|---|
| 70 |
|
|---|
| 71 | G4ReactionProduct *
|
|---|
| 72 | G4NeutronHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /*targetMass*/,
|
|---|
| 73 | G4int angularRep, G4int interpolE )
|
|---|
| 74 | {
|
|---|
| 75 | G4ReactionProduct * result = new G4ReactionProduct;
|
|---|
| 76 | G4int Z = static_cast<G4int>(massCode/1000);
|
|---|
| 77 | G4int A = static_cast<G4int>(massCode-1000*Z);
|
|---|
| 78 | if(massCode==0)
|
|---|
| 79 | {
|
|---|
| 80 | result->SetDefinition(G4Gamma::Gamma());
|
|---|
| 81 | }
|
|---|
| 82 | else if(A==0)
|
|---|
| 83 | {
|
|---|
| 84 | result->SetDefinition(G4Electron::Electron());
|
|---|
| 85 | if(Z==1) result->SetDefinition(G4Positron::Positron());
|
|---|
| 86 | }
|
|---|
| 87 | else if(A==1)
|
|---|
| 88 | {
|
|---|
| 89 | result->SetDefinition(G4Neutron::Neutron());
|
|---|
| 90 | if(Z==1) result->SetDefinition(G4Proton::Proton());
|
|---|
| 91 | }
|
|---|
| 92 | else if(A==2)
|
|---|
| 93 | {
|
|---|
| 94 | result->SetDefinition(G4Deuteron::Deuteron());
|
|---|
| 95 | }
|
|---|
| 96 | else if(A==3)
|
|---|
| 97 | {
|
|---|
| 98 | result->SetDefinition(G4Triton::Triton());
|
|---|
| 99 | if(Z==2) result->SetDefinition(G4He3::He3());
|
|---|
| 100 | }
|
|---|
| 101 | else if(A==4)
|
|---|
| 102 | {
|
|---|
| 103 | result->SetDefinition(G4Alpha::Alpha());
|
|---|
| 104 | if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unknown ion case 1");
|
|---|
| 105 | }
|
|---|
| 106 | else
|
|---|
| 107 | {
|
|---|
| 108 | result->SetDefinition(G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z));
|
|---|
| 109 | }
|
|---|
| 110 | G4int i(0);
|
|---|
| 111 | G4int it(0);
|
|---|
| 112 | G4double fsEnergy(0);
|
|---|
| 113 | G4double cosTh(0);
|
|---|
| 114 | if(angularRep==1)
|
|---|
| 115 | {
|
|---|
| 116 | // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
|
|---|
| 117 | if (interpolE == 2)
|
|---|
| 118 | {
|
|---|
| 119 |
|
|---|
| 120 | //TK080711
|
|---|
| 121 | if ( fresh == true )
|
|---|
| 122 | {
|
|---|
| 123 | remaining_energy = theAngular[0].GetLabel();
|
|---|
| 124 | fresh = false;
|
|---|
| 125 | }
|
|---|
| 126 | //TK080711
|
|---|
| 127 |
|
|---|
| 128 | G4double random = G4UniformRand();
|
|---|
| 129 | G4double * running = new G4double[nEnergies+1];
|
|---|
| 130 | running[0]=0;
|
|---|
| 131 |
|
|---|
| 132 | for(i=1; i<nEnergies+1; i++)
|
|---|
| 133 | {
|
|---|
| 134 | //TK080711
|
|---|
| 135 | if ( remaining_energy >= theAngular[ i-1 ].GetLabel() )
|
|---|
| 136 | running[i] = running[i-1] + theAngular[i-1].GetValue(0);
|
|---|
| 137 | else
|
|---|
| 138 | running[i] = running[i-1];
|
|---|
| 139 | //TK080711
|
|---|
| 140 | }
|
|---|
| 141 |
|
|---|
| 142 | //080730
|
|---|
| 143 | if ( running[ nEnergies ] != 0 )
|
|---|
| 144 | {
|
|---|
| 145 |
|
|---|
| 146 | for ( i = 1 ; i < nEnergies+1 ; i++ )
|
|---|
| 147 | {
|
|---|
| 148 | it = i-1;
|
|---|
| 149 | if ( random > running[ i-1 ]/running[ nEnergies ] && random <= running[ i ] / running[ nEnergies ] ) break;
|
|---|
| 150 | }
|
|---|
| 151 | fsEnergy = theAngular[ it ].GetLabel();
|
|---|
| 152 |
|
|---|
| 153 | }
|
|---|
| 154 |
|
|---|
| 155 | //TK080711
|
|---|
| 156 | if ( i == nEnergies+1 || running[ nEnergies ] == 0 ) fsEnergy = remaining_energy;
|
|---|
| 157 | //TK080711 //080730
|
|---|
| 158 |
|
|---|
| 159 | G4NeutronHPLegendreStore theStore(1);
|
|---|
| 160 | theStore.Init(0,fsEnergy,nAngularParameters);
|
|---|
| 161 | for(i=0;i<nAngularParameters;i++)
|
|---|
| 162 | {
|
|---|
| 163 | theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
|
|---|
| 164 | }
|
|---|
| 165 | // use it to sample.
|
|---|
| 166 | cosTh = theStore.SampleMax(fsEnergy);
|
|---|
| 167 |
|
|---|
| 168 | //TK080711
|
|---|
| 169 | remaining_energy -= fsEnergy;
|
|---|
| 170 | //TK080711
|
|---|
| 171 |
|
|---|
| 172 | //080801b
|
|---|
| 173 | delete[] running;
|
|---|
| 174 | //080801b
|
|---|
| 175 | }
|
|---|
| 176 | else
|
|---|
| 177 | {
|
|---|
| 178 |
|
|---|
| 179 | //080714
|
|---|
| 180 | if ( fresh == true )
|
|---|
| 181 | {
|
|---|
| 182 | remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
|
|---|
| 183 | fresh = false;
|
|---|
| 184 | }
|
|---|
| 185 | //080714
|
|---|
| 186 |
|
|---|
| 187 |
|
|---|
| 188 | G4double random = G4UniformRand();
|
|---|
| 189 | G4double * running = new G4double[nEnergies];
|
|---|
| 190 | running[0]=0;
|
|---|
| 191 | G4double weighted = 0;
|
|---|
| 192 | for(i=1; i<nEnergies; i++)
|
|---|
| 193 | {
|
|---|
| 194 | /*
|
|---|
| 195 | if(i!=0)
|
|---|
| 196 | {
|
|---|
| 197 | running[i]=running[i-1];
|
|---|
| 198 | }
|
|---|
| 199 | running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 200 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 201 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 202 | weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 203 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 204 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 205 | */
|
|---|
| 206 |
|
|---|
| 207 | running[i]=running[i-1];
|
|---|
| 208 | if ( remaining_energy >= theAngular[i].GetLabel() )
|
|---|
| 209 | {
|
|---|
| 210 | running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 211 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 212 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 213 | weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 214 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 215 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 216 | }
|
|---|
| 217 | }
|
|---|
| 218 | // cash the mean energy in this distribution
|
|---|
| 219 | //080409 TKDB
|
|---|
| 220 | if ( nEnergies == 1 || running[nEnergies-1] == 0 )
|
|---|
| 221 | currentMeanEnergy = 0.0;
|
|---|
| 222 | else
|
|---|
| 223 | {
|
|---|
| 224 | currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 225 | }
|
|---|
| 226 |
|
|---|
| 227 | //080409 TKDB
|
|---|
| 228 | if ( nEnergies == 1 ) it = 0;
|
|---|
| 229 |
|
|---|
| 230 | //080729
|
|---|
| 231 | if ( running[nEnergies-1] != 0 )
|
|---|
| 232 | {
|
|---|
| 233 | for ( i = 1 ; i < nEnergies ; i++ )
|
|---|
| 234 | {
|
|---|
| 235 | it = i;
|
|---|
| 236 | if ( random < running [ i ] / running [ nEnergies-1 ] ) break;
|
|---|
| 237 | }
|
|---|
| 238 | }
|
|---|
| 239 |
|
|---|
| 240 | //080714
|
|---|
| 241 | if ( running [ nEnergies-1 ] == 0 ) it = 0;
|
|---|
| 242 | //080714
|
|---|
| 243 |
|
|---|
| 244 | if(it<nDiscreteEnergies||it==0)
|
|---|
| 245 | {
|
|---|
| 246 | if(it == 0)
|
|---|
| 247 | {
|
|---|
| 248 | fsEnergy = theAngular[it].GetLabel();
|
|---|
| 249 | G4NeutronHPLegendreStore theStore(1);
|
|---|
| 250 | theStore.Init(0,fsEnergy,nAngularParameters);
|
|---|
| 251 | for(i=0;i<nAngularParameters;i++)
|
|---|
| 252 | {
|
|---|
| 253 | theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
|
|---|
| 254 | }
|
|---|
| 255 | // use it to sample.
|
|---|
| 256 | cosTh = theStore.SampleMax(fsEnergy);
|
|---|
| 257 | }
|
|---|
| 258 | else
|
|---|
| 259 | {
|
|---|
| 260 | G4double e1, e2;
|
|---|
| 261 | e1 = theAngular[it-1].GetLabel();
|
|---|
| 262 | e2 = theAngular[it].GetLabel();
|
|---|
| 263 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
|
|---|
| 264 | random,
|
|---|
| 265 | running[it-1]/running[nEnergies-1],
|
|---|
| 266 | running[it]/running[nEnergies-1],
|
|---|
| 267 | e1, e2);
|
|---|
| 268 | // fill a Legendrestore
|
|---|
| 269 | G4NeutronHPLegendreStore theStore(2);
|
|---|
| 270 | theStore.Init(0,e1,nAngularParameters);
|
|---|
| 271 | theStore.Init(1,e2,nAngularParameters);
|
|---|
| 272 | for(i=0;i<nAngularParameters;i++)
|
|---|
| 273 | {
|
|---|
| 274 | theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
|
|---|
| 275 | theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
|
|---|
| 276 | }
|
|---|
| 277 | // use it to sample.
|
|---|
| 278 | theStore.SetManager(theManager);
|
|---|
| 279 | cosTh = theStore.SampleMax(fsEnergy);
|
|---|
| 280 | }
|
|---|
| 281 | }
|
|---|
| 282 | else // continuum contribution
|
|---|
| 283 | {
|
|---|
| 284 | G4double x1 = running[it-1]/running[nEnergies-1];
|
|---|
| 285 | G4double x2 = running[it]/running[nEnergies-1];
|
|---|
| 286 | G4double y1 = theAngular[it-1].GetLabel();
|
|---|
| 287 | G4double y2 = theAngular[it].GetLabel();
|
|---|
| 288 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
|
|---|
| 289 | random,x1,x2,y1,y2);
|
|---|
| 290 | G4NeutronHPLegendreStore theStore(2);
|
|---|
| 291 | theStore.Init(0,y1,nAngularParameters);
|
|---|
| 292 | theStore.Init(1,y2,nAngularParameters);
|
|---|
| 293 | theStore.SetManager(theManager);
|
|---|
| 294 | for(i=0;i<nAngularParameters;i++)
|
|---|
| 295 | {
|
|---|
| 296 | theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
|
|---|
| 297 | theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
|
|---|
| 298 | }
|
|---|
| 299 | // use it to sample.
|
|---|
| 300 | cosTh = theStore.SampleMax(fsEnergy);
|
|---|
| 301 | }
|
|---|
| 302 | delete [] running;
|
|---|
| 303 |
|
|---|
| 304 | //080714
|
|---|
| 305 | remaining_energy -= fsEnergy;
|
|---|
| 306 | //080714
|
|---|
| 307 |
|
|---|
| 308 | }
|
|---|
| 309 |
|
|---|
| 310 | }
|
|---|
| 311 | else if(angularRep==2)
|
|---|
| 312 | {
|
|---|
| 313 | // first get the energy (already the right for this incoming energy)
|
|---|
| 314 | G4int i;
|
|---|
| 315 | G4double * running = new G4double[nEnergies];
|
|---|
| 316 | running[0]=0;
|
|---|
| 317 | G4double weighted = 0;
|
|---|
| 318 | for(i=1; i<nEnergies; i++)
|
|---|
| 319 | {
|
|---|
| 320 | if(i!=0) running[i]=running[i-1];
|
|---|
| 321 | running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 322 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 323 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 324 | weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 325 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 326 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 327 | }
|
|---|
| 328 | // cash the mean energy in this distribution
|
|---|
| 329 | //080409 TKDB
|
|---|
| 330 | //currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 331 | if ( nEnergies == 1 )
|
|---|
| 332 | currentMeanEnergy = 0.0;
|
|---|
| 333 | else
|
|---|
| 334 | currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 335 |
|
|---|
| 336 | G4int it(0);
|
|---|
| 337 | G4double randkal = G4UniformRand();
|
|---|
| 338 | //080409 TKDB
|
|---|
| 339 | //for(i=0; i<nEnergies; i++)
|
|---|
| 340 | for(i=1; i<nEnergies; i++)
|
|---|
| 341 | {
|
|---|
| 342 | it = i;
|
|---|
| 343 | if(randkal<running[i]/running[nEnergies-1]) break;
|
|---|
| 344 | }
|
|---|
| 345 |
|
|---|
| 346 | // interpolate the secondary energy.
|
|---|
| 347 | G4double x, x1,x2,y1,y2;
|
|---|
| 348 | if(it==0) it=1;
|
|---|
| 349 | x = randkal*running[nEnergies-1];
|
|---|
| 350 | x1 = running[it-1];
|
|---|
| 351 | x2 = running[it];
|
|---|
| 352 | G4double compoundFraction;
|
|---|
| 353 | // interpolate energy
|
|---|
| 354 | y1 = theAngular[it-1].GetLabel();
|
|---|
| 355 | y2 = theAngular[it].GetLabel();
|
|---|
| 356 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it-1),
|
|---|
| 357 | x, x1,x2,y1,y2);
|
|---|
| 358 | // for theta interpolate the compoundFractions
|
|---|
| 359 | G4double cLow = theAngular[it-1].GetValue(1);
|
|---|
| 360 | G4double cHigh = theAngular[it].GetValue(1);
|
|---|
| 361 | compoundFraction = theInt.Interpolate(theManager.GetScheme(it),
|
|---|
| 362 | fsEnergy, y1, y2, cLow,cHigh);
|
|---|
| 363 | delete [] running;
|
|---|
| 364 |
|
|---|
| 365 | // get cosTh
|
|---|
| 366 | G4double incidentEnergy = anEnergy;
|
|---|
| 367 | G4double incidentMass = G4Neutron::Neutron()->GetPDGMass();
|
|---|
| 368 | G4double productEnergy = fsEnergy;
|
|---|
| 369 | G4double productMass = result->GetMass();
|
|---|
| 370 | G4int targetZ = G4int(theTargetCode/1000);
|
|---|
| 371 | G4int targetA = G4int(theTargetCode-1000*targetZ);
|
|---|
| 372 | // To correspond to natural composition (-nat-) data files.
|
|---|
| 373 | if ( targetA == 0 )
|
|---|
| 374 | targetA = int ( theTarget->GetMass()/amu_c2 + 0.5 );
|
|---|
| 375 | G4double targetMass = theTarget->GetMass();
|
|---|
| 376 | G4int residualA = targetA+1-A;
|
|---|
| 377 | G4int residualZ = targetZ-Z;
|
|---|
| 378 | G4double residualMass = residualZ*G4Proton::Proton()->GetPDGMass();
|
|---|
| 379 | residualMass +=(residualA-residualZ)*G4Neutron::Neutron()->GetPDGMass();
|
|---|
| 380 | residualMass -= G4NucleiProperties::GetBindingEnergy( residualA , residualZ );
|
|---|
| 381 | G4NeutronHPKallbachMannSyst theKallbach(compoundFraction,
|
|---|
| 382 | incidentEnergy, incidentMass,
|
|---|
| 383 | productEnergy, productMass,
|
|---|
| 384 | residualMass, residualA, residualZ,
|
|---|
| 385 | targetMass, targetA, targetZ);
|
|---|
| 386 | cosTh = theKallbach.Sample(anEnergy);
|
|---|
| 387 | }
|
|---|
| 388 | else if(angularRep>10&&angularRep<16)
|
|---|
| 389 | {
|
|---|
| 390 | G4double random = G4UniformRand();
|
|---|
| 391 | G4double * running = new G4double[nEnergies];
|
|---|
| 392 | running[0]=0;
|
|---|
| 393 | G4double weighted = 0;
|
|---|
| 394 | for(i=1; i<nEnergies; i++)
|
|---|
| 395 | {
|
|---|
| 396 | if(i!=0) running[i]=running[i-1];
|
|---|
| 397 | running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 398 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 399 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 400 | weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
|
|---|
| 401 | theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
|
|---|
| 402 | theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
|
|---|
| 403 | }
|
|---|
| 404 | // cash the mean energy in this distribution
|
|---|
| 405 | //currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 406 | if ( nEnergies == 1 )
|
|---|
| 407 | currentMeanEnergy = 0.0;
|
|---|
| 408 | else
|
|---|
| 409 | currentMeanEnergy = weighted/running[nEnergies-1];
|
|---|
| 410 |
|
|---|
| 411 | //080409 TKDB
|
|---|
| 412 | if ( nEnergies == 1 ) it = 0;
|
|---|
| 413 | //for(i=0; i<nEnergies; i++)
|
|---|
| 414 | for(i=1; i<nEnergies; i++)
|
|---|
| 415 | {
|
|---|
| 416 | it = i;
|
|---|
| 417 | if(random<running[i]/running[nEnergies-1]) break;
|
|---|
| 418 | }
|
|---|
| 419 | if(it<nDiscreteEnergies||it==0)
|
|---|
| 420 | {
|
|---|
| 421 | if(it==0)
|
|---|
| 422 | {
|
|---|
| 423 | fsEnergy = theAngular[0].GetLabel();
|
|---|
| 424 | G4NeutronHPVector theStore;
|
|---|
| 425 | G4int aCounter = 0;
|
|---|
| 426 | for(G4int i=1; i<nAngularParameters; i+=2)
|
|---|
| 427 | {
|
|---|
| 428 | theStore.SetX(aCounter, theAngular[0].GetValue(i));
|
|---|
| 429 | theStore.SetY(aCounter, theAngular[0].GetValue(i+1));
|
|---|
| 430 | aCounter++;
|
|---|
| 431 | }
|
|---|
| 432 | G4InterpolationManager aMan;
|
|---|
| 433 | aMan.Init(angularRep-10, nAngularParameters-1);
|
|---|
| 434 | theStore.SetInterpolationManager(aMan);
|
|---|
| 435 | cosTh = theStore.Sample();
|
|---|
| 436 | }
|
|---|
| 437 | else
|
|---|
| 438 | {
|
|---|
| 439 | fsEnergy = theAngular[it].GetLabel();
|
|---|
| 440 | G4NeutronHPVector theStore;
|
|---|
| 441 | G4InterpolationManager aMan;
|
|---|
| 442 | aMan.Init(angularRep-10, nAngularParameters-1);
|
|---|
| 443 | theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 444 | G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
|
|---|
| 445 | G4int aCounter = 0;
|
|---|
| 446 | for(G4int i=1; i<nAngularParameters; i+=2)
|
|---|
| 447 | {
|
|---|
| 448 | theStore.SetX(aCounter, theAngular[it].GetValue(i));
|
|---|
| 449 | theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
|
|---|
| 450 | random,
|
|---|
| 451 | running[it-1]/running[nEnergies-1],
|
|---|
| 452 | running[it]/running[nEnergies-1],
|
|---|
| 453 | theAngular[it-1].GetValue(i+1),
|
|---|
| 454 | theAngular[it].GetValue(i+1)));
|
|---|
| 455 | aCounter++;
|
|---|
| 456 | }
|
|---|
| 457 | cosTh = theStore.Sample();
|
|---|
| 458 | }
|
|---|
| 459 | }
|
|---|
| 460 | else
|
|---|
| 461 | {
|
|---|
| 462 | G4double x1 = running[it-1]/running[nEnergies-1];
|
|---|
| 463 | G4double x2 = running[it]/running[nEnergies-1];
|
|---|
| 464 | G4double y1 = theAngular[it-1].GetLabel();
|
|---|
| 465 | G4double y2 = theAngular[it].GetLabel();
|
|---|
| 466 | fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
|
|---|
| 467 | random,x1,x2,y1,y2);
|
|---|
| 468 | G4NeutronHPVector theBuff1;
|
|---|
| 469 | G4NeutronHPVector theBuff2;
|
|---|
| 470 | G4InterpolationManager aMan;
|
|---|
| 471 | aMan.Init(angularRep-10, nAngularParameters-1);
|
|---|
| 472 | // theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 473 | // theBuff2.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 474 | for(i=0; i<nAngularParameters; i++) // i=1 ist wichtig!
|
|---|
| 475 | {
|
|---|
| 476 | theBuff1.SetX(i, theAngular[it-1].GetValue(i));
|
|---|
| 477 | theBuff1.SetY(i, theAngular[it-1].GetValue(i+1));
|
|---|
| 478 | theBuff2.SetX(i, theAngular[it].GetValue(i));
|
|---|
| 479 | theBuff2.SetY(i, theAngular[it].GetValue(i+1));
|
|---|
| 480 | i++;
|
|---|
| 481 | }
|
|---|
| 482 | G4NeutronHPVector theStore;
|
|---|
| 483 | theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
|
|---|
| 484 | x1 = y1;
|
|---|
| 485 | x2 = y2;
|
|---|
| 486 | G4double x, y;
|
|---|
| 487 | //for(i=0;i<theBuff1.GetVectorLength(); i++);
|
|---|
| 488 | for(i=0;i<theBuff1.GetVectorLength(); i++)
|
|---|
| 489 | {
|
|---|
| 490 | x = theBuff1.GetX(i); // costh binning identical
|
|---|
| 491 | y1 = theBuff1.GetY(i);
|
|---|
| 492 | y2 = theBuff2.GetY(i);
|
|---|
| 493 | y = theInt.Interpolate(theManager.GetScheme(it),
|
|---|
| 494 | fsEnergy, theAngular[it-1].GetLabel(),
|
|---|
| 495 | theAngular[it].GetLabel(), y1, y2);
|
|---|
| 496 | theStore.SetX(i, x);
|
|---|
| 497 | theStore.SetY(i, y);
|
|---|
| 498 | }
|
|---|
| 499 | cosTh = theStore.Sample();
|
|---|
| 500 | }
|
|---|
| 501 | delete [] running;
|
|---|
| 502 | }
|
|---|
| 503 | else
|
|---|
| 504 | {
|
|---|
| 505 | throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar::Sample: Unknown angular representation");
|
|---|
| 506 | }
|
|---|
| 507 | result->SetKineticEnergy(fsEnergy);
|
|---|
| 508 | G4double phi = twopi*G4UniformRand();
|
|---|
| 509 | G4double theta = std::acos(cosTh);
|
|---|
| 510 | G4double sinth = std::sin(theta);
|
|---|
| 511 | G4double mtot = result->GetTotalMomentum();
|
|---|
| 512 | G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
|
|---|
| 513 | result->SetMomentum(tempVector);
|
|---|
| 514 | // return the result.
|
|---|
| 515 | return result;
|
|---|
| 516 | }
|
|---|
