Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ParticleHPContAngularPar.cc
이 파일의 문서화 페이지로 가기
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 
71  theProjectile = projectile;
72 
73  theEnergy = 0.0;
74  nEnergies = 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 
87  /*if( getenv("G4PHPTEST") )*/
88  theEnergy *= eV;
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 
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));
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();
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();
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));
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));
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)
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();
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);
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++) {
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;
736  theManager = angpar1.theManager;
737  theEnergy = anEnergy;
738 
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;
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
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++) {
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);
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 }
Float_t x
Definition: compare.C:6
void SetInterpolationManager(const G4InterpolationManager &aManager)
static G4He3 * He3()
Definition: G4He3.cc:94
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4InterpolationScheme
std::map< G4double, G4int > theDiscreteEnergiesOwn
Float_t y1[n_points_granero]
Definition: compare.C:5
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double GetX(G4int i) const
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetMomentum(const G4double x, const G4double y, const G4double z)
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
void SetCoeff(G4int i, G4int l, G4double coeff)
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double GetPDGMass() const
G4InterpolationScheme GetScheme(G4int index) const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4double Interpolate2(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
Float_t y2[n_points_geant4]
Definition: compare.C:26
G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass, G4int angularRep, G4int interpol)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
Float_t Z
void Init(std::istream &aDataFile, G4ParticleDefinition *projectile)
static G4double GetBindingEnergy(const G4int A, const G4int Z)
double G4double
Definition: G4Types.hh:76
static constexpr double amu_c2
void SetKineticEnergy(const G4double en)
G4double GetY(G4double x)
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
const XML_Char int const XML_Char * value
Definition: expat.h:331
void SetManager(G4InterpolationManager &aManager)
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double GetTotalMomentum() const
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
double A(double temperature)
void Init(G4int aScheme, G4int aRange)
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static G4Positron * Positron()
Definition: G4Positron.cc:94
static constexpr double eV
Definition: G4SIunits.hh:215
void Init(std::istream &aDataFile, G4int nPar, G4double unit=1.)
static G4Electron * Electron()
Definition: G4Electron.cc:94
void BuildByInterpolation(G4double anEnergy, G4InterpolationScheme aScheme, G4ParticleHPContAngularPar &store1, G4ParticleHPContAngularPar &store2)
G4double G4ParticleHPJENDLHEData::G4double result
void SetX(G4int i, G4double e)
G4InterpolationScheme GetInverseScheme(G4int index)
G4double GetMass() const
G4double GetValue(G4int i)
int G4int
Definition: G4Types.hh:78
void SetValue(G4int i, G4double y)
void SetY(G4int i, G4double x)
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double SampleMax(G4double energy)
G4GLOB_DLL std::ostream G4cout
G4int GetVectorLength() const
void PrepareTableInterpolation(const G4ParticleHPContAngularPar *angularPrev)
void SetLabel(G4double aLabel)
Float_t x2[n_points_geant4]
Definition: compare.C:26
std::set< G4double > GetEnergiesTransformed() const
#define DBL_MAX
Definition: templates.hh:83
std::map< G4double, G4int > GetDiscreteEnergiesOwn() const
void Init(G4int i, G4double e, G4int n)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments