source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelastic.cc@ 1350

Last change on this file since 1350 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 7.0 KB
Line 
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: G4NeutronHPInelastic.cc,v 1.24 2008/12/03 22:28:48 tkoi Exp $
36// GEANT4 tag $Name: geant4-09-03-ref-09 $
37//
38// 070523 bug fix for G4FPE_DEBUG on by A. Howard (and T. Koi)
39// 081203 limit maximum trial for creating final states add protection for 1H isotope case by T. Koi
40//
41#include "G4NeutronHPInelastic.hh"
42
43 G4NeutronHPInelastic::G4NeutronHPInelastic()
44 :G4HadronicInteraction("NeutronHPInelastic")
45 {
46 SetMinEnergy( 0.0 );
47 SetMaxEnergy( 20.*MeV );
48 system("echo $G4NEUTRONHPDATA");
49// G4cout << " entering G4NeutronHPInelastic constructor"<<G4endl;
50 if(!getenv("G4NEUTRONHPDATA"))
51 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
52 dirName = getenv("G4NEUTRONHPDATA");
53 G4String tString = "/Inelastic/";
54 dirName = dirName + tString;
55 numEle = G4Element::GetNumberOfElements();
56 theInelastic = new G4NeutronHPChannelList[numEle];
57
58 for (G4int i=0; i<numEle; i++)
59 {
60 theInelastic[i].Init((*(G4Element::GetElementTable()))[i], dirName);
61 G4int itry = 0;
62 do
63 {
64 theInelastic[i].Register(&theNFS, "F01"); // has
65 theInelastic[i].Register(&theNXFS, "F02");
66 theInelastic[i].Register(&the2NDFS, "F03");
67 theInelastic[i].Register(&the2NFS, "F04"); // has, E Done
68 theInelastic[i].Register(&the3NFS, "F05"); // has, E Done
69 theInelastic[i].Register(&theNAFS, "F06");
70 theInelastic[i].Register(&theN3AFS, "F07");
71 theInelastic[i].Register(&the2NAFS, "F08");
72 theInelastic[i].Register(&the3NAFS, "F09");
73 theInelastic[i].Register(&theNPFS, "F10");
74 theInelastic[i].Register(&theN2AFS, "F11");
75 theInelastic[i].Register(&the2N2AFS, "F12");
76 theInelastic[i].Register(&theNDFS, "F13");
77 theInelastic[i].Register(&theNTFS, "F14");
78 theInelastic[i].Register(&theNHe3FS, "F15");
79 theInelastic[i].Register(&theND2AFS, "F16");
80 theInelastic[i].Register(&theNT2AFS, "F17");
81 theInelastic[i].Register(&the4NFS, "F18"); // has, E Done
82 theInelastic[i].Register(&the2NPFS, "F19");
83 theInelastic[i].Register(&the3NPFS, "F20");
84 theInelastic[i].Register(&theN2PFS, "F21");
85 theInelastic[i].Register(&theNPAFS, "F22");
86 theInelastic[i].Register(&thePFS, "F23");
87 theInelastic[i].Register(&theDFS, "F24");
88 theInelastic[i].Register(&theTFS, "F25");
89 theInelastic[i].Register(&theHe3FS, "F26");
90 theInelastic[i].Register(&theAFS, "F27");
91 theInelastic[i].Register(&the2AFS, "F28");
92 theInelastic[i].Register(&the3AFS, "F29");
93 theInelastic[i].Register(&the2PFS, "F30");
94 theInelastic[i].Register(&thePAFS, "F31");
95 theInelastic[i].Register(&theD2AFS, "F32");
96 theInelastic[i].Register(&theT2AFS, "F33");
97 theInelastic[i].Register(&thePDFS, "F34");
98 theInelastic[i].Register(&thePTFS, "F35");
99 theInelastic[i].Register(&theDAFS, "F36");
100 theInelastic[i].RestartRegistration();
101 itry++;
102 }
103 //while(!theInelastic[i].HasDataInAnyFinalState());
104 while( !theInelastic[i].HasDataInAnyFinalState() && itry < 6 );
105 // 6 is corresponding to the value(5) of G4NeutronHPChannel. TK
106 if ( itry == 6 )
107 {
108 // No Final State at all.
109 G4bool exceptional = false;
110 if ( (*(G4Element::GetElementTable()))[i]->GetNumberOfIsotopes() == 1 )
111 {
112 if ( (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetZ() == 1 && (*(G4Element::GetElementTable()))[i]->GetIsotope( 0 )->GetN() == 1 ) exceptional = true; //1H
113 }
114 if ( !exceptional ) throw G4HadronicException(__FILE__, __LINE__, "Channel: Do not know what to do with this element");
115 }
116 }
117 }
118
119 G4NeutronHPInelastic::~G4NeutronHPInelastic()
120 {
121 delete [] theInelastic;
122 }
123
124 #include "G4NeutronHPThermalBoost.hh"
125
126 G4HadFinalState * G4NeutronHPInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& )
127 {
128 const G4Material * theMaterial = aTrack.GetMaterial();
129 G4int n = theMaterial->GetNumberOfElements();
130 G4int index = theMaterial->GetElement(0)->GetIndex();
131 G4int it=0;
132 if(n!=1)
133 {
134 xSec = new G4double[n];
135 G4double sum=0;
136 G4int i;
137 const G4double * NumAtomsPerVolume = theMaterial->GetVecNbOfAtomsPerVolume();
138 G4double rWeight;
139 G4NeutronHPThermalBoost aThermalE;
140 for (i=0; i<n; i++)
141 {
142 index = theMaterial->GetElement(i)->GetIndex();
143 rWeight = NumAtomsPerVolume[i];
144 xSec[i] = theInelastic[index].GetXsec(aThermalE.GetThermalEnergy(aTrack,
145 theMaterial->GetElement(i),
146 theMaterial->GetTemperature()));
147 xSec[i] *= rWeight;
148 sum+=xSec[i];
149 }
150 G4double random = G4UniformRand();
151 G4double running = 0;
152 for (i=0; i<n; i++)
153 {
154 running += xSec[i];
155 index = theMaterial->GetElement(i)->GetIndex();
156 it = i;
157 //if(random<=running/sum) break;
158 if( sum == 0 || random<=running/sum) break;
159 }
160 delete [] xSec;
161 }
162 return theInelastic[index].ApplyYourself(theMaterial->GetElement(it), aTrack);
163 }
Note: See TracBrowser for help on using the repository browser.