Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
G4ParticleHPDiscreteTwoBody Class Reference

#include <G4ParticleHPDiscreteTwoBody.hh>

Inheritance diagram for G4ParticleHPDiscreteTwoBody:
Inheritance graph
[legend]
Collaboration diagram for G4ParticleHPDiscreteTwoBody:
Collaboration graph
[legend]

Public Member Functions

 G4ParticleHPDiscreteTwoBody ()
 
 ~G4ParticleHPDiscreteTwoBody ()
 
void Init (std::istream &aDataFile)
 
G4ReactionProduct * Sample (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)
 
G4ReactionProduct * GetTarget ()
 
G4ReactionProduct * GetProjectileRP ()
 
G4ReactionProduct * GetCMS ()
 
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  {
65  G4double energy;
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  }
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition: G4InterpolationManager.hh:84
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleHPLegendreTable::SetCoeff
void SetCoeff(G4int l, G4double coeff)
Definition: G4ParticleHPLegendreTable.hh:95
CLHEP::eV
static constexpr double eV
Definition: SystemOfUnits.h:175
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleHPLegendreTable::Init
void Init(std::istream &aDataFile)
Definition: G4ParticleHPLegendreTable.hh:70

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 @@@@@
51  G4ReactionProduct * result = new G4ReactionProduct;
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  {
61  result->SetDefinition(G4Electron::Electron());
62  if(Z==1) result->SetDefinition(G4Positron::Positron());
63  }
64  else if(A==1)
65  {
66  result->SetDefinition(G4Neutron::Neutron());
67  if(Z==1) result->SetDefinition(G4Proton::Proton());
68  }
69  else if(A==2)
70  {
71  result->SetDefinition(G4Deuteron::Deuteron());
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 @@
302  G4double A1 = GetTarget()->GetMass()/GetProjectileRP()->GetMass();
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 @@
310  G4double phi = CLHEP::twopi*G4UniformRand();
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 }
result
G4double G4ParticleHPJENDLHEData::G4double result
Definition: G4ParticleHPJENDLHEData.cc:257
G4ParticleHPVector::GetVectorLength
G4int GetVectorLength() const
Definition: G4ParticleHPVector.hh:199
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition: G4InterpolationManager.hh:84
G4ReactionProduct::SetKineticEnergy
void SetKineticEnergy(const G4double en)
Definition: G4ReactionProduct.hh:132
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:171
G4ParticleHPVector::SetData
void SetData(G4int i, G4double x, G4double y)
Definition: G4ParticleHPVector.hh:96
LINLIN
Definition: G4InterpolationScheme.hh:39
G4ParticleHPVector::SetInterpolationManager
void SetInterpolationManager(const G4InterpolationManager &aManager)
Definition: G4ParticleHPVector.hh:483
test.x
tuple x
Definition: test.py:50
G4int
int G4int
Definition: G4Types.hh:78
G4ReactionProduct::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:160
G4ParticleHPVector::SetY
void SetY(G4int i, G4double x)
Definition: G4ParticleHPVector.hh:113
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
G4InterpolationManager::GetScheme
G4InterpolationScheme GetScheme(G4int index) const
Definition: G4InterpolationManager.hh:131
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4ParticleHPVector::Merge
void Merge(G4ParticleHPVector *active, G4ParticleHPVector *passive)
Definition: G4ParticleHPVector.hh:273
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4ParticleHPVector::GetX
G4double GetX(G4int i) const
Definition: G4ParticleHPVector.hh:127
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4ParticleHPVector::GetY
G4double GetY(G4double x)
Definition: G4ParticleHPVector.hh:198
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4ParticleHPInterpolator::Interpolate
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
Definition: G4ParticleHPInterpolator.hh:91
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double
double G4double
Definition: G4Types.hh:76
LOGLIN
Definition: G4InterpolationScheme.hh:39
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:94
G4ParticleHPVector::SetX
void SetX(G4int i, G4double e)
Definition: G4ParticleHPVector.hh:103
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
CLHEP::twopi
static constexpr double twopi
Definition: SystemOfUnits.h:55
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262

Here is the call graph for this function:

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