ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPElasticFS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPElasticFS.cc
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 // 25-08-06 New Final State type (refFlag==3 , Legendre (Low Energy) + Probability (High Energy) )
31 // is added by T. KOI
32 // 080904 Add Protection for negative energy results in very low energy ( 1E-6 eV ) scattering by T. Koi
33 //
34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
35 //
36 #include "G4ParticleHPElasticFS.hh"
37 #include "G4ParticleHPManager.hh"
38 
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ReactionProduct.hh"
42 #include "G4Nucleus.hh"
43 #include "G4Proton.hh"
44 #include "G4Deuteron.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"
51 #include "G4Pow.hh"
52 #include "zlib.h"
53 
55  G4String& dirName, G4String&,
57 {
58  G4String tString = "/FS";
59  G4bool dbool;
60  G4ParticleHPDataUsed aFile =
61  theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
62  G4String filename = aFile.GetName();
63  SetAZMs( A, Z, M, aFile );
64  //theBaseA = aFile.GetA();
65  //theBaseZ = aFile.GetZ();
66  if (!dbool) {
67  hasAnyData = false;
68  hasFSData = false;
69  hasXsec = false;
70  return;
71  }
72 
73  //130205 For compressed data files
74  std::istringstream theData(std::ios::in);
76  //130205 END
77  theData >> repFlag >> targetMass >> frameFlag;
78 
79  if (repFlag == 1) {
80  G4int nEnergy;
81  theData >> nEnergy;
84  G4double temp, energy;
85  G4int tempdep, nLegendre;
86  G4int i, ii;
87  for (i=0; i < nEnergy; i++) {
88  theData >> temp >> energy >> tempdep >> nLegendre;
89  energy *=eV;
90  theCoefficients->Init(i, energy, nLegendre);
92  G4double coeff = 0;
93  for (ii = 0; ii < nLegendre; ii++) {
94  // load legendre coefficients.
95  theData >> coeff;
96  theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
97  }
98  }
99 
100  } else if (repFlag == 2) {
101  G4int nEnergy;
102  theData >> nEnergy;
103  theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
105  G4double temp, energy;
106  G4int tempdep, nPoints;
107  for (G4int i = 0; i < nEnergy; i++) {
108  theData >> temp >> energy >> tempdep >> nPoints;
109  energy *= eV;
110  theProbArray->InitInterpolation(i, theData);
111  theProbArray->SetT(i, temp);
112  theProbArray->SetX(i, energy);
113  G4double prob, costh;
114  for (G4int ii = 0; ii < nPoints; ii++) {
115  // fill probability arrays.
116  theData >> costh >> prob;
117  theProbArray->SetX(i, ii, costh);
118  theProbArray->SetY(i, ii, prob);
119  }
120  theProbArray->DoneSetXY( i );
121  }
122 
123  } else if (repFlag == 3) {
124  G4int nEnergy_Legendre;
125  theData >> nEnergy_Legendre;
126  if (nEnergy_Legendre <= 0 ) {
127  std::stringstream iss;
128  iss << "G4ParticleHPElasticFS::Init Data Error repFlag is 3 but nEnergy_Legendre <= 0";
129  iss << "Z, A and M of problematic file is " << theNDLDataZ << ", "
130  << theNDLDataA << " and " << theNDLDataM << " respectively.";
131  throw G4HadronicException(__FILE__, __LINE__, iss.str() );
132  }
133  theCoefficients = new G4ParticleHPLegendreStore( nEnergy_Legendre );
135  G4double temp, energy;
136  G4int tempdep, nLegendre;
137 
138  for (G4int i = 0; i < nEnergy_Legendre; i++) {
139  theData >> temp >> energy >> tempdep >> nLegendre;
140  energy *=eV;
141  theCoefficients->Init( i , energy , nLegendre );
142  theCoefficients->SetTemperature( i , temp );
143  G4double coeff = 0;
144  for (G4int ii = 0; ii < nLegendre; ii++) {
145  // load legendre coefficients.
146  theData >> coeff;
147  theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
148  }
149  }
150 
152 
153  G4int nEnergy_Prob;
154  theData >> nEnergy_Prob;
155  theProbArray = new G4ParticleHPPartial( nEnergy_Prob , nEnergy_Prob );
156  theProbArray->InitInterpolation( theData );
157  G4int nPoints;
158  for (G4int i = 0; i < nEnergy_Prob; i++) {
159  theData >> temp >> energy >> tempdep >> nPoints;
160  energy *= eV;
161 
162  // consistency check
163  if (i == 0)
164  //if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines
165  if (std::abs(energy - tE_of_repFlag3) / tE_of_repFlag3 > 1.0e-15)
166  G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl;
167 
168  theProbArray->InitInterpolation( i , theData );
169  theProbArray->SetT( i , temp );
170  theProbArray->SetX( i , energy );
171  G4double prob, costh;
172  for (G4int ii = 0; ii < nPoints; ii++) {
173  // fill probability arrays.
174  theData >> costh >> prob;
175  theProbArray->SetX( i , ii , costh );
176  theProbArray->SetY( i , ii , prob );
177  }
178  theProbArray->DoneSetXY( i );
179  }
180 
181  } else if (repFlag==0) {
182  theData >> frameFlag;
183 
184  } else {
185  G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
186  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
187  }
188  //130205 For compressed data files(theData changed from ifstream to istringstream)
189  //theData.close();
190 }
191 
192 
195 {
196  if (theResult.Get() == NULL) theResult.Put(new G4HadFinalState);
197  theResult.Get()->Clear();
198  G4double eKinetic = theTrack.GetKineticEnergy();
199  const G4HadProjectile *incidentParticle = &theTrack;
200  G4ReactionProduct theNeutron(const_cast<G4ParticleDefinition*>(incidentParticle->GetDefinition() ));
201  theNeutron.SetMomentum(incidentParticle->Get4Momentum().vect() );
202  theNeutron.SetKineticEnergy(eKinetic);
203 
205  G4Nucleus aNucleus;
206  G4ThreeVector neuVelo =
207  (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
208  theTarget =
209  aNucleus.GetBiasedThermalNucleus(targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
210 
211  // Neutron and target defined as G4ReactionProducts
212  // Prepare Lorentz transformation to lab
213 
214  G4ThreeVector the3Neutron = theNeutron.GetMomentum();
215  G4double nEnergy = theNeutron.GetTotalEnergy();
216  G4ThreeVector the3Target = theTarget.GetMomentum();
217  G4double tEnergy = theTarget.GetTotalEnergy();
218  G4ReactionProduct theCMS;
219  G4double totE = nEnergy+tEnergy;
220  G4ThreeVector the3CMS = the3Target+the3Neutron;
221  theCMS.SetMomentum(the3CMS);
222  G4double cmsMom = std::sqrt(the3CMS*the3CMS);
223  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
224  theCMS.SetMass(sqrts);
225  theCMS.SetTotalEnergy(totE);
226 
227  // Data come as function of n-energy in nuclear rest frame
228  G4ReactionProduct boosted;
229  boosted.Lorentz(theNeutron, theTarget);
230  eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
231  G4double cosTh = -2;
232 
233  if (repFlag == 1) {
234  cosTh = theCoefficients->SampleElastic(eKinetic);
235 
236  } else if (repFlag == 2) {
237  cosTh = theProbArray->Sample(eKinetic);
238 
239  } else if (repFlag == 3) {
240  if (eKinetic <= tE_of_repFlag3) {
241  cosTh = theCoefficients->SampleElastic(eKinetic);
242  } else {
243  cosTh = theProbArray->Sample(eKinetic);
244  }
245 
246  } else if (repFlag == 0) {
247  cosTh = 2.*G4UniformRand() - 1.;
248 
249  } else {
250  G4cout << "Unusable number for repFlag: repFlag=" << repFlag << G4endl;
251  throw G4HadronicException(__FILE__, __LINE__,
252  "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
253  }
254 
255  if (cosTh < -1.1) { return 0; }
256 
258  G4double cosPhi = std::cos(phi);
259  G4double sinPhi = std::sin(phi);
260  G4double theta = std::acos(cosTh);
261  G4double sinth = std::sin(theta);
262 
263  if (frameFlag == 1) {
264  // Projectile scattering values cosTh are in target rest frame
265  // In this frame, do relativistic calculation of scattered projectile and
266  // target 4-momenta
267 
268  theNeutron.Lorentz(theNeutron, theTarget);
269  G4double mN = theNeutron.GetMass();
270  G4double Pinit = theNeutron.GetTotalMomentum(); // Incident momentum
271  G4double Einit = theNeutron.GetTotalEnergy(); // Incident energy
272  G4double mT = theTarget.GetMass();
273 
274  G4double ratio = mT/mN;
275  G4double sqt = std::sqrt(ratio*ratio - 1.0 + cosTh*cosTh);
276  G4double beta = Pinit/(mT + Einit); // CMS beta
277  G4double denom = 1. - beta*beta*cosTh*cosTh;
278  G4double term1 = cosTh*(Einit*ratio + mN)/(mN*ratio + Einit);
279  G4double pN = beta*mN*(term1 + sqt)/denom;
280 
281  // Get the scattered momentum and rotate it in theta and phi
282  G4ThreeVector pDir = theNeutron.GetMomentum()/Pinit;
283  G4double px = pN*pDir.x();
284  G4double py = pN*pDir.y();
285  G4double pz = pN*pDir.z();
286 
287  G4ThreeVector pcmRot;
288  pcmRot.setX(px*cosTh*cosPhi - py*sinPhi + pz*sinth*cosPhi);
289  pcmRot.setY(px*cosTh*sinPhi + py*cosPhi + pz*sinth*sinPhi);
290  pcmRot.setZ(-px*sinth + pz*cosTh);
291  theNeutron.SetMomentum(pcmRot);
292  G4double eN = std::sqrt(pN*pN + mN*mN); // Scattered neutron energy
293  theNeutron.SetTotalEnergy(eN);
294 
295  // Get the scattered target momentum
296  G4ReactionProduct toLab(-1.*theTarget);
297  theTarget.SetMomentum(pDir*Pinit - pcmRot);
298  G4double eT = Einit - eN + mT;
299  theTarget.SetTotalEnergy(eT);
300 
301  // Now back to lab frame
302  theNeutron.Lorentz(theNeutron, toLab);
303  theTarget.Lorentz(theTarget, toLab);
304 
305  //111005 Protection for not producing 0 kinetic energy target
306  if (theNeutron.GetKineticEnergy() <= 0)
307  theNeutron.SetTotalEnergy(theNeutron.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
308  if (theTarget.GetKineticEnergy() <= 0)
309  theTarget.SetTotalEnergy(theTarget.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
310 
311  } else if (frameFlag == 2) {
312  // Projectile scattering values cosTh taken from center of mass tabulation
313 
314  G4LorentzVector proj(nEnergy, the3Neutron);
315  G4LorentzVector targ(tEnergy, the3Target);
316  G4ThreeVector boostToCM = proj.findBoostToCM(targ);
317  proj.boost(boostToCM);
318  targ.boost(boostToCM);
319 
320  // Rotate projectile and target momenta by CM scattering angle
321  // Note: at this point collision axis is not along z axis, due to
322  // momentum given target nucleus by thermal process
323  G4double px = proj.px();
324  G4double py = proj.py();
325  G4double pz = proj.pz();
326 
327  G4ThreeVector pcmRot;
328  pcmRot.setX(px*cosTh*cosPhi - py*sinPhi + pz*sinth*cosPhi);
329  pcmRot.setY(px*cosTh*sinPhi + py*cosPhi + pz*sinth*sinPhi);
330  pcmRot.setZ(-px*sinth + pz*cosTh);
331  proj.setVect(pcmRot);
332  targ.setVect(-pcmRot);
333 
334  // Back to lab frame
335  proj.boost(-boostToCM);
336  targ.boost(-boostToCM);
337 
338  theNeutron.SetMomentum(proj.vect() );
339  theNeutron.SetTotalEnergy(proj.e() );
340 
341  theTarget.SetMomentum(targ.vect() );
342  theTarget.SetTotalEnergy(targ.e() );
343 
344  //080904 Add Protection for very low energy (1e-6eV) scattering
345  if (theNeutron.GetKineticEnergy() <= 0) {
346  theNeutron.SetTotalEnergy(theNeutron.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
347  }
348 
349  //080904 Add Protection for very low energy (1e-6eV) scattering
350  if (theTarget.GetKineticEnergy() <= 0) {
351  theTarget.SetTotalEnergy(theTarget.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
352  }
353 
354  } else {
355  G4cout << "Value of frameFlag (1=LAB, 2=CMS): " << frameFlag;
356  throw G4HadronicException(__FILE__, __LINE__,
357  "G4ParticleHPElasticFS::ApplyYourSelf frameflag incorrect");
358  }
359 
360  // Everything is now in the lab frame
361  // Set energy change and momentum change
363  theResult.Get()->SetMomentumChange(theNeutron.GetMomentum().unit());
364 
365  // Make recoil a G4DynamicParticle
366  G4DynamicParticle* theRecoil = new G4DynamicParticle;
367  theRecoil->SetDefinition(G4IonTable::GetIonTable()->GetIon(static_cast<G4int>(theBaseZ),
368  static_cast<G4int>(theBaseA), 0) );
369  theRecoil->SetMomentum(theTarget.GetMomentum());
370  theResult.Get()->AddSecondary(theRecoil);
371 
372  // Postpone the tracking of the primary neutron
374  return theResult.Get();
375 }
376