Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4NeutronHPContAngularPar Class Reference

#include <G4NeutronHPContAngularPar.hh>

Public Member Functions

 G4NeutronHPContAngularPar ()
 
 ~G4NeutronHPContAngularPar ()
 
void Init (std::ifstream &aDataFile)
 
G4ReactionProduct * Sample (G4double anEnergy, G4double massCode, G4double mass, G4int angularRep, G4int interpol)
 
G4double GetEnergy ()
 
void SetPrimary (G4ReactionProduct *aPrimary)
 
void SetTarget (G4ReactionProduct *aTarget)
 
void SetTargetCode (G4double aTargetCode)
 
void SetInterpolation (G4int theInterpolation)
 
void Merge (G4double anEnergy, G4InterpolationScheme &aScheme, G4NeutronHPContAngularPar &store1, G4NeutronHPContAngularPar &store2)
 
G4double MeanEnergyOfThisInteraction ()
 
void ClearHistories ()
 

Detailed Description

Definition at line 42 of file G4NeutronHPContAngularPar.hh.

Constructor & Destructor Documentation

◆ G4NeutronHPContAngularPar()

G4NeutronHPContAngularPar::G4NeutronHPContAngularPar ( )
inline

Definition at line 46 of file G4NeutronHPContAngularPar.hh.

47 {
48 theAngular = 0;
49 currentMeanEnergy = -2;
50 fresh = true;
51 }

◆ ~G4NeutronHPContAngularPar()

G4NeutronHPContAngularPar::~G4NeutronHPContAngularPar ( )
inline

Definition at line 52 of file G4NeutronHPContAngularPar.hh.

53 {
54 if(theAngular!=0) delete [] theAngular;
55 }

Member Function Documentation

◆ ClearHistories()

void G4NeutronHPContAngularPar::ClearHistories ( )
inline

Definition at line 155 of file G4NeutronHPContAngularPar.hh.

155{ fresh = true; };

◆ GetEnergy()

G4double G4NeutronHPContAngularPar::GetEnergy ( )
inline

Definition at line 62 of file G4NeutronHPContAngularPar.hh.

62{ return theEnergy; }

◆ Init()

void G4NeutronHPContAngularPar::Init ( std::ifstream &  aDataFile)

Definition at line 58 of file G4NeutronHPContAngularPar.cc.

59 {
60 aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
61 theEnergy *= eV;
62 theAngular = new G4NeutronHPList [nEnergies];
63 for(G4int i=0; i<nEnergies; i++)
64 {
65 G4double sEnergy;
66 aDataFile >> sEnergy;
67 sEnergy*=eV;
68 theAngular[i].SetLabel(sEnergy);
69 theAngular[i].Init(aDataFile, nAngularParameters, 1.);
70 }
71 }
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4NeutronHPList::SetLabel
void SetLabel(G4double aLabel)
Definition: G4NeutronHPList.hh:65
G4NeutronHPList::Init
void Init(std::ifstream &aDataFile, G4int nPar, G4double unit=1.)
Definition: G4NeutronHPList.cc:51

Referenced by G4NeutronHPContEnergyAngular::Init().

◆ MeanEnergyOfThisInteraction()

G4double G4NeutronHPContAngularPar::MeanEnergyOfThisInteraction ( )
inline

Definition at line 111 of file G4NeutronHPContAngularPar.hh.

112 {
113 G4double result;
114 if(currentMeanEnergy<-1)
115 {
116 return 0;
117 // throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Logical error in Product class");
118 }
119 else
120 {
121 result = currentMeanEnergy;
122 }
123 currentMeanEnergy = -2;
124 return result;
125 }

Referenced by G4NeutronHPContEnergyAngular::Sample().

◆ Merge()

void G4NeutronHPContAngularPar::Merge ( G4double  anEnergy,
G4InterpolationScheme &  aScheme,
G4NeutronHPContAngularPar &  store1,
G4NeutronHPContAngularPar &  store2 
)
inline

Definition at line 81 of file G4NeutronHPContAngularPar.hh.

84 {
85 nDiscreteEnergies = store1.nDiscreteEnergies;
86 nAngularParameters = store1.nAngularParameters;
87 nEnergies = store1.nEnergies;
88 theManager = store1.theManager;
89 theEnergy = anEnergy;
90 if(theAngular != 0) delete [] theAngular;
91 theAngular = new G4NeutronHPList[nEnergies];
92 G4int i, ii;
93 G4double value;
94 for(i=0; i<nEnergies; i++)
95 {
96 theAngular[i].SetLabel(store1.theAngular[i].GetLabel());
97 for(ii=0; ii<nAngularParameters; ii++)
98 {
99// G4cout <<"test "<<i<<" "<<store1.theEnergy<<" "<<store2.theEnergy<<" "
100// << store1.theAngular[i].GetValue(ii)<<" "<<
101// store2.theAngular[i].GetValue(ii)<<G4endl;
102 value = theInt.Interpolate(aScheme, anEnergy,
103 store1.theEnergy, store2.theEnergy,
104 store1.theAngular[i].GetValue(ii),
105 store2.theAngular[i].GetValue(ii));
106 theAngular[i].SetValue(ii, value);
107 }
108 }
109 };
G4NeutronHPInterpolator::Interpolate
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
Definition: G4NeutronHPInterpolator.hh:86
G4NeutronHPList::SetValue
void SetValue(G4int i, G4double y)
Definition: G4NeutronHPList.hh:50
G4NeutronHPList::GetLabel
G4double GetLabel()
Definition: G4NeutronHPList.hh:67
G4NeutronHPList::GetValue
G4double GetValue(G4int i)
Definition: G4NeutronHPList.cc:74

◆ Sample()

G4ReactionProduct * G4NeutronHPContAngularPar::Sample ( G4double  anEnergy,
G4double  massCode,
G4double  mass,
G4int  angularRep,
G4int  interpol 
)

Definition at line 74 of file G4NeutronHPContAngularPar.cc.

76 {
77 G4ReactionProduct * result = new G4ReactionProduct;
78 G4int Z = static_cast<G4int>(massCode/1000);
79 G4int A = static_cast<G4int>(massCode-1000*Z);
80 if(massCode==0)
81 {
82 result->SetDefinition(G4Gamma::Gamma());
83 }
84 else if(A==0)
85 {
86 result->SetDefinition(G4Electron::Electron());
87 if(Z==1) result->SetDefinition(G4Positron::Positron());
88 }
89 else if(A==1)
90 {
91 result->SetDefinition(G4Neutron::Neutron());
92 if(Z==1) result->SetDefinition(G4Proton::Proton());
93 }
94 else if(A==2)
95 {
96 result->SetDefinition(G4Deuteron::Deuteron());
97 }
98 else if(A==3)
99 {
100 result->SetDefinition(G4Triton::Triton());
101 if(Z==2) result->SetDefinition(G4He3::He3());
102 }
103 else if(A==4)
104 {
105 result->SetDefinition(G4Alpha::Alpha());
106 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unknown ion case 1");
107 }
108 else
109 {
110 result->SetDefinition(G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z));
111 }
112 G4int i(0);
113 G4int it(0);
114 G4double fsEnergy(0);
115 G4double cosTh(0);
116
117 if( angularRep == 1 )
118 {
119// 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
120 //if (interpolE == 2)
121//110609 above was wrong interupition, pointed out by E.Mendoza and D.Cano (CIMAT)
122//Following are reviesd version written by T.Koi (SLAC)
123 if ( nDiscreteEnergies != 0 )
124 {
125
126//1st check remaining_energy
127// if this is the first set it. (How?)
128 if ( fresh == true )
129 {
130 //Discrete Lines, larger energies come first
131 //Continues Emssions, low to high LAST
132 remaining_energy = std::max ( theAngular[0].GetLabel() , theAngular[nEnergies-1].GetLabel() );
133 fresh = false;
134 }
135
136 //Cheating for small remaining_energy
137 //TEMPORAL SOLUTION
138 if ( nDiscreteEnergies == nEnergies )
139 {
140 remaining_energy = std::max ( remaining_energy , theAngular[nDiscreteEnergies-1].GetLabel() ); //Minimum Line
141 }
142 else
143 {
144 //G4double cont_min = theAngular[nDiscreteEnergies].GetLabel();
145 //if ( theAngular[nDiscreteEnergies].GetLabel() == 0.0 ) cont_min = theAngular[nDiscreteEnergies+1].GetLabel();
146 G4double cont_min=0.0;
147 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
148 {
149 cont_min = theAngular[j].GetLabel();
150 if ( theAngular[j].GetValue(0) != 0.0 ) break;
151 }
152 remaining_energy = std::max ( remaining_energy , std::min ( theAngular[nDiscreteEnergies-1].GetLabel() , cont_min ) ); //Minimum Line or grid
153 }
154//
155 G4double random = G4UniformRand();
156
157 G4double * running = new G4double[nEnergies+1];
158 running[0] = 0.0;
159
160 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
161 {
162 G4double delta = 0.0;
163 if ( theAngular[j].GetLabel() <= remaining_energy ) delta = theAngular[i].GetValue(0);
164 running[j+1] = running[j] + delta;
165 }
166 G4double tot_prob_DIS = running[ nDiscreteEnergies ];
167
168 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
169 {
170 G4double delta = 0.0;
171 G4double e_low = 0.0;
172 G4double e_high = 0.0;
173 if ( theAngular[j].GetLabel() <= remaining_energy ) delta = theAngular[j].GetValue(0);
174
175 //To calculate Prob. e_low and e_high should be in eV
176 //There are two case
177 //1:theAngular[nDiscreteEnergies].GetLabel() != 0.0
178 // delta should be used between j-1 and j
179 // At j = nDiscreteEnergies (the first) e_low should be set explicitly
180 if ( theAngular[j].GetLabel() != 0 )
181 {
182 if ( j == nDiscreteEnergies ) {
183 e_low = 0.0/eV;
184 } else {
185 e_low = theAngular[j-1].GetLabel()/eV;
186 }
187 e_high = theAngular[j].GetLabel()/eV;
188 }
189 //2:theAngular[nDiscreteEnergies].GetLabel() == 0.0
190 // delta should be used between j and j+1
191 if ( theAngular[j].GetLabel() == 0.0 ) {
192 e_low = theAngular[j].GetLabel()/eV;
193 if ( j != nEnergies-1 ) {
194 e_high = theAngular[j+1].GetLabel()/eV;
195 } else {
196 e_high = theAngular[j].GetLabel()/eV;
197 if ( theAngular[j].GetValue(0) != 0.0 ) {
198 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar: Unexpected non zero value of theAngular[nEnergies-1].GetValue(0)");
199 }
200 }
201 }
202
203 running[j+1] = running[j] + ( ( e_high - e_low ) * delta );
204 }
205 G4double tot_prob_CON = running[ nEnergies ] - running[ nDiscreteEnergies ];
206
207/*
208 For FPE debugging
209 if (tot_prob_DIS + tot_prob_CON == 0 ) {
210 G4cout << "TKDB tot_prob_DIS + tot_prob_CON " << tot_prob_DIS + tot_prob_CON << G4endl;
211 G4cout << "massCode " << massCode << G4endl;
212 G4cout << "nDiscreteEnergies " << nDiscreteEnergies << " nEnergies " << nEnergies << G4endl;
213 for ( int j = nDiscreteEnergies ; j < nEnergies ; j++ ) {
214 G4cout << j << " " << theAngular[j].GetLabel() << " " << theAngular[j].GetValue(0) << G4endl;
215 }
216 }
217*/
218 // Normalize random
219 random *= (tot_prob_DIS + tot_prob_CON);
220//2nd Judge Discrete or not This shoudl be relatively close to 1 For safty
221 if ( random <= ( tot_prob_DIS / ( tot_prob_DIS + tot_prob_CON ) ) || nDiscreteEnergies == nEnergies )
222 {
223// Discrete Emission
224 for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
225 {
226 //Here we should use i+1
227 if ( random < running[ j+1 ] )
228 {
229 it = j;
230 break;
231 }
232 }
233 fsEnergy = theAngular[ it ].GetLabel();
234
235 G4NeutronHPLegendreStore theStore(1);
236 theStore.Init(0,fsEnergy,nAngularParameters);
237 for (G4int j=0;j<nAngularParameters;j++)
238 {
239 theStore.SetCoeff(0,j,theAngular[it].GetValue(j));
240 }
241 // use it to sample.
242 cosTh = theStore.SampleMax(fsEnergy);
243 //Done
244 }
245 else
246 {
247// Continuous Emission
248 for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
249 {
250 //Here we should use i
251 if ( random < running[ j ] )
252 {
253 it = j;
254 break;
255 }
256 }
257
258 G4double x1 = running[it-1];
259 G4double x2 = running[it];
260
261 G4double y1 = 0.0;
262 if ( it != nDiscreteEnergies )
263 y1 = theAngular[it-1].GetLabel();
264 G4double y2 = theAngular[it].GetLabel();
265
266 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
267 random,x1,x2,y1,y2);
268
269 G4NeutronHPLegendreStore theStore(2);
270 theStore.Init(0,y1,nAngularParameters);
271 theStore.Init(1,y2,nAngularParameters);
272 theStore.SetManager(theManager);
273 for (G4int j=0;j<nAngularParameters;j++)
274 {
275 G4int itt = it;
276 if ( it == nDiscreteEnergies ) itt = it+1; //"This case "it-1" has data for Discrete, so we will use an extrpolate values it and it+1
277 if ( it == 0 )
278 {
279 //Safty for unexpected it = 0;
280 //G4cout << "110611 G4NeutronHPContAngularPar::Sample it = 0; invetigation required " << G4endl;
281 itt = it+1;
282 }
283 theStore.SetCoeff(0,j,theAngular[itt-1].GetValue(j));
284 theStore.SetCoeff(1,j,theAngular[itt].GetValue(j));
285 }
286 // use it to sample.
287 cosTh = theStore.SampleMax(fsEnergy);
288
289 //Done
290 }
291
292 //TK080711
293 remaining_energy -= fsEnergy;
294 //TK080711
295
296 //080801b
297 delete[] running;
298 //080801b
299 }
300 else
301 {
302 // Only continue, TK will clean up
303
304 //080714
305 if ( fresh == true )
306 {
307 remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
308 fresh = false;
309 }
310 //080714
311 G4double random = G4UniformRand();
312 G4double * running = new G4double[nEnergies];
313 running[0]=0;
314 G4double weighted = 0;
315 for(i=1; i<nEnergies; i++)
316 {
317/*
318 if(i!=0)
319 {
320 running[i]=running[i-1];
321 }
322 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
323 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
324 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
325 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
326 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
327 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
328*/
329
330 running[i]=running[i-1];
331 if ( remaining_energy >= theAngular[i].GetLabel() )
332 {
333 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
334 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
335 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
336 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
337 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
338 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
339 }
340 }
341 // cash the mean energy in this distribution
342 //080409 TKDB
343 if ( nEnergies == 1 || running[nEnergies-1] == 0 )
344 currentMeanEnergy = 0.0;
345 else
346 {
347 currentMeanEnergy = weighted/running[nEnergies-1];
348 }
349
350 //080409 TKDB
351 if ( nEnergies == 1 ) it = 0;
352
353 //080729
354 if ( running[nEnergies-1] != 0 )
355 {
356 for ( i = 1 ; i < nEnergies ; i++ )
357 {
358 it = i;
359 if ( random < running [ i ] / running [ nEnergies-1 ] ) break;
360 }
361 }
362
363 //080714
364 if ( running [ nEnergies-1 ] == 0 ) it = 0;
365 //080714
366
367 if (it<nDiscreteEnergies||it==0)
368 {
369 if(it == 0)
370 {
371 fsEnergy = theAngular[it].GetLabel();
372 G4NeutronHPLegendreStore theStore(1);
373 theStore.Init(0,fsEnergy,nAngularParameters);
374 for(i=0;i<nAngularParameters;i++)
375 {
376 theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
377 }
378 // use it to sample.
379 cosTh = theStore.SampleMax(fsEnergy);
380 }
381 else
382 {
383 G4double e1, e2;
384 e1 = theAngular[it-1].GetLabel();
385 e2 = theAngular[it].GetLabel();
386 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
387 random,
388 running[it-1]/running[nEnergies-1],
389 running[it]/running[nEnergies-1],
390 e1, e2);
391 // fill a Legendrestore
392 G4NeutronHPLegendreStore theStore(2);
393 theStore.Init(0,e1,nAngularParameters);
394 theStore.Init(1,e2,nAngularParameters);
395 for(i=0;i<nAngularParameters;i++)
396 {
397 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
398 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
399 }
400 // use it to sample.
401 theStore.SetManager(theManager);
402 cosTh = theStore.SampleMax(fsEnergy);
403 }
404 }
405 else // continuum contribution
406 {
407 G4double x1 = running[it-1]/running[nEnergies-1];
408 G4double x2 = running[it]/running[nEnergies-1];
409 G4double y1 = theAngular[it-1].GetLabel();
410 G4double y2 = theAngular[it].GetLabel();
411 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
412 random,x1,x2,y1,y2);
413 G4NeutronHPLegendreStore theStore(2);
414 theStore.Init(0,y1,nAngularParameters);
415 theStore.Init(1,y2,nAngularParameters);
416 theStore.SetManager(theManager);
417 for(i=0;i<nAngularParameters;i++)
418 {
419 theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
420 theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
421 }
422 // use it to sample.
423 cosTh = theStore.SampleMax(fsEnergy);
424 }
425 delete [] running;
426
427 //080714
428 remaining_energy -= fsEnergy;
429 //080714
430 }
431 }
432 else if(angularRep==2)
433 {
434 // first get the energy (already the right for this incoming energy)
435 G4int j;
436 G4double * running = new G4double[nEnergies];
437 running[0]=0;
438 G4double weighted = 0;
439 for(j=1; j<nEnergies; j++)
440 {
441 if(j!=0) running[j]=running[j-1];
442 running[j] += theInt.GetBinIntegral(theManager.GetScheme(j-1),
443 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
444 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
445 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(j-1),
446 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
447 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
448 }
449 // cash the mean energy in this distribution
450 //080409 TKDB
451 //currentMeanEnergy = weighted/running[nEnergies-1];
452 if ( nEnergies == 1 )
453 currentMeanEnergy = 0.0;
454 else
455 currentMeanEnergy = weighted/running[nEnergies-1];
456
457 G4int itt(0);
458 G4double randkal = G4UniformRand();
459 //080409 TKDB
460 //for(i=0; i<nEnergies; i++)
461 for(j=1; j<nEnergies; j++)
462 {
463 itt = j;
464 if(randkal<running[j]/running[nEnergies-1]) break;
465 }
466
467 // interpolate the secondary energy.
468 G4double x, x1,x2,y1,y2;
469 if(itt==0) itt=1;
470 x = randkal*running[nEnergies-1];
471 x1 = running[itt-1];
472 x2 = running[itt];
473 G4double compoundFraction;
474 // interpolate energy
475 y1 = theAngular[itt-1].GetLabel();
476 y2 = theAngular[itt].GetLabel();
477 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(itt-1),
478 x, x1,x2,y1,y2);
479 // for theta interpolate the compoundFractions
480 G4double cLow = theAngular[itt-1].GetValue(1);
481 G4double cHigh = theAngular[itt].GetValue(1);
482 compoundFraction = theInt.Interpolate(theManager.GetScheme(itt),
483 fsEnergy, y1, y2, cLow,cHigh);
484 delete [] running;
485
486 // get cosTh
487 G4double incidentEnergy = anEnergy;
488 G4double incidentMass = G4Neutron::Neutron()->GetPDGMass();
489 G4double productEnergy = fsEnergy;
490 G4double productMass = result->GetMass();
491 G4int targetZ = G4int(theTargetCode/1000);
492 G4int targetA = G4int(theTargetCode-1000*targetZ);
493 // To correspond to natural composition (-nat-) data files.
494 if ( targetA == 0 )
495 targetA = G4int ( theTarget->GetMass()/amu_c2 + 0.5 );
496 G4double targetMass = theTarget->GetMass();
497 G4int residualA = targetA+1-A;
498 G4int residualZ = targetZ-Z;
499 G4double residualMass = residualZ*G4Proton::Proton()->GetPDGMass();
500 residualMass +=(residualA-residualZ)*G4Neutron::Neutron()->GetPDGMass();
501 residualMass -= G4NucleiProperties::GetBindingEnergy( residualA , residualZ );
502 G4NeutronHPKallbachMannSyst theKallbach(compoundFraction,
503 incidentEnergy, incidentMass,
504 productEnergy, productMass,
505 residualMass, residualA, residualZ,
506 targetMass, targetA, targetZ);
507 cosTh = theKallbach.Sample(anEnergy);
508 }
509 else if(angularRep>10&&angularRep<16)
510 {
511 G4double random = G4UniformRand();
512 G4double * running = new G4double[nEnergies];
513 running[0]=0;
514 G4double weighted = 0;
515 for(i=1; i<nEnergies; i++)
516 {
517 if(i!=0) running[i]=running[i-1];
518 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
519 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
520 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
521 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
522 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
523 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
524 }
525 // cash the mean energy in this distribution
526 //currentMeanEnergy = weighted/running[nEnergies-1];
527 if ( nEnergies == 1 )
528 currentMeanEnergy = 0.0;
529 else
530 currentMeanEnergy = weighted/running[nEnergies-1];
531
532 //080409 TKDB
533 if ( nEnergies == 1 ) it = 0;
534 //for(i=0; i<nEnergies; i++)
535 for(i=1; i<nEnergies; i++)
536 {
537 it = i;
538 if(random<running[i]/running[nEnergies-1]) break;
539 }
540 if(it<nDiscreteEnergies||it==0)
541 {
542 if(it==0)
543 {
544 fsEnergy = theAngular[0].GetLabel();
545 G4NeutronHPVector theStore;
546 G4int aCounter = 0;
547 for(G4int j=1; j<nAngularParameters; j+=2)
548 {
549 theStore.SetX(aCounter, theAngular[0].GetValue(j));
550 theStore.SetY(aCounter, theAngular[0].GetValue(j+1));
551 aCounter++;
552 }
553 G4InterpolationManager aMan;
554 aMan.Init(angularRep-10, nAngularParameters-1);
555 theStore.SetInterpolationManager(aMan);
556 cosTh = theStore.Sample();
557 }
558 else
559 {
560 fsEnergy = theAngular[it].GetLabel();
561 G4NeutronHPVector theStore;
562 G4InterpolationManager aMan;
563 aMan.Init(angularRep-10, nAngularParameters-1);
564 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
565 G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
566 G4int aCounter = 0;
567 for(G4int j=1; j<nAngularParameters; j+=2)
568 {
569 theStore.SetX(aCounter, theAngular[it].GetValue(j));
570 theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
571 random,
572 running[it-1]/running[nEnergies-1],
573 running[it]/running[nEnergies-1],
574 theAngular[it-1].GetValue(j+1),
575 theAngular[it].GetValue(j+1)));
576 aCounter++;
577 }
578 cosTh = theStore.Sample();
579 }
580 }
581 else
582 {
583 G4double x1 = running[it-1]/running[nEnergies-1];
584 G4double x2 = running[it]/running[nEnergies-1];
585 G4double y1 = theAngular[it-1].GetLabel();
586 G4double y2 = theAngular[it].GetLabel();
587 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
588 random,x1,x2,y1,y2);
589 G4NeutronHPVector theBuff1;
590 G4NeutronHPVector theBuff2;
591 G4InterpolationManager aMan;
592 aMan.Init(angularRep-10, nAngularParameters-1);
593// theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
594// theBuff2.SetInterpolationManager(aMan); // Store interpolates f(costh)
595// Bug Report #1366 from L. Russell
596 //for(i=0; i<nAngularParameters; i++) // i=1 ist wichtig!
597 //{
598 // theBuff1.SetX(i, theAngular[it-1].GetValue(i));
599 // theBuff1.SetY(i, theAngular[it-1].GetValue(i+1));
600 // theBuff2.SetX(i, theAngular[it].GetValue(i));
601 // theBuff2.SetY(i, theAngular[it].GetValue(i+1));
602 // i++;
603 //}
604 {
605 G4int j;
606 for(i=0,j=1; i<nAngularParameters; i++,j+=2)
607 {
608 theBuff1.SetX(i, theAngular[it-1].GetValue(j));
609 theBuff1.SetY(i, theAngular[it-1].GetValue(j+1));
610 theBuff2.SetX(i, theAngular[it].GetValue(j));
611 theBuff2.SetY(i, theAngular[it].GetValue(j+1));
612 }
613 }
614 G4NeutronHPVector theStore;
615 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
616 x1 = y1;
617 x2 = y2;
618 G4double x, y;
619 //for(i=0;i<theBuff1.GetVectorLength(); i++);
620 for(i=0;i<theBuff1.GetVectorLength(); i++)
621 {
622 x = theBuff1.GetX(i); // costh binning identical
623 y1 = theBuff1.GetY(i);
624 y2 = theBuff2.GetY(i);
625 y = theInt.Interpolate(theManager.GetScheme(it),
626 fsEnergy, theAngular[it-1].GetLabel(),
627 theAngular[it].GetLabel(), y1, y2);
628 theStore.SetX(i, x);
629 theStore.SetY(i, y);
630 }
631 cosTh = theStore.Sample();
632 }
633 delete [] running;
634 }
635 else
636 {
637 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPContAngularPar::Sample: Unknown angular representation");
638 }
639 result->SetKineticEnergy(fsEnergy);
640 G4double phi = twopi*G4UniformRand();
641 G4double theta = std::acos(cosTh);
642 G4double sinth = std::sin(theta);
643 G4double mtot = result->GetTotalMomentum();
644 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
645 result->SetMomentum(tempVector);
646// return the result.
647 return result;
648 }
G4InterpolationScheme
G4InterpolationScheme
Definition: G4InterpolationScheme.hh:39
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition: G4InterpolationManager.hh:85
G4InterpolationManager::GetScheme
G4InterpolationScheme GetScheme(G4int index) const
Definition: G4InterpolationManager.hh:132
G4InterpolationManager::GetInverseScheme
G4InterpolationScheme GetInverseScheme(G4int index)
Definition: G4InterpolationManager.hh:149
G4NeutronHPInterpolator::GetWeightedBinIntegral
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
Definition: G4NeutronHPInterpolator.cc:95
G4NeutronHPInterpolator::GetBinIntegral
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
Definition: G4NeutronHPInterpolator.cc:33
G4NeutronHPVector::SetX
void SetX(G4int i, G4double e)
Definition: G4NeutronHPVector.hh:97
G4NeutronHPVector::GetVectorLength
G4int GetVectorLength() const
Definition: G4NeutronHPVector.hh:193
G4NeutronHPVector::GetX
G4double GetX(G4int i) const
Definition: G4NeutronHPVector.hh:121
G4NeutronHPVector::GetY
G4double GetY(G4double x)
Definition: G4NeutronHPVector.hh:192
G4NeutronHPVector::SetInterpolationManager
void SetInterpolationManager(const G4InterpolationManager &aManager)
Definition: G4NeutronHPVector.hh:477
G4NeutronHPVector::SetY
void SetY(G4int i, G4double x)
Definition: G4NeutronHPVector.hh:107
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4NucleiProperties::GetBindingEnergy
static G4double GetBindingEnergy(const G4int A, const G4int Z)
Definition: G4NucleiProperties.cc:223
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:63
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:166
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:123
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:129
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:155
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:147
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95

Referenced by G4NeutronHPContEnergyAngular::Sample().

◆ SetInterpolation()

void G4NeutronHPContAngularPar::SetInterpolation ( G4int  theInterpolation)
inline

Definition at line 76 of file G4NeutronHPContAngularPar.hh.

77 {
78 theManager.Init(theInterpolation, nEnergies); // one range only
79 }

Referenced by G4NeutronHPContEnergyAngular::Init().

◆ SetPrimary()

void G4NeutronHPContAngularPar::SetPrimary ( G4ReactionProduct *  aPrimary)
inline

Definition at line 64 of file G4NeutronHPContAngularPar.hh.

65 {
66 thePrimary = aPrimary;
67 }

Referenced by G4NeutronHPContEnergyAngular::Sample().

◆ SetTarget()

void G4NeutronHPContAngularPar::SetTarget ( G4ReactionProduct *  aTarget)
inline

Definition at line 69 of file G4NeutronHPContAngularPar.hh.

70 {
71 theTarget = aTarget;
72 }

Referenced by G4NeutronHPContEnergyAngular::Sample().

◆ SetTargetCode()

void G4NeutronHPContAngularPar::SetTargetCode ( G4double  aTargetCode)
inline

Definition at line 74 of file G4NeutronHPContAngularPar.hh.

74{ theTargetCode = aTargetCode; }

Referenced by G4NeutronHPContEnergyAngular::Sample().

The documentation for this class was generated from the following files: