Geant4  10.01
G4NeutronHPAngular.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 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
32 // 110505 protection for object is created but not initialized
33 // 110510 delete above protection with more coordinated work to other classes
34 //
35 #include "G4NeutronHPAngular.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 
39 void G4NeutronHPAngular::Init(std::istream & aDataFile)
40 {
41 // G4cout << "here we are entering the Angular Init"<<G4endl;
42  aDataFile >> theAngularDistributionType >> targetMass;
43  aDataFile >> frameFlag;
45  {
46  theIsoFlag = true;
47  }
48  else if(theAngularDistributionType==1)
49  {
50  theIsoFlag = false;
51  G4int nEnergy;
52  aDataFile >> nEnergy;
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);
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 }
97 
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( *fCache.Get().theNeutron, *fCache.Get().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( *fCache.Get().theNeutron, *fCache.Get().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 = fCache.Get().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(*fCache.Get().theTarget, *fCache.Get().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.*(*fCache.Get().theTarget));
287 // G4cout << aHadron.GetMomentum()<<" ";
288 // G4cout << aHadron.GetTotalMomentum()<<G4endl;
289 }
G4double SampleMax(G4double energy)
G4double GetTotalMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
CLHEP::Hep3Vector G4ThreeVector
void SetKineticEnergy(const G4double en)
void Init(G4int i, G4double e, G4int n)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetCoeff(G4int i, G4int l, G4double coeff)
int G4int
Definition: G4Types.hh:78
void SetT(G4int i, G4double x)
void SetMass(const G4double mas)
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
void SetX(G4int i, G4double x)
void SetTotalEnergy(const G4double en)
void InitData(G4int i, std::istream &aDataFile, G4double unit=1.)
void Init(std::istream &aDataFile)
void InitInterpolation(std::istream &aDataFile)
G4NeutronHPPartial * theProbArray
G4double GetKineticEnergy() const
void SampleAndUpdate(G4ReactionProduct &aNeutron)
static const double eV
Definition: G4SIunits.hh:194
G4double GetTotalEnergy() const
G4Cache< toBeCached > fCache
G4double energy(const ThreeVector &p, const G4double m)
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
void InitInterpolation(G4int i, std::istream &aDataFile)
G4NeutronHPLegendreStore * theCoefficients
double G4double
Definition: G4Types.hh:76
void SetTemperature(G4int i, G4double temp)
G4double GetMass() const
G4double Sample(G4double x)