source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPAngular.cc@ 1036

Last change on this file since 1036 was 962, checked in by garnier, 17 years ago

update processes
File size: 10.0 KB
Line 
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//
33#include "G4NeutronHPAngular.hh"
34
35void G4NeutronHPAngular::Init(std::ifstream & aDataFile)
36{
37// G4cout << "here we are entering the Angular Init"<<G4endl;
38 aDataFile >> theAngularDistributionType >> targetMass;
39 aDataFile >> frameFlag;
40 if(theAngularDistributionType == 0)
41 {
42 theIsoFlag = true;
43 }
44 else if(theAngularDistributionType==1)
45 {
46 G4int nEnergy;
47 aDataFile >> nEnergy;
48 theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
49 theCoefficients->InitInterpolation(aDataFile);
50 G4double temp, energy;
51 G4int tempdep, nLegendre;
52 G4int i, ii;
53 for (i=0; i<nEnergy; i++)
54 {
55 aDataFile >> temp >> energy >> tempdep >> nLegendre;
56 energy *=eV;
57 theCoefficients->Init(i, energy, nLegendre);
58 theCoefficients->SetTemperature(i, temp);
59 G4double coeff=0;
60 for(ii=0; ii<nLegendre; ii++)
61 {
62 aDataFile >> coeff;
63 theCoefficients->SetCoeff(i, ii+1, coeff);
64 }
65 }
66 }
67 else if (theAngularDistributionType==2)
68 {
69 G4int nEnergy;
70 aDataFile >> nEnergy;
71 theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
72 theProbArray->InitInterpolation(aDataFile);
73 G4double temp, energy;
74 G4int tempdep;
75 for(G4int i=0; i<nEnergy; i++)
76 {
77 aDataFile >> temp >> energy >> tempdep;
78 energy *= eV;
79 theProbArray->SetT(i, temp);
80 theProbArray->SetX(i, energy);
81 theProbArray->InitData(i, aDataFile);
82 }
83 }
84 else
85 {
86 theIsoFlag = false;
87 G4cout << "unknown distribution found for Angular"<<G4endl;
88 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
89 }
90}
91
92void G4NeutronHPAngular::SampleAndUpdate(G4ReactionProduct & aHadron)
93{
94 if(theIsoFlag)
95 {
96// G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
97// @@@ add code for isotropic emission in CMS.
98 G4double costheta = 2.*G4UniformRand()-1;
99 G4double theta = std::acos(costheta);
100 G4double phi = twopi*G4UniformRand();
101 G4double sinth = std::sin(theta);
102 G4double en = aHadron.GetTotalMomentum();
103 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
104 aHadron.SetMomentum( temp );
105 aHadron.Lorentz(aHadron, -1.*theTarget);
106 }
107 else
108 {
109 if(frameFlag == 1) // LAB
110 {
111 G4double en = aHadron.GetTotalMomentum();
112 G4ReactionProduct boosted;
113 boosted.Lorentz(theNeutron, theTarget);
114 G4double kineticEnergy = boosted.GetKineticEnergy();
115 G4double cosTh = 0.0;
116 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
117 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
118 G4double theta = std::acos(cosTh);
119 G4double phi = twopi*G4UniformRand();
120 G4double sinth = std::sin(theta);
121 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
122 aHadron.SetMomentum( temp );
123 }
124 else if(frameFlag == 2) // costh in CMS
125 {
126 G4ReactionProduct boostedN;
127 boostedN.Lorentz(theNeutron, theTarget);
128 G4double kineticEnergy = boostedN.GetKineticEnergy();
129
130 G4double cosTh = 0.0;
131 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
132 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
133
134//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern)
135/*
136 if(theAngularDistributionType == 1) // LAB
137 {
138 G4double en = aHadron.GetTotalMomentum();
139 G4ReactionProduct boosted;
140 boosted.Lorentz(theNeutron, theTarget);
141 G4double kineticEnergy = boosted.GetKineticEnergy();
142 G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
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(theAngularDistributionType == 2) // costh in CMS
150 {
151*/
152
153// G4ReactionProduct boostedN;
154// boostedN.Lorentz(theNeutron, theTarget);
155// G4double kineticEnergy = boostedN.GetKineticEnergy();
156// G4double cosTh = theProbArray->Sample(kineticEnergy);
157
158 G4double theta = std::acos(cosTh);
159 G4double phi = twopi*G4UniformRand();
160 G4double sinth = std::sin(theta);
161
162 G4ThreeVector temp(sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); //CMS
163
164//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
165/*
166 G4double en = aHadron.GetTotalEnergy(); // Target rest
167
168 // get trafo from Target rest frame to CMS
169 G4ReactionProduct boostedT;
170 boostedT.Lorentz(theTarget, theTarget);
171
172 G4ThreeVector the3Neutron = boostedN.GetMomentum();
173 G4double nEnergy = boostedN.GetTotalEnergy();
174 G4ThreeVector the3Target = boostedT.GetMomentum();
175 G4double tEnergy = boostedT.GetTotalEnergy();
176 G4double totE = nEnergy+tEnergy;
177 G4ThreeVector the3trafo = -the3Target-the3Neutron;
178 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
179 trafo.SetMomentum(the3trafo);
180 G4double cmsMom = std::sqrt(the3trafo*the3trafo);
181 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
182 trafo.SetMass(sqrts);
183 trafo.SetTotalEnergy(totE);
184
185 G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
186 G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
187 G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
188 fac*=gamma;
189
190 G4double mom;
191// For G4FPE_DEBUG ON
192 G4double mom2 = ( en*fac*en*fac -
193 (fac*fac - gamma*gamma)*
194 (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
195 );
196 if ( mom2 > 0.0 )
197 mom = std::sqrt( mom2 );
198 else
199 mom = 0.0;
200
201 mom = -en*fac - mom;
202 mom /= (fac*fac-gamma*gamma);
203 temp = mom*temp;
204
205 aHadron.SetMomentum( temp ); // now all in CMS
206 aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
207 aHadron.Lorentz(aHadron, trafo); // now in target rest frame
208*/
209 // Determination of the hadron kinetic energy in CMS
210 // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
211 // kineticEnergy is incident neutron kinetic energy in CMS (or target frame)
212 G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
213 G4double A1 = theTarget.GetMass()/boostedN.GetMass();
214 G4double A1prim = aHadron.GetMass()/ boostedN.GetMass();
215 G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*kineticEnergy+(1+A1)*QValue);
216 G4double totalE = kinE + aHadron.GetMass();
217 G4double mom2 = totalE*totalE - aHadron.GetMass()*aHadron.GetMass();
218 G4double mom;
219 if ( mom2 > 0.0 ) mom = std::sqrt( mom2 );
220 else mom = 0.0;
221
222 aHadron.SetMomentum( mom*temp ); // Set momentum in CMS
223 aHadron.SetKineticEnergy(kinE); // Set kinetic energy in CMS
224
225 // get trafo from Target rest frame to CMS
226 G4ReactionProduct boostedT;
227 boostedT.Lorentz(theTarget, theTarget);
228
229 G4ThreeVector the3Neutron = boostedN.GetMomentum();
230 G4double nEnergy = boostedN.GetTotalEnergy();
231 G4ThreeVector the3Target = boostedT.GetMomentum();
232 G4double tEnergy = boostedT.GetTotalEnergy();
233 G4double totE = nEnergy+tEnergy;
234 G4ThreeVector the3trafo = -the3Target-the3Neutron;
235 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
236 trafo.SetMomentum(the3trafo);
237 G4double cmsMom = std::sqrt(the3trafo*the3trafo);
238 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
239 trafo.SetMass(sqrts);
240 trafo.SetTotalEnergy(totE);
241
242 aHadron.Lorentz(aHadron, trafo);
243
244 }
245 else
246 {
247 throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
248 }
249 }
250 aHadron.Lorentz(aHadron, -1.*theTarget);
251// G4cout << aHadron.GetMomentum()<<" ";
252// G4cout << aHadron.GetTotalMomentum()<<G4endl;
253}
Note: See TracBrowser for help on using the repository browser.