Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPContAngularPar.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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
40 //
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
45 #include "G4Gamma.hh"
46 #include "G4Electron.hh"
47 #include "G4Positron.hh"
48 #include "G4Neutron.hh"
49 #include "G4Proton.hh"
50 #include "G4Deuteron.hh"
51 #include "G4Triton.hh"
52 #include "G4He3.hh"
53 #include "G4Alpha.hh"
54 #include "G4ParticleHPVector.hh"
55 #include "G4NucleiProperties.hh"
57 #include "G4IonTable.hh"
58 #include <set>
59 
61 {
62  theAngular = 0;
63  if ( fCache.Get() == NULL ) cacheInit();
64  fCache.Get()->currentMeanEnergy = -2;
65  fCache.Get()->fresh = true;
66  adjustResult = true;
67  if ( getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) adjustResult = false;
68 
69  theMinEner = DBL_MAX;
70  theMaxEner = -DBL_MAX;
71  theProjectile = projectile;
72 
73  theEnergy = 0.0;
74  nEnergies = 0;
75  nDiscreteEnergies = 0;
76  nAngularParameters = 0;
77 }
78 
79  void G4ParticleHPContAngularPar::Init(std::istream & aDataFile, G4ParticleDefinition* projectile)
80  {
81  adjustResult = true;
82  if ( getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) adjustResult = false;
83 
84  theProjectile = projectile;
85 
86  aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
87  /*if( getenv("G4PHPTEST") )*/
88  theEnergy *= eV;
89  theAngular = new G4ParticleHPList [nEnergies];
90  for(G4int i=0; i<nEnergies; i++)
91  {
92  G4double sEnergy;
93  aDataFile >> sEnergy;
94  sEnergy*=eV;
95  theAngular[i].SetLabel(sEnergy);
96  theAngular[i].Init(aDataFile, nAngularParameters, 1.);
97  theMinEner = std::min(theMinEner,sEnergy);
98  theMaxEner = std::max(theMaxEner,sEnergy);
99  }
100  }
101 
103  G4ParticleHPContAngularPar::Sample(G4double anEnergy, G4double massCode, G4double /*targetMass*/,
104  G4int angularRep, G4int /*interpolE*/ )
105  {
106  if( getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample " << anEnergy << " " << massCode << " " << angularRep << G4endl; //GDEB
107  if ( fCache.Get() == NULL ) cacheInit();
109  G4int Z = static_cast<G4int>(massCode/1000);
110  G4int A = static_cast<G4int>(massCode-1000*Z);
111  if(massCode==0)
112  {
113  result->SetDefinition(G4Gamma::Gamma());
114  }
115  else if(A==0)
116  {
118  if(Z==1) result->SetDefinition(G4Positron::Positron());
119  }
120  else if(A==1)
121  {
123  if(Z==1) result->SetDefinition(G4Proton::Proton());
124  }
125  else if(A==2)
126  {
128  }
129  else if(A==3)
130  {
131  result->SetDefinition(G4Triton::Triton());
132  if(Z==2) result->SetDefinition(G4He3::He3());
133  }
134  else if(A==4)
135  {
136  result->SetDefinition(G4Alpha::Alpha());
137  if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar: Unknown ion case 1");
138  }
139  else
140  {
141  result->SetDefinition(G4IonTable::GetIonTable()->GetIon(Z,A,0));
142  }
143  G4int i(0);
144  G4int it(0);
145  G4double fsEnergy(0);
146  G4double cosTh(0);
147 
148  if( angularRep == 1 )
149  {
150 // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
151  //if (interpolE == 2)
152 //110609 above was wrong interupition, pointed out by E.Mendoza and D.Cano (CIMAT)
153 //Following are reviesd version written by T.Koi (SLAC)
154  if ( nDiscreteEnergies != 0 )
155  {
156 
157 //1st check remaining_energy
158 // if this is the first set it. (How?)
159  if ( fCache.Get()->fresh == true )
160  {
161  //Discrete Lines, larger energies come first
162  //Continues Emssions, low to high LAST
163  fCache.Get()->remaining_energy = std::max ( theAngular[0].GetLabel() , theAngular[nEnergies-1].GetLabel() );
164  fCache.Get()->fresh = false;
165  }
166 
167  //Cheating for small remaining_energy
168  //TEMPORAL SOLUTION
169  if ( nDiscreteEnergies == nEnergies )
170  {
171  fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , theAngular[nDiscreteEnergies-1].GetLabel() ); //Minimum Line
172  }
173  else
174  {
175  //G4double cont_min = theAngular[nDiscreteEnergies].GetLabel();
176  //if ( theAngular[nDiscreteEnergies].GetLabel() == 0.0 ) cont_min = theAngular[nDiscreteEnergies+1].GetLabel();
177  G4double cont_min=0.0;
178  for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
179  {
180  cont_min = theAngular[j].GetLabel();
181  if ( theAngular[j].GetValue(0) != 0.0 ) break;
182  }
183  fCache.Get()->remaining_energy = std::max ( fCache.Get()->remaining_energy , std::min ( theAngular[nDiscreteEnergies-1].GetLabel() , cont_min ) ); //Minimum Line or grid
184  }
185 //
186  G4double random = G4UniformRand();
187 
188  G4double * running = new G4double[nEnergies+1];
189  running[0] = 0.0;
190 
191  for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
192  {
193  G4double delta = 0.0;
194  if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[i].GetValue(0);
195  running[j+1] = running[j] + delta;
196  }
197  G4double tot_prob_DIS = running[ nDiscreteEnergies ];
198 
199  for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
200  {
201  G4double delta = 0.0;
202  G4double e_low = 0.0;
203  G4double e_high = 0.0;
204  if ( theAngular[j].GetLabel() <= fCache.Get()->remaining_energy ) delta = theAngular[j].GetValue(0);
205 
206  //To calculate Prob. e_low and e_high should be in eV
207  //There are two case
208  //1:theAngular[nDiscreteEnergies].GetLabel() != 0.0
209  // delta should be used between j-1 and j
210  // At j = nDiscreteEnergies (the first) e_low should be set explicitly
211  if ( theAngular[j].GetLabel() != 0 )
212  {
213  if ( j == nDiscreteEnergies ) {
214  e_low = 0.0/eV;
215  } else {
216  e_low = theAngular[j-1].GetLabel()/eV;
217  }
218  e_high = theAngular[j].GetLabel()/eV;
219  }
220  //2:theAngular[nDiscreteEnergies].GetLabel() == 0.0
221  // delta should be used between j and j+1
222  if ( theAngular[j].GetLabel() == 0.0 ) {
223  e_low = theAngular[j].GetLabel()/eV;
224  if ( j != nEnergies-1 ) {
225  e_high = theAngular[j+1].GetLabel()/eV;
226  } else {
227  e_high = theAngular[j].GetLabel()/eV;
228  if ( theAngular[j].GetValue(0) != 0.0 ) {
229  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar: Unexpected non zero value of theAngular[nEnergies-1].GetValue(0)");
230  }
231  }
232  }
233 
234  running[j+1] = running[j] + ( ( e_high - e_low ) * delta );
235  }
236  G4double tot_prob_CON = running[ nEnergies ] - running[ nDiscreteEnergies ];
237 
238 /*
239  For FPE debugging
240  if (tot_prob_DIS + tot_prob_CON == 0 ) {
241  G4cout << "TKDB tot_prob_DIS + tot_prob_CON " << tot_prob_DIS + tot_prob_CON << G4endl;
242  G4cout << "massCode " << massCode << G4endl;
243  G4cout << "nDiscreteEnergies " << nDiscreteEnergies << " nEnergies " << nEnergies << G4endl;
244  for ( int j = nDiscreteEnergies ; j < nEnergies ; j++ ) {
245  G4cout << j << " " << theAngular[j].GetLabel() << " " << theAngular[j].GetValue(0) << G4endl;
246  }
247  }
248 */
249  // Normalize random
250  random *= (tot_prob_DIS + tot_prob_CON);
251 //2nd Judge Discrete or not This shoudl be relatively close to 1 For safty
252  if ( random <= ( tot_prob_DIS / ( tot_prob_DIS + tot_prob_CON ) ) || nDiscreteEnergies == nEnergies )
253  {
254 // Discrete Emission
255  for ( G4int j = 0 ; j < nDiscreteEnergies ; j++ )
256  {
257  //Here we should use i+1
258  if ( random < running[ j+1 ] )
259  {
260  it = j;
261  break;
262  }
263  }
264  fsEnergy = theAngular[ it ].GetLabel();
265 
266  G4ParticleHPLegendreStore theStore(1);
267  theStore.Init(0,fsEnergy,nAngularParameters);
268  for (G4int j=0;j<nAngularParameters;j++)
269  {
270  theStore.SetCoeff(0,j,theAngular[it].GetValue(j));
271  }
272  // use it to sample.
273  cosTh = theStore.SampleMax(fsEnergy);
274  //Done
275  }
276  else
277  {
278 // Continuous Emission
279  for ( G4int j = nDiscreteEnergies ; j < nEnergies ; j++ )
280  {
281  //Here we should use i
282  if ( random < running[ j ] )
283  {
284  it = j;
285  break;
286  }
287  }
288 
289  G4double x1 = running[it-1];
290  G4double x2 = running[it];
291 
292  G4double y1 = 0.0;
293  if ( it != nDiscreteEnergies )
294  y1 = theAngular[it-1].GetLabel();
295  G4double y2 = theAngular[it].GetLabel();
296 
297  fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
298  random,x1,x2,y1,y2);
299 
300  G4ParticleHPLegendreStore theStore(2);
301  theStore.Init(0,y1,nAngularParameters);
302  theStore.Init(1,y2,nAngularParameters);
303  theStore.SetManager(theManager);
304  for (G4int j=0;j<nAngularParameters;j++)
305  {
306  G4int itt = it;
307  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
308  if ( it == 0 )
309  {
310  //Safty for unexpected it = 0;
311  //G4cout << "110611 G4ParticleHPContAngularPar::Sample it = 0; invetigation required " << G4endl;
312  itt = it+1;
313  }
314  theStore.SetCoeff(0,j,theAngular[itt-1].GetValue(j));
315  theStore.SetCoeff(1,j,theAngular[itt].GetValue(j));
316  }
317  // use it to sample.
318  cosTh = theStore.SampleMax(fsEnergy);
319 
320  //Done
321  }
322 
323  //TK080711
324  if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
325  //TK080711
326 
327  //080801b
328  delete[] running;
329  //080801b
330  }
331  else
332  {
333  // Only continue, TK will clean up
334 
335  //080714
336  if ( fCache.Get()->fresh == true )
337  {
338  fCache.Get()->remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
339  fCache.Get()->fresh = false;
340  }
341  //080714
342  G4double random = G4UniformRand();
343  G4double * running = new G4double[nEnergies];
344  running[0]=0;
345  G4double weighted = 0;
346  for(i=1; i<nEnergies; i++)
347  {
348 /*
349  if(i!=0)
350  {
351  running[i]=running[i-1];
352  }
353  running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
354  theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
355  theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
356  weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
357  theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
358  theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
359 */
360 
361  running[i]=running[i-1];
362  if ( fCache.Get()->remaining_energy >= theAngular[i].GetLabel() )
363  {
364  running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
365  theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
366  theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
367  weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
368  theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
369  theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
370  }
371  }
372  // cash the mean energy in this distribution
373  //080409 TKDB
374  if ( nEnergies == 1 || running[nEnergies-1] == 0 )
375  fCache.Get()->currentMeanEnergy = 0.0;
376  else
377  {
378  fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
379  }
380 
381  //080409 TKDB
382  if ( nEnergies == 1 ) it = 0;
383 
384  //080729
385  if ( running[nEnergies-1] != 0 )
386  {
387  for ( i = 1 ; i < nEnergies ; i++ )
388  {
389  it = i;
390  if ( random < running [ i ] / running [ nEnergies-1 ] ) break;
391  }
392  }
393 
394  //080714
395  if ( running [ nEnergies-1 ] == 0 ) it = 0;
396  //080714
397 
398  if (it<nDiscreteEnergies||it==0)
399  {
400  if(it == 0)
401  {
402  fsEnergy = theAngular[it].GetLabel();
403  G4ParticleHPLegendreStore theStore(1);
404  theStore.Init(0,fsEnergy,nAngularParameters);
405  for(i=0;i<nAngularParameters;i++)
406  {
407  theStore.SetCoeff(0,i,theAngular[it].GetValue(i));
408  }
409  // use it to sample.
410  cosTh = theStore.SampleMax(fsEnergy);
411  }
412  else
413  {
414  G4double e1, e2;
415  e1 = theAngular[it-1].GetLabel();
416  e2 = theAngular[it].GetLabel();
417  fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
418  random,
419  running[it-1]/running[nEnergies-1],
420  running[it]/running[nEnergies-1],
421  e1, e2);
422  // fill a Legendrestore
423  G4ParticleHPLegendreStore theStore(2);
424  theStore.Init(0,e1,nAngularParameters);
425  theStore.Init(1,e2,nAngularParameters);
426  for(i=0;i<nAngularParameters;i++)
427  {
428  theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
429  theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
430  }
431  // use it to sample.
432  theStore.SetManager(theManager);
433  cosTh = theStore.SampleMax(fsEnergy);
434  }
435  }
436  else // continuum contribution
437  {
438  G4double x1 = running[it-1]/running[nEnergies-1];
439  G4double x2 = running[it]/running[nEnergies-1];
440  G4double y1 = theAngular[it-1].GetLabel();
441  G4double y2 = theAngular[it].GetLabel();
442  fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
443  random,x1,x2,y1,y2);
444  G4ParticleHPLegendreStore theStore(2);
445  theStore.Init(0,y1,nAngularParameters);
446  theStore.Init(1,y2,nAngularParameters);
447  theStore.SetManager(theManager);
448  for(i=0;i<nAngularParameters;i++)
449  {
450  theStore.SetCoeff(0,i,theAngular[it-1].GetValue(i));
451  theStore.SetCoeff(1,i,theAngular[it].GetValue(i));
452  }
453  // use it to sample.
454  cosTh = theStore.SampleMax(fsEnergy);
455  }
456  delete [] running;
457 
458  //080714
459  if( adjustResult ) fCache.Get()->remaining_energy -= fsEnergy;
460  //080714
461  }
462  }
463  else if(angularRep==2)
464  {
465  // first get the energy (already the right for this incoming energy)
466  G4int j;
467  G4double * running = new G4double[nEnergies];
468  running[0]=0;
469  G4double weighted = 0;
470  if( getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample nEnergies " << nEnergies << G4endl;
471  for(j=1; j<nEnergies; j++)
472  {
473  if(j!=0) running[j]=running[j-1];
474  running[j] += theInt.GetBinIntegral(theManager.GetScheme(j-1),
475  theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
476  theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
477  weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(j-1),
478  theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
479  theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
480  if( getenv("G4PHPTEST") ) G4cout << " G4ParticleHPContAngularPar::Sample " << j << " running " << running[j]
481  << " " << theManager.GetScheme(j-1) << " " << theAngular[j-1].GetLabel() << " " << theAngular[j].GetLabel() << " " << theAngular[j-1].GetValue(0) << " " << theAngular[j].GetValue(0) << G4endl; //GDEB
482  }
483  // cash the mean energy in this distribution
484  //080409 TKDB
485  //currentMeanEnergy = weighted/running[nEnergies-1];
486  if ( nEnergies == 1 )
487  fCache.Get()->currentMeanEnergy = 0.0;
488  else
489  fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
490 
491  G4int itt(0);
492  G4double randkal = G4UniformRand();
493  //080409 TKDB
494  //for(i=0; i<nEnergies; i++)
495  for(j=1; j<nEnergies; j++)
496  {
497  itt = j;
498  if(randkal<running[j]/running[nEnergies-1]) break;
499  }
500 
501  // interpolate the secondary energy.
502  G4double x, x1,x2,y1,y2;
503  if(itt==0) itt=1;
504  x = randkal*running[nEnergies-1];
505  x1 = running[itt-1];
506  x2 = running[itt];
507  G4double compoundFraction;
508  // interpolate energy
509  y1 = theAngular[itt-1].GetLabel();
510  y2 = theAngular[itt].GetLabel();
511  fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(itt-1),
512  x, x1,x2,y1,y2);
513  if( getenv("G4PHPTEST") ) G4cout << itt << " G4particleHPContAngularPar fsEnergy " << fsEnergy << " " << theManager.GetInverseScheme(itt-1) << " x " << x << " " << x1 << " " << x2 << " y " << y1 << " " << y2 << G4endl; //GDEB
514  // for theta interpolate the compoundFractions
515  G4double cLow = theAngular[itt-1].GetValue(1);
516  G4double cHigh = theAngular[itt].GetValue(1);
517  compoundFraction = theInt.Interpolate(theManager.GetScheme(itt),
518  fsEnergy, y1, y2, cLow,cHigh);
519  if( getenv("G4PHPTEST") ) G4cout << itt << " G4particleHPContAngularPar compoundFraction " << compoundFraction << " E " << fsEnergy << " " << theManager.GetScheme(itt) << " ener " << fsEnergy << " y " << y1 << " " << y2 << " cLH " << cLow << " " << cHigh << G4endl; //GDEB
520  delete [] running;
521 
522  // get cosTh
523  G4double incidentEnergy = anEnergy;
524  G4double incidentMass = theProjectile->GetPDGMass();
525  G4double productEnergy = fsEnergy;
526  G4double productMass = result->GetMass();
527  G4int targetZ = G4int(fCache.Get()->theTargetCode/1000);
528  G4int targetA = G4int(fCache.Get()->theTargetCode-1000*targetZ);
529  // To correspond to natural composition (-nat-) data files.
530  if ( targetA == 0 )
531  targetA = G4int ( fCache.Get()->theTarget->GetMass()/amu_c2 + 0.5 );
532  G4double targetMass = fCache.Get()->theTarget->GetMass();
533  G4int residualA = targetA+1-A;
534  G4int residualZ = targetZ-Z;
535  G4double residualMass = residualZ*G4Proton::Proton()->GetPDGMass();
536  residualMass +=(residualA-residualZ)*theProjectile->GetPDGMass();
537  residualMass -= G4NucleiProperties::GetBindingEnergy( residualA , residualZ );
538  G4ParticleHPKallbachMannSyst theKallbach(compoundFraction,
539  incidentEnergy, incidentMass,
540  productEnergy, productMass,
541  residualMass, residualA, residualZ,
542  targetMass, targetA, targetZ);
543  cosTh = theKallbach.Sample(anEnergy);
544  if( getenv("G4PHPTEST") ) G4cout << " G4ParticleHPKallbachMannSyst::Sample resulttest " << cosTh << G4endl; //GDEB
545  }
546  else if(angularRep>10&&angularRep<16)
547  {
548  G4double random = G4UniformRand();
549  G4double * running = new G4double[nEnergies];
550  running[0]=0;
551  G4double weighted = 0;
552  for(i=1; i<nEnergies; i++)
553  {
554  if(i!=0) running[i]=running[i-1];
555  running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
556  theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
557  theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
558  weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
559  theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
560  theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
561  }
562  // cash the mean energy in this distribution
563  //currentMeanEnergy = weighted/running[nEnergies-1];
564  if ( nEnergies == 1 )
565  fCache.Get()->currentMeanEnergy = 0.0;
566  else
567  fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
568 
569  //080409 TKDB
570  if ( nEnergies == 1 ) it = 0;
571  //for(i=0; i<nEnergies; i++)
572  for(i=1; i<nEnergies; i++)
573  {
574  it = i;
575  if(random<running[i]/running[nEnergies-1]) break;
576  }
577  if(it<nDiscreteEnergies||it==0)
578  {
579  if(it==0)
580  {
581  fsEnergy = theAngular[0].GetLabel();
582  G4ParticleHPVector theStore;
583  G4int aCounter = 0;
584  for(G4int j=1; j<nAngularParameters; j+=2)
585  {
586  theStore.SetX(aCounter, theAngular[0].GetValue(j));
587  theStore.SetY(aCounter, theAngular[0].GetValue(j+1));
588  aCounter++;
589  }
591  aMan.Init(angularRep-10, nAngularParameters-1);
592  theStore.SetInterpolationManager(aMan);
593  cosTh = theStore.Sample();
594  }
595  else
596  {
597  fsEnergy = theAngular[it].GetLabel();
598  G4ParticleHPVector theStore;
600  aMan.Init(angularRep-10, nAngularParameters-1);
601  theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
602  G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
603  G4int aCounter = 0;
604  for(G4int j=1; j<nAngularParameters; j+=2)
605  {
606  theStore.SetX(aCounter, theAngular[it].GetValue(j));
607  theStore.SetY(aCounter, theInt.Interpolate(currentScheme,
608  random,
609  running[it-1]/running[nEnergies-1],
610  running[it]/running[nEnergies-1],
611  theAngular[it-1].GetValue(j+1),
612  theAngular[it].GetValue(j+1)));
613  aCounter++;
614  }
615  cosTh = theStore.Sample();
616  }
617  }
618  else
619  {
620  G4double x1 = running[it-1]/running[nEnergies-1];
621  G4double x2 = running[it]/running[nEnergies-1];
622  G4double y1 = theAngular[it-1].GetLabel();
623  G4double y2 = theAngular[it].GetLabel();
624  fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
625  random,x1,x2,y1,y2);
626  G4ParticleHPVector theBuff1;
627  G4ParticleHPVector theBuff2;
629  aMan.Init(angularRep-10, nAngularParameters-1);
630 // theBuff1.SetInterpolationManager(aMan); // Store interpolates f(costh)
631 // theBuff2.SetInterpolationManager(aMan); // Store interpolates f(costh)
632 // Bug Report #1366 from L. Russell
633  //for(i=0; i<nAngularParameters; i++) // i=1 ist wichtig!
634  //{
635  // theBuff1.SetX(i, theAngular[it-1].GetValue(i));
636  // theBuff1.SetY(i, theAngular[it-1].GetValue(i+1));
637  // theBuff2.SetX(i, theAngular[it].GetValue(i));
638  // theBuff2.SetY(i, theAngular[it].GetValue(i+1));
639  // i++;
640  //}
641  {
642  G4int j;
643  for(i=0,j=1; i<nAngularParameters; i++,j+=2)
644  {
645  theBuff1.SetX(i, theAngular[it-1].GetValue(j));
646  theBuff1.SetY(i, theAngular[it-1].GetValue(j+1));
647  theBuff2.SetX(i, theAngular[it].GetValue(j));
648  theBuff2.SetY(i, theAngular[it].GetValue(j+1));
649  }
650  }
651  G4ParticleHPVector theStore;
652  theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
653  x1 = y1;
654  x2 = y2;
655  G4double x, y;
656  //for(i=0;i<theBuff1.GetVectorLength(); i++);
657  for(i=0;i<theBuff1.GetVectorLength(); i++)
658  {
659  x = theBuff1.GetX(i); // costh binning identical
660  y1 = theBuff1.GetY(i);
661  y2 = theBuff2.GetY(i);
662  y = theInt.Interpolate(theManager.GetScheme(it),
663  fsEnergy, theAngular[it-1].GetLabel(),
664  theAngular[it].GetLabel(), y1, y2);
665  theStore.SetX(i, x);
666  theStore.SetY(i, y);
667  }
668  cosTh = theStore.Sample();
669  }
670  delete [] running;
671  }
672  else
673  {
674  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::Sample: Unknown angular representation");
675  }
676  result->SetKineticEnergy(fsEnergy);
677  G4double phi = twopi*G4UniformRand();
678  G4double theta = std::acos(cosTh);
679  G4double sinth = std::sin(theta);
680  G4double mtot = result->GetTotalMomentum();
681  G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
682  result->SetMomentum(tempVector);
683 // return the result.
684  return result;
685  }
686 
687 
688 #define MERGE_NEW
689 
691 {
692 
693  //----- Discrete energies: store own energies in a map for faster searching
694  G4int ie;
695  for(ie=0; ie<nDiscreteEnergies; ie++) {
696  theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
697  }
698  if( !angParPrev ) return;
699 
700  //----- Discrete energies: use energies that appear in one or another
701  for(ie=0; ie<nDiscreteEnergies; ie++) {
702  theDiscreteEnergies.insert(theAngular[ie].GetLabel());
703  }
704  G4int nDiscreteEnergiesPrev = angParPrev->GetNDiscreteEnergies();
705  for(ie=0; ie<nDiscreteEnergiesPrev; ie++) {
706  theDiscreteEnergies.insert(angParPrev->theAngular[ie].GetLabel());
707  }
708 
709  //--- Get the values for which interpolation will be done : all energies of this and previous ContAngularPar
710  for(ie=nDiscreteEnergies; ie<nEnergies; ie++) {
711  G4double ener = theAngular[ie].GetLabel();
712  G4double enerT = (ener-theMinEner)/(theMaxEner-theMinEner);
713  theEnergiesTransformed.insert(enerT);
714  //- if( getenv("G4PHPTEST2") ) G4cout <<this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed1 " << enerT << G4endl; //GDEB
715  }
716  G4int nEnergiesPrev = angParPrev->GetNEnergies();
717  G4double minEnerPrev = angParPrev->GetMinEner();
718  G4double maxEnerPrev = angParPrev->GetMaxEner();
719  for(ie=nDiscreteEnergiesPrev; ie<nEnergiesPrev; ie++) {
720  G4double ener = angParPrev->theAngular[ie].GetLabel();
721  G4double enerT = (ener-minEnerPrev)/(maxEnerPrev-minEnerPrev);
722  theEnergiesTransformed.insert(enerT);
723  //- if( getenv("G4PHPTEST2") ) G4cout << this << " G4ParticleHPContAngularPar::PrepareTableInterpolation theEnergiesTransformed2 " << enerT << G4endl; //GDEB
724  }
725  // add the maximum energy
726  theEnergiesTransformed.insert(1.);
727 
728 }
729 
731  G4ParticleHPContAngularPar & angpar1,
732  G4ParticleHPContAngularPar & angpar2)
733 {
734  G4int ie,ie1,ie2, ie1Prev, ie2Prev;
735  nAngularParameters = angpar1.nAngularParameters;
736  theManager = angpar1.theManager;
737  theEnergy = anEnergy;
738 
739  nDiscreteEnergies = theDiscreteEnergies.size();
740  std::set<G4double>::const_iterator itede;
741  std::map<G4double,G4int> discEnerOwn1 = angpar1.GetDiscreteEnergiesOwn();
742  std::map<G4double,G4int> discEnerOwn2 = angpar2.GetDiscreteEnergiesOwn();
743  std::map<G4double,G4int>::const_iterator itedeo;
744  ie = 0;
745  for( itede = theDiscreteEnergies.begin(); itede != theDiscreteEnergies.end(); itede++, ie++ ) {
746  G4double discEner = *itede;
747  itedeo = discEnerOwn1.find(discEner);
748  if( itedeo == discEnerOwn1.end() ) {
749  ie1 = 0;
750  } else {
751  ie1 = -1;
752  }
753  itedeo = discEnerOwn2.find(discEner);
754  if( itedeo == discEnerOwn2.end() ) {
755  ie2 = 0;
756  } else {
757  ie2 = -1;
758  }
759 
760  theAngular[ie].SetLabel(discEner);
761  G4double val1, val2;
762  for(G4int ip=0; ip<nAngularParameters; ip++) {
763  if( ie1 != -1 ) {
764  val1 = angpar1.theAngular[ie1].GetValue(ip);
765  } else {
766  val1 = 0.;
767  }
768  if( ie2 != -1 ) {
769  val2 = angpar2.theAngular[ie2].GetValue(ip);
770  } else {
771  val2 = 0.;
772  }
773 
774  G4double value = theInt.Interpolate(aScheme, anEnergy,
775  angpar1.theEnergy, angpar2.theEnergy,
776  val1,
777  val2);
778  if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge DiscreteEnergies val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
779 
780  theAngular[ie].SetValue(ip, value);
781  }
782  }
783 
784  if(theAngular != 0) delete [] theAngular;
785  nEnergies = nDiscreteEnergies + angpar2.GetNEnergiesTransformed();
786  theAngular = new G4ParticleHPList [nEnergies];
787 
788  //---- Get minimum and maximum energy interpolating
789  theMinEner = angpar1.GetMinEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMinEner()-angpar1.GetMinEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());
790  theMaxEner = angpar1.GetMaxEner() + (theEnergy-angpar1.GetEnergy()) * (angpar2.GetMaxEner()-angpar1.GetMaxEner())/(angpar2.GetEnergy()-angpar1.GetEnergy());
791 
792  if( getenv("G4PHPTEST2") ) G4cout << " G4ParticleHPContAngularPar::Merge E " << anEnergy << " minmax " << theMinEner << " " << theMaxEner << G4endl; //GDEB
793 
794  //--- Loop to energies of new set
795  std::set<G4double> energiesTransformed = angpar2.GetEnergiesTransformed();
796  std::set<G4double>::const_iterator iteet = energiesTransformed.begin();
797  G4int nEnergies1 = angpar1.GetNEnergies();
798  G4int nDiscreteEnergies1 = angpar1.GetNDiscreteEnergies();
799  G4double minEner1 = angpar1.GetMinEner();
800  G4double maxEner1 = angpar1.GetMaxEner();
801  G4int nEnergies2 = angpar2.GetNEnergies();
802  G4int nDiscreteEnergies2 = angpar2.GetNDiscreteEnergies();
803  G4double minEner2 = angpar2.GetMinEner();
804  G4double maxEner2 = angpar2.GetMaxEner();
805  for(ie=nDiscreteEnergies; ie<nEnergies; ie++,iteet++) {
806  G4double eT = (*iteet);
807 
808  //--- Use eT1 = eT: Get energy and parameters of angpar1 for this eT
809  G4double e1 = (maxEner1-minEner1) * eT + minEner1;
810  //----- Get parameter value corresponding to this e1
811  for(ie1=nDiscreteEnergies1; ie1<nEnergies1; ie1++) {
812  if( (angpar1.theAngular[ie1].GetLabel() - e1) > 1.E-10*e1 ) break;
813  }
814  ie1Prev = ie1 - 1;
815  if( ie1 == 0 ) ie1Prev++;
816  if( ie1 == nEnergies1 ) {
817  ie1--;
818  ie1Prev = ie1;
819  }
820  //--- Use eT2 = eT: Get energy and parameters of angpar2 for this eT
821  G4double e2 = (maxEner2-minEner2) * eT + minEner2;
822  //----- Get parameter value corresponding to this e2
823  for(ie2=nDiscreteEnergies2; ie2<nEnergies2; ie2++) {
824  // G4cout << " GET IE2 " << ie2 << " - " << angpar2.theAngular[ie2].GetLabel() - e2 << " " << angpar2.theAngular[ie2].GetLabel() << " " << e2 <<G4endl;
825  if( (angpar2.theAngular[ie2].GetLabel() - e2) > 1.E-10*e2 ) break;
826  }
827  ie2Prev = ie2 - 1;
828  if( ie2 == 0 ) ie2Prev++;
829  if( ie2 == nEnergies2 ) {
830  ie2--;
831  ie2Prev = ie2;
832  }
833 
834  //---- Energy corresponding to energy transformed
835  G4double eN = (theMaxEner-theMinEner) * eT + theMinEner;
836  if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ie1 << " " << ie2 << " G4ParticleHPContAngularPar::loop eT " << eT << " -> eN " << eN << " e1 " << e1 << " e2 " << e2 << G4endl; //GDEB
837 
838  theAngular[ie].SetLabel(eN);
839 
840  for(G4int ip=0; ip<nAngularParameters; ip++) {
841  G4double val1 = theInt.Interpolate2(theManager.GetScheme(ie),
842  e1,
843  angpar1.theAngular[ie1Prev].GetLabel(),
844  angpar1.theAngular[ie1].GetLabel(),
845  angpar1.theAngular[ie1Prev].GetValue(ip),
846  angpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);
847  G4double val2 = theInt.Interpolate2(theManager.GetScheme(ie),
848  e2,
849  angpar2.theAngular[ie2Prev].GetLabel(),
850  angpar2.theAngular[ie2].GetLabel(),
851  angpar2.theAngular[ie2Prev].GetValue(ip),
852  angpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);
853 
854  G4double value = theInt.Interpolate(aScheme, anEnergy,
855  angpar1.theEnergy, angpar2.theEnergy,
856  val1,
857  val2);
858  //value /= (theMaxEner-theMinEner);
859  if ( theMaxEner != theMinEner ) {
860  value /= (theMaxEner-theMinEner);
861  } else if ( value != 0 ) {
862  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPContAngularPar::PrepareTableInterpolation theMaxEner == theMinEner and value != 0.");
863  }
864  if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::Merge val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
865  //- val1 = angpar1.theAngular[ie1-1].GetValue(ip) * (maxEner1-minEner1);
866  //- val2 = angpar2.theAngular[ie2-1].GetValue(ip) * (maxEner2-minEner2);
867  //- if( getenv("G4PHPTEST2") ) G4cout << ie << " " << ip << " G4ParticleHPContAngularPar::MergeOLD val1 " << val1 << " val2 " << val2 << " value " << value << G4endl; //GDEB
868 
869  theAngular[ie].SetValue(ip, value);
870  }
871  }
872 
873  if( getenv("G4PHPTEST2") ) {
874  G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR1 " << G4endl; //GDEB
875  angpar1.Dump();
876  G4cout << " G4ParticleHPContAngularPar::Merge ANGPAR2 " << G4endl;
877  angpar2.Dump();
878  G4cout << " G4ParticleHPContAngularPar::Merge ANGPARNEW " << G4endl;
879  Dump();
880  }
881 }
882 
884 {
885  G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies << " " << nAngularParameters << G4endl;
886 
887  for(G4int ii=0; ii<nEnergies; ii++) {
888  theAngular[ii].Dump();
889  }
890 
891 }
G4double G4ParticleHPJENDLHEData::G4double result
std::map< G4double, G4int > GetDiscreteEnergiesOwn() const
void Init(G4int i, G4double e, G4int n)
G4int GetVectorLength() const
G4double GetTotalMomentum() const
void Init(std::istream &aDataFile, G4int nPar, G4double unit=1.)
G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass, G4int angularRep, G4int interpol)
void BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme, G4ParticleHPContAngularPar &store1, G4ParticleHPContAngularPar &store2)
G4double Interpolate2(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
void Init(G4int aScheme, G4int aRange)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4double SampleMax(G4double energy)
value_type & Get() const
Definition: G4Cache.hh:282
void SetInterpolationManager(const G4InterpolationManager &aManager)
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static constexpr double twopi
Definition: G4SIunits.hh:76
std::set< G4double > GetEnergiesTransformed() const
void SetY(G4int i, G4double x)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
void Init(std::istream &aDataFile, G4ParticleDefinition *projectile)
const XML_Char int const XML_Char * value
Definition: expat.h:331
void SetValue(G4int i, G4double y)
G4InterpolationScheme GetScheme(G4int index) const
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static constexpr double eV
Definition: G4SIunits.hh:215
G4double GetX(G4int i) const
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetY(G4double x)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
void SetLabel(G4double aLabel)
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
void SetManager(G4InterpolationManager &aManager)
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4InterpolationScheme
G4double GetPDGMass() const
static G4double GetBindingEnergy(const G4int A, const G4int Z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetCoeff(G4int i, G4int l, G4double coeff)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4InterpolationScheme GetInverseScheme(G4int index)
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4double GetValue(G4int i)
#define G4endl
Definition: G4ios.hh:61
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
double G4double
Definition: G4Types.hh:76
void PrepareTableInterpolation(const G4ParticleHPContAngularPar *angularPrev)
float amu_c2
Definition: hepunit.py:277
static G4He3 * He3()
Definition: G4He3.cc:94
void SetX(G4int i, G4double e)
#define DBL_MAX
Definition: templates.hh:83
G4double GetMass() const