Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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::ifstream & aDataFile)
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 }
97 
98 void G4NeutronHPAngular::SampleAndUpdate(G4ReactionProduct & aHadron)
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 }
Generated on Sat May 25 2013 14:34:03 for Geant4 by   doxygen 1.8.4