source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPAngular.cc @ 1340

Last change on this file since 1340 was 962, checked in by garnier, 15 years ago

update processes
File size: 10.0 KB
Line 
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//
33#include "G4NeutronHPAngular.hh"
34
35void G4NeutronHPAngular::Init(std::ifstream & aDataFile)
36{
37//  G4cout << "here we are entering the Angular Init"<<G4endl;
38  aDataFile >> theAngularDistributionType >> targetMass;
39  aDataFile >> frameFlag;
40  if(theAngularDistributionType == 0)
41  {
42    theIsoFlag = true;
43  }
44  else if(theAngularDistributionType==1)
45  {
46    G4int nEnergy;
47    aDataFile >> nEnergy; 
48    theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
49    theCoefficients->InitInterpolation(aDataFile);
50    G4double temp, energy;
51    G4int tempdep, nLegendre;
52    G4int i, ii;
53    for (i=0; i<nEnergy; i++)
54    {
55      aDataFile >> temp >> energy >> tempdep >> nLegendre;
56      energy *=eV;
57      theCoefficients->Init(i, energy, nLegendre);
58      theCoefficients->SetTemperature(i, temp);
59      G4double coeff=0;
60      for(ii=0; ii<nLegendre; ii++)
61      {
62        aDataFile >> coeff;
63        theCoefficients->SetCoeff(i, ii+1, coeff);
64      }
65    }
66  }
67  else if (theAngularDistributionType==2)
68  {
69    G4int nEnergy;
70    aDataFile >> nEnergy;
71    theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
72    theProbArray->InitInterpolation(aDataFile);
73    G4double temp, energy;
74    G4int tempdep;
75    for(G4int i=0; i<nEnergy; i++)
76    {
77      aDataFile >> temp >> energy >> tempdep;
78      energy *= eV;
79      theProbArray->SetT(i, temp);
80      theProbArray->SetX(i, energy);
81      theProbArray->InitData(i, aDataFile);
82    }
83  }
84  else
85  {
86    theIsoFlag = false;
87    G4cout << "unknown distribution found for Angular"<<G4endl;
88    throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
89  }   
90}
91
92void G4NeutronHPAngular::SampleAndUpdate(G4ReactionProduct & aHadron)
93{
94  if(theIsoFlag)
95  {
96//  G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
97// @@@ add code for isotropic emission in CMS.
98      G4double costheta = 2.*G4UniformRand()-1;
99      G4double theta = std::acos(costheta);
100      G4double phi = twopi*G4UniformRand();
101      G4double sinth = std::sin(theta);
102      G4double en = aHadron.GetTotalMomentum();
103      G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
104      aHadron.SetMomentum( temp );
105      aHadron.Lorentz(aHadron, -1.*theTarget);
106  }
107  else
108  {
109    if(frameFlag == 1) // LAB
110    {
111      G4double en = aHadron.GetTotalMomentum();
112      G4ReactionProduct boosted;
113      boosted.Lorentz(theNeutron, theTarget);
114      G4double kineticEnergy = boosted.GetKineticEnergy();
115      G4double cosTh = 0.0; 
116      if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy); 
117      if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy); 
118      G4double theta = std::acos(cosTh);
119      G4double phi = twopi*G4UniformRand();
120      G4double sinth = std::sin(theta);
121      G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
122      aHadron.SetMomentum( temp );
123    }
124    else if(frameFlag == 2) // costh in CMS
125    {
126      G4ReactionProduct boostedN;
127      boostedN.Lorentz(theNeutron, theTarget);
128      G4double kineticEnergy = boostedN.GetKineticEnergy();
129
130      G4double cosTh = 0.0; 
131      if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy); 
132      if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy); 
133
134//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern)
135/*
136    if(theAngularDistributionType == 1) // LAB
137    {
138      G4double en = aHadron.GetTotalMomentum();
139      G4ReactionProduct boosted;
140      boosted.Lorentz(theNeutron, theTarget);
141      G4double kineticEnergy = boosted.GetKineticEnergy();
142      G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
143      G4double theta = std::acos(cosTh);
144      G4double phi = twopi*G4UniformRand();
145      G4double sinth = std::sin(theta);
146      G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
147      aHadron.SetMomentum( temp );
148    }
149    else if(theAngularDistributionType == 2) // costh in CMS
150    {
151*/
152
153//      G4ReactionProduct boostedN;
154//      boostedN.Lorentz(theNeutron, theTarget);
155//      G4double kineticEnergy = boostedN.GetKineticEnergy();
156//      G4double cosTh = theProbArray->Sample(kineticEnergy);
157
158      G4double theta = std::acos(cosTh);
159      G4double phi = twopi*G4UniformRand();
160      G4double sinth = std::sin(theta);     
161     
162      G4ThreeVector temp(sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); //CMS
163
164//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
165/*
166      G4double en = aHadron.GetTotalEnergy(); // Target rest
167     
168      // get trafo from Target rest frame to CMS
169      G4ReactionProduct boostedT;
170      boostedT.Lorentz(theTarget, theTarget);
171     
172      G4ThreeVector the3Neutron = boostedN.GetMomentum();
173      G4double nEnergy = boostedN.GetTotalEnergy();
174      G4ThreeVector the3Target = boostedT.GetMomentum();
175      G4double tEnergy = boostedT.GetTotalEnergy();
176      G4double totE = nEnergy+tEnergy;
177      G4ThreeVector the3trafo = -the3Target-the3Neutron;
178      G4ReactionProduct trafo; // for transformation from CMS to target rest frame
179      trafo.SetMomentum(the3trafo);
180      G4double cmsMom = std::sqrt(the3trafo*the3trafo);
181      G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
182      trafo.SetMass(sqrts);
183      trafo.SetTotalEnergy(totE);
184     
185      G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
186      G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
187      G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
188      fac*=gamma;
189     
190      G4double mom;
191//    For G4FPE_DEBUG ON
192      G4double mom2 = ( en*fac*en*fac -
193                   (fac*fac - gamma*gamma)*
194                   (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
195                );
196      if ( mom2 > 0.0 )
197        mom = std::sqrt( mom2 );
198      else
199        mom = 0.0;
200
201      mom = -en*fac - mom;
202      mom /= (fac*fac-gamma*gamma);
203      temp = mom*temp;
204     
205      aHadron.SetMomentum( temp ); // now all in CMS
206      aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
207      aHadron.Lorentz(aHadron, trafo); // now in target rest frame
208*/
209      // Determination of the hadron kinetic energy in CMS
210      // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
211      // kineticEnergy is incident neutron kinetic energy  in CMS (or target frame) 
212      G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
213      G4double A1     =   theTarget.GetMass()/boostedN.GetMass(); 
214      G4double A1prim =   aHadron.GetMass()/ boostedN.GetMass();
215      G4double kinE   = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*kineticEnergy+(1+A1)*QValue);
216      G4double totalE = kinE + aHadron.GetMass();
217      G4double mom2   = totalE*totalE - aHadron.GetMass()*aHadron.GetMass();
218      G4double mom;
219      if ( mom2 > 0.0 ) mom = std::sqrt( mom2 );
220      else mom = 0.0;     
221
222      aHadron.SetMomentum( mom*temp ); // Set momentum in CMS
223      aHadron.SetKineticEnergy(kinE);  // Set kinetic energy in CMS
224
225      // get trafo from Target rest frame to CMS
226      G4ReactionProduct boostedT;
227      boostedT.Lorentz(theTarget, theTarget);
228     
229      G4ThreeVector the3Neutron = boostedN.GetMomentum();
230      G4double nEnergy = boostedN.GetTotalEnergy();
231      G4ThreeVector the3Target = boostedT.GetMomentum();
232      G4double tEnergy = boostedT.GetTotalEnergy();
233      G4double totE = nEnergy+tEnergy;
234      G4ThreeVector the3trafo = -the3Target-the3Neutron;
235      G4ReactionProduct trafo; // for transformation from CMS to target rest frame
236      trafo.SetMomentum(the3trafo);
237      G4double cmsMom = std::sqrt(the3trafo*the3trafo);
238      G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
239      trafo.SetMass(sqrts);
240      trafo.SetTotalEnergy(totE);
241
242      aHadron.Lorentz(aHadron, trafo);
243
244    }
245    else
246    {
247      throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
248    }
249  }
250  aHadron.Lorentz(aHadron, -1.*theTarget); 
251//  G4cout << aHadron.GetMomentum()<<" ";
252//  G4cout << aHadron.GetTotalMomentum()<<G4endl;
253}
Note: See TracBrowser for help on using the repository browser.