Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ParticleHPDiscreteTwoBody.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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 //080612 Bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2,3
31 //080709 Bug fix Sampling Legendre expansion by T. Koi
32 //101110 Bug fix in MF=6, LAW=2 case; contribution from E. Mendoza, D. Cano-Ott (CIEMAT)
33 //
34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
35 //
37 #include "G4Gamma.hh"
38 #include "G4Electron.hh"
39 #include "G4Positron.hh"
40 #include "G4Neutron.hh"
41 #include "G4Proton.hh"
42 #include "G4Deuteron.hh"
43 #include "G4Triton.hh"
44 #include "G4He3.hh"
45 #include "G4Alpha.hh"
46 #include "G4ParticleHPVector.hh"
48 
50 { // Interpolation still only for the most used parts; rest to be Done @@@@@
52  G4int Z = static_cast<G4int>(massCode/1000);
53  G4int A = static_cast<G4int>(massCode-1000*Z);
54 
55  if(massCode==0)
56  {
57  result->SetDefinition(G4Gamma::Gamma());
58  }
59  else if(A==0)
60  {
62  if(Z==1) result->SetDefinition(G4Positron::Positron());
63  }
64  else if(A==1)
65  {
67  if(Z==1) result->SetDefinition(G4Proton::Proton());
68  }
69  else if(A==2)
70  {
72  }
73  else if(A==3)
74  {
75  result->SetDefinition(G4Triton::Triton());
76  if(Z==2) result->SetDefinition(G4He3::He3());
77  }
78  else if(A==4)
79  {
80  result->SetDefinition(G4Alpha::Alpha());
81  if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
82  }
83  else
84  {
85  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPDiscreteTwoBody: Unknown ion case 2");
86  }
87 
88 // get cosine(theta)
89  G4int i(0), it(0);
90  G4double cosTh(0);
91  for(i=0; i<nEnergy; i++)
92  {
93  it = i;
94  if(theCoeff[i].GetEnergy()>anEnergy) break;
95  }
96  if(it==0||it==nEnergy-1)
97  {
98  if(theCoeff[it].GetRepresentation()==0)
99  {
100 //TK Legendre expansion
101  G4ParticleHPLegendreStore theStore(1);
102  theStore.SetCoeff(0, theCoeff);
103  theStore.SetManager(theManager);
104  //cosTh = theStore.SampleMax(anEnergy);
105  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
106  cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
107  }
108  else if(theCoeff[it].GetRepresentation()==12) // means LINLIN
109  {
110  G4ParticleHPVector theStore;
111  G4InterpolationManager aManager;
112  aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
113  theStore.SetInterpolationManager(aManager);
114  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
115  {
116  //101110
117  //theStore.SetX(i, theCoeff[it].GetCoeff(i));
118  //theStore.SetY(i, theCoeff[it].GetCoeff(i));
119  theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
120  theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
121  }
122  cosTh = theStore.Sample();
123  }
124  else if(theCoeff[it].GetRepresentation()==14) //this is LOGLIN
125  {
126  G4ParticleHPVector theStore;
127  G4InterpolationManager aManager;
128  aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
129  theStore.SetInterpolationManager(aManager);
130  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
131  {
132  //101110
133  //theStore.SetX(i, theCoeff[it].GetCoeff(i));
134  //theStore.SetY(i, theCoeff[it].GetCoeff(i));
135  theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
136  theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
137  }
138  cosTh = theStore.Sample();
139  }
140  else
141  {
142  throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering");
143  }
144  }
145  else
146  {
147  if(!bCheckDiffCoeffRepr || theCoeff[it].GetRepresentation() == theCoeff[it-1].GetRepresentation())
148  {
149  if(theCoeff[it].GetRepresentation()==0)
150  {
151 //TK Legendre expansion
152  G4ParticleHPLegendreStore theStore(2);
153  theStore.SetCoeff(0, &(theCoeff[it-1]));
154  theStore.SetCoeff(1, &(theCoeff[it]));
155  G4InterpolationManager aManager;
156  aManager.Init(theManager.GetScheme(it), 2);
157  theStore.SetManager(aManager);
158  //cosTh = theStore.SampleMax(anEnergy);
159 //080709 TKDB
160  cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
161  }
162  else if(theCoeff[it].GetRepresentation()==12) // LINLIN
163  {
164  G4ParticleHPVector theBuff1;
165  G4InterpolationManager aManager1;
166  aManager1.Init(LINLIN, theCoeff[it-1].GetNumberOfPoly()/2);
167  theBuff1.SetInterpolationManager(aManager1);
168  for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i+=2)
169  {
170  //101110
171  //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
172  //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
173  theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
174  theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
175  }
176  G4ParticleHPVector theBuff2;
177  G4InterpolationManager aManager2;
178  aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly()/2);
179  theBuff2.SetInterpolationManager(aManager2);
180  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
181  {
182  //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
183  //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
184  theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
185  theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
186  }
187 
188  G4double x1 = theCoeff[it-1].GetEnergy();
189  G4double x2 = theCoeff[it].GetEnergy();
190  G4double x = anEnergy;
191  G4double y1, y2, y, mu;
192 
193  G4ParticleHPVector theStore1;
194  theStore1.SetInterpolationManager(aManager1);
195  G4ParticleHPVector theStore2;
196  theStore2.SetInterpolationManager(aManager2);
197  G4ParticleHPVector theStore;
198 
199  // for fixed mu get p1, p2 and interpolate according to x
200  for(i=0; i<theBuff1.GetVectorLength(); i++)
201  {
202  mu = theBuff1.GetX(i);
203  y1 = theBuff1.GetY(i);
204  y2 = theBuff2.GetY(mu);
206  theStore1.SetData(i, mu, y);
207  }
208  for(i=0; i<theBuff2.GetVectorLength(); i++)
209  {
210  mu = theBuff2.GetX(i);
211  y1 = theBuff2.GetY(i);
212  y2 = theBuff1.GetY(mu);
214  theStore2.SetData(i, mu, y);
215  }
216  theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
217  cosTh = theStore.Sample();
218  }
219  else if(theCoeff[it].GetRepresentation()==14) //TK LOG_LIN
220  {
221  G4ParticleHPVector theBuff1;
222  G4InterpolationManager aManager1;
223  aManager1.Init(LOGLIN, theCoeff[it-1].GetNumberOfPoly()/2);
224  theBuff1.SetInterpolationManager(aManager1);
225  for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i+=2)
226  {
227  //101110
228  //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
229  //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
230  theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
231  theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
232  }
233 
234  G4ParticleHPVector theBuff2;
235  G4InterpolationManager aManager2;
236  aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly()/2);
237  theBuff2.SetInterpolationManager(aManager2);
238  for(i=0;i<theCoeff[it].GetNumberOfPoly(); i+=2)
239  {
240  //101110
241  //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
242  //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
243  theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
244  theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
245  }
246 
247  G4double x1 = theCoeff[it-1].GetEnergy();
248  G4double x2 = theCoeff[it].GetEnergy();
249  G4double x = anEnergy;
250  G4double y1, y2, y, mu;
251 
252  G4ParticleHPVector theStore1;
253  theStore1.SetInterpolationManager(aManager1);
254  G4ParticleHPVector theStore2;
255  theStore2.SetInterpolationManager(aManager2);
256  G4ParticleHPVector theStore;
257 
258  // for fixed mu get p1, p2 and interpolate according to x
259  for(i=0; i<theBuff1.GetVectorLength(); i++)
260  {
261  mu = theBuff1.GetX(i);
262  y1 = theBuff1.GetY(i);
263  y2 = theBuff2.GetY(mu);
265  theStore1.SetData(i, mu, y);
266  }
267  for(i=0; i<theBuff2.GetVectorLength(); i++)
268  {
269  mu = theBuff2.GetX(i);
270  y1 = theBuff2.GetY(i);
271  y2 = theBuff1.GetY(mu);
273  theStore2.SetData(i, mu, y);
274  }
275  theStore.Merge(&theStore1, &theStore2);
276  cosTh = theStore.Sample();
277  }
278  else
279  {
280  throw G4HadronicException(__FILE__, __LINE__, "Two neighbouring distributions with different interpolation");
281  }
282  }
283  else
284  {
285  G4cout << " theCoeff[it].GetRepresent MEM " << it << " " << &theCoeff[it] << " " << &theCoeff[it-1] << G4endl;
286  G4cout << " theCoeff[it].GetRepresent " << it << " " << theCoeff[it].GetRepresentation() << " != " << theCoeff[it-1].GetRepresentation() << G4endl;
287 
288  throw G4HadronicException(__FILE__, __LINE__, "unknown representation type in Two-body scattering, case 2");
289  }
290  }
291 
292 // now get the energy from kinematics and Q-value.
293 
294  //G4double restEnergy = anEnergy+GetQValue();
295 
296 // assumed to be in CMS @@@@@@@@@@@@@@@@@
297 
298  //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
299  //G4double residualMass = GetTarget()->GetMass() + GetNeutron()->GetMass()
300  // - result->GetMass() - GetQValue();
301  //G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
303  G4double A1prim = result->GetMass()/GetProjectileRP()->GetMass();
304  //G4double E1 = (A1+1)*(A1+1)/A1/A1*anEnergy;
305  //Bug fix Bugzilla #1815
306  G4double E1 = anEnergy;
307  G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue());
308 
309  result->SetKineticEnergy(kinE); // non relativistic @@
311  G4double theta = std::acos(cosTh);
312  G4double sinth = std::sin(theta);
313  G4double mtot = result->GetTotalMomentum();
314  G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
315  result->SetMomentum(tempVector);
316 
317 // some garbage collection
318 
319 // return the result
320  return result;
321 }
Float_t x
Definition: compare.C:6
void SetInterpolationManager(const G4InterpolationManager &aManager)
static G4He3 * He3()
Definition: G4He3.cc:94
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
G4ReactionProduct * Sample(G4double anEnergy, G4double massCode, G4double mass)
void SetCoeff(G4int i, G4int l, G4double coeff)
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 SetManager(G4InterpolationManager &aManager)
G4double SampleDiscreteTwoBody(G4double anEnergy)
#define G4UniformRand()
Definition: Randomize.hh:53
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
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)
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4double GetMass() const
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
G4GLOB_DLL std::ostream G4cout
G4int GetVectorLength() const
Float_t x2[n_points_geant4]
Definition: compare.C:26
G4ParticleHPLegendreTable * theCoeff
void Merge(G4ParticleHPVector *active, G4ParticleHPVector *passive)