Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronHPorLEInelasticData.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 //
27 // 05-11-21 NeutronHP or Low Energy Parameterization Models
28 // Implemented by T. Koi (SLAC/SCCS)
29 // If NeutronHP data do not available for an element, then Low Energy
30 // Parameterization models handle the interactions of the element.
31 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
32 //
33 
35 #include "G4SystemOfUnits.hh"
36 #include "G4Neutron.hh"
37 #include "G4ElementTable.hh"
38 #include "G4NeutronHPData.hh"
39 
40 #include "G4PhysicsVector.hh"
41 
43 {
44  SetMinKinEnergy( 0*MeV );
45  SetMaxKinEnergy( 20*MeV );
46 
47  ke_cache = 0.0;
48  xs_cache = 0.0;
49  element_cache = NULL;
50  material_cache = NULL;
51 // BuildPhysicsTable(*G4Neutron::Neutron());
52 }
53 
55 {
56 // delete theCrossSections;
57 }
58 
60  G4int /*Z*/ , G4int /*A*/ ,
61  const G4Element* element ,
62  const G4Material* /*mat*/ )
63 {
64  G4double eKin = dp->GetKineticEnergy();
65  if ( eKin > GetMaxKinEnergy()
66  || eKin < GetMinKinEnergy()
67  || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
68  if ( unavailable_elements->find( element->GetName() ) != unavailable_elements->end() ) return false;
69 
70  return true;
71 }
72 
74  G4int /*Z*/ , G4int /*A*/ ,
75  const G4Isotope* /*iso*/ ,
76  const G4Element* element ,
77  const G4Material* material )
78 {
79  if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
80 
81  ke_cache = dp->GetKineticEnergy();
82  element_cache = element;
83  material_cache = material;
84  G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
85  xs_cache = xs;
86  return xs;
87  //return GetCrossSection( dp , element , material->GetTemperature() );
88 }
89 
91 :G4VCrossSectionDataSet("NeutronHPorLEInelasticXS")
92 {
93  theInelasticChannel = pChannel;
94  unavailable_elements = pSet;
95  //BuildPhysicsTable(*G4Neutron::Neutron());
96 
97  SetMinKinEnergy( 0*MeV );
98  SetMaxKinEnergy( 20*MeV );
99 
100  ke_cache = 0.0;
101  xs_cache = 0.0;
102  element_cache = NULL;
103  material_cache = NULL;
104 }
105 
106 /*
107 G4bool G4NeutronHPorLEInelasticData::IsApplicable(const G4DynamicParticle*aP, const G4Element* anElement)
108 {
109  G4bool result = true;
110  G4double eKin = aP->GetKineticEnergy();
111  if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
112  if ( unavailable_elements->find( anElement->GetName() ) != unavailable_elements->end() ) result = false;
113  return result;
114 }
115 */
116 
118 #include "G4LPhysicsFreeVector.hh"
119 //#include "G4NeutronHPElementData.hh"
120 
122 {
123  if( &aP!=G4Neutron::Neutron() )
124  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
125 
126  size_t numberOfElements = G4Element::GetNumberOfElements();
127  theCrossSections = new G4PhysicsTable( numberOfElements );
128 
129  static const G4ElementTable *theElementTable = G4Element::GetElementTable();
130  for ( size_t i=0 ; i < numberOfElements; ++i )
131  {
132  G4PhysicsVector* thePhysVec = new G4LPhysicsFreeVector(0, 0, 0);
133 
134  if ( unavailable_elements->find( (*theElementTable)[i]->GetName() ) == unavailable_elements->end() )
135  {
136 
137  G4NeutronHPElementData* theElementData = new G4NeutronHPElementData();
138  theElementData->Init( (*theElementTable)[i] );
139 
140  G4NeutronHPVector* theHPVector = theElementData->GetData( (G4NeutronHPInelasticData*)this );
141 
142  G4int len = theHPVector->GetVectorLength();
143 
144  if ( len!=0 )
145  {
146  G4double emin = theHPVector->GetX(0);
147  G4double emax = theHPVector->GetX(len-1);
148 
149  G4LPhysicsFreeVector* aPhysVector= new G4LPhysicsFreeVector ( len , emin , emax );
150  for ( G4int ii=0; ii<len; ii++ )
151  {
152  aPhysVector->PutValues( ii , theHPVector->GetX(ii) , theHPVector->GetY(ii) );
153  }
154  delete thePhysVec;
155  thePhysVec = aPhysVector;
156  }
157 
158  //G4PhysicsVector* physVec = G4NeutronHPData::
159  //Instance()->MakePhysicsVector((*theElementTable)[i], this);
160  //theCrossSections->push_back(physVec);
161  }
162 
163  theCrossSections->push_back(thePhysVec);
164  }
165 }
166 
167 
168 
170 {
171  if(&aP!=G4Neutron::Neutron())
172  throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
173 // G4cout << "G4NeutronHPorLEInelasticData::DumpPhysicsTable still to be implemented"<<G4endl;
174 }
175 
176 
177 
178 #include "G4Nucleus.hh"
179 #include "G4NucleiProperties.hh"
180 #include "G4Neutron.hh"
181 #include "G4Electron.hh"
182 
185 {
186 
187  // G4cout << "Choice G4NeutronHPorLEInelasticData for element " << anE->GetName() << G4endl;
188  G4double result = 0;
189  //G4bool outOfRange;
190  G4int index = anE->GetIndex();
191 
192  // prepare neutron
193  G4double eKinetic = aP->GetKineticEnergy();
194  G4ReactionProduct theNeutron( aP->GetDefinition() );
195  theNeutron.SetMomentum( aP->GetMomentum() );
196  theNeutron.SetKineticEnergy( eKinetic );
197 
198  // prepare thermal nucleus
199  G4Nucleus aNuc;
200  G4double eps = 0.0001;
201  G4double theA = anE->GetN();
202  G4double theZ = anE->GetZ();
203  G4double eleMass;
204  eleMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA+eps), static_cast<G4int>(theZ+eps))
206 
207  G4ReactionProduct boosted;
208  G4double aXsection;
209 
210  // MC integration loop
211  G4int counter = 0;
212  G4double buffer = 0;
213  G4int size = G4int(std::max(10., aT/60*kelvin));
214  G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
215  G4double neutronVMag = neutronVelocity.mag();
216  while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.03*buffer)
217  {
218  if(counter) buffer = result/counter;
219  while (counter<size)
220  {
221  counter ++;
222  G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
223  boosted.Lorentz(theNeutron, aThermalNuc);
224  G4double theEkin = boosted.GetKineticEnergy();
225  //aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
226  aXsection = theInelasticChannel[index].GetXsec( theEkin );
227  // velocity correction.
228  G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
229  aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
230  result += aXsection;
231  }
232  size += size;
233  }
234  result /= counter;
235  //return result;
236  return result*barn;
237 }