Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ParticleHPLabAngularEnergy.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 // 080808 Bug fix in serching mu bin and index for theBuff2b by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "Randomize.hh"
38 #include "G4Gamma.hh"
39 #include "G4Electron.hh"
40 #include "G4Positron.hh"
41 #include "G4Neutron.hh"
42 #include "G4Proton.hh"
43 #include "G4Deuteron.hh"
44 #include "G4Triton.hh"
45 #include "G4He3.hh"
46 #include "G4Alpha.hh"
47 
48 void G4ParticleHPLabAngularEnergy::Init(std::istream & aDataFile)
49 {
50  aDataFile >> nEnergies;
51  theManager.Init(aDataFile);
53  nCosTh = new G4int[nEnergies];
56  for(G4int i=0; i<nEnergies; i++)
57  {
58  aDataFile >> theEnergies[i];
59  theEnergies[i]*=eV;
60  aDataFile >> nCosTh[i];
61  theSecondManager[i].Init(aDataFile);
62  theData[i] = new G4ParticleHPVector[nCosTh[i]];
63  G4double label;
64  for(G4int ii=0; ii<nCosTh[i]; ii++)
65  {
66  aDataFile >> label;
67  theData[i][ii].SetLabel(label);
68  theData[i][ii].Init(aDataFile, eV);
69  }
70  }
71 }
72 
74 {
76  G4int Z = static_cast<G4int>(massCode/1000);
77  G4int A = static_cast<G4int>(massCode-1000*Z);
78 
79  if(massCode==0)
80  {
81  result->SetDefinition(G4Gamma::Gamma());
82  }
83  else if(A==0)
84  {
86  if(Z==1) result->SetDefinition(G4Positron::Positron());
87  }
88  else if(A==1)
89  {
91  if(Z==1) result->SetDefinition(G4Proton::Proton());
92  }
93  else if(A==2)
94  {
96  }
97  else if(A==3)
98  {
99  result->SetDefinition(G4Triton::Triton());
100  if(Z==2) result->SetDefinition(G4He3::He3());
101  }
102  else if(A==4)
103  {
104  result->SetDefinition(G4Alpha::Alpha());
105  if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
106  }
107  else
108  {
109  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPLabAngularEnergy: Unknown ion case 2");
110  }
111 
112  // get theta, E
113  G4double cosTh, secEnergy;
114  G4int i, it(0);
115  // find the energy bin
116  for(i=0; i<nEnergies; i++)
117  {
118  it = i;
119  if ( anEnergy < theEnergies[i] ) break;
120  }
121  //080808
122  //if ( it == 0 || it == nEnergies-1 ) // it marks the energy bin
123  if ( it == 0 ) // it marks the energy bin
124  {
125  G4cout << "080808 Something unexpected is happen in G4ParticleHPLabAngularEnergy " << G4endl;
126  // integrate the prob for each costh, and select theta.
127  G4double * running = new G4double [nCosTh[it]];
128  running[0]=0;
129  for(i=0;i<nCosTh[it]; i++)
130  {
131  if(i!=0) running[i] = running[i-1];
132  running[i]+=theData[it][i].GetIntegral(); // Does interpolated integral.
133  }
134  G4double random = running[nCosTh[it]-1]*G4UniformRand();
135  G4int ith(0);
136  for(i=0;i<nCosTh[it]; i++)
137  {
138  ith = i;
139  if(random<running[i]) break;
140  }
141  //080807
142  //if ( ith == 0 || ith == nCosTh[it]-1 ) //ith marks the angluar bin
143  if ( ith == 0 ) //ith marks the angluar bin
144  {
145  cosTh = theData[it][ith].GetLabel();
146  secEnergy = theData[it][ith].Sample();
147  currentMeanEnergy = theData[it][ith].GetMeanX();
148  }
149  else
150  {
151  //080808
152  //G4double x1 = theData[it][ith-1].GetIntegral();
153  //G4double x2 = theData[it][ith].GetIntegral();
154  G4double x1 = running [ ith-1 ];
155  G4double x2 = running [ ith ];
156  G4double x = random;
157  G4double y1 = theData[it][ith-1].GetLabel();
158  G4double y2 = theData[it][ith].GetLabel();
159  cosTh = theInt.Interpolate(theSecondManager[it].GetInverseScheme(ith),
160  x, x1, x2, y1, y2);
161  G4ParticleHPVector theBuff1;
162  theBuff1.SetInterpolationManager(theData[it][ith-1].GetInterpolationManager());
163  G4ParticleHPVector theBuff2;
164  theBuff2.SetInterpolationManager(theData[it][ith].GetInterpolationManager());
165  x1=y1;
166  x2=y2;
167  G4double y, mu;
168  for(i=0;i<theData[it][ith-1].GetVectorLength(); i++)
169  {
170  mu = theData[it][ith-1].GetX(i);
171  y1 = theData[it][ith-1].GetY(i);
172  y2 = theData[it][ith].GetY(mu);
173 
174  y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
175  cosTh, x1,x2,y1,y2);
176  theBuff1.SetData(i, mu, y);
177  }
178  for(i=0;i<theData[it][ith].GetVectorLength(); i++)
179  {
180  mu = theData[it][ith].GetX(i);
181  y1 = theData[it][ith-1].GetY(mu);
182  y2 = theData[it][ith].GetY(i);
183  y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
184  cosTh, x1,x2,y1,y2);
185  theBuff2.SetData(i, mu, y);
186  }
187  G4ParticleHPVector theStore;
188  theStore.Merge(&theBuff1, &theBuff2);
189  secEnergy = theStore.Sample();
190  currentMeanEnergy = theStore.GetMeanX();
191  }
192  delete [] running;
193  }
194  else // this is the small big else.
195  {
196  G4double x, x1, x2, y1, y2, y, tmp, E;
197  // integrate the prob for each costh, and select theta.
198  G4ParticleHPVector run1;
199  run1.SetY(0, 0.);
200  for(i=0;i<nCosTh[it-1]; i++)
201  {
202  if(i!=0) run1.SetY(i, run1.GetY(i-1));
203  run1.SetX(i, theData[it-1][i].GetLabel());
204  run1.SetY(i, run1.GetY(i)+theData[it-1][i].GetIntegral());
205  }
206  G4ParticleHPVector run2;
207  run2.SetY(0, 0.);
208  for(i=0;i<nCosTh[it]; i++)
209  {
210  if(i!=0) run2.SetY(i, run2.GetY(i-1));
211  run2.SetX(i, theData[it][i].GetLabel());
212  run2.SetY(i, run2.GetY(i)+theData[it][i].GetIntegral());
213  }
214  // get the distributions for the correct neutron energy
215  x = anEnergy;
216  x1 = theEnergies[it-1];
217  x2 = theEnergies[it];
218  G4ParticleHPVector thBuff1; // to be interpolated as run1.
220  for(i=0; i<run1.GetVectorLength(); i++)
221  {
222  tmp = run1.GetX(i); //theta
223  y1 = run1.GetY(i); // integral
224  y2 = run2.GetY(tmp);
226  thBuff1.SetData(i, tmp, y);
227  }
228  G4ParticleHPVector thBuff2;
230  for(i=0; i<run2.GetVectorLength(); i++)
231  {
232  tmp = run2.GetX(i); //theta
233  y1 = run1.GetY(tmp); // integral
234  y2 = run2.GetY(i);
235  y = theInt.Lin(x, x1,x2,y1,y2);
236  thBuff2.SetData(i, tmp, y);
237  }
238  G4ParticleHPVector theThVec;
239  theThVec.Merge(&thBuff1 ,&thBuff2); // takes care of interpolation
240  G4double random = (theThVec.GetY(theThVec.GetVectorLength()-1)
241  -theThVec.GetY(0)) *G4UniformRand();
242  G4int ith(0);
243  for(i=1;i<theThVec.GetVectorLength(); i++)
244  {
245  ith = i;
246  if(random<theThVec.GetY(i)-theThVec.GetY(0)) break;
247  }
248  {
249  // calculate theta
250  G4double xx, xx1, xx2, yy1, yy2;
251  xx = random;
252  xx1 = theThVec.GetY(ith-1)-theThVec.GetY(0); // integrals
253  xx2 = theThVec.GetY(ith)-theThVec.GetY(0);
254  yy1 = theThVec.GetX(ith-1); // std::cos(theta)
255  yy2 = theThVec.GetX(ith);
256  cosTh = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
257  xx, xx1,xx2,yy1,yy2);
258  }
259  G4int i1(0), i2(0);
260  // get the indixes of the vectors close to theta for low energy
261  // first it-1 !!!! i.e. low in energy
262  for(i=0; i<nCosTh[it-1]; i++)
263  {
264  i1 = i;
265  if(cosTh<theData[it-1][i].GetLabel()) break;
266  }
267  // now get the prob at this energy for the right theta value
268  x = cosTh;
269  x1 = theData[it-1][i1-1].GetLabel();
270  x2 = theData[it-1][i1].GetLabel();
271  G4ParticleHPVector theBuff1a;
272  theBuff1a.SetInterpolationManager(theData[it-1][i1-1].GetInterpolationManager());
273  for(i=0;i<theData[it-1][i1-1].GetVectorLength(); i++)
274  {
275  E = theData[it-1][i1-1].GetX(i);
276  y1 = theData[it-1][i1-1].GetY(i);
277  y2 = theData[it-1][i1].GetY(E);
278  y = theInt.Lin(x, x1,x2,y1,y2);
279  theBuff1a.SetData(i, E, y); // wrong E, right theta.
280  }
281  G4ParticleHPVector theBuff2a;
282  theBuff2a.SetInterpolationManager(theData[it-1][i1].GetInterpolationManager());
283  for(i=0;i<theData[it-1][i1].GetVectorLength(); i++)
284  {
285  E = theData[it-1][i1].GetX(i);
286  y1 = theData[it-1][i1-1].GetY(E);
287  y2 = theData[it-1][i1].GetY(i);
288  y = theInt.Lin(x, x1,x2,y1,y2);
289  theBuff2a.SetData(i, E, y); // wrong E, right theta.
290  }
291  G4ParticleHPVector theStore1;
292  theStore1.Merge(&theBuff1a, &theBuff2a); // wrong E, right theta, complete binning
293 
294  // get the indixes of the vectors close to theta for high energy
295  // then it !!!! i.e. high in energy
296  for(i=0; i<nCosTh[it]; i++)
297  {
298  i2 = i;
299  if(cosTh<theData[it][i2].GetLabel()) break;
300  } // sonderfaelle mit i1 oder i2 head on fehlen. @@@@@
301  x1 = theData[it][i2-1].GetLabel();
302  x2 = theData[it][i2].GetLabel();
303  G4ParticleHPVector theBuff1b;
304  theBuff1b.SetInterpolationManager(theData[it][i2-1].GetInterpolationManager());
305  for(i=0;i<theData[it][i2-1].GetVectorLength(); i++)
306  {
307  E = theData[it][i2-1].GetX(i);
308  y1 = theData[it][i2-1].GetY(i);
309  y2 = theData[it][i2].GetY(E);
310  y = theInt.Lin(x, x1,x2,y1,y2);
311  theBuff1b.SetData(i, E, y); // wrong E, right theta.
312  }
313  G4ParticleHPVector theBuff2b;
314  theBuff2b.SetInterpolationManager(theData[it][i2].GetInterpolationManager());
315  //080808 i1 -> i2
316  //for(i=0;i<theData[it][i1].GetVectorLength(); i++)
317  for(i=0;i<theData[it][i2].GetVectorLength(); i++)
318  {
319  //E = theData[it][i1].GetX(i);
320  //y1 = theData[it][i1-1].GetY(E);
321  //y2 = theData[it][i1].GetY(i);
322  E = theData[it][i2].GetX(i);
323  y1 = theData[it][i2-1].GetY(E);
324  y2 = theData[it][i2].GetY(i);
325  y = theInt.Lin(x, x1,x2,y1,y2);
326  theBuff2b.SetData(i, E, y); // wrong E, right theta.
327  }
328  G4ParticleHPVector theStore2;
329  theStore2.Merge(&theBuff1b, &theBuff2b); // wrong E, right theta, complete binning
330  // now get to the right energy.
331 
332  x = anEnergy;
333  x1 = theEnergies[it-1];
334  x2 = theEnergies[it];
335  G4ParticleHPVector theOne1;
337  for(i=0; i<theStore1.GetVectorLength(); i++)
338  {
339  E = theStore1.GetX(i);
340  y1 = theStore1.GetY(i);
341  y2 = theStore2.GetY(E);
343  theOne1.SetData(i, E, y); // both correct
344  }
345  G4ParticleHPVector theOne2;
347  for(i=0; i<theStore2.GetVectorLength(); i++)
348  {
349  E = theStore2.GetX(i);
350  y1 = theStore1.GetY(E);
351  y2 = theStore2.GetY(i);
353  theOne2.SetData(i, E, y); // both correct
354  }
355  G4ParticleHPVector theOne;
356  theOne.Merge(&theOne1, &theOne2); // both correct, complete binning
357 
358  secEnergy = theOne.Sample();
359  currentMeanEnergy = theOne.GetMeanX();
360  }
361 
362 // now do random direction in phi, and fill the result.
363 
364  result->SetKineticEnergy(secEnergy);
365 
366  G4double phi = twopi*G4UniformRand();
367  G4double theta = std::acos(cosTh);
368  G4double sinth = std::sin(theta);
369  G4double mtot = result->GetTotalMomentum();
370  G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
371  result->SetMomentum(tempVector);
372 
373  return result;
374 }
Float_t x
Definition: compare.C:6
void SetInterpolationManager(const G4InterpolationManager &aManager)
static G4He3 * He3()
Definition: G4He3.cc:94
Double_t xx
Float_t y1[n_points_granero]
Definition: compare.C:5
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 Init(std::istream &aDataFile)
Float_t tmp
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4InterpolationScheme GetScheme(G4int index) const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
Float_t y2[n_points_geant4]
Definition: compare.C:26
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
Float_t Z
double G4double
Definition: G4Types.hh:76
void SetKineticEnergy(const G4double en)
G4double GetY(G4double x)
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
const G4InterpolationManager & GetInterpolationManager() const
G4double GetTotalMomentum() const
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 SetData(G4int i, G4double x, G4double y)
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4double G4ParticleHPJENDLHEData::G4double result
void SetX(G4int i, G4double e)
G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass)
void SetLabel(G4double aLabel)
int G4int
Definition: G4Types.hh:78
void SetY(G4int i, G4double x)
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
G4GLOB_DLL std::ostream G4cout
G4int GetVectorLength() const
Float_t x2[n_points_geant4]
Definition: compare.C:26
void Merge(G4ParticleHPVector *active, G4ParticleHPVector *passive)