Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
252 const 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 
261 void 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 }
Generated on Sat May 25 2013 14:34:04 for Geant4 by   doxygen 1.8.4