| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. How
31 // 080612 bug fix contribution from Benoit Pir
32 // 110505 protection for object is created but
33 // 110510 delete above protection with more co
34 //
35 // P. Arce, June-2014 Conversion neutron_hp to
36 //
37 #include "G4ParticleHPAngular.hh"
38
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41
42 void G4ParticleHPAngular::Init(std::istream& a
43 {
44 // G4cout << "here we are entering the Angu
45 aDataFile >> theAngularDistributionType >> t
46 aDataFile >> frameFlag;
47 if (theAngularDistributionType == 0) {
48 theIsoFlag = true;
49 }
50 else if (theAngularDistributionType == 1) {
51 theIsoFlag = false;
52 G4int nEnergy;
53 aDataFile >> nEnergy;
54 theCoefficients = new G4ParticleHPLegendre
55 theCoefficients->InitInterpolation(aDataFi
56 G4double temp, energy;
57 G4int tempdep, nLegendre;
58 G4int i, ii;
59 for (i = 0; i < nEnergy; i++) {
60 aDataFile >> temp >> energy >> tempdep >
61 energy *= eV;
62 theCoefficients->Init(i, energy, nLegend
63 theCoefficients->SetTemperature(i, temp)
64 G4double coeff = 0;
65 for (ii = 0; ii < nLegendre; ii++) {
66 aDataFile >> coeff;
67 theCoefficients->SetCoeff(i, ii + 1, c
68 }
69 }
70 }
71 else if (theAngularDistributionType == 2) {
72 theIsoFlag = false;
73 G4int nEnergy;
74 aDataFile >> nEnergy;
75 theProbArray = new G4ParticleHPPartial(nEn
76 theProbArray->InitInterpolation(aDataFile)
77 G4double temp, energy;
78 G4int tempdep;
79 for (G4int i = 0; i < nEnergy; i++) {
80 aDataFile >> temp >> energy >> tempdep;
81 energy *= eV;
82 theProbArray->SetT(i, temp);
83 theProbArray->SetX(i, energy);
84 theProbArray->InitData(i, aDataFile);
85 }
86 }
87 else {
88 theIsoFlag = false;
89 G4cout << "unknown distribution found for
90 throw G4HadronicException(__FILE__, __LINE
91 }
92 }
93
94 void G4ParticleHPAngular::SampleAndUpdate(G4Re
95 {
96 //******************************************
97 // EMendoza -> sampling can be isotropic in
98 /*
99 if(theIsoFlag)
100 {
101 // G4cout << "Angular result "<<aHadron.GetTo
102 // @@@ add code for isotropic emission in CMS.
103 G4double costheta = 2.*G4UniformRand()-1
104 G4double theta = std::acos(costheta);
105 G4double phi = twopi*G4UniformRand();
106 G4double sinth = std::sin(theta);
107 G4double en = aHadron.GetTotalMomentum()
108 G4ThreeVector temp(en*sinth*std::cos(phi
109 aHadron.SetMomentum( temp );
110 aHadron.Lorentz(aHadron, -1.*theTarget);
111 }
112 else
113 {
114 */
115 //******************************************
116 if (frameFlag == 1) // LAB
117 {
118 G4double en = aHadron.GetTotalMomentum();
119 G4ReactionProduct boosted;
120 boosted.Lorentz(*fCache.Get().theProjectil
121 G4double kineticEnergy = boosted.GetKineti
122 G4double cosTh = 0.0;
123 //****************************************
124 // EMendoza --> sampling can be also isotr
125 /*
126 if(theAngularDistributionType == 1) cosTh
127 if(theAngularDistributionType == 2) cosTh
128 */
129 //****************************************
130 if (theIsoFlag) {
131 cosTh = 2. * G4UniformRand() - 1;
132 }
133 else if (theAngularDistributionType == 1)
134 cosTh = theCoefficients->SampleMax(kinet
135 }
136 else if (theAngularDistributionType == 2)
137 cosTh = theProbArray->Sample(kineticEner
138 }
139 else {
140 G4cout << "unknown distribution found fo
141 throw G4HadronicException(__FILE__, __LI
142 }
143 //****************************************
144 G4double theta = std::acos(cosTh);
145 G4double phi = twopi * G4UniformRand();
146 G4double sinth = std::sin(theta);
147 G4ThreeVector temp(en * sinth * std::cos(p
148 en * std::cos(theta));
149 aHadron.SetMomentum(temp);
150 }
151 else if (frameFlag == 2) // costh in CMS
152 {
153 G4ReactionProduct boostedN;
154 boostedN.Lorentz(*fCache.Get().theProjecti
155 G4double kineticEnergy = boostedN.GetKinet
156
157 G4double cosTh = 0.0;
158 //****************************************
159 // EMendoza --> sampling can be also isotr
160 /*
161 if(theAngularDistributionType == 1) cosTh
162 if(theAngularDistributionType == 2) cosTh
163 */
164 //****************************************
165 if (theIsoFlag) {
166 cosTh = 2. * G4UniformRand() - 1;
167 }
168 else if (theAngularDistributionType == 1)
169 cosTh = theCoefficients->SampleMax(kinet
170 }
171 else if (theAngularDistributionType == 2)
172 cosTh = theProbArray->Sample(kineticEner
173 }
174 else {
175 G4cout << "unknown distribution found fo
176 throw G4HadronicException(__FILE__, __LI
177 }
178 //****************************************
179 // 080612TK bug fix contribution from Beno
180 /*
181 if(theAngularDistributionType == 1) //
182 {
183 G4double en = aHadron.GetTotalMoment
184 G4ReactionProduct boosted;
185 boosted.Lorentz(theProjectile, theTa
186 G4double kineticEnergy = boosted.Get
187 G4double cosTh = theCoefficients->Sa
188 G4double theta = std::acos(cosTh);
189 G4double phi = twopi*G4UniformRand()
190 G4double sinth = std::sin(theta);
191 G4ThreeVector temp(en*sinth*std::cos
192 aHadron.SetMomentum( temp );
193 }
194 else if(theAngularDistributionType ==
195 }
196 */
197
198 // G4ReactionProduct boostedN;
199 // boostedN.Lorentz(theProjectile, th
200 // G4double kineticEnergy = boostedN.
201 // G4double cosTh = theProbArray->Sam
202
203 G4double theta = std::acos(cosTh);
204 G4double phi = twopi * G4UniformRand();
205 G4double sinth = std::sin(theta);
206
207 G4ThreeVector temp(sinth * std::cos(phi),
208
209 // 080612TK bug fix contribution from Beno
210 /*
211 G4double en = aHadron.GetTotalEnergy
212
213 // get trafo from Target rest frame
214 G4ReactionProduct boostedT;
215 boostedT.Lorentz(theTarget, theTarge
216
217 G4ThreeVector the3IncidentParticle =
218 G4double nEnergy = boostedN.GetTotal
219 G4ThreeVector the3Target = boostedT.
220 G4double tEnergy = boostedT.GetTotal
221 G4double totE = nEnergy+tEnergy;
222 G4ThreeVector the3trafo = -the3Targe
223 G4ReactionProduct trafo; // for tran
224 trafo.SetMomentum(the3trafo);
225 G4double cmsMom = std::sqrt(the3traf
226 G4double sqrts = std::sqrt((totE-cms
227 trafo.SetMass(sqrts);
228 trafo.SetTotalEnergy(totE);
229
230 G4double gamma = trafo.GetTotalEnerg
231 G4double cosalpha = temp*trafo.GetMo
232 G4double fac = cosalpha*trafo.GetTot
233 fac*=gamma;
234
235 G4double mom;
236 // For G4FPE_DEBUG ON
237 G4double mom2 = ( en*fac*en*fac -
238 (fac*fac - gamma*gamma)
239 (en*en - gamma*gamma*aH
240 );
241 if ( mom2 > 0.0 )
242 mom = std::sqrt( mom2 );
243 else
244 mom = 0.0;
245
246 mom = -en*fac - mom;
247 mom /= (fac*fac-gamma*gamma);
248 temp = mom*temp;
249
250 aHadron.SetMomentum( temp ); // now
251 aHadron.SetTotalEnergy( std::sqrt( m
252 aHadron.Lorentz(aHadron, trafo); //
253 */
254 // Determination of the hadron kinetic ene
255 // aHadron.GetKineticEnergy() is actually
256 // kineticEnergy is incident neutron kinet
257 G4double QValue = aHadron.GetKineticEnergy
258 G4double A1 = fCache.Get().theTarget->GetM
259 G4double A1prim = aHadron.GetMass() / boos
260 G4double kinE =
261 (A1 + 1 - A1prim) / (A1 + 1) / (A1 + 1)
262 G4double totalE = kinE + aHadron.GetMass()
263 G4double mom2 = totalE * totalE - aHadron.
264 G4double mom;
265 if (mom2 > 0.0)
266 mom = std::sqrt(mom2);
267 else
268 mom = 0.0;
269
270 aHadron.SetMomentum(mom * temp); // Set m
271 aHadron.SetKineticEnergy(kinE); // Set ki
272
273 // get trafo from Target rest frame to CMS
274 G4ReactionProduct boostedT;
275 boostedT.Lorentz(*fCache.Get().theTarget,
276
277 G4ThreeVector the3IncidentParticle = boost
278 G4double nEnergy = boostedN.GetTotalEnergy
279 G4ThreeVector the3Target = boostedT.GetMom
280 G4double tEnergy = boostedT.GetTotalEnergy
281 G4double totE = nEnergy + tEnergy;
282 G4ThreeVector the3trafo = -the3Target - th
283 G4ReactionProduct trafo; // for transform
284 trafo.SetMomentum(the3trafo);
285 G4double cmsMom = std::sqrt(the3trafo * th
286 G4double sqrts = std::sqrt((totE - cmsMom)
287 trafo.SetMass(sqrts);
288 trafo.SetTotalEnergy(totE);
289
290 aHadron.Lorentz(aHadron, trafo);
291 }
292 else {
293 throw G4HadronicException(__FILE__, __LINE
294 }
295 aHadron.Lorentz(aHadron, -1. * (*fCache.Get(
296 // G4cout << aHadron.GetMomentum()<<" ";
297 // G4cout << aHadron.GetTotalMomentum()<<G4
298 }
299