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