| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Class Description
27 // Final state production model for a LEND (Lo
28 // LEND is Geant4 interface for GIDI (General
29 // which gives a discription of nuclear and at
30 // Binary collision cross sections
31 // Particle number multiplicity distributio
32 // Energy and angular distributions of reac
33 // Derived calculational constants
34 // GIDI is developped at Lawrence Livermore Na
35 // Class Description - End
36
37 // 071025 First implementation done by T. Koi
38 // 101118 Name modifications for release T. Ko
39
40 #include "G4LENDModel.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4NistManager.hh"
44 #include "G4PhysicsModelCatalog.hh"
45
46 double MyRNG(void*) { return G4Random::getThe
47
48 G4LENDModel::G4LENDModel( G4String name )
49 :G4HadronicInteraction( name ), secID(-1)
50 {
51
52 proj = NULL; //will be set in an inherited
53
54 SetMinEnergy( 0.*eV );
55 SetMaxEnergy( 20.*MeV );
56
57 //default_evaluation = "endl99";
58 //default_evaluation = "ENDF.B-VII.0";
59 default_evaluation = "ENDF/BVII.1";
60
61 allow_nat = false;
62 allow_any = false;
63
64 lend_manager = G4LENDManager::GetInstance()
65
66 secID = G4PhysicsModelCatalog::GetModelID(
67 }
68
69 G4LENDModel::~G4LENDModel()
70 {
71 for ( std::map< G4int , G4LENDUsedTarget* >
72 it = usedTarget_map.begin() ; it != u
73 {
74 delete it->second;
75 }
76 }
77
78
79 void G4LENDModel::recreate_used_target_map()
80 {
81
82 for ( std::map< G4int , G4LENDUsedTarget* >
83 it = usedTarget_map.begin() ; it != u
84 {
85 delete it->second;
86 }
87 usedTarget_map.clear();
88
89 create_used_target_map();
90
91 }
92
93
94
95 void G4LENDModel::create_used_target_map()
96 {
97
98 lend_manager->RequestChangeOfVerboseLevel(
99
100 std::size_t numberOfElements = G4Element::G
101 static const G4ElementTable* theElementTabl
102
103 for ( std::size_t i = 0 ; i < numberOfEleme
104 {
105
106 const G4Element* anElement = (*theElemen
107 G4int numberOfIsotope = (G4int)anElement
108
109 if ( numberOfIsotope > 0 )
110 {
111 // User Defined Abundances
112 for ( G4int i_iso = 0 ; i_iso < numbe
113 {
114 G4int iZ = anElement->GetIsotope(
115 G4int iA = anElement->GetIsotope(
116 G4int iIsomer = anElement->GetIsot
117
118 G4LENDUsedTarget* aTarget = new G4
119 if ( allow_nat == true ) aTarget->
120 if ( allow_any == true ) aTarget->
121 usedTarget_map.insert( std::pair<
122 }
123 }
124 else
125 {
126 // Natural Abundances
127 G4NistElementBuilder* nistElementBuil
128 G4int iZ = int ( anElement->GetZ() );
129 //G4cout << nistElementBuild->GetNumb
130 G4int numberOfNistIso = nistElementBu
131
132 for ( G4int ii = 0 ; ii < numberOfNis
133 {
134 //G4cout << nistElementBuild->GetI
135 if ( nistElementBuild->GetIsotopeA
136 {
137 G4int iMass = nistElementBuild-
138 //G4cout << iZ << " " << nistEl
139 G4int iIsomer = 0;
140
141 G4LENDUsedTarget* aTarget = new
142 if ( allow_nat == true ) aTarge
143 if ( allow_any == true ) aTarge
144 usedTarget_map.insert( std::pai
145
146 }
147
148 }
149
150 }
151 }
152
153 DumpLENDTargetInfo();
154 }
155
156
157
158 #include "G4IonTable.hh"
159
160 G4HadFinalState * G4LENDModel::ApplyYourself(c
161 {
162
163 G4double temp = aTrack.GetMaterial()->GetTe
164
165 //G4int iZ = int ( aTarg.GetZ() );
166 //G4int iA = int ( aTarg.GetN() );
167 //migrate to integer A and Z (GetN_asInt re
168 G4int iZ = aTarg.GetZ_asInt();
169 G4int iA = aTarg.GetA_asInt();
170 G4int iM = 0;
171 if ( aTarg.GetIsotope() != NULL ) {
172 iM = aTarg.GetIsotope()->Getm();
173 }
174
175 G4double ke = aTrack.GetKineticEnergy();
176
177 G4HadFinalState* theResult = new G4HadFinal
178
179 G4GIDI_target* aTarget = usedTarget_map.fin
180
181 G4double aMu = aTarget->getElasticFinalStat
182
183 G4double phi = twopi*G4UniformRand();
184 G4double theta = std::acos( aMu );
185 //G4double sinth = std::sin( theta );
186
187 G4ReactionProduct theNeutron( aTrack.GetDef
188 theNeutron.SetMomentum( aTrack.Get4Momentum
189 theNeutron.SetKineticEnergy( ke );
190
191 G4ParticleDefinition* pd = G4IonTable::GetI
192 G4ReactionProduct theTarget( pd );
193
194 G4double mass = pd->GetPDGMass();
195
196 // add Thermal motion
197 G4double kT = k_Boltzmann*temp;
198 G4ThreeVector v ( G4RandGauss::shoot() * st
199 , G4RandGauss::shoot() * st
200 , G4RandGauss::shoot() * st
201
202 theTarget.SetMomentum( v );
203
204
205 G4ThreeVector the3Neutron = theNeutron.Ge
206 G4double nEnergy = theNeutron.GetTotalEne
207 G4ThreeVector the3Target = theTarget.GetM
208 G4double tEnergy = theTarget.GetTotalEner
209 G4ReactionProduct theCMS;
210 G4double totE = nEnergy+tEnergy;
211 G4ThreeVector the3CMS = the3Target+the3Ne
212 theCMS.SetMomentum(the3CMS);
213 G4double cmsMom = std::sqrt(the3CMS*the3C
214 G4double sqrts = std::sqrt((totE-cmsMom)*
215 theCMS.SetMass(sqrts);
216 theCMS.SetTotalEnergy(totE);
217
218 theNeutron.Lorentz(theNeutron, theCMS);
219 theTarget.Lorentz(theTarget, theCMS);
220 G4double en = theNeutron.GetTotalMoment
221 G4ThreeVector cms3Mom=theNeutron.GetMom
222 G4double cms_theta=cms3Mom.theta();
223 G4double cms_phi=cms3Mom.phi();
224 G4ThreeVector tempVector;
225 tempVector.setX(std::cos(theta)*std::si
226 +std::sin(theta)*std::c
227 -std::sin(theta)*std::s
228 tempVector.setY(std::cos(theta)*std::si
229 +std::sin(theta)*std::c
230 +std::sin(theta)*std::s
231 tempVector.setZ(std::cos(theta)*std::co
232 -std::sin(theta)*std::c
233 tempVector *= en;
234 theNeutron.SetMomentum(tempVector);
235 theTarget.SetMomentum(-tempVector);
236 G4double tP = theTarget.GetTotalMomentu
237 G4double tM = theTarget.GetMass();
238 theTarget.SetTotalEnergy(std::sqrt((tP+
239 theNeutron.Lorentz(theNeutron, -1.*theC
240 theTarget.Lorentz(theTarget, -1.*theCMS
241
242 theResult->SetEnergyChange(theNeutron.Get
243 theResult->SetMomentumChange(theNeutron.G
244 G4DynamicParticle* theRecoil = new G4Dyna
245
246 theRecoil->SetDefinition( G4IonTable::Get
247 theRecoil->SetMomentum( theTarget.GetMome
248
249 theResult->AddSecondary( theRecoil );
250
251 return theResult;
252
253 }
254
255 G4HadFinalState* G4LENDModel::returnUnchanged(
256 if ( lend_manager->GetVerboseLevel() >= 1 )
257 G4String message;
258 message = "Produce unchanged final state
259 message += this->GetModelName();
260 message += ". Cross section and model li
261 G4Exception( "G4LENDModel::returnUnchang
262 message );
263 }
264 theResult->SetEnergyChange( aTrack.GetKinet
265 theResult->SetMomentumChange( aTrack.Get4Mo
266 return theResult;
267 }
268
269 G4GIDI_target* G4LENDModel::get_target_from_ma
270 G4GIDI_target* target = NULL;
271 if ( usedTarget_map.find( nuclear_code ) !=
272 target = usedTarget_map.find( nuclear_co
273 }
274 return target;
275 }
276
277 void G4LENDModel::DumpLENDTargetInfo( G4bool f
278
279 if ( lend_manager->GetVerboseLevel() >= 1 |
280 if ( usedTarget_map.size() == 0 ) create
281 G4cout << "Dumping UsedTarget of " << Ge
282 G4cout << "Requested Evaluation, Z , A -
283 for ( std::map< G4int , G4LENDUsedTarget
284 it = usedTarget_map.begin() ; it != u
285 G4cout
286 << " " << it->second->GetWantedEvalua
287 << ", " << it->second->GetWantedZ()
288 << ", " << it->second->GetWantedA()
289 << " -> " << it->second->GetActualEva
290 << ", " << it->second->GetActualZ()
291 << ", " << it->second->GetActualA()
292 << G4endl;
293 }
294 }
295 }
296