Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4NeutronHPAngular Class Reference

#include <G4NeutronHPAngular.hh>

Public Member Functions

 G4NeutronHPAngular ()
 
 ~G4NeutronHPAngular ()
 
void Init (std::istream &aDataFile)
 
void SampleAndUpdate (G4ReactionProduct &aNeutron)
 
void SetTarget (const G4ReactionProduct &aTarget)
 
void SetNeutron (const G4ReactionProduct &aNeutron)
 
G4double GetTargetMass ()
 

Detailed Description

Definition at line 41 of file G4NeutronHPAngular.hh.

Constructor & Destructor Documentation

G4NeutronHPAngular::G4NeutronHPAngular ( )
inline

Definition at line 45 of file G4NeutronHPAngular.hh.

46  {
47 //TKDB110511
48  //theAngularDistributionType = 0;
49  //theIsoFlag = false;
50  theIsoFlag = true;
51 // TKDB
52  theCoefficients = 0;
53  theProbArray = 0;
54  }
G4NeutronHPAngular::~G4NeutronHPAngular ( )
inline

Definition at line 56 of file G4NeutronHPAngular.hh.

57  {
58 // TKDB
59  delete theCoefficients;
60  delete theProbArray;
61  }

Member Function Documentation

G4double G4NeutronHPAngular::GetTargetMass ( )
inline

Definition at line 71 of file G4NeutronHPAngular.hh.

Referenced by G4NeutronHPInelasticBaseFS::BaseApply().

71 { return targetMass; }
void G4NeutronHPAngular::Init ( std::istream &  aDataFile)

Definition at line 39 of file G4NeutronHPAngular.cc.

References energy(), python.hepunit::eV, G4cout, G4endl, G4NeutronHPLegendreStore::Init(), G4NeutronHPPartial::InitData(), G4NeutronHPPartial::InitInterpolation(), G4NeutronHPLegendreStore::InitInterpolation(), G4NeutronHPLegendreStore::SetCoeff(), G4NeutronHPPartial::SetT(), G4NeutronHPLegendreStore::SetTemperature(), and G4NeutronHPPartial::SetX().

Referenced by G4NeutronHPFSFissionFS::Init(), G4NeutronHPFissionBaseFS::Init(), G4NeutronHPInelasticBaseFS::Init(), G4NeutronHPInelasticCompFS::Init(), and G4FissionLibrary::Init().

40 {
41 // G4cout << "here we are entering the Angular Init"<<G4endl;
42  aDataFile >> theAngularDistributionType >> targetMass;
43  aDataFile >> frameFlag;
44  if(theAngularDistributionType == 0 )
45  {
46  theIsoFlag = true;
47  }
48  else if(theAngularDistributionType==1)
49  {
50  theIsoFlag = false;
51  G4int nEnergy;
52  aDataFile >> nEnergy;
53  theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
54  theCoefficients->InitInterpolation(aDataFile);
55  G4double temp, energy;
56  G4int tempdep, nLegendre;
57  G4int i, ii;
58  for (i=0; i<nEnergy; i++)
59  {
60  aDataFile >> temp >> energy >> tempdep >> nLegendre;
61  energy *=eV;
62  theCoefficients->Init(i, energy, nLegendre);
63  theCoefficients->SetTemperature(i, temp);
64  G4double coeff=0;
65  for(ii=0; ii<nLegendre; ii++)
66  {
67  aDataFile >> coeff;
68  theCoefficients->SetCoeff(i, ii+1, coeff);
69  }
70  }
71  }
72  else if (theAngularDistributionType==2)
73  {
74  theIsoFlag = false;
75  G4int nEnergy;
76  aDataFile >> nEnergy;
77  theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
78  theProbArray->InitInterpolation(aDataFile);
79  G4double temp, energy;
80  G4int tempdep;
81  for(G4int i=0; i<nEnergy; i++)
82  {
83  aDataFile >> temp >> energy >> tempdep;
84  energy *= eV;
85  theProbArray->SetT(i, temp);
86  theProbArray->SetX(i, energy);
87  theProbArray->InitData(i, aDataFile);
88  }
89  }
90  else
91  {
92  theIsoFlag = false;
93  G4cout << "unknown distribution found for Angular"<<G4endl;
94  throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
95  }
96 }
G4NeutronHPLegendreStore::Init
void Init(G4int i, G4double e, G4int n)
Definition: G4NeutronHPLegendreStore.hh:51
G4NeutronHPLegendreStore::SetCoeff
void SetCoeff(G4int i, G4int l, G4double coeff)
Definition: G4NeutronHPLegendreStore.hh:58
G4int
int G4int
Definition: G4Types.hh:78
G4NeutronHPPartial::SetT
void SetT(G4int i, G4double x)
Definition: G4NeutronHPPartial.hh:108
energy
double precision function energy(A, Z)
Definition: dpm25nuc6.f:4106
G4cout
G4GLOB_DLL std::ostream G4cout
G4NeutronHPPartial::SetX
void SetX(G4int i, G4double x)
Definition: G4NeutronHPPartial.hh:107
G4NeutronHPPartial::InitData
void InitData(G4int i, std::istream &aDataFile, G4double unit=1.)
Definition: G4NeutronHPPartial.hh:83
G4NeutronHPLegendreStore::InitInterpolation
void InitInterpolation(std::istream &aDataFile)
Definition: G4NeutronHPLegendreStore.hh:78
G4endl
#define G4endl
Definition: G4ios.hh:61
G4NeutronHPPartial::InitInterpolation
void InitInterpolation(G4int i, std::istream &aDataFile)
Definition: G4NeutronHPPartial.hh:60
G4double
double G4double
Definition: G4Types.hh:76
G4NeutronHPLegendreStore::SetTemperature
void SetTemperature(G4int i, G4double temp)
Definition: G4NeutronHPLegendreStore.hh:57
void G4NeutronHPAngular::SampleAndUpdate ( G4ReactionProduct &  aNeutron)

Definition at line 98 of file G4NeutronHPAngular.cc.

References G4cout, G4endl, G4UniformRand, G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4ReactionProduct::Lorentz(), G4NeutronHPPartial::Sample(), G4NeutronHPLegendreStore::SampleMax(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTotalEnergy(), and python.hepunit::twopi.

Referenced by G4NeutronHPFSFissionFS::ApplyYourself(), G4NeutronHPFissionBaseFS::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), and G4NeutronHPInelasticCompFS::CompositeApply().

99 {
100 
101  //********************************************************************
102  //EMendoza -> sampling can be isotropic in LAB or in CMS
103  /*
104  if(theIsoFlag)
105  {
106 // G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
107 // @@@ add code for isotropic emission in CMS.
108  G4double costheta = 2.*G4UniformRand()-1;
109  G4double theta = std::acos(costheta);
110  G4double phi = twopi*G4UniformRand();
111  G4double sinth = std::sin(theta);
112  G4double en = aHadron.GetTotalMomentum();
113  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
114  aHadron.SetMomentum( temp );
115  aHadron.Lorentz(aHadron, -1.*theTarget);
116  }
117  else
118  {
119  */
120  //********************************************************************
121  if(frameFlag == 1) // LAB
122  {
123  G4double en = aHadron.GetTotalMomentum();
124  G4ReactionProduct boosted;
125  boosted.Lorentz(theNeutron, theTarget);
126  G4double kineticEnergy = boosted.GetKineticEnergy();
127  G4double cosTh = 0.0;
128  //********************************************************************
129  //EMendoza --> sampling can be also isotropic
130  /*
131  if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
132  if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
133  */
134  //********************************************************************
135  if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
136  else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
137  else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
138  else{
139  G4cout << "unknown distribution found for Angular"<<G4endl;
140  throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
141  }
142  //********************************************************************
143  G4double theta = std::acos(cosTh);
144  G4double phi = twopi*G4UniformRand();
145  G4double sinth = std::sin(theta);
146  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
147  aHadron.SetMomentum( temp );
148  }
149  else if(frameFlag == 2) // costh in CMS
150  {
151  G4ReactionProduct boostedN;
152  boostedN.Lorentz(theNeutron, theTarget);
153  G4double kineticEnergy = boostedN.GetKineticEnergy();
154 
155  G4double cosTh = 0.0;
156  //********************************************************************
157  //EMendoza --> sampling can be also isotropic
158  /*
159  if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
160  if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
161  */
162  //********************************************************************
163  if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
164  else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
165  else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
166  else{
167  G4cout << "unknown distribution found for Angular"<<G4endl;
168  throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
169  }
170  //********************************************************************
171 //080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern)
172 /*
173  if(theAngularDistributionType == 1) // LAB
174  {
175  G4double en = aHadron.GetTotalMomentum();
176  G4ReactionProduct boosted;
177  boosted.Lorentz(theNeutron, theTarget);
178  G4double kineticEnergy = boosted.GetKineticEnergy();
179  G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
180  G4double theta = std::acos(cosTh);
181  G4double phi = twopi*G4UniformRand();
182  G4double sinth = std::sin(theta);
183  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
184  aHadron.SetMomentum( temp );
185  }
186  else if(theAngularDistributionType == 2) // costh in CMS {
187  }
188 */
189 
190 // G4ReactionProduct boostedN;
191 // boostedN.Lorentz(theNeutron, theTarget);
192 // G4double kineticEnergy = boostedN.GetKineticEnergy();
193 // G4double cosTh = theProbArray->Sample(kineticEnergy);
194 
195  G4double theta = std::acos(cosTh);
196  G4double phi = twopi*G4UniformRand();
197  G4double sinth = std::sin(theta);
198 
199  G4ThreeVector temp(sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); //CMS
200 
201 //080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
202 /*
203  G4double en = aHadron.GetTotalEnergy(); // Target rest
204 
205  // get trafo from Target rest frame to CMS
206  G4ReactionProduct boostedT;
207  boostedT.Lorentz(theTarget, theTarget);
208 
209  G4ThreeVector the3Neutron = boostedN.GetMomentum();
210  G4double nEnergy = boostedN.GetTotalEnergy();
211  G4ThreeVector the3Target = boostedT.GetMomentum();
212  G4double tEnergy = boostedT.GetTotalEnergy();
213  G4double totE = nEnergy+tEnergy;
214  G4ThreeVector the3trafo = -the3Target-the3Neutron;
215  G4ReactionProduct trafo; // for transformation from CMS to target rest frame
216  trafo.SetMomentum(the3trafo);
217  G4double cmsMom = std::sqrt(the3trafo*the3trafo);
218  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
219  trafo.SetMass(sqrts);
220  trafo.SetTotalEnergy(totE);
221 
222  G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
223  G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
224  G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
225  fac*=gamma;
226 
227  G4double mom;
228 // For G4FPE_DEBUG ON
229  G4double mom2 = ( en*fac*en*fac -
230  (fac*fac - gamma*gamma)*
231  (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
232  );
233  if ( mom2 > 0.0 )
234  mom = std::sqrt( mom2 );
235  else
236  mom = 0.0;
237 
238  mom = -en*fac - mom;
239  mom /= (fac*fac-gamma*gamma);
240  temp = mom*temp;
241 
242  aHadron.SetMomentum( temp ); // now all in CMS
243  aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
244  aHadron.Lorentz(aHadron, trafo); // now in target rest frame
245 */
246  // Determination of the hadron kinetic energy in CMS
247  // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
248  // kineticEnergy is incident neutron kinetic energy in CMS (or target frame)
249  G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
250  G4double A1 = theTarget.GetMass()/boostedN.GetMass();
251  G4double A1prim = aHadron.GetMass()/ boostedN.GetMass();
252  G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*kineticEnergy+(1+A1)*QValue);
253  G4double totalE = kinE + aHadron.GetMass();
254  G4double mom2 = totalE*totalE - aHadron.GetMass()*aHadron.GetMass();
255  G4double mom;
256  if ( mom2 > 0.0 ) mom = std::sqrt( mom2 );
257  else mom = 0.0;
258 
259  aHadron.SetMomentum( mom*temp ); // Set momentum in CMS
260  aHadron.SetKineticEnergy(kinE); // Set kinetic energy in CMS
261 
262  // get trafo from Target rest frame to CMS
263  G4ReactionProduct boostedT;
264  boostedT.Lorentz(theTarget, theTarget);
265 
266  G4ThreeVector the3Neutron = boostedN.GetMomentum();
267  G4double nEnergy = boostedN.GetTotalEnergy();
268  G4ThreeVector the3Target = boostedT.GetMomentum();
269  G4double tEnergy = boostedT.GetTotalEnergy();
270  G4double totE = nEnergy+tEnergy;
271  G4ThreeVector the3trafo = -the3Target-the3Neutron;
272  G4ReactionProduct trafo; // for transformation from CMS to target rest frame
273  trafo.SetMomentum(the3trafo);
274  G4double cmsMom = std::sqrt(the3trafo*the3trafo);
275  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
276  trafo.SetMass(sqrts);
277  trafo.SetTotalEnergy(totE);
278 
279  aHadron.Lorentz(aHadron, trafo);
280 
281  }
282  else
283  {
284  throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
285  }
286  aHadron.Lorentz(aHadron, -1.*theTarget);
287 // G4cout << aHadron.GetMomentum()<<" ";
288 // G4cout << aHadron.GetTotalMomentum()<<G4endl;
289 }
G4NeutronHPLegendreStore::SampleMax
G4double SampleMax(G4double energy)
Definition: G4NeutronHPLegendreStore.cc:117
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:199
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:165
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:147
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4NeutronHPPartial::Sample
G4double Sample(G4double x)
Definition: G4NeutronHPPartial.cc:103
void G4NeutronHPAngular::SetNeutron ( const G4ReactionProduct &  aNeutron)
inline

Definition at line 69 of file G4NeutronHPAngular.hh.

Referenced by G4FissionLibrary::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::InitDistributionInitialState(), G4NeutronHPFissionBaseFS::SetNeutron(), and G4NeutronHPFSFissionFS::SetNeutron().

69 { theNeutron = aNeutron; }
void G4NeutronHPAngular::SetTarget ( const G4ReactionProduct &  aTarget)
inline

Definition at line 67 of file G4NeutronHPAngular.hh.

Referenced by G4FissionLibrary::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::InitDistributionInitialState(), G4NeutronHPFissionBaseFS::SetTarget(), and G4NeutronHPFSFissionFS::SetTarget().

67 { theTarget = aTarget; }
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:46 for Geant4.10 by   doxygen 1.8.7