Geant4  10.02.p01
G4ParticleHPAngular.cc
Go to the documentation of this file.
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 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 
41 void G4ParticleHPAngular::Init(std::istream & aDataFile)
42 {
43 // G4cout << "here we are entering the Angular Init"<<G4endl;
44  aDataFile >> theAngularDistributionType >> targetMass;
45  aDataFile >> frameFlag;
47  {
48  theIsoFlag = true;
49  }
50  else if(theAngularDistributionType==1)
51  {
52  theIsoFlag = false;
53  G4int nEnergy;
54  aDataFile >> nEnergy;
57  G4double temp, energy;
58  G4int tempdep, nLegendre;
59  G4int i, ii;
60  for (i=0; i<nEnergy; i++)
61  {
62  aDataFile >> temp >> energy >> tempdep >> nLegendre;
63  energy *=eV;
64  theCoefficients->Init(i, energy, nLegendre);
66  G4double coeff=0;
67  for(ii=0; ii<nLegendre; ii++)
68  {
69  aDataFile >> coeff;
70  theCoefficients->SetCoeff(i, ii+1, coeff);
71  }
72  }
73  }
74  else if (theAngularDistributionType==2)
75  {
76  theIsoFlag = false;
77  G4int nEnergy;
78  aDataFile >> nEnergy;
79  theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
80  theProbArray->InitInterpolation(aDataFile);
81  G4double temp, energy;
82  G4int tempdep;
83  for(G4int i=0; i<nEnergy; i++)
84  {
85  aDataFile >> temp >> energy >> tempdep;
86  energy *= eV;
87  theProbArray->SetT(i, temp);
88  theProbArray->SetX(i, energy);
89  theProbArray->InitData(i, aDataFile);
90  }
91  }
92  else
93  {
94  theIsoFlag = false;
95  G4cout << "unknown distribution found for Angular: "<< theAngularDistributionType << G4endl;
96  throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
97  }
98 }
99 
101 {
102 
103  //********************************************************************
104  //EMendoza -> sampling can be isotropic in LAB or in CMS
105  /*
106  if(theIsoFlag)
107  {
108 // G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
109 // @@@ add code for isotropic emission in CMS.
110  G4double costheta = 2.*G4UniformRand()-1;
111  G4double theta = std::acos(costheta);
112  G4double phi = twopi*G4UniformRand();
113  G4double sinth = std::sin(theta);
114  G4double en = aHadron.GetTotalMomentum();
115  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
116  aHadron.SetMomentum( temp );
117  aHadron.Lorentz(aHadron, -1.*theTarget);
118  }
119  else
120  {
121  */
122  //********************************************************************
123  if(frameFlag == 1) // LAB
124  {
125  G4double en = aHadron.GetTotalMomentum();
126  G4ReactionProduct boosted;
127  boosted.Lorentz(*fCache.Get().theProjectileRP, *fCache.Get().theTarget);
128  G4double kineticEnergy = boosted.GetKineticEnergy();
129  G4double cosTh = 0.0;
130  //********************************************************************
131  //EMendoza --> sampling can be also isotropic
132  /*
133  if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
134  if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
135  */
136  //********************************************************************
137  if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
138  else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
139  else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
140  else{
141  G4cout << "unknown distribution found for Angular: "<< theAngularDistributionType << G4endl;
142  throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
143  }
144  //********************************************************************
145  G4double theta = std::acos(cosTh);
146  G4double phi = twopi*G4UniformRand();
147  G4double sinth = std::sin(theta);
148  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), 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){cosTh =2.*G4UniformRand()-1;}
166  else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
167  else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
168  else{
169  G4cout << "unknown distribution found for Angular: "<< theAngularDistributionType << G4endl;
170  throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
171  }
172  //********************************************************************
173 //080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern)
174 /*
175  if(theAngularDistributionType == 1) // LAB
176  {
177  G4double en = aHadron.GetTotalMomentum();
178  G4ReactionProduct boosted;
179  boosted.Lorentz(theProjectile, theTarget);
180  G4double kineticEnergy = boosted.GetKineticEnergy();
181  G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
182  G4double theta = std::acos(cosTh);
183  G4double phi = twopi*G4UniformRand();
184  G4double sinth = std::sin(theta);
185  G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
186  aHadron.SetMomentum( temp );
187  }
188  else if(theAngularDistributionType == 2) // costh in CMS {
189  }
190 */
191 
192 // G4ReactionProduct boostedN;
193 // boostedN.Lorentz(theProjectile, theTarget);
194 // G4double kineticEnergy = boostedN.GetKineticEnergy();
195 // G4double cosTh = theProbArray->Sample(kineticEnergy);
196 
197  G4double theta = std::acos(cosTh);
198  G4double phi = twopi*G4UniformRand();
199  G4double sinth = std::sin(theta);
200 
201  G4ThreeVector temp(sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); //CMS
202 
203 //080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
204 /*
205  G4double en = aHadron.GetTotalEnergy(); // Target rest
206 
207  // get trafo from Target rest frame to CMS
208  G4ReactionProduct boostedT;
209  boostedT.Lorentz(theTarget, theTarget);
210 
211  G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
212  G4double nEnergy = boostedN.GetTotalEnergy();
213  G4ThreeVector the3Target = boostedT.GetMomentum();
214  G4double tEnergy = boostedT.GetTotalEnergy();
215  G4double totE = nEnergy+tEnergy;
216  G4ThreeVector the3trafo = -the3Target-the3IncidentParticle;
217  G4ReactionProduct trafo; // for transformation from CMS to target rest frame
218  trafo.SetMomentum(the3trafo);
219  G4double cmsMom = std::sqrt(the3trafo*the3trafo);
220  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
221  trafo.SetMass(sqrts);
222  trafo.SetTotalEnergy(totE);
223 
224  G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
225  G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
226  G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
227  fac*=gamma;
228 
229  G4double mom;
230 // For G4FPE_DEBUG ON
231  G4double mom2 = ( en*fac*en*fac -
232  (fac*fac - gamma*gamma)*
233  (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
234  );
235  if ( mom2 > 0.0 )
236  mom = std::sqrt( mom2 );
237  else
238  mom = 0.0;
239 
240  mom = -en*fac - mom;
241  mom /= (fac*fac-gamma*gamma);
242  temp = mom*temp;
243 
244  aHadron.SetMomentum( temp ); // now all in CMS
245  aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
246  aHadron.Lorentz(aHadron, trafo); // now in target rest frame
247 */
248  // Determination of the hadron kinetic energy in CMS
249  // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
250  // kineticEnergy is incident neutron kinetic energy in CMS (or target frame)
251  G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
252  G4double A1 = fCache.Get().theTarget->GetMass()/boostedN.GetMass();
253  G4double A1prim = aHadron.GetMass()/ boostedN.GetMass();
254  G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*kineticEnergy+(1+A1)*QValue);
255  G4double totalE = kinE + aHadron.GetMass();
256  G4double mom2 = totalE*totalE - aHadron.GetMass()*aHadron.GetMass();
257  G4double mom;
258  if ( mom2 > 0.0 ) mom = std::sqrt( mom2 );
259  else mom = 0.0;
260 
261  aHadron.SetMomentum( mom*temp ); // Set momentum in CMS
262  aHadron.SetKineticEnergy(kinE); // Set kinetic energy in CMS
263 
264  // get trafo from Target rest frame to CMS
265  G4ReactionProduct boostedT;
266  boostedT.Lorentz(*fCache.Get().theTarget, *fCache.Get().theTarget);
267 
268  G4ThreeVector the3IncidentParticle = boostedN.GetMomentum();
269  G4double nEnergy = boostedN.GetTotalEnergy();
270  G4ThreeVector the3Target = boostedT.GetMomentum();
271  G4double tEnergy = boostedT.GetTotalEnergy();
272  G4double totE = nEnergy+tEnergy;
273  G4ThreeVector the3trafo = -the3Target-the3IncidentParticle;
274  G4ReactionProduct trafo; // for transformation from CMS to target rest frame
275  trafo.SetMomentum(the3trafo);
276  G4double cmsMom = std::sqrt(the3trafo*the3trafo);
277  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
278  trafo.SetMass(sqrts);
279  trafo.SetTotalEnergy(totE);
280 
281  aHadron.Lorentz(aHadron, trafo);
282 
283  }
284  else
285  {
286  throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
287  }
288  aHadron.Lorentz(aHadron, -1.*(*fCache.Get().theTarget));
289 // G4cout << aHadron.GetMomentum()<<" ";
290 // G4cout << aHadron.GetTotalMomentum()<<G4endl;
291 }
void Init(G4int i, G4double e, G4int n)
G4ParticleHPPartial * theProbArray
G4double GetTotalMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
CLHEP::Hep3Vector G4ThreeVector
void InitInterpolation(G4int i, std::istream &aDataFile)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4double SampleMax(G4double energy)
int G4int
Definition: G4Types.hh:78
void InitInterpolation(std::istream &aDataFile)
G4Cache< toBeCached > fCache
void SetMass(const G4double mas)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void SampleAndUpdate(G4ReactionProduct &anIncidentParticle)
G4double Sample(G4double x)
void SetTotalEnergy(const G4double en)
static const double twopi
Definition: G4SIunits.hh:75
G4double GetKineticEnergy() const
void SetTemperature(G4int i, G4double temp)
void SetX(G4int i, G4double x)
static const double eV
Definition: G4SIunits.hh:212
G4double GetTotalEnergy() const
G4double energy(const ThreeVector &p, const G4double m)
void SetCoeff(G4int i, G4int l, G4double coeff)
void InitData(G4int i, std::istream &aDataFile, G4double unit=1.)
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
void SetT(G4int i, G4double x)
double G4double
Definition: G4Types.hh:76
void Init(std::istream &aDataFile)
G4ParticleHPLegendreStore * theCoefficients
G4double GetMass() const