Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPElasticFS.cc
Go to the documentation of this file.
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 // 25-08-06 New Final State type (refFlag==3 , Legendre (Low Energy) + Probability (High Energy) )
31 // is added by T. KOI
32 // 080904 Add Protection for negative energy results in very low energy ( 1E-6 eV ) scattering by T. Koi
33 //
34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
35 //
36 #include "G4ParticleHPElasticFS.hh"
37 #include "G4ParticleHPManager.hh"
38 
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ReactionProduct.hh"
42 #include "G4Nucleus.hh"
43 #include "G4Proton.hh"
44 #include "G4Deuteron.hh"
45 #include "G4Triton.hh"
46 #include "G4Alpha.hh"
47 #include "G4ThreeVector.hh"
48 #include "G4LorentzVector.hh"
49 #include "G4IonTable.hh"
50 #include "G4ParticleHPDataUsed.hh"
51 #include "G4Pow.hh"
52 #include "zlib.h"
53 
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  }
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.
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 G4ParticleHPManager * GetInstance()
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
void Init(G4int i, G4double e, G4int n)
G4Cache< G4HadFinalState * > theResult
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void InitInterpolation(G4int i, std::istream &aDataFile)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
value_type & Get() const
Definition: G4Cache.hh:282
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 setY(double)
void InitInterpolation(std::istream &aDataFile)
void setZ(double)
void setX(double)
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetStatusChange(G4HadFinalStateStatus aS)
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)
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
double A(double temperature)
const G4ParticleDefinition * GetDefinition() const
G4double Sample(G4double x)
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
void SetTotalEnergy(const G4double en)
G4HadFinalState * ApplyYourself(const G4HadProjectile &theTrack)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
static constexpr double eV
Definition: G4SIunits.hh:215
double phi() const
const G4LorentzVector & Get4Momentum() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
double theta() const
G4double GetKineticEnergy() const
void SetTemperature(G4int i, G4double temp)
void SetX(G4int i, G4double x)
void SetEnergyChange(G4double anEnergy)
G4double SampleElastic(G4double anEnergy)
G4double GetTotalEnergy() const
G4double GetPDGMass() const
G4double energy(const ThreeVector &p, const G4double m)
Hep3Vector unit() const
void SetCoeff(G4int i, G4int l, G4double coeff)
G4ThreeVector GetMomentum() const
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
G4ParticleHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
void SetT(G4int i, G4double x)
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