| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 25-08-06 New Final State type (refFlag==3 , << 30 // 25-08-06 New Final State type (refFlag==3 , Legendre (Low Energy) + Probability (High Energy) )
31 // is added by T. KOI 31 // is added by T. KOI
32 // 080904 Add Protection for negative energy r << 32 // 080904 Add Protection for negative energy results in very low energy ( 1E-6 eV ) scattering by T. Koi
33 // Koi <<
34 // 33 //
35 // P. Arce, June-2014 Conversion neutron_hp to 34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 // 35 //
37 #include "G4ParticleHPElasticFS.hh" 36 #include "G4ParticleHPElasticFS.hh"
38 <<
39 #include "G4Alpha.hh" <<
40 #include "G4Deuteron.hh" <<
41 #include "G4HadronicParameters.hh" <<
42 #include "G4IonTable.hh" <<
43 #include "G4LorentzVector.hh" <<
44 #include "G4Nucleus.hh" <<
45 #include "G4ParticleHPDataUsed.hh" <<
46 #include "G4ParticleHPManager.hh" 37 #include "G4ParticleHPManager.hh"
>> 38
47 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
48 #include "G4PhysicsModelCatalog.hh" <<
49 #include "G4Pow.hh" <<
50 #include "G4Proton.hh" <<
51 #include "G4ReactionProduct.hh" <<
52 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
53 #include "G4ThreeVector.hh" << 41 #include "G4ReactionProduct.hh"
>> 42 #include "G4Nucleus.hh"
>> 43 #include "G4Proton.hh"
>> 44 #include "G4Deuteron.hh"
54 #include "G4Triton.hh" 45 #include "G4Triton.hh"
>> 46 #include "G4Alpha.hh"
>> 47 #include "G4ThreeVector.hh"
>> 48 #include "G4LorentzVector.hh"
>> 49 #include "G4IonTable.hh"
>> 50 #include "G4ParticleHPDataUsed.hh"
55 51
56 #include "zlib.h" 52 #include "zlib.h"
57 53
58 G4ParticleHPElasticFS::G4ParticleHPElasticFS() << 54 void G4ParticleHPElasticFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String &, G4ParticleDefinition* )
59 { << 55 {
60 svtEmax = 0.0; << 56 G4String tString = "/FS";
61 dbrcEmax = 0.0; << 57 G4bool dbool;
62 dbrcEmin = 0.0; << 58 G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
63 dbrcAmin = 0.0; << 59 G4String filename = aFile.GetName();
64 dbrcUse = false; << 60 SetAZMs( A, Z, M, aFile );
65 xsForDBRC = nullptr; << 61 //theBaseA = aFile.GetA();
66 << 62 //theBaseZ = aFile.GetZ();
67 secID = G4PhysicsModelCatalog::GetModelID("m << 63 if(!dbool)
68 << 64 {
69 hasXsec = false; << 65 hasAnyData = false;
70 theCoefficients = nullptr; << 66 hasFSData = false;
71 theProbArray = nullptr; << 67 hasXsec = false;
72 << 68 return;
73 repFlag = 0; <<
74 tE_of_repFlag3 = 0.0; <<
75 targetMass = 0.0; <<
76 frameFlag = 0; <<
77 } <<
78 <<
79 void G4ParticleHPElasticFS::Init(G4double A, G <<
80 const G4Strin <<
81 G4ParticleDef <<
82 { <<
83 G4String tString = "/FS"; <<
84 G4bool dbool = true; <<
85 SetA_Z(A, Z, M); <<
86 const G4ParticleHPDataUsed& aFile = <<
87 theNames.GetName(theBaseA, theBaseZ, M, di <<
88 const G4String& filename = aFile.GetName(); <<
89 SetAZMs(aFile); <<
90 if (!dbool) { <<
91 hasAnyData = false; <<
92 hasFSData = false; <<
93 hasXsec = false; <<
94 return; <<
95 } <<
96 <<
97 // 130205 For compressed data files <<
98 std::istringstream theData(std::ios::in); <<
99 G4ParticleHPManager::GetInstance()->GetDataS <<
100 // 130205 END <<
101 theData >> repFlag >> targetMass >> frameFla <<
102 <<
103 if (repFlag == 1) { <<
104 G4int nEnergy; <<
105 theData >> nEnergy; <<
106 theCoefficients = new G4ParticleHPLegendre <<
107 theCoefficients->InitInterpolation(theData <<
108 G4double temp, energy; <<
109 G4int tempdep, nLegendre; <<
110 G4int i, ii; <<
111 for (i = 0; i < nEnergy; i++) { <<
112 theData >> temp >> energy >> tempdep >> <<
113 energy *= eV; <<
114 theCoefficients->Init(i, energy, nLegend <<
115 theCoefficients->SetTemperature(i, temp) <<
116 G4double coeff = 0; <<
117 for (ii = 0; ii < nLegendre; ii++) { <<
118 // load legendre coefficients. <<
119 theData >> coeff; <<
120 theCoefficients->SetCoeff(i, ii + 1, c <<
121 } <<
122 } 69 }
123 } << 70 //130205 For compressed data files
124 else if (repFlag == 2) { << 71 std::istringstream theData(std::ios::in);
125 G4int nEnergy; << 72 G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
126 theData >> nEnergy; << 73 //130205 END
127 theProbArray = new G4ParticleHPPartial(nEn << 74 theData >> repFlag >> targetMass >> frameFlag;
128 theProbArray->InitInterpolation(theData); << 75 if(repFlag==1)
129 G4double temp, energy; << 76 {
130 G4int tempdep, nPoints; << 77 G4int nEnergy;
131 for (G4int i = 0; i < nEnergy; i++) { << 78 theData >> nEnergy;
132 theData >> temp >> energy >> tempdep >> << 79 theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
133 energy *= eV; << 80 theCoefficients->InitInterpolation(theData);
134 theProbArray->InitInterpolation(i, theDa << 81 G4double temp, energy;
135 theProbArray->SetT(i, temp); << 82 G4int tempdep, nLegendre;
136 theProbArray->SetX(i, energy); << 83 G4int i, ii;
137 G4double prob, costh; << 84 for (i=0; i<nEnergy; i++)
138 for (G4int ii = 0; ii < nPoints; ii++) { << 85 {
139 // fill probability arrays. << 86 theData >> temp >> energy >> tempdep >> nLegendre;
140 theData >> costh >> prob; << 87 energy *=eV;
141 theProbArray->SetX(i, ii, costh); << 88 theCoefficients->Init(i, energy, nLegendre);
142 theProbArray->SetY(i, ii, prob); << 89 theCoefficients->SetTemperature(i, temp);
>> 90 G4double coeff=0;
>> 91 for(ii=0; ii<nLegendre; ii++)
>> 92 {
>> 93 // load legendre coefficients.
>> 94 theData >> coeff;
>> 95 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
>> 96 }
143 } 97 }
144 theProbArray->DoneSetXY(i); <<
145 } 98 }
146 } << 99 else if (repFlag==2)
147 else if (repFlag == 3) { << 100 {
148 G4int nEnergy_Legendre; << 101 G4int nEnergy;
149 theData >> nEnergy_Legendre; << 102 theData >> nEnergy;
150 if (nEnergy_Legendre <= 0) { << 103 theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
151 std::stringstream iss; << 104 theProbArray->InitInterpolation(theData);
152 iss << "G4ParticleHPElasticFS::Init Data << 105 G4double temp, energy;
153 iss << "Z, A and M of problematic file i << 106 G4int tempdep, nPoints;
154 << theNDLDataM << " respectively."; << 107 for(G4int i=0; i<nEnergy; i++)
155 throw G4HadronicException(__FILE__, __LI << 108 {
156 } << 109 theData >> temp >> energy >> tempdep >> nPoints;
157 theCoefficients = new G4ParticleHPLegendre << 110 energy *= eV;
158 theCoefficients->InitInterpolation(theData << 111 theProbArray->InitInterpolation(i, theData);
159 G4double temp, energy; << 112 theProbArray->SetT(i, temp);
160 G4int tempdep, nLegendre; << 113 theProbArray->SetX(i, energy);
161 << 114 G4double prob, costh;
162 for (G4int i = 0; i < nEnergy_Legendre; i+ << 115 for(G4int ii=0; ii<nPoints; ii++)
163 theData >> temp >> energy >> tempdep >> << 116 {
164 energy *= eV; << 117 // fill probability arrays.
165 theCoefficients->Init(i, energy, nLegend << 118 theData >> costh >> prob;
166 theCoefficients->SetTemperature(i, temp) << 119 theProbArray->SetX(i, ii, costh);
167 G4double coeff = 0; << 120 theProbArray->SetY(i, ii, prob);
168 for (G4int ii = 0; ii < nLegendre; ii++) << 121 }
169 // load legendre coefficients. <<
170 theData >> coeff; <<
171 theCoefficients->SetCoeff(i, ii + 1, c <<
172 } 122 }
173 } 123 }
174 << 124 else if ( repFlag==3 )
175 tE_of_repFlag3 = energy; << 125 {
176 << 126 G4int nEnergy_Legendre;
177 G4int nEnergy_Prob; << 127 theData >> nEnergy_Legendre;
178 theData >> nEnergy_Prob; << 128 theCoefficients = new G4ParticleHPLegendreStore( nEnergy_Legendre );
179 theProbArray = new G4ParticleHPPartial(nEn << 129 theCoefficients->InitInterpolation( theData );
180 theProbArray->InitInterpolation(theData); << 130 G4double temp, energy;
181 G4int nPoints; << 131 G4int tempdep, nLegendre;
182 for (G4int i = 0; i < nEnergy_Prob; i++) { << 132 //G4int i, ii;
183 theData >> temp >> energy >> tempdep >> << 133 for ( G4int i = 0 ; i < nEnergy_Legendre ; i++ )
184 energy *= eV; << 134 {
185 << 135 theData >> temp >> energy >> tempdep >> nLegendre;
186 // consistency check << 136 energy *=eV;
187 if (i == 0) << 137 theCoefficients->Init( i , energy , nLegendre );
188 // if ( energy != tE_of_repFlag3 ) //1 << 138 theCoefficients->SetTemperature( i , temp );
189 if (std::abs(energy - tE_of_repFlag3) << 139 G4double coeff = 0;
190 G4cout << "Warning Transition Energy << 140 for (G4int ii = 0 ; ii < nLegendre ; ii++ )
191 << 141 {
192 theProbArray->InitInterpolation(i, theDa << 142 // load legendre coefficients.
193 theProbArray->SetT(i, temp); << 143 theData >> coeff;
194 theProbArray->SetX(i, energy); << 144 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
195 G4double prob, costh; << 145 }
196 for (G4int ii = 0; ii < nPoints; ii++) { << 146 }
197 // fill probability arrays. << 147
198 theData >> costh >> prob; << 148 tE_of_repFlag3 = energy;
199 theProbArray->SetX(i, ii, costh); << 149
200 theProbArray->SetY(i, ii, prob); << 150 G4int nEnergy_Prob;
201 } << 151 theData >> nEnergy_Prob;
202 theProbArray->DoneSetXY(i); << 152 theProbArray = new G4ParticleHPPartial( nEnergy_Prob , nEnergy_Prob );
>> 153 theProbArray->InitInterpolation( theData );
>> 154 G4int nPoints;
>> 155 for ( G4int i=0 ; i < nEnergy_Prob ; i++ )
>> 156 {
>> 157 theData >> temp >> energy >> tempdep >> nPoints;
>> 158
>> 159 energy *= eV;
>> 160
>> 161 // consistency check
>> 162 if ( i == 0 )
>> 163 //if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines
>> 164 if ( std::abs( energy - tE_of_repFlag3 ) / tE_of_repFlag3 > 1.0e-15 )
>> 165 G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl;
>> 166
>> 167 theProbArray->InitInterpolation( i , theData );
>> 168 theProbArray->SetT( i , temp );
>> 169 theProbArray->SetX( i , energy );
>> 170 G4double prob, costh;
>> 171 for( G4int ii = 0 ; ii < nPoints ; ii++ )
>> 172 {
>> 173 // fill probability arrays.
>> 174 theData >> costh >> prob;
>> 175 theProbArray->SetX( i , ii , costh );
>> 176 theProbArray->SetY( i , ii , prob );
>> 177 }
>> 178 }
203 } 179 }
>> 180 else if (repFlag==0)
>> 181 {
>> 182 theData >> frameFlag;
>> 183 }
>> 184 else
>> 185 {
>> 186 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
>> 187 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
>> 188 }
>> 189 //130205 For compressed data files(theData changed from ifstream to istringstream)
>> 190 //theData.close();
204 } 191 }
205 else if (repFlag == 0) { << 192 G4HadFinalState * G4ParticleHPElasticFS::ApplyYourself(const G4HadProjectile & theTrack)
206 theData >> frameFlag; << 193 {
207 } << 194 // G4cout << "G4ParticleHPElasticFS::ApplyYourself+"<<G4endl;
208 else { << 195 theResult.Clear();
209 G4cout << "unusable number for repFlag: re << 196 G4double eKinetic = theTrack.GetKineticEnergy();
210 throw G4HadronicException(__FILE__, __LINE << 197 const G4HadProjectile *incidentParticle = &theTrack;
211 "G4ParticleHPEla << 198 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition() ));
212 } << 199 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
213 // 130205 For compressed data files(theData << 200 theNeutron.SetKineticEnergy( eKinetic );
214 // theData.close(); << 201 // G4cout << "G4ParticleHPElasticFS::ApplyYourself++"<<eKinetic<<" "<<G4endl;
215 } << 202 // G4cout << "CMSVALUES 0 "<<theNeutron.GetTotalMomentum()<<G4endl;
216 << 203
217 G4HadFinalState* G4ParticleHPElasticFS::ApplyY << 204 G4ReactionProduct theTarget;
218 { << 205 G4Nucleus aNucleus;
219 if (theResult.Get() == nullptr) theResult.Pu << 206 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
220 theResult.Get()->Clear(); << 207 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
221 G4double eKinetic = theTrack.GetKineticEnerg << 208 //t theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) ); //TESTPHP
222 const G4HadProjectile* incidentParticle = &t << 209 // G4cout << "Nucleus-test"<<" "<<targetMass<<" ";
223 G4ReactionProduct theNeutron( << 210 // G4cout << theTarget.GetMomentum().x()<<" ";
224 const_cast<G4ParticleDefinition*>(incident << 211 // G4cout << theTarget.GetMomentum().y()<<" ";
225 theNeutron.SetMomentum(incidentParticle->Get << 212 // G4cout << theTarget.GetMomentum().z()<<G4endl;
226 theNeutron.SetKineticEnergy(eKinetic); << 213
227 << 214 // neutron and target defined as reaction products.
228 G4ThreeVector neuVelo = << 215
229 (1. / incidentParticle->GetDefinition()->G << 216 // prepare lorentz-transformation to Lab.
230 G4ReactionProduct theTarget = << 217
231 GetBiasedThermalNucleus(targetMass, neuVel << 218 G4ThreeVector the3Neutron = theNeutron.GetMomentum();
232 << 219 G4double nEnergy = theNeutron.GetTotalEnergy();
233 // Neutron and target defined as G4ReactionP << 220 G4ThreeVector the3Target = theTarget.GetMomentum();
234 // Prepare Lorentz transformation to lab << 221 // cout << "@@@" << the3Target<<G4endl;
235 << 222 G4double tEnergy = theTarget.GetTotalEnergy();
236 G4ThreeVector the3Neutron = theNeutron.GetMo << 223 G4ReactionProduct theCMS;
237 G4double nEnergy = theNeutron.GetTotalEnergy << 224 G4double totE = nEnergy+tEnergy;
238 G4ThreeVector the3Target = theTarget.GetMome << 225 G4ThreeVector the3CMS = the3Target+the3Neutron;
239 G4double tEnergy = theTarget.GetTotalEnergy( << 226 theCMS.SetMomentum(the3CMS);
240 G4ReactionProduct theCMS; << 227 G4double cmsMom = std::sqrt(the3CMS*the3CMS);
241 G4double totE = nEnergy + tEnergy; << 228 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
242 G4ThreeVector the3CMS = the3Target + the3Neu << 229 theCMS.SetMass(sqrts);
243 theCMS.SetMomentum(the3CMS); << 230 theCMS.SetTotalEnergy(totE);
244 G4double cmsMom = std::sqrt(the3CMS * the3CM << 231
245 G4double sqrts = std::sqrt((totE - cmsMom) * << 232 // data come as fcn of n-energy in nuclear rest frame
246 theCMS.SetMass(sqrts); << 233 G4ReactionProduct boosted;
247 theCMS.SetTotalEnergy(totE); << 234 boosted.Lorentz(theNeutron, theTarget);
248 << 235 eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
249 // Data come as function of n-energy in nucl << 236 G4double cosTh = -2;
250 G4ReactionProduct boosted; << 237 if(repFlag == 1)
251 boosted.Lorentz(theNeutron, theTarget); << 238 {
252 eKinetic = boosted.GetKineticEnergy(); // g <<
253 G4double cosTh = -2; <<
254 <<
255 if (repFlag == 1) { <<
256 cosTh = theCoefficients->SampleElastic(eKi <<
257 } <<
258 else if (repFlag == 2) { <<
259 cosTh = theProbArray->Sample(eKinetic); <<
260 } <<
261 else if (repFlag == 3) { <<
262 if (eKinetic <= tE_of_repFlag3) { <<
263 cosTh = theCoefficients->SampleElastic(e 239 cosTh = theCoefficients->SampleElastic(eKinetic);
264 } 240 }
265 else { << 241
>> 242 else if (repFlag==2)
>> 243 {
266 cosTh = theProbArray->Sample(eKinetic); 244 cosTh = theProbArray->Sample(eKinetic);
267 } 245 }
268 } << 246 else if (repFlag==3)
269 else if (repFlag == 0) { << 247 {
270 cosTh = 2. * G4UniformRand() - 1.; << 248 if ( eKinetic <= tE_of_repFlag3 )
271 } << 249 {
272 else { << 250 cosTh = theCoefficients->SampleElastic(eKinetic);
273 G4cout << "Unusable number for repFlag: re << 251 }
274 throw G4HadronicException(__FILE__, __LINE << 252 else
275 "G4ParticleHPEla << 253 {
276 } << 254 cosTh = theProbArray->Sample(eKinetic);
277 << 255 }
278 if (cosTh < -1.1) { <<
279 return nullptr; <<
280 } <<
281 <<
282 G4double phi = twopi * G4UniformRand(); <<
283 G4double cosPhi = std::cos(phi); <<
284 G4double sinPhi = std::sin(phi); <<
285 G4double theta = std::acos(cosTh); <<
286 G4double sinth = std::sin(theta); <<
287 <<
288 if (frameFlag == 1) { <<
289 // Projectile scattering values cosTh are <<
290 // In this frame, do relativistic calculat <<
291 // target 4-momenta <<
292 <<
293 theNeutron.Lorentz(theNeutron, theTarget); <<
294 G4double mN = theNeutron.GetMass(); <<
295 G4double Pinit = theNeutron.GetTotalMoment <<
296 G4double Einit = theNeutron.GetTotalEnergy <<
297 G4double mT = theTarget.GetMass(); <<
298 <<
299 G4double ratio = mT / mN; <<
300 G4double sqt = std::sqrt(ratio * ratio - 1 <<
301 G4double beta = Pinit / (mT + Einit); // <<
302 G4double denom = 1. - beta * beta * cosTh <<
303 G4double term1 = cosTh * (Einit * ratio + <<
304 G4double pN = beta * mN * (term1 + sqt) / <<
305 <<
306 // Get the scattered momentum and rotate i <<
307 G4ThreeVector pDir = theNeutron.GetMomentu <<
308 G4double px = pN * pDir.x(); <<
309 G4double py = pN * pDir.y(); <<
310 G4double pz = pN * pDir.z(); <<
311 <<
312 G4ThreeVector pcmRot; <<
313 pcmRot.setX(px * cosTh * cosPhi - py * sin <<
314 pcmRot.setY(px * cosTh * sinPhi + py * cos <<
315 pcmRot.setZ(-px * sinth + pz * cosTh); <<
316 theNeutron.SetMomentum(pcmRot); <<
317 G4double eN = std::sqrt(pN * pN + mN * mN) <<
318 theNeutron.SetTotalEnergy(eN); <<
319 <<
320 // Get the scattered target momentum <<
321 G4ReactionProduct toLab(-1. * theTarget); <<
322 theTarget.SetMomentum(pDir * Pinit - pcmRo <<
323 G4double eT = Einit - eN + mT; <<
324 theTarget.SetTotalEnergy(eT); <<
325 <<
326 // Now back to lab frame <<
327 theNeutron.Lorentz(theNeutron, toLab); <<
328 theTarget.Lorentz(theTarget, toLab); <<
329 <<
330 // 111005 Protection for not producing 0 k <<
331 if (theNeutron.GetKineticEnergy() <= 0) <<
332 theNeutron.SetTotalEnergy(theNeutron.Get <<
333 * (1. + G4Pow: <<
334 if (theTarget.GetKineticEnergy() <= 0) <<
335 theTarget.SetTotalEnergy(theTarget.GetMa <<
336 } <<
337 else if (frameFlag == 2) { <<
338 // Projectile scattering values cosTh take <<
339 <<
340 G4LorentzVector proj(nEnergy, the3Neutron) <<
341 G4LorentzVector targ(tEnergy, the3Target); <<
342 G4ThreeVector boostToCM = proj.findBoostTo <<
343 proj.boost(boostToCM); <<
344 targ.boost(boostToCM); <<
345 <<
346 // Rotate projectile and target momenta by <<
347 // Note: at this point collision axis is n <<
348 // momentum given target nucleus by <<
349 G4double px = proj.px(); <<
350 G4double py = proj.py(); <<
351 G4double pz = proj.pz(); <<
352 <<
353 G4ThreeVector pcmRot; <<
354 pcmRot.setX(px * cosTh * cosPhi - py * sin <<
355 pcmRot.setY(px * cosTh * sinPhi + py * cos <<
356 pcmRot.setZ(-px * sinth + pz * cosTh); <<
357 proj.setVect(pcmRot); <<
358 targ.setVect(-pcmRot); <<
359 <<
360 // Back to lab frame <<
361 proj.boost(-boostToCM); <<
362 targ.boost(-boostToCM); <<
363 <<
364 theNeutron.SetMomentum(proj.vect()); <<
365 theNeutron.SetTotalEnergy(proj.e()); <<
366 <<
367 theTarget.SetMomentum(targ.vect()); <<
368 theTarget.SetTotalEnergy(targ.e()); <<
369 <<
370 // 080904 Add Protection for very low ener <<
371 if (theNeutron.GetKineticEnergy() <= 0) { <<
372 theNeutron.SetTotalEnergy(theNeutron.Get <<
373 * (1. + G4Pow: <<
374 } <<
375 <<
376 // 080904 Add Protection for very low ener <<
377 if (theTarget.GetKineticEnergy() <= 0) { <<
378 theTarget.SetTotalEnergy(theTarget.GetMa <<
379 } 256 }
380 } << 257 else if (repFlag==0)
381 else { << 258 {
382 G4cout << "Value of frameFlag (1=LAB, 2=CM << 259 cosTh = 2.*G4UniformRand()-1.;
383 throw G4HadronicException(__FILE__, __LINE << 260 }
384 "G4ParticleHPEla << 261 else
385 } << 262 {
386 << 263 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
387 // Everything is now in the lab frame << 264 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
388 // Set energy change and momentum change << 265 }
389 theResult.Get()->SetEnergyChange(theNeutron. << 266 if(cosTh<-1.1) { return 0; }
390 theResult.Get()->SetMomentumChange(theNeutro << 267 G4double phi = twopi*G4UniformRand();
391 << 268 G4double theta = std::acos(cosTh);
392 // Make recoil a G4DynamicParticle << 269 G4double sinth = std::sin(theta);
393 auto theRecoil = new G4DynamicParticle; << 270 if (frameFlag == 1) // final state data given in target rest frame.
394 theRecoil->SetDefinition(G4IonTable::GetIonT << 271 {
395 << 272 // we have the scattering angle, now we need the energy, then do the
396 theRecoil->SetMomentum(theTarget.GetMomentum << 273 // boosting.
397 theResult.Get()->AddSecondary(theRecoil, sec << 274 // relativistic elastic scattering energy angular correlation:
398 << 275 theNeutron.Lorentz(theNeutron, theTarget);
399 // Postpone the tracking of the primary neut << 276 G4double e0 = theNeutron.GetTotalEnergy();
400 theResult.Get()->SetStatusChange(suspend); << 277 G4double p0 = theNeutron.GetTotalMomentum();
401 return theResult.Get(); << 278 G4double mN = theNeutron.GetMass();
402 } << 279 G4double mT = theTarget.GetMass();
403 << 280 G4double eE = e0+mT;
404 void G4ParticleHPElasticFS::InitializeScatteri << 281 G4double ap = (mT+eE)*(mT-eE) + (p0+mN)*(p0-mN);
>> 282 G4double a = 4*(eE+p0*cosTh)*(eE-p0*cosTh);
>> 283 G4double b = 4*ap*p0*cosTh;
>> 284 G4double c = (2.*eE*mN-ap)*(2.*eE*mN+ap);
>> 285 G4double en = (-b+std::sqrt(b*b - 4*a*c) )/(2*a);
>> 286 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
>> 287 theNeutron.SetMomentum(tempVector);
>> 288 theNeutron.SetTotalEnergy(std::sqrt(en*en+theNeutron.GetMass()*theNeutron.GetMass()));
>> 289 // first to lab
>> 290 theNeutron.Lorentz(theNeutron, -1.*theTarget);
>> 291 // now to CMS
>> 292 theNeutron.Lorentz(theNeutron, theCMS);
>> 293 theTarget.SetMomentum(-theNeutron.GetMomentum());
>> 294 theTarget.SetTotalEnergy(theNeutron.GetTotalEnergy());
>> 295 // and back to lab
>> 296 theNeutron.Lorentz(theNeutron, -1.*theCMS);
>> 297 theTarget.Lorentz(theTarget, -1.*theCMS);
>> 298 //111005 Protection for not producing 0 kinetic energy target
>> 299 if ( theNeutron.GetKineticEnergy() <= 0 ) theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
>> 300 if ( theTarget.GetKineticEnergy() <= 0 ) theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
>> 301 }
>> 302 else if (frameFlag == 2) // CMS
>> 303 {
>> 304 theNeutron.Lorentz(theNeutron, theCMS);
>> 305 theTarget.Lorentz(theTarget, theCMS);
>> 306 G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
>> 307 G4ThreeVector cmsMom_tmp=theNeutron.GetMomentum(); // for neutron direction in CMS
>> 308 G4double cms_theta=cmsMom_tmp.theta();
>> 309 G4double cms_phi=cmsMom_tmp.phi();
>> 310 G4ThreeVector tempVector;
>> 311 tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
>> 312 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
>> 313 -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
>> 314 tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
>> 315 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
>> 316 +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
>> 317 tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
>> 318 -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
>> 319 tempVector *= en;
>> 320 theNeutron.SetMomentum(tempVector);
>> 321 theTarget.SetMomentum(-tempVector);
>> 322 G4double tP = theTarget.GetTotalMomentum();
>> 323 G4double tM = theTarget.GetMass();
>> 324 theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
>> 325
>> 326 /*
>> 327 For debug purpose.
>> 328 Same transformation G4ReactionProduct.Lorentz() by 4vectors
405 { 329 {
406 // Initialize DBRC variables << 330 G4LorentzVector n4p = G4LorentzVector ( theNeutron.GetMomentum() , theNeutron.GetKineticEnergy() + theNeutron.GetMass() );
407 svtEmax = G4HadronicParameters::Instance()-> << 331 G4cout << "before " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
408 G4ParticleHPManager* manager = G4ParticleHPM << 332 G4LorentzVector cm4p = G4LorentzVector ( theCMS.GetMomentum() , theCMS.GetKineticEnergy() + theCMS.GetMass() );
409 dbrcUse = manager->GetUseDBRC(); << 333 n4p.boost( cm4p.boostVector() );
410 dbrcEmax = manager->GetMaxEnergyDBRC(); << 334 G4cout << cm4p/eV << G4endl;
411 dbrcEmin = manager->GetMinEnergyDBRC(); << 335 G4cout << "after " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
412 dbrcAmin = manager->GetMinADBRC(); <<
413 } 336 }
>> 337 */
414 338
415 G4ReactionProduct G4ParticleHPElasticFS::GetBi << 339 theNeutron.Lorentz(theNeutron, -1.*theCMS);
416 << 340 //080904 Add Protection for very low energy (1e-6eV) scattering
417 << 341 if ( theNeutron.GetKineticEnergy() <= 0 )
418 { << 342 {
419 // This new method implements the DBRC (Dopp << 343 //theNeutron.SetMomentum( G4ThreeVector(0) );
420 // on top of the SVT (Sampling of the Veloci << 344 //theNeutron.SetTotalEnergy ( theNeutron.GetMass() );
421 // The SVT algorithm was written by Loic Thu << 345 //110822 Protection for not producing 0 kinetic energy neutron
422 // the method G4Nucleus::GetBiasedThermalNuc << 346 theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
423 // implemented the DBRC algorithm on top of << 347 }
424 // While the SVT algorithm is still present <<
425 // the DBRC algorithm on top of the SVT one <<
426 // order to avoid a cycle dependency between <<
427 <<
428 InitializeScatteringKernelParameters(); <<
429 <<
430 // Set threshold for SVT algorithm <<
431 G4double E_threshold = svtEmax; <<
432 if (svtEmax == -1.) { <<
433 // If E_neutron <= 400*kB*T (400 is a comm <<
434 // then apply the Sampling ot the Velocity <<
435 // else consider the target nucleus being <<
436 E_threshold = 400.0 * 8.617333262E-11 * te <<
437 } <<
438 <<
439 // If DBRC is enabled and the nucleus is hea <<
440 if (dbrcUse && aMass >= dbrcAmin) { <<
441 E_threshold = std::max(svtEmax, dbrcEmax); <<
442 } <<
443 <<
444 G4double E_neutron = 0.5 * aVelocity.mag2() <<
445 <<
446 G4bool dbrcIsOn = dbrcUse && E_neutron >= db <<
447 <<
448 G4Nucleus aNucleus; <<
449 if (E_neutron > E_threshold || !dbrcIsOn) { <<
450 // Apply only the SVT algorithm, not the D <<
451 return aNucleus.GetBiasedThermalNucleus(ta <<
452 } <<
453 <<
454 G4ReactionProduct result; <<
455 result.SetMass(aMass * G4Neutron::Neutron()- <<
456 348
457 // Beta = sqrt(m/2kT) << 349 theTarget.Lorentz(theTarget, -1.*theCMS);
458 G4double beta = << 350 //080904 Add Protection for very low energy (1e-6eV) scattering
459 std::sqrt(result.GetMass() << 351 if ( theTarget.GetKineticEnergy() < 0 )
460 / (2. * 8.617333262E-11 * temp)) << 352 {
461 << 353 //theTarget.SetMomentum( G4ThreeVector(0) );
462 // Neutron speed vn << 354 //theTarget.SetTotalEnergy ( theTarget.GetMass() );
463 G4double vN_norm = aVelocity.mag(); << 355 //110822 Protection for not producing 0 kinetic energy target
464 G4double vN_norm2 = vN_norm * vN_norm; << 356 theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
465 G4double y = beta * vN_norm; <<
466 <<
467 // Normalize neutron velocity <<
468 aVelocity = (1. / vN_norm) * aVelocity; <<
469 <<
470 // Variables for sampling of target speed an <<
471 G4double x2; <<
472 G4double randThresholdSVT; <<
473 G4double vT_norm, vT_norm2, mu; <<
474 G4double acceptThresholdSVT; <<
475 G4double vRelativeSpeed; <<
476 G4double cdf0 = 2. / (2. + std::sqrt(CLHEP:: <<
477 <<
478 // DBRC variables <<
479 G4double xsRelative = -99.; <<
480 G4double randThresholdDBRC = 0.; <<
481 // Calculate max cross-section in interval f <<
482 G4double eMin = <<
483 0.5 * G4Neutron::Neutron()->GetPDGMass() * <<
484 G4double eMax = <<
485 0.5 * G4Neutron::Neutron()->GetPDGMass() * <<
486 G4double xsMax = xsForDBRC->GetMaxY(eMin, eM <<
487 <<
488 do { <<
489 do { <<
490 // Sample the target velocity vT in the <<
491 if (G4UniformRand() < cdf0) { <<
492 // Sample in C45 from https://laws.lan <<
493 x2 = -std::log(G4UniformRand() * G4Uni <<
494 } 357 }
495 else { << 358 }
496 // Sample in C61 from https://laws.lan << 359 else
497 G4double ampl = std::cos(CLHEP::pi / 2 << 360 {
498 x2 = -std::log(G4UniformRand()) - std: << 361 G4cout <<"Value of frameFlag (1=LAB, 2=CMS): "<<frameFlag;
>> 362 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::ApplyYourSelf frameflag incorrect");
>> 363 }
>> 364 // now all in Lab
>> 365 // nun den recoil generieren...und energy change, momentum change angeben.
>> 366 theResult.SetEnergyChange(theNeutron.GetKineticEnergy());
>> 367 theResult.SetMomentumChange(theNeutron.GetMomentum().unit());
>> 368 G4DynamicParticle* theRecoil = new G4DynamicParticle;
>> 369 if(targetMass<4.5)
>> 370 {
>> 371 if(targetMass<1)
>> 372 {
>> 373 // proton
>> 374 theRecoil->SetDefinition(G4Proton::Proton());
499 } 375 }
500 << 376 else if(targetMass<2 )
501 vT_norm = std::sqrt(x2) / beta; << 377 {
502 vT_norm2 = vT_norm * vT_norm; << 378 // deuteron
503 << 379 theRecoil->SetDefinition(G4Deuteron::Deuteron());
504 // Sample cosine between the incident ne << 380 }
505 mu = 2. * G4UniformRand() - 1.; << 381 else if(targetMass<2.999 )
506 << 382 {
507 // Define acceptance threshold for SVT << 383 // 3He
508 vRelativeSpeed = std::sqrt(vN_norm2 + vT << 384 theRecoil->SetDefinition(G4He3::He3());
509 acceptThresholdSVT = vRelativeSpeed / (v << 385 }
510 randThresholdSVT = G4UniformRand(); << 386 else if(targetMass<3 )
511 } while (randThresholdSVT >= acceptThresho << 387 {
512 << 388 // Triton
513 // Apply DBRC rejection << 389 theRecoil->SetDefinition(G4Triton::Triton());
514 xsRelative = xsForDBRC->GetXsec(0.5 * G4Ne << 390 }
515 * vRelativ << 391 else
516 randThresholdDBRC = G4UniformRand(); << 392 {
517 << 393 // alpha
518 } while (randThresholdDBRC >= xsRelative / x << 394 theRecoil->SetDefinition(G4Alpha::Alpha());
519 << 395 }
520 aNucleus.DoKinematicsOfThermalNucleus(mu, vT << 396 }
521 << 397 else
522 return result; << 398 {
523 } << 399 theRecoil->SetDefinition(G4IonTable::GetIonTable()
>> 400 ->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0 ));
>> 401 }
>> 402 theRecoil->SetMomentum(theTarget.GetMomentum());
>> 403 theResult.AddSecondary(theRecoil);
>> 404 // G4cout << "G4ParticleHPElasticFS::ApplyYourself 10+"<<G4endl;
>> 405 // postpone the tracking of the primary neutron
>> 406 theResult.SetStatusChange(suspend);
>> 407 return &theResult;
>> 408 }
524 409