Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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