| Geant4 Cross Reference |
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 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 //
37 #include "G4ParticleHPAngular.hh"
38
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41
42 void G4ParticleHPAngular::Init(std::istream& aDataFile)
43 {
44 // G4cout << "here we are entering the Angular Init"<<G4endl;
45 aDataFile >> theAngularDistributionType >> targetMass;
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 G4ParticleHPLegendreStore(nEnergy);
55 theCoefficients->InitInterpolation(aDataFile);
56 G4double temp, energy;
57 G4int tempdep, nLegendre;
58 G4int i, ii;
59 for (i = 0; i < nEnergy; i++) {
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 aDataFile >> coeff;
67 theCoefficients->SetCoeff(i, ii + 1, coeff);
68 }
69 }
70 }
71 else if (theAngularDistributionType == 2) {
72 theIsoFlag = false;
73 G4int nEnergy;
74 aDataFile >> nEnergy;
75 theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
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 Angular: " << theAngularDistributionType << G4endl;
90 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
91 }
92 }
93
94 void G4ParticleHPAngular::SampleAndUpdate(G4ReactionProduct& aHadron)
95 {
96 //********************************************************************
97 // EMendoza -> sampling can be isotropic in LAB or in CMS
98 /*
99 if(theIsoFlag)
100 {
101 // G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
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), en*sinth*std::sin(phi), en*std::cos(theta) );
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().theProjectileRP, *fCache.Get().theTarget);
121 G4double kineticEnergy = boosted.GetKineticEnergy();
122 G4double cosTh = 0.0;
123 //********************************************************************
124 // EMendoza --> sampling can be also isotropic
125 /*
126 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
127 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
128 */
129 //********************************************************************
130 if (theIsoFlag) {
131 cosTh = 2. * G4UniformRand() - 1;
132 }
133 else if (theAngularDistributionType == 1) {
134 cosTh = theCoefficients->SampleMax(kineticEnergy);
135 }
136 else if (theAngularDistributionType == 2) {
137 cosTh = theProbArray->Sample(kineticEnergy);
138 }
139 else {
140 G4cout << "unknown distribution found for Angular: " << theAngularDistributionType << G4endl;
141 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
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(phi), en * sinth * std::sin(phi),
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().theProjectileRP, *fCache.Get().theTarget);
155 G4double kineticEnergy = boostedN.GetKineticEnergy();
156
157 G4double cosTh = 0.0;
158 //********************************************************************
159 // EMendoza --> sampling can be also isotropic
160 /*
161 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
162 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
163 */
164 //********************************************************************
165 if (theIsoFlag) {
166 cosTh = 2. * G4UniformRand() - 1;
167 }
168 else if (theAngularDistributionType == 1) {
169 cosTh = theCoefficients->SampleMax(kineticEnergy);
170 }
171 else if (theAngularDistributionType == 2) {
172 cosTh = theProbArray->Sample(kineticEnergy);
173 }
174 else {
175 G4cout << "unknown distribution found for Angular: " << theAngularDistributionType << G4endl;
176 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
177 }
178 //********************************************************************
179 // 080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern)
180 /*
181 if(theAngularDistributionType == 1) // LAB
182 {
183 G4double en = aHadron.GetTotalMomentum();
184 G4ReactionProduct boosted;
185 boosted.Lorentz(theProjectile, theTarget);
186 G4double kineticEnergy = boosted.GetKineticEnergy();
187 G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
188 G4double theta = std::acos(cosTh);
189 G4double phi = twopi*G4UniformRand();
190 G4double sinth = std::sin(theta);
191 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
192 aHadron.SetMomentum( temp );
193 }
194 else if(theAngularDistributionType == 2) // costh in CMS {
195 }
196 */
197
198 // G4ReactionProduct boostedN;
199 // boostedN.Lorentz(theProjectile, theTarget);
200 // G4double kineticEnergy = boostedN.GetKineticEnergy();
201 // G4double cosTh = theProbArray->Sample(kineticEnergy);
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), sinth * std::sin(phi), std::cos(theta)); // CMS
208
209 // 080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
210 /*
211 G4double en = aHadron.GetTotalEnergy(); // Target rest
212
213 // get trafo from Target rest frame to CMS
214 G4ReactionProduct boostedT;
215 boostedT.Lorentz(theTarget, theTarget);
216
217 G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
218 G4double nEnergy = boostedN.GetTotalEnergy();
219 G4ThreeVector the3Target = boostedT.GetMomentum();
220 G4double tEnergy = boostedT.GetTotalEnergy();
221 G4double totE = nEnergy+tEnergy;
222 G4ThreeVector the3trafo = -the3Target-the3IncidentParticle;
223 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
224 trafo.SetMomentum(the3trafo);
225 G4double cmsMom = std::sqrt(the3trafo*the3trafo);
226 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
227 trafo.SetMass(sqrts);
228 trafo.SetTotalEnergy(totE);
229
230 G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
231 G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
232 G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
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*aHadron.GetMass()*aHadron.GetMass())
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 all in CMS
251 aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
252 aHadron.Lorentz(aHadron, trafo); // now in target rest frame
253 */
254 // Determination of the hadron kinetic energy in CMS
255 // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
256 // kineticEnergy is incident neutron kinetic energy in CMS (or target frame)
257 G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
258 G4double A1 = fCache.Get().theTarget->GetMass() / boostedN.GetMass();
259 G4double A1prim = aHadron.GetMass() / boostedN.GetMass();
260 G4double kinE =
261 (A1 + 1 - A1prim) / (A1 + 1) / (A1 + 1) * (A1 * kineticEnergy + (1 + A1) * QValue);
262 G4double totalE = kinE + aHadron.GetMass();
263 G4double mom2 = totalE * totalE - aHadron.GetMass() * aHadron.GetMass();
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 momentum in CMS
271 aHadron.SetKineticEnergy(kinE); // Set kinetic energy in CMS
272
273 // get trafo from Target rest frame to CMS
274 G4ReactionProduct boostedT;
275 boostedT.Lorentz(*fCache.Get().theTarget, *fCache.Get().theTarget);
276
277 G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
278 G4double nEnergy = boostedN.GetTotalEnergy();
279 G4ThreeVector the3Target = boostedT.GetMomentum();
280 G4double tEnergy = boostedT.GetTotalEnergy();
281 G4double totE = nEnergy + tEnergy;
282 G4ThreeVector the3trafo = -the3Target - the3IncidentParticle;
283 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
284 trafo.SetMomentum(the3trafo);
285 G4double cmsMom = std::sqrt(the3trafo * the3trafo);
286 G4double sqrts = std::sqrt((totE - cmsMom) * (totE + cmsMom));
287 trafo.SetMass(sqrts);
288 trafo.SetTotalEnergy(totE);
289
290 aHadron.Lorentz(aHadron, trafo);
291 }
292 else {
293 throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
294 }
295 aHadron.Lorentz(aHadron, -1. * (*fCache.Get().theTarget));
296 // G4cout << aHadron.GetMomentum()<<" ";
297 // G4cout << aHadron.GetTotalMomentum()<<G4endl;
298 }
299