Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPDiscreteTwoBody Class Reference

#include <G4ParticleHPDiscreteTwoBody.hh>

Inheritance diagram for G4ParticleHPDiscreteTwoBody:
Collaboration diagram for G4ParticleHPDiscreteTwoBody:

Public Member Functions

 G4ParticleHPDiscreteTwoBody ()
 
 ~G4ParticleHPDiscreteTwoBody ()
 
void Init (std::istream &aDataFile)
 
G4ReactionProductSample (G4double anEnergy, G4double massCode, G4double mass)
 
G4double MeanEnergyOfThisInteraction ()
 
- Public Member Functions inherited from G4VParticleHPEnergyAngular
 G4VParticleHPEnergyAngular ()
 
virtual ~G4VParticleHPEnergyAngular ()
 
void SetProjectileRP (G4ReactionProduct *aIncidentParticleRP)
 
void SetTarget (G4ReactionProduct *aTarget)
 
G4ReactionProductGetTarget ()
 
G4ReactionProductGetProjectileRP ()
 
G4ReactionProductGetCMS ()
 
void SetQValue (G4double aValue)
 
virtual void ClearHistories ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4VParticleHPEnergyAngular
G4double GetQValue ()
 

Detailed Description

Definition at line 42 of file G4ParticleHPDiscreteTwoBody.hh.

Constructor & Destructor Documentation

G4ParticleHPDiscreteTwoBody::G4ParticleHPDiscreteTwoBody ( )
inline

Definition at line 46 of file G4ParticleHPDiscreteTwoBody.hh.

47  {
48  theCoeff = 0;
49  bCheckDiffCoeffRepr = true;
50  if ( getenv( "G4PHP_DO_NOT_CHECK_DIFF_COEFF_REPR" ) ) bCheckDiffCoeffRepr = false;
51  nEnergy = 0;
52  }
G4ParticleHPDiscreteTwoBody::~G4ParticleHPDiscreteTwoBody ( )
inline

Definition at line 53 of file G4ParticleHPDiscreteTwoBody.hh.

54  {
55  if(theCoeff!=0) delete [] theCoeff;
56  }

Member Function Documentation

void G4ParticleHPDiscreteTwoBody::Init ( std::istream &  aDataFile)
inlinevirtual

Implements G4VParticleHPEnergyAngular.

Definition at line 58 of file G4ParticleHPDiscreteTwoBody.hh.

59  {
60  aDataFile >> nEnergy;
61  theManager.Init(aDataFile);
62  theCoeff = new G4ParticleHPLegendreTable[nEnergy];
63  for(G4int i=0; i<nEnergy; i++)
64  {
66  G4int aRep, nCoeff;
67  aDataFile >> energy >> aRep >> nCoeff;
68  //- G4cout << this << " " << i << " G4ParticleHPDiscreteTwoBody READ DATA " << energy << " " << aRep << " " << nCoeff << G4endl;
69  energy*=CLHEP::eV;
70  G4int nPoints=nCoeff;
71  if(aRep>0) nPoints*=2;
72  //theCoeff[i].Init(energy, nPoints);
73 
74  theCoeff[i].Init(energy, nPoints-1);
75  theCoeff[i].SetRepresentation(aRep);
76  for(G4int ii=0; ii<nPoints; ii++)
77  {
78  G4double y;
79  aDataFile >> y;
80  theCoeff[i].SetCoeff(ii, y);
81  }
82  }
83  }
void Init(G4int aScheme, G4int aRange)
int G4int
Definition: G4Types.hh:78
void SetCoeff(G4int l, G4double coeff)
static constexpr double eV
G4double energy(const ThreeVector &p, const G4double m)
double G4double
Definition: G4Types.hh:76
void Init(std::istream &aDataFile)

Here is the call graph for this function:

G4double G4ParticleHPDiscreteTwoBody::MeanEnergyOfThisInteraction ( )
inlinevirtual

Implements G4VParticleHPEnergyAngular.

Definition at line 86 of file G4ParticleHPDiscreteTwoBody.hh.

86 { return -1; }
G4ReactionProduct * G4ParticleHPDiscreteTwoBody::Sample ( G4double  anEnergy,
G4double  massCode,
G4double  mass 
)
virtual

Implements G4VParticleHPEnergyAngular.

Definition at line 49 of file G4ParticleHPDiscreteTwoBody.cc.

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);
205  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
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);
213  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
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);
264  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
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);
272  y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
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 }
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetVectorLength() const
G4double GetTotalMomentum() const
void Init(G4int aScheme, G4int aRange)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetData(G4int i, G4double x, G4double y)
void SetInterpolationManager(const G4InterpolationManager &aManager)
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetY(G4int i, G4double x)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4InterpolationScheme GetScheme(G4int index) const
static G4Triton * Triton()
Definition: G4Triton.cc:95
void Merge(G4ParticleHPVector *active, G4ParticleHPVector *passive)
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double GetX(G4int i) const
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetY(G4double x)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
static G4Positron * Positron()
Definition: G4Positron.cc:94
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
void SetX(G4int i, G4double e)
G4double GetMass() const
static constexpr double twopi
Definition: SystemOfUnits.h:55

Here is the call graph for this function:


The documentation for this class was generated from the following files: