Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronHPContAngularPar.cc
Go to the documentation of this file.
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 "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4NeutronHPLegendreStore.hh"
44 #include "G4Gamma.hh"
45 #include "G4Electron.hh"
46 #include "G4Positron.hh"
47 #include "G4Neutron.hh"
48 #include "G4Proton.hh"
49 #include "G4Deuteron.hh"
50 #include "G4Triton.hh"
51 #include "G4He3.hh"
52 #include "G4Alpha.hh"
53 #include "G4NeutronHPVector.hh"
54 #include "G4NucleiProperties.hh"
55 #include "G4NeutronHPKallbachMannSyst.hh"
56 #include "G4ParticleTable.hh"
57 
58  void G4NeutronHPContAngularPar::Init(std::ifstream & aDataFile)
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  }
72 
73  G4ReactionProduct *
74  G4NeutronHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /*targetMass*/,
75  G4int angularRep, G4int /*interpolE*/ )
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  }
Generated on Sat May 25 2013 14:34:03 for Geant4 by   doxygen 1.8.4