Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4NeutronHPInelastic.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// this code implementation is the intellectual property of
27// neutron_hp -- source file
28// J.P. Wellisch, Nov-1996
29// A prototype of the low energy neutron transport model.
30//
31// By copying, distributing or modifying the Program (or any work
32// based on the Program) you indicate your acceptance of this statement,
33// and all its terms.
34//
35// $Id$
36//
37// 070523 bug fix for G4FPE_DEBUG on by A. Howard (and T. Koi)
38// 081203 limit maximum trial for creating final states add protection for 1H isotope case by T. Koi
39//
40#include "G4NeutronHPInelastic.hh"
41#include "G4SystemOfUnits.hh"
42
43#include "G4NeutronHPManager.hh"
44
45 G4NeutronHPInelastic::G4NeutronHPInelastic()
46 :G4HadronicInteraction("NeutronHPInelastic")
47 {
48 SetMinEnergy( 0.0 );
49 SetMaxEnergy( 20.*MeV );
50
51 G4int istatus = system("echo $G4NEUTRONHPDATA");
52 if ( istatus < 0 )
53 {
54 G4cout << "Warning! system(\"echo $G4NEUTRONHPDATA\") returns error value at G4NeutronHPInelastic" << G4endl;
55 }
56
57// G4cout << " entering G4NeutronHPInelastic constructor"<<G4endl;
58 if(!getenv("G4NEUTRONHPDATA"))
59 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
60 dirName = getenv("G4NEUTRONHPDATA");
61 G4String tString = "/Inelastic";
62 dirName = dirName + tString;
63 numEle = G4Element::GetNumberOfElements();
64/*
65 theInelastic = new G4NeutronHPChannelList[numEle];
66 for (G4int i=0; i<numEle; i++)
67 {
68 theInelastic[i].Init((*(G4Element::GetElementTable()))[i], dirName);
69 G4int itry = 0;
70 do
71 {
72 theInelastic[i].Register(&theNFS, "F01"); // has
73 theInelastic[i].Register(&theNXFS, "F02");
74 theInelastic[i].Register(&the2NDFS, "F03");
75 theInelastic[i].Register(&the2NFS, "F04"); // has, E Done
76 theInelastic[i].Register(&the3NFS, "F05"); // has, E Done
77 theInelastic[i].Register(&theNAFS, "F06");
78 theInelastic[i].Register(&theN3AFS, "F07");
79 theInelastic[i].Register(&the2NAFS, "F08");
80 theInelastic[i].Register(&the3NAFS, "F09");
81 theInelastic[i].Register(&theNPFS, "F10");
82 theInelastic[i].Register(&theN2AFS, "F11");
83 theInelastic[i].Register(&the2N2AFS, "F12");
84 theInelastic[i].Register(&theNDFS, "F13");
85 theInelastic[i].Register(&theNTFS, "F14");
86 theInelastic[i].Register(&theNHe3FS, "F15");
87 theInelastic[i].Register(&theND2AFS, "F16");
88 theInelastic[i].Register(&theNT2AFS, "F17");
89 theInelastic[i].Register(&the4NFS, "F18"); // has, E Done
90 theInelastic[i].Register(&the2NPFS, "F19");
91 theInelastic[i].Register(&the3NPFS, "F20");
92 theInelastic[i].Register(&theN2PFS, "F21");
93 theInelastic[i].Register(&theNPAFS, "F22");
94 theInelastic[i].Register(&thePFS, "F23");
95 theInelastic[i].Register(&theDFS, "F24");
96 theInelastic[i].Register(&theTFS, "F25");
97 theInelastic[i].Register(&theHe3FS, "F26");
98 theInelastic[i].Register(&theAFS, "F27");
99 theInelastic[i].Register(&the2AFS, "F28");
100 theInelastic[i].Register(&the3AFS, "F29");
101 theInelastic[i].Register(&the2PFS, "F30");
102 theInelastic[i].Register(&thePAFS, "F31");
103 theInelastic[i].Register(&theD2AFS, "F32");
104 theInelastic[i].Register(&theT2AFS, "F33");
105 theInelastic[i].Register(&thePDFS, "F34");
106 theInelastic[i].Register(&thePTFS, "F35");
107 theInelastic[i].Register(&theDAFS, "F36");
108 theInelastic[i].RestartRegistration();
109 itry++;
110 }
111 //while(!theInelastic[i].HasDataInAnyFinalState());
112 while( !theInelastic[i].HasDataInAnyFinalState() && itry < 6 );
113 // 6 is corresponding to the value(5) of G4NeutronHPChannel. TK
114
115 if ( itry == 6 )
116 {
117 // No Final State at all.
118 G4bool exceptional = false;
119 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
120 {
121 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
122 }
123 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
124 }
125 }
126*/
127
128 for (G4int i=0; i<numEle; i++)
129 {
130 theInelastic.push_back( new G4NeutronHPChannelList );
131 (*theInelastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
132 G4int itry = 0;
133 do
134 {
135 (*theInelastic[i]).Register(&theNFS, "F01"); // has
136 (*theInelastic[i]).Register(&theNXFS, "F02");
137 (*theInelastic[i]).Register(&the2NDFS, "F03");
138 (*theInelastic[i]).Register(&the2NFS, "F04"); // has, E Done
139 (*theInelastic[i]).Register(&the3NFS, "F05"); // has, E Done
140 (*theInelastic[i]).Register(&theNAFS, "F06");
141 (*theInelastic[i]).Register(&theN3AFS, "F07");
142 (*theInelastic[i]).Register(&the2NAFS, "F08");
143 (*theInelastic[i]).Register(&the3NAFS, "F09");
144 (*theInelastic[i]).Register(&theNPFS, "F10");
145 (*theInelastic[i]).Register(&theN2AFS, "F11");
146 (*theInelastic[i]).Register(&the2N2AFS, "F12");
147 (*theInelastic[i]).Register(&theNDFS, "F13");
148 (*theInelastic[i]).Register(&theNTFS, "F14");
149 (*theInelastic[i]).Register(&theNHe3FS, "F15");
150 (*theInelastic[i]).Register(&theND2AFS, "F16");
151 (*theInelastic[i]).Register(&theNT2AFS, "F17");
152 (*theInelastic[i]).Register(&the4NFS, "F18"); // has, E Done
153 (*theInelastic[i]).Register(&the2NPFS, "F19");
154 (*theInelastic[i]).Register(&the3NPFS, "F20");
155 (*theInelastic[i]).Register(&theN2PFS, "F21");
156 (*theInelastic[i]).Register(&theNPAFS, "F22");
157 (*theInelastic[i]).Register(&thePFS, "F23");
158 (*theInelastic[i]).Register(&theDFS, "F24");
159 (*theInelastic[i]).Register(&theTFS, "F25");
160 (*theInelastic[i]).Register(&theHe3FS, "F26");
161 (*theInelastic[i]).Register(&theAFS, "F27");
162 (*theInelastic[i]).Register(&the2AFS, "F28");
163 (*theInelastic[i]).Register(&the3AFS, "F29");
164 (*theInelastic[i]).Register(&the2PFS, "F30");
165 (*theInelastic[i]).Register(&thePAFS, "F31");
166 (*theInelastic[i]).Register(&theD2AFS, "F32");
167 (*theInelastic[i]).Register(&theT2AFS, "F33");
168 (*theInelastic[i]).Register(&thePDFS, "F34");
169 (*theInelastic[i]).Register(&thePTFS, "F35");
170 (*theInelastic[i]).Register(&theDAFS, "F36");
171 (*theInelastic[i]).RestartRegistration();
172 itry++;
173 }
174 while( !(*theInelastic[i]).HasDataInAnyFinalState() && itry < 6 );
175 // 6 is corresponding to the value(5) of G4NeutronHPChannel. TK
176
177 if ( itry == 6 )
178 {
179 // No Final State at all.
180 G4bool exceptional = false;
181 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
182 {
183 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
184 }
185 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
186 }
187
188 }
189 }
190
191 G4NeutronHPInelastic::~G4NeutronHPInelastic()
192 {
193// delete [] theInelastic;
194 for ( std::vector<G4NeutronHPChannelList*>::iterator
195 it = theInelastic.begin() ; it != theInelastic.end() ; it++ )
196 {
197 delete *it;
198 }
199 theInelastic.clear();
200 }
201
202 #include "G4NeutronHPThermalBoost.hh"
203
204 G4HadFinalState * G4NeutronHPInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
205 {
206 if ( numEle < (G4int)G4Element::GetNumberOfElements() ) addChannelForNewElement();
207 G4NeutronHPManager::GetInstance()->OpenReactionWhiteBoard();
208 const G4Material * theMaterial = aTrack.GetMaterial();
209 G4int n = theMaterial->GetNumberOfElements();
210 G4int index = theMaterial->GetElement(0)->GetIndex();
211 G4int it=0;
212 if(n!=1)
213 {
214 xSec = new G4double[n];
215 G4double sum=0;
216 G4int i;
217 const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
218 G4double rWeight;
219 G4NeutronHPThermalBoost aThermalE;
220 for (i=0; i<n; i++)
221 {
222 index = theMaterial->GetElement(i)->GetIndex();
223 rWeight = NumAtomsPerVolume[i];
224 //xSec[i] = theInelastic[index].GetXsec(aThermalE.GetThermalEnergy(aTrack,
225 xSec[i] = (*theInelastic[index]).GetXsec(aThermalE.GetThermalEnergy(aTrack,
226 theMaterial->GetElement(i),
227 theMaterial->GetTemperature()));
228 xSec[i] *= rWeight;
229 sum+=xSec[i];
230 }
231 G4double random = G4UniformRand();
232 G4double running = 0;
233 for (i=0; i<n; i++)
234 {
235 running += xSec[i];
236 index = theMaterial->GetElement(i)->GetIndex();
237 it = i;
238 //if(random<=running/sum) break;
239 if( sum == 0 || random<=running/sum) break;
240 }
241 delete [] xSec;
242 }
243
244 //return theInelastic[index].ApplyYourself(theMaterial->GetElement(it), aTrack);
245 G4HadFinalState* result = (*theInelastic[index]).ApplyYourself(theMaterial->GetElement(it), aTrack);
246 //
247 aNucleus.SetParameters(G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargA(),G4NeutronHPManager::GetInstance()->GetReactionWhiteBoard()->GetTargZ());
248 G4NeutronHPManager::GetInstance()->CloseReactionWhiteBoard();
249 return result;
250 }
251
252const std::pair<G4double, G4double> G4NeutronHPInelastic::GetFatalEnergyCheckLevels() const
253{
254 // max energy non-conservation is mass of heavy nucleus
255// if ( getenv("G4NEUTRONHP_DO_NOT_ADJUST_FINAL_STATE") ) return std::pair<G4double, G4double>(5*perCent,250*GeV);
256 // This should be same to the hadron default value
257// return std::pair<G4double, G4double>(10*perCent,10*GeV);
258 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
259}
260
261void G4NeutronHPInelastic::addChannelForNewElement()
262{
263 for ( G4int i = numEle ; i < (G4int)G4Element::GetNumberOfElements() ; i++ )
264 {
265 G4cout << "G4NeutronHPInelastic Prepairing Data for the new element of " << (*(G4Element::GetElementTable()))[i]->GetName() << G4endl;
266
267 theInelastic.push_back( new G4NeutronHPChannelList );
268 (*theInelastic[i]).Init((*(G4Element::GetElementTable()))[i], dirName);
269 G4int itry = 0;
270 do
271 {
272 (*theInelastic[i]).Register(&theNFS, "F01"); // has
273 (*theInelastic[i]).Register(&theNXFS, "F02");
274 (*theInelastic[i]).Register(&the2NDFS, "F03");
275 (*theInelastic[i]).Register(&the2NFS, "F04"); // has, E Done
276 (*theInelastic[i]).Register(&the3NFS, "F05"); // has, E Done
277 (*theInelastic[i]).Register(&theNAFS, "F06");
278 (*theInelastic[i]).Register(&theN3AFS, "F07");
279 (*theInelastic[i]).Register(&the2NAFS, "F08");
280 (*theInelastic[i]).Register(&the3NAFS, "F09");
281 (*theInelastic[i]).Register(&theNPFS, "F10");
282 (*theInelastic[i]).Register(&theN2AFS, "F11");
283 (*theInelastic[i]).Register(&the2N2AFS, "F12");
284 (*theInelastic[i]).Register(&theNDFS, "F13");
285 (*theInelastic[i]).Register(&theNTFS, "F14");
286 (*theInelastic[i]).Register(&theNHe3FS, "F15");
287 (*theInelastic[i]).Register(&theND2AFS, "F16");
288 (*theInelastic[i]).Register(&theNT2AFS, "F17");
289 (*theInelastic[i]).Register(&the4NFS, "F18"); // has, E Done
290 (*theInelastic[i]).Register(&the2NPFS, "F19");
291 (*theInelastic[i]).Register(&the3NPFS, "F20");
292 (*theInelastic[i]).Register(&theN2PFS, "F21");
293 (*theInelastic[i]).Register(&theNPAFS, "F22");
294 (*theInelastic[i]).Register(&thePFS, "F23");
295 (*theInelastic[i]).Register(&theDFS, "F24");
296 (*theInelastic[i]).Register(&theTFS, "F25");
297 (*theInelastic[i]).Register(&theHe3FS, "F26");
298 (*theInelastic[i]).Register(&theAFS, "F27");
299 (*theInelastic[i]).Register(&the2AFS, "F28");
300 (*theInelastic[i]).Register(&the3AFS, "F29");
301 (*theInelastic[i]).Register(&the2PFS, "F30");
302 (*theInelastic[i]).Register(&thePAFS, "F31");
303 (*theInelastic[i]).Register(&theD2AFS, "F32");
304 (*theInelastic[i]).Register(&theT2AFS, "F33");
305 (*theInelastic[i]).Register(&thePDFS, "F34");
306 (*theInelastic[i]).Register(&thePTFS, "F35");
307 (*theInelastic[i]).Register(&theDAFS, "F36");
308 (*theInelastic[i]).RestartRegistration();
309 itry++;
310 }
311 while( !(*theInelastic[i]).HasDataInAnyFinalState() && itry < 6 );
312 // 6 is corresponding to the value(5) of G4NeutronHPChannel. TK
313
314 if ( itry == 6 )
315 {
316 // No Final State at all.
317 G4bool exceptional = false;
318 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
319 {
320 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
321 }
322 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
323 }
324 }
325
326 numEle = (G4int)G4Element::GetNumberOfElements();
327}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4Element::GetNumberOfElements
static size_t GetNumberOfElements()
Definition: G4Element.cc:406
G4Element::GetIndex
size_t GetIndex() const
Definition: G4Element.hh:182
G4Element::GetElementTable
static const G4ElementTable * GetElementTable()
Definition: G4Element.cc:399
G4HadProjectile::GetMaterial
const G4Material * GetMaterial() const
Definition: G4HadProjectile.hh:76
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103
G4Material::GetTemperature
G4double GetTemperature() const
Definition: G4Material.hh:181
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:201
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
G4Material::GetVecNbOfAtomsPerVolume
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:205
G4NeutronHPInelastic::ApplyYourself
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4NeutronHPInelastic.cc:204
G4NeutronHPInelastic::GetFatalEnergyCheckLevels
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
Definition: G4NeutronHPInelastic.cc:252
G4NeutronHPManager::GetInstance
static G4NeutronHPManager * GetInstance()
Definition: G4NeutronHPManager.hh:42
G4NeutronHPManager::GetReactionWhiteBoard
G4NeutronHPReactionWhiteBoard * GetReactionWhiteBoard()
Definition: G4NeutronHPManager.cc:53
G4NeutronHPThermalBoost::GetThermalEnergy
G4double GetThermalEnergy(const G4HadProjectile &aP, const G4Element *anE, G4double aT)
Definition: G4NeutronHPThermalBoost.hh:42
G4Nucleus::SetParameters
void SetParameters(const G4double A, const G4double Z)
Definition: G4Nucleus.cc:198
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83