ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPAngular.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPAngular.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 // 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;
46  if(theAngularDistributionType == 0 )
47  {
48  theIsoFlag = true;
49  }
50  else if(theAngularDistributionType==1)
51  {
52  theIsoFlag = false;
53  G4int nEnergy;
54  aDataFile >> nEnergy;
55  theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
56  theCoefficients->InitInterpolation(aDataFile);
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);
65  theCoefficients->SetTemperature(i, temp);
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 
100 void G4ParticleHPAngular::SampleAndUpdate(G4ReactionProduct & aHadron)
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 }