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

#include <G4ParticleHPElasticFS.hh>

Inheritance diagram for G4ParticleHPElasticFS:
Collaboration diagram for G4ParticleHPElasticFS:

Public Member Functions

 G4ParticleHPElasticFS ()
 
 ~G4ParticleHPElasticFS ()
 
void Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)
 
G4HadFinalStateApplyYourself (const G4HadProjectile &theTrack)
 
G4ParticleHPFinalStateNew ()
 
- Public Member Functions inherited from G4ParticleHPFinalState
 G4ParticleHPFinalState ()
 
virtual ~G4ParticleHPFinalState ()
 
void Init (G4double A, G4double Z, G4String &dirName, G4String &aFSType, G4ParticleDefinition *projectile)
 
G4bool HasXsec ()
 
G4bool HasFSData ()
 
G4bool HasAnyData ()
 
virtual G4double GetXsec (G4double)
 
virtual G4ParticleHPVectorGetXsec ()
 
void SetA_Z (G4double anA, G4double aZ, G4int aM=0)
 
G4double GetZ ()
 
G4double GetN ()
 
G4double GetA ()
 
G4int GetM ()
 
void SetAZMs (G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
 
void SetProjectile (G4ParticleDefinition *projectile)
 

Additional Inherited Members

- Protected Member Functions inherited from G4ParticleHPFinalState
void adjust_final_state (G4LorentzVector)
 
G4bool DoNotAdjustFinalState ()
 
- Protected Attributes inherited from G4ParticleHPFinalState
G4bool hasXsec
 
G4bool hasFSData
 
G4bool hasAnyData
 
G4ParticleHPNames theNames
 
G4Cache< G4HadFinalState * > theResult
 
G4double theBaseA
 
G4double theBaseZ
 
G4int theBaseM
 
G4int theNDLDataZ
 
G4int theNDLDataA
 
G4int theNDLDataM
 
G4ParticleDefinitiontheProjectile
 

Detailed Description

Definition at line 46 of file G4ParticleHPElasticFS.hh.

Constructor & Destructor Documentation

G4ParticleHPElasticFS::G4ParticleHPElasticFS ( )
inline

Definition at line 50 of file G4ParticleHPElasticFS.hh.

51  {
52  hasXsec = false;
53  theCoefficients = 0;
54  theProbArray = 0;
55 
56  repFlag = 0;
57  tE_of_repFlag3 = 0.0;
58  targetMass = 0.0;
59  frameFlag = 0;
60  }

Here is the caller graph for this function:

G4ParticleHPElasticFS::~G4ParticleHPElasticFS ( )
inline

Definition at line 61 of file G4ParticleHPElasticFS.hh.

62  {
63  if(theCoefficients!=0) delete theCoefficients;
64  if(theProbArray!=0) delete theProbArray;
65  }

Member Function Documentation

G4HadFinalState * G4ParticleHPElasticFS::ApplyYourself ( const G4HadProjectile theTrack)
virtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 200 of file G4ParticleHPElasticFS.cc.

201  {
202 // G4cout << "G4ParticleHPElasticFS::ApplyYourself+"<<G4endl;
203  if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
204  theResult.Get()->Clear();
205  G4double eKinetic = theTrack.GetKineticEnergy();
206  const G4HadProjectile *incidentParticle = &theTrack;
207  G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition() ));
208  theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
209  theNeutron.SetKineticEnergy( eKinetic );
210 // G4cout << "G4ParticleHPElasticFS::ApplyYourself++"<<eKinetic<<" "<<G4endl;
211 // G4cout << "CMSVALUES 0 "<<theNeutron.GetTotalMomentum()<<G4endl;
212 
214  G4Nucleus aNucleus;
215  G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
216  theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
217  //t theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) ); //TESTPHP
218 // G4cout << "Nucleus-test"<<" "<<targetMass<<" ";
219 // G4cout << theTarget.GetMomentum().x()<<" ";
220 // G4cout << theTarget.GetMomentum().y()<<" ";
221 // G4cout << theTarget.GetMomentum().z()<<G4endl;
222 
223  // neutron and target defined as reaction products.
224 
225 // prepare lorentz-transformation to Lab.
226 
227  G4ThreeVector the3Neutron = theNeutron.GetMomentum();
228  G4double nEnergy = theNeutron.GetTotalEnergy();
229  G4ThreeVector the3Target = theTarget.GetMomentum();
230 // cout << "@@@" << the3Target<<G4endl;
231  G4double tEnergy = theTarget.GetTotalEnergy();
232  G4ReactionProduct theCMS;
233  G4double totE = nEnergy+tEnergy;
234  G4ThreeVector the3CMS = the3Target+the3Neutron;
235  theCMS.SetMomentum(the3CMS);
236  G4double cmsMom = std::sqrt(the3CMS*the3CMS);
237  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
238  theCMS.SetMass(sqrts);
239  theCMS.SetTotalEnergy(totE);
240 
241  // data come as fcn of n-energy in nuclear rest frame
242  G4ReactionProduct boosted;
243  boosted.Lorentz(theNeutron, theTarget);
244  eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
245  G4double cosTh = -2;
246  if(repFlag == 1)
247  {
248  cosTh = theCoefficients->SampleElastic(eKinetic);
249  }
250 
251  else if (repFlag==2)
252  {
253  cosTh = theProbArray->Sample(eKinetic);
254  }
255  else if (repFlag==3)
256  {
257  if ( eKinetic <= tE_of_repFlag3 )
258  {
259  cosTh = theCoefficients->SampleElastic(eKinetic);
260  }
261  else
262  {
263  cosTh = theProbArray->Sample(eKinetic);
264  }
265  }
266  else if (repFlag==0)
267  {
268  cosTh = 2.*G4UniformRand()-1.;
269  }
270  else
271  {
272  G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
273  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
274  }
275  if(cosTh<-1.1) { return 0; }
276  G4double phi = twopi*G4UniformRand();
277  G4double theta = std::acos(cosTh);
278  G4double sinth = std::sin(theta);
279  if (frameFlag == 1) // final state data given in target rest frame.
280  {
281  // we have the scattering angle, now we need the energy, then do the
282  // boosting.
283  // relativistic elastic scattering energy angular correlation:
284  theNeutron.Lorentz(theNeutron, theTarget);
285  G4double e0 = theNeutron.GetTotalEnergy();
286  G4double p0 = theNeutron.GetTotalMomentum();
287  G4double mN = theNeutron.GetMass();
288  G4double mT = theTarget.GetMass();
289  G4double eE = e0+mT;
290  G4double ap = (mT+eE)*(mT-eE) + (p0+mN)*(p0-mN);
291  G4double a = 4*(eE+p0*cosTh)*(eE-p0*cosTh);
292  G4double b = 4*ap*p0*cosTh;
293  G4double c = (2.*eE*mN-ap)*(2.*eE*mN+ap);
294  G4double en = (-b+std::sqrt(b*b - 4*a*c) )/(2*a);
295  G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
296  theNeutron.SetMomentum(tempVector);
297  theNeutron.SetTotalEnergy(std::sqrt(en*en+theNeutron.GetMass()*theNeutron.GetMass()));
298  // first to lab
299  theNeutron.Lorentz(theNeutron, -1.*theTarget);
300  // now to CMS
301  theNeutron.Lorentz(theNeutron, theCMS);
302  theTarget.SetMomentum(-theNeutron.GetMomentum());
303  theTarget.SetTotalEnergy(theNeutron.GetTotalEnergy());
304  // and back to lab
305  theNeutron.Lorentz(theNeutron, -1.*theCMS);
306  theTarget.Lorentz(theTarget, -1.*theCMS);
307 //111005 Protection for not producing 0 kinetic energy target
308  if ( theNeutron.GetKineticEnergy() <= 0 ) theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
309  if ( theTarget.GetKineticEnergy() <= 0 ) theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
310  }
311  else if (frameFlag == 2) // CMS
312  {
313  theNeutron.Lorentz(theNeutron, theCMS);
314  theTarget.Lorentz(theTarget, theCMS);
315  G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
316  G4ThreeVector cmsMom_tmp=theNeutron.GetMomentum(); // for neutron direction in CMS
317  G4double cms_theta=cmsMom_tmp.theta();
318  G4double cms_phi=cmsMom_tmp.phi();
319  G4ThreeVector tempVector;
320  tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
321  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
322  -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
323  tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
324  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
325  +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
326  tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
327  -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
328  tempVector *= en;
329  theNeutron.SetMomentum(tempVector);
330  theTarget.SetMomentum(-tempVector);
331  G4double tP = theTarget.GetTotalMomentum();
332  G4double tM = theTarget.GetMass();
333  theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
334 
335 /*
336 For debug purpose.
337 Same transformation G4ReactionProduct.Lorentz() by 4vectors
338 {
339 G4LorentzVector n4p = G4LorentzVector ( theNeutron.GetMomentum() , theNeutron.GetKineticEnergy() + theNeutron.GetMass() );
340 G4cout << "before " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
341 G4LorentzVector cm4p = G4LorentzVector ( theCMS.GetMomentum() , theCMS.GetKineticEnergy() + theCMS.GetMass() );
342 n4p.boost( cm4p.boostVector() );
343 G4cout << cm4p/eV << G4endl;
344 G4cout << "after " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
345 }
346 */
347 
348  theNeutron.Lorentz(theNeutron, -1.*theCMS);
349 //080904 Add Protection for very low energy (1e-6eV) scattering
350  if ( theNeutron.GetKineticEnergy() <= 0 )
351  {
352  //theNeutron.SetMomentum( G4ThreeVector(0) );
353  //theNeutron.SetTotalEnergy ( theNeutron.GetMass() );
354 //110822 Protection for not producing 0 kinetic energy neutron
355  theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
356  }
357 
358  theTarget.Lorentz(theTarget, -1.*theCMS);
359 //080904 Add Protection for very low energy (1e-6eV) scattering
360  if ( theTarget.GetKineticEnergy() < 0 )
361  {
362  //theTarget.SetMomentum( G4ThreeVector(0) );
363  //theTarget.SetTotalEnergy ( theTarget.GetMass() );
364 //110822 Protection for not producing 0 kinetic energy target
365  theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
366  }
367  }
368  else
369  {
370  G4cout <<"Value of frameFlag (1=LAB, 2=CMS): "<<frameFlag;
371  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::ApplyYourSelf frameflag incorrect");
372  }
373  // now all in Lab
374 // nun den recoil generieren...und energy change, momentum change angeben.
375  theResult.Get()->SetEnergyChange(theNeutron.GetKineticEnergy());
376  theResult.Get()->SetMomentumChange(theNeutron.GetMomentum().unit());
377  G4DynamicParticle* theRecoil = new G4DynamicParticle;
378  theRecoil->SetDefinition( G4IonTable::GetIonTable()->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0 ) );
379  theRecoil->SetMomentum(theTarget.GetMomentum());
380  theResult.Get()->AddSecondary(theRecoil);
381 // G4cout << "G4ParticleHPElasticFS::ApplyYourself 10+"<<G4endl;
382  // postpone the tracking of the primary neutron
384  return theResult.Get();
385  }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4Cache< G4HadFinalState * > theResult
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void SetMomentum(const G4double x, const G4double y, const G4double z)
value_type & Get() const
Definition: G4Cache.hh:282
void setY(double)
void setZ(double)
void setX(double)
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetStatusChange(G4HadFinalStateStatus aS)
tuple b
Definition: test.py:12
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:113
void SetMass(const G4double mas)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
G4double Sample(G4double x)
G4double GetKineticEnergy() const
void SetTotalEnergy(const G4double en)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
double phi() const
const G4LorentzVector & Get4Momentum() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
double theta() const
G4double GetKineticEnergy() const
void SetEnergyChange(G4double anEnergy)
G4double SampleElastic(G4double anEnergy)
G4double GetTotalEnergy() const
G4double GetPDGMass() const
G4ThreeVector GetMomentum() const
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
const G4Material * GetMaterial() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
tuple c
Definition: test.py:13
void Put(const value_type &val) const
Definition: G4Cache.hh:286
void SetMomentumChange(const G4ThreeVector &aV)
G4double GetMass() const

Here is the call graph for this function:

void G4ParticleHPElasticFS::Init ( G4double  A,
G4double  Z,
G4int  M,
G4String dirName,
G4String aFSType,
G4ParticleDefinition  
)
virtual

Implements G4ParticleHPFinalState.

Definition at line 54 of file G4ParticleHPElasticFS.cc.

55  {
56  G4String tString = "/FS";
57  G4bool dbool;
58  G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
59  G4String filename = aFile.GetName();
60  SetAZMs( A, Z, M, aFile );
61  //theBaseA = aFile.GetA();
62  //theBaseZ = aFile.GetZ();
63  if(!dbool)
64  {
65  hasAnyData = false;
66  hasFSData = false;
67  hasXsec = false;
68  return;
69  }
70  //130205 For compressed data files
71  std::istringstream theData(std::ios::in);
73  //130205 END
74  theData >> repFlag >> targetMass >> frameFlag;
75  if(repFlag==1)
76  {
77  G4int nEnergy;
78  theData >> nEnergy;
79  theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
80  theCoefficients->InitInterpolation(theData);
81  G4double temp, energy;
82  G4int tempdep, nLegendre;
83  G4int i, ii;
84  for (i=0; i<nEnergy; i++)
85  {
86  theData >> temp >> energy >> tempdep >> nLegendre;
87  energy *=eV;
88  theCoefficients->Init(i, energy, nLegendre);
89  theCoefficients->SetTemperature(i, temp);
90  G4double coeff=0;
91  for(ii=0; ii<nLegendre; ii++)
92  {
93  // load legendre coefficients.
94  theData >> coeff;
95  theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
96  }
97  }
98  }
99  else if (repFlag==2)
100  {
101  G4int nEnergy;
102  theData >> nEnergy;
103  theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
104  theProbArray->InitInterpolation(theData);
105  G4double temp, energy;
106  G4int tempdep, nPoints;
107  for(G4int i=0; i<nEnergy; i++)
108  {
109  theData >> temp >> energy >> tempdep >> nPoints;
110  energy *= eV;
111  theProbArray->InitInterpolation(i, theData);
112  theProbArray->SetT(i, temp);
113  theProbArray->SetX(i, energy);
114  G4double prob, costh;
115  for(G4int ii=0; ii<nPoints; ii++)
116  {
117  // fill probability arrays.
118  theData >> costh >> prob;
119  theProbArray->SetX(i, ii, costh);
120  theProbArray->SetY(i, ii, prob);
121  }
122  theProbArray->DoneSetXY( i );
123  }
124  }
125  else if ( repFlag==3 )
126  {
127  G4int nEnergy_Legendre;
128  theData >> nEnergy_Legendre;
129  if ( nEnergy_Legendre <= 0 ) {
130  std::stringstream iss;
131  iss << "G4ParticleHPElasticFS::Init Data Error repFlag is 3 but nEnergy_Legendre <= 0";
132  iss << "Z, A and M of problematic file is " << theNDLDataZ << ", " << theNDLDataA << " and " << theNDLDataM << " respectively.";
133  throw G4HadronicException(__FILE__, __LINE__, iss.str() );
134  }
135  theCoefficients = new G4ParticleHPLegendreStore( nEnergy_Legendre );
136  theCoefficients->InitInterpolation( theData );
137  G4double temp, energy;
138  G4int tempdep, nLegendre;
139  //G4int i, ii;
140  for ( G4int i = 0 ; i < nEnergy_Legendre ; i++ )
141  {
142  theData >> temp >> energy >> tempdep >> nLegendre;
143  energy *=eV;
144  theCoefficients->Init( i , energy , nLegendre );
145  theCoefficients->SetTemperature( i , temp );
146  G4double coeff = 0;
147  for (G4int ii = 0 ; ii < nLegendre ; ii++ )
148  {
149  // load legendre coefficients.
150  theData >> coeff;
151  theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
152  }
153  }
154 
155  tE_of_repFlag3 = energy;
156 
157  G4int nEnergy_Prob;
158  theData >> nEnergy_Prob;
159  theProbArray = new G4ParticleHPPartial( nEnergy_Prob , nEnergy_Prob );
160  theProbArray->InitInterpolation( theData );
161  G4int nPoints;
162  for ( G4int i=0 ; i < nEnergy_Prob ; i++ )
163  {
164  theData >> temp >> energy >> tempdep >> nPoints;
165 
166  energy *= eV;
167 
168 // consistency check
169  if ( i == 0 )
170  //if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines
171  if ( std::abs( energy - tE_of_repFlag3 ) / tE_of_repFlag3 > 1.0e-15 )
172  G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl;
173 
174  theProbArray->InitInterpolation( i , theData );
175  theProbArray->SetT( i , temp );
176  theProbArray->SetX( i , energy );
177  G4double prob, costh;
178  for( G4int ii = 0 ; ii < nPoints ; ii++ )
179  {
180  // fill probability arrays.
181  theData >> costh >> prob;
182  theProbArray->SetX( i , ii , costh );
183  theProbArray->SetY( i , ii , prob );
184  }
185  theProbArray->DoneSetXY( i );
186  }
187  }
188  else if (repFlag==0)
189  {
190  theData >> frameFlag;
191  }
192  else
193  {
194  G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
195  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
196  }
197  //130205 For compressed data files(theData changed from ifstream to istringstream)
198  //theData.close();
199  }
static G4ParticleHPManager * GetInstance()
void Init(G4int i, G4double e, G4int n)
void InitInterpolation(G4int i, std::istream &aDataFile)
void SetY(G4int i, G4int j, G4double y)
void SetAZMs(G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
void GetDataStream(G4String, std::istringstream &iss)
int G4int
Definition: G4Types.hh:78
void InitInterpolation(std::istream &aDataFile)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
static constexpr double eV
Definition: G4SIunits.hh:215
void SetTemperature(G4int i, G4double temp)
void SetX(G4int i, G4double x)
G4double energy(const ThreeVector &p, const G4double m)
void SetCoeff(G4int i, G4int l, G4double coeff)
#define G4endl
Definition: G4ios.hh:61
G4ParticleHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
void SetT(G4int i, G4double x)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4ParticleHPFinalState* G4ParticleHPElasticFS::New ( )
inlinevirtual

Implements G4ParticleHPFinalState.

Definition at line 68 of file G4ParticleHPElasticFS.hh.

Here is the call graph for this function:


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