source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPElasticFS.cc @ 1196

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

update processes
File size: 14.5 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// 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#include "G4NeutronHPElasticFS.hh"
35#include "G4ReactionProduct.hh"
36#include "G4Nucleus.hh"
37#include "G4Proton.hh"
38#include "G4Deuteron.hh"
39#include "G4Triton.hh"
40#include "G4Alpha.hh"
41#include "G4ThreeVector.hh"
42#include "G4LorentzVector.hh"
43#include "G4ParticleTable.hh"
44#include "G4NeutronHPDataUsed.hh"
45
46  void G4NeutronHPElasticFS::Init (G4double A, G4double Z, G4String & dirName, G4String & )
47  {
48    G4String tString = "/FS/";
49    G4bool dbool;
50    G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), dirName, tString, dbool);
51    G4String filename = aFile.GetName();
52    theBaseA = aFile.GetA();
53    theBaseZ = aFile.GetZ();
54    if(!dbool)
55    {
56      hasAnyData = false;
57      hasFSData = false; 
58      hasXsec = false;
59      return;
60    }
61    std::ifstream theData(filename, std::ios::in);
62    theData >> repFlag >> targetMass >> frameFlag;
63    if(repFlag==1)
64    {
65      G4int nEnergy;
66      theData >> nEnergy; 
67      theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
68      theCoefficients->InitInterpolation(theData);
69      G4double temp, energy;
70      G4int tempdep, nLegendre;
71      G4int i, ii;
72      for (i=0; i<nEnergy; i++)
73      {
74        theData >> temp >> energy >> tempdep >> nLegendre;
75        energy *=eV;
76        theCoefficients->Init(i, energy, nLegendre);
77        theCoefficients->SetTemperature(i, temp);
78        G4double coeff=0;
79        for(ii=0; ii<nLegendre; ii++)
80        {
81          // load legendre coefficients.
82          theData >> coeff;
83          theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
84        }
85      }
86    }
87    else if (repFlag==2)
88    {
89      G4int nEnergy;
90      theData >> nEnergy;
91      theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
92      theProbArray->InitInterpolation(theData);
93      G4double temp, energy;
94      G4int tempdep, nPoints;
95      for(G4int i=0; i<nEnergy; i++)
96      {
97        theData >> temp >> energy >> tempdep >> nPoints;
98        energy *= eV;
99        theProbArray->InitInterpolation(i, theData);
100        theProbArray->SetT(i, temp);
101        theProbArray->SetX(i, energy);
102        G4double prob, costh;
103        for(G4int ii=0; ii<nPoints; ii++)
104        {
105          // fill probability arrays.
106          theData >> costh >> prob;
107          theProbArray->SetX(i, ii, costh);
108          theProbArray->SetY(i, ii, prob);
109        }
110      }
111    }
112    else if ( repFlag==3 )
113    {
114       G4int nEnergy_Legendre;
115       theData >> nEnergy_Legendre; 
116       theCoefficients = new G4NeutronHPLegendreStore( nEnergy_Legendre );
117       theCoefficients->InitInterpolation( theData );
118       G4double temp, energy;
119       G4int tempdep, nLegendre;
120       G4int i, ii;
121       for ( i = 0 ; i < nEnergy_Legendre ; i++ )
122       {
123          theData >> temp >> energy >> tempdep >> nLegendre;
124          energy *=eV;
125          theCoefficients->Init( i , energy , nLegendre );
126          theCoefficients->SetTemperature( i , temp );
127          G4double coeff = 0;
128          for ( ii = 0 ; ii < nLegendre ; ii++ )
129          {
130             // load legendre coefficients.
131             theData >> coeff;
132             theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
133          }
134       } 
135
136       tE_of_repFlag3 = energy; 
137
138       G4int nEnergy_Prob;
139       theData >> nEnergy_Prob;
140       theProbArray = new G4NeutronHPPartial( nEnergy_Prob , nEnergy_Prob );
141       theProbArray->InitInterpolation( theData );
142       G4int nPoints;
143       for ( G4int i=0 ; i < nEnergy_Prob ; i++ )
144       {
145          theData >> temp >> energy >> tempdep >> nPoints;
146
147          energy *= eV;
148
149//        consistensy check 
150          if ( i == 0 )
151             if ( energy !=  tE_of_repFlag3 )
152                G4cout << "Warning Trangition Energy of repFlag3 is not consistent." << G4endl; 
153
154          theProbArray->InitInterpolation( i , theData );
155          theProbArray->SetT( i , temp );
156          theProbArray->SetX( i , energy );
157          G4double prob, costh;
158          for( G4int ii = 0 ; ii < nPoints ; ii++ )
159          {
160             // fill probability arrays.
161             theData >> costh >> prob;
162             theProbArray->SetX( i , ii , costh );
163             theProbArray->SetY( i , ii , prob );
164          }
165       }
166    }
167    else if (repFlag==0)
168    {
169      theData >> frameFlag;
170    }
171    else
172    {
173      G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
174      throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
175    }
176    theData.close();
177  }
178  G4HadFinalState * G4NeutronHPElasticFS::ApplyYourself(const G4HadProjectile & theTrack)
179  { 
180//    G4cout << "G4NeutronHPElasticFS::ApplyYourself+"<<G4endl;
181    theResult.Clear();
182    G4double eKinetic = theTrack.GetKineticEnergy();
183    const G4HadProjectile *incidentParticle = &theTrack;
184    G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
185    theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
186    theNeutron.SetKineticEnergy( eKinetic );
187//    G4cout << "G4NeutronHPElasticFS::ApplyYourself++"<<eKinetic<<" "<<G4endl;
188//    G4cout << "CMSVALUES 0 "<<theNeutron.GetTotalMomentum()<<G4endl;
189   
190    G4ReactionProduct theTarget; 
191    G4Nucleus aNucleus;
192    G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
193    theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
194//     G4cout << "Nucleus-test"<<" "<<targetMass<<" ";
195//     G4cout << theTarget.GetMomentum().x()<<" ";
196//     G4cout << theTarget.GetMomentum().y()<<" ";
197//     G4cout << theTarget.GetMomentum().z()<<G4endl;
198   
199    // neutron and target defined as reaction products.
200
201// prepare lorentz-transformation to Lab.
202
203    G4ThreeVector the3Neutron = theNeutron.GetMomentum();
204    G4double nEnergy = theNeutron.GetTotalEnergy();
205    G4ThreeVector the3Target = theTarget.GetMomentum();
206//    cout << "@@@" << the3Target<<G4endl;
207    G4double tEnergy = theTarget.GetTotalEnergy();
208    G4ReactionProduct theCMS;
209    G4double totE = nEnergy+tEnergy;
210    G4ThreeVector the3CMS = the3Target+the3Neutron;
211    theCMS.SetMomentum(the3CMS);
212    G4double cmsMom = std::sqrt(the3CMS*the3CMS);
213    G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
214    theCMS.SetMass(sqrts);
215    theCMS.SetTotalEnergy(totE);
216   
217    // data come as fcn of n-energy in nuclear rest frame
218    G4ReactionProduct boosted;
219    boosted.Lorentz(theNeutron, theTarget);
220    eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
221    G4double cosTh = -2;
222    if(repFlag == 1)
223    {
224      cosTh = theCoefficients->SampleElastic(eKinetic);
225    }
226   
227    else if (repFlag==2)
228    {
229      cosTh = theProbArray->Sample(eKinetic);
230    }
231    else if (repFlag==3)
232    {
233       if ( eKinetic <= tE_of_repFlag3 )
234       {
235          cosTh = theCoefficients->SampleElastic(eKinetic);
236       }
237       else
238       {
239          cosTh = theProbArray->Sample(eKinetic);
240       }
241    }
242    else if (repFlag==0)
243    {
244      cosTh = 2.*G4UniformRand()-1.;
245    }
246    else
247    {
248      G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
249      throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
250    }
251    if(cosTh<-1.1) { return 0; }
252    G4double phi = twopi*G4UniformRand();
253    G4double theta = std::acos(cosTh);
254    G4double sinth = std::sin(theta);
255    if (frameFlag == 1) // final state data given in target rest frame.
256    {
257      // we have the scattering angle, now we need the energy, then do the
258      // boosting.
259      // relativistic elastic scattering energy angular correlation:
260      theNeutron.Lorentz(theNeutron, theTarget);
261      G4double e0 = theNeutron.GetTotalEnergy();
262      G4double p0 = theNeutron.GetTotalMomentum();
263      G4double mN = theNeutron.GetMass();
264      G4double mT = theTarget.GetMass();
265      G4double eE = e0+mT;
266      G4double ap = (mT+eE)*(mT-eE) + (p0+mN)*(p0-mN);
267      G4double a = 4*(eE+p0*cosTh)*(eE-p0*cosTh);
268      G4double b = 4*ap*p0*cosTh;
269      G4double c = (2.*eE*mN-ap)*(2.*eE*mN+ap);
270      G4double en = (-b+std::sqrt(b*b - 4*a*c) )/(2*a);
271      G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
272      theNeutron.SetMomentum(tempVector);
273      theNeutron.SetTotalEnergy(std::sqrt(en*en+theNeutron.GetMass()*theNeutron.GetMass()));
274      // first to lab     
275      theNeutron.Lorentz(theNeutron, -1.*theTarget);
276      // now to CMS
277      theNeutron.Lorentz(theNeutron, theCMS);
278      theTarget.SetMomentum(-theNeutron.GetMomentum());
279      theTarget.SetTotalEnergy(theNeutron.GetTotalEnergy());
280      // and back to lab
281      theNeutron.Lorentz(theNeutron, -1.*theCMS);
282      theTarget.Lorentz(theTarget, -1.*theCMS);     
283    }
284    else if (frameFlag == 2) // CMS
285    {
286      theNeutron.Lorentz(theNeutron, theCMS);
287      theTarget.Lorentz(theTarget, theCMS);
288      G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
289      G4ThreeVector cmsMom=theNeutron.GetMomentum(); // for neutron direction in CMS
290      G4double cms_theta=cmsMom.theta();
291      G4double cms_phi=cmsMom.phi();
292      G4ThreeVector tempVector;
293      tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
294                      +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
295                      -std::sin(theta)*std::sin(phi)*std::sin(cms_phi)  );
296      tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
297                      +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
298                      +std::sin(theta)*std::sin(phi)*std::cos(cms_phi)  );
299      tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
300                      -std::sin(theta)*std::cos(phi)*std::sin(cms_theta)  );
301      tempVector *= en;
302      theNeutron.SetMomentum(tempVector);
303      theTarget.SetMomentum(-tempVector);
304      G4double tP = theTarget.GetTotalMomentum();
305      G4double tM = theTarget.GetMass();
306      theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
307
308/*
309For debug purpose.
310Same transformation G4ReactionProduct.Lorentz() by 4vectors
311{
312G4LorentzVector n4p = G4LorentzVector ( theNeutron.GetMomentum() , theNeutron.GetKineticEnergy() + theNeutron.GetMass() );   
313G4cout << "before " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
314G4LorentzVector cm4p = G4LorentzVector ( theCMS.GetMomentum() , theCMS.GetKineticEnergy() + theCMS.GetMass() );   
315n4p.boost( cm4p.boostVector() );
316G4cout << cm4p/eV << G4endl;
317G4cout << "after " <<  ( n4p.e() - n4p.m() ) / eV<< G4endl;
318}
319*/
320
321      theNeutron.Lorentz(theNeutron, -1.*theCMS);
322//080904 Add Protection for very low energy (1e-6eV) scattering
323      if ( theNeutron.GetKineticEnergy() < 0 )
324      {
325         theNeutron.SetMomentum( G4ThreeVector(0) ); 
326         theNeutron.SetTotalEnergy ( theNeutron.GetMass() );
327      }
328
329      theTarget.Lorentz(theTarget, -1.*theCMS);
330//080904 Add Protection for very low energy (1e-6eV) scattering
331      if ( theTarget.GetKineticEnergy() < 0 )
332      {
333         theTarget.SetMomentum( G4ThreeVector(0) ); 
334         theTarget.SetTotalEnergy ( theTarget.GetMass() );
335      }
336    }
337    else
338    {
339      G4cout <<"Value of frameFlag (1=LAB, 2=CMS): "<<frameFlag;
340      throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::ApplyYourSelf frameflag incorrect");
341    }
342    // now all in Lab
343// nun den recoil generieren...und energy change, momentum change angeben.
344    theResult.SetEnergyChange(theNeutron.GetKineticEnergy());
345    theResult.SetMomentumChange(theNeutron.GetMomentum().unit());
346    G4DynamicParticle* theRecoil = new G4DynamicParticle;
347    if(targetMass<4.5)
348    {
349      if(targetMass<1)
350      {
351        // proton
352        theRecoil->SetDefinition(G4Proton::Proton());
353      }
354      else if(targetMass<2 )
355      {
356        // deuteron
357        theRecoil->SetDefinition(G4Deuteron::Deuteron());
358      }
359      else if(targetMass<2.999 )
360      {
361        // 3He
362        theRecoil->SetDefinition(G4He3::He3());
363      }
364      else if(targetMass<3 )
365      {
366        // Triton
367        theRecoil->SetDefinition(G4Triton::Triton());
368      }
369      else
370      {
371        // alpha
372        theRecoil->SetDefinition(G4Alpha::Alpha());
373      }
374    }
375    else
376    {
377      theRecoil->SetDefinition(G4ParticleTable::GetParticleTable()
378                               ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ)));
379    }
380    theRecoil->SetMomentum(theTarget.GetMomentum());
381    theResult.AddSecondary(theRecoil);
382//    G4cout << "G4NeutronHPElasticFS::ApplyYourself 10+"<<G4endl;
383    // postpone the tracking of the primary neutron
384     theResult.SetStatusChange(suspend);
385    return &theResult;
386  }
Note: See TracBrowser for help on using the repository browser.