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