Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4LENDModel.cc
이 파일의 문서화 페이지로 가기
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 // Class Description
27 // Final state production model for a LEND (Low Energy Nuclear Data)
28 // LEND is Geant4 interface for GIDI (General Interaction Data Interface)
29 // which gives a discription of nuclear and atomic reactions, such as
30 // Binary collision cross sections
31 // Particle number multiplicity distributions of reaction products
32 // Energy and angular distributions of reaction products
33 // Derived calculational constants
34 // GIDI is developped at Lawrence Livermore National Laboratory
35 // Class Description - End
36 
37 // 071025 First implementation done by T. Koi (SLAC/SCCS)
38 // 101118 Name modifications for release T. Koi (SLAC/PPA)
39 
40 #include "G4LENDModel.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4NistManager.hh"
44 
45 double MyRNG(void*) { return G4Random::getTheEngine()->flat(); }
46 
48 :G4HadronicInteraction( name )
49 {
50 
51  proj = NULL; //will be set in an inherited class
52 
53  SetMinEnergy( 0.*eV );
54  SetMaxEnergy( 20.*MeV );
55 
56  //default_evaluation = "endl99";
57  //default_evaluation = "ENDF.B-VII.0";
58  default_evaluation = "ENDF/BVII.1";
59 
60  allow_nat = false;
61  allow_any = false;
62 
64 
65 }
66 
68 {
69  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
70  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
71  {
72  delete it->second;
73  }
74 }
75 
76 
78 {
79 
80  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
81  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
82  {
83  delete it->second;
84  }
85  usedTarget_map.clear();
86 
88 
89 }
90 
91 
92 
94 {
95 
97 
98  size_t numberOfElements = G4Element::GetNumberOfElements();
99  static const G4ElementTable* theElementTable = G4Element::GetElementTable();
100 
101  for ( size_t i = 0 ; i < numberOfElements ; ++i )
102  {
103 
104  const G4Element* anElement = (*theElementTable)[i];
105  G4int numberOfIsotope = anElement->GetNumberOfIsotopes();
106 
107  if ( numberOfIsotope > 0 )
108  {
109  // User Defined Abundances
110  for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
111  {
112  G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
113  G4int iA = anElement->GetIsotope( i_iso )->GetN();
114  G4int iIsomer = anElement->GetIsotope( i_iso )->Getm();
115 
116  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iA , iIsomer );
117  if ( allow_nat == true ) aTarget->AllowNat();
118  if ( allow_any == true ) aTarget->AllowAny();
119  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iA , iIsomer ) , aTarget ) );
120  }
121  }
122  else
123  {
124  // Natural Abundances
126  G4int iZ = int ( anElement->GetZ() );
127  //G4cout << nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ) << G4endl;
128  G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) );
129 
130  for ( G4int ii = 0 ; ii < numberOfNistIso ; ii++ )
131  {
132  //G4cout << nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + i ) << G4endl;
133  if ( nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii ) > 0 )
134  {
135  G4int iMass = nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii;
136  //G4cout << iZ << " " << nistElementBuild->GetNistFirstIsotopeN( iZ ) + i << " " << nistElementBuild->GetIsotopeAbundance ( iZ , iMass ) << G4endl;
137  G4int iIsomer = 0;
138 
139  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iMass );
140  if ( allow_nat == true ) aTarget->AllowNat();
141  if ( allow_any == true ) aTarget->AllowAny();
142  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iMass , iIsomer ) , aTarget ) );
143 
144  }
145 
146  }
147 
148  }
149  }
150 
152 }
153 
154 
155 
156 #include "G4IonTable.hh"
157 
159 {
160 
161  G4double temp = aTrack.GetMaterial()->GetTemperature();
162 
163  //G4int iZ = int ( aTarg.GetZ() );
164  //G4int iA = int ( aTarg.GetN() );
165  //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
166  G4int iZ = aTarg.GetZ_asInt();
167  G4int iA = aTarg.GetA_asInt();
168  G4int iM = 0;
169  if ( aTarg.GetIsotope() != NULL ) {
170  iM = aTarg.GetIsotope()->Getm();
171  }
172 
173  G4double ke = aTrack.GetKineticEnergy();
174 
175  G4HadFinalState* theResult = new G4HadFinalState();
176 
177  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA , iM ) )->second->GetTarget();
178 
179  G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, NULL, NULL );
180 
181  G4double phi = twopi*G4UniformRand();
182  G4double theta = std::acos( aMu );
183  //G4double sinth = std::sin( theta );
184 
185  G4ReactionProduct theNeutron( aTrack.GetDefinition() );
186  theNeutron.SetMomentum( aTrack.Get4Momentum().vect() );
187  theNeutron.SetKineticEnergy( ke );
188 
189  G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( iZ , iA , iM );
191 
192  G4double mass = pd->GetPDGMass();
193 
194 // add Thermal motion
195  G4double kT = k_Boltzmann*temp;
196  G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass )
197  , G4RandGauss::shoot() * std::sqrt( kT*mass )
198  , G4RandGauss::shoot() * std::sqrt( kT*mass ) );
199 
200  theTarget.SetMomentum( v );
201 
202 
203  G4ThreeVector the3Neutron = theNeutron.GetMomentum();
204  G4double nEnergy = theNeutron.GetTotalEnergy();
205  G4ThreeVector the3Target = theTarget.GetMomentum();
206  G4double tEnergy = theTarget.GetTotalEnergy();
207  G4ReactionProduct theCMS;
208  G4double totE = nEnergy+tEnergy;
209  G4ThreeVector the3CMS = the3Target+the3Neutron;
210  theCMS.SetMomentum(the3CMS);
211  G4double cmsMom = std::sqrt(the3CMS*the3CMS);
212  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
213  theCMS.SetMass(sqrts);
214  theCMS.SetTotalEnergy(totE);
215 
216  theNeutron.Lorentz(theNeutron, theCMS);
217  theTarget.Lorentz(theTarget, theCMS);
218  G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
219  G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS
220  G4double cms_theta=cms3Mom.theta();
221  G4double cms_phi=cms3Mom.phi();
222  G4ThreeVector tempVector;
223  tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
224  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
225  -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
226  tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
227  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
228  +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
229  tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
230  -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
231  tempVector *= en;
232  theNeutron.SetMomentum(tempVector);
233  theTarget.SetMomentum(-tempVector);
234  G4double tP = theTarget.GetTotalMomentum();
235  G4double tM = theTarget.GetMass();
236  theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
237  theNeutron.Lorentz(theNeutron, -1.*theCMS);
238  theTarget.Lorentz(theTarget, -1.*theCMS);
239 
240  theResult->SetEnergyChange(theNeutron.GetKineticEnergy());
241  theResult->SetMomentumChange(theNeutron.GetMomentum().unit());
242  G4DynamicParticle* theRecoil = new G4DynamicParticle;
243 
244  theRecoil->SetDefinition( G4IonTable::GetIonTable()->GetIon( iZ , iA , iM , iZ ) );
245  theRecoil->SetMomentum( theTarget.GetMomentum() );
246 
247  theResult->AddSecondary( theRecoil );
248 
249  return theResult;
250 
251 }
252 
254  if ( lend_manager->GetVerboseLevel() >= 1 ) {
256  message = "Produce unchanged final state is requested in ";
257  message += this->GetModelName();
258  message += ". Cross section and model likely have an inconsistency.";
259  G4Exception( "G4LENDModel::returnUnchanged(,)" , "LENDModel-01" , JustWarning ,
260  message );
261  }
262  theResult->SetEnergyChange( aTrack.GetKineticEnergy() );
263  theResult->SetMomentumChange( aTrack.Get4Momentum().getV().unit() );
264  return theResult;
265 }
266 
268  G4GIDI_target* target = NULL;
269  if ( usedTarget_map.find( nuclear_code ) != usedTarget_map.end() ) {
270  target = usedTarget_map.find( nuclear_code )->second->GetTarget();
271  }
272  return target;
273 }
274 
276 
277  if ( lend_manager->GetVerboseLevel() >= 1 || force ) {
278  if ( usedTarget_map.size() == 0 ) create_used_target_map();
279  G4cout << "Dumping UsedTarget of " << GetModelName() << " for " << proj->GetParticleName() << G4endl;
280  G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) " << G4endl;
281  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
282  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ ) {
283  G4cout
284  << " " << it->second->GetWantedEvaluation()
285  << ", " << it->second->GetWantedZ()
286  << ", " << it->second->GetWantedA()
287  << " -> " << it->second->GetActualEvaluation()
288  << ", " << it->second->GetActualZ()
289  << ", " << it->second->GetActualA()
290  << G4endl;
291  }
292  }
293 }
const XML_Char * name
Definition: expat.h:151
G4int GetNistFirstIsotopeN(G4int Z) const
void SetMinEnergy(G4double anEnergy)
G4LENDModel(G4String name="LENDModel")
Definition: G4LENDModel.cc:47
void SetMaxEnergy(const G4double anEnergy)
const XML_Char * target
Definition: expat.h:268
static constexpr double MeV
Definition: G4SIunits.hh:214
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
void SetMomentumChange(const G4ThreeVector &aV)
void SetMomentum(const G4double x, const G4double y, const G4double z)
#define G4endl
Definition: G4ios.hh:61
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
const G4String & GetParticleName() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
static constexpr double second
Definition: G4SIunits.hh:157
G4bool RequestChangeOfVerboseLevel(G4int)
G4int GetZ() const
Definition: G4Isotope.hh:91
G4int GetVerboseLevel()
G4double GetPDGMass() const
void setX(double)
double theta() const
void SetTotalEnergy(const G4double en)
G4LENDManager * lend_manager
Definition: G4LENDModel.hh:84
void SetMass(const G4double mas)
void SetEnergyChange(G4double anEnergy)
G4NistElementBuilder * GetNistElementBuilder()
G4double GetTotalEnergy() const
G4bool allow_nat
Definition: G4LENDModel.hh:72
G4int GetNucleusEncoding(G4int iZ, G4int iA, G4int iM)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
const G4Material * GetMaterial() const
double MyRNG(void *)
Definition: G4LENDModel.cc:45
void setZ(double)
G4double GetIsotopeAbundance(G4int Z, G4int N) const
const G4Isotope * GetIsotope()
Definition: G4Nucleus.hh:119
G4HadFinalState * returnUnchanged(const G4HadProjectile &aTrack, G4HadFinalState *theResult)
Definition: G4LENDModel.cc:253
G4int GetNumberOfNistIsotopes(G4int Z) const
std::vector< G4Element * > G4ElementTable
double getElasticFinalState(double e_in, double temperature, double(*rng)(void *), void *rngState)
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
G4ErrorTarget * theTarget
Definition: errprop.cc:59
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4double GetTotalMomentum() const
std::map< G4int, G4LENDUsedTarget * > usedTarget_map
Definition: G4LENDModel.hh:85
static constexpr double eV
Definition: G4SIunits.hh:215
G4ThreeVector GetMomentum() const
G4int Getm() const
Definition: G4Isotope.hh:100
G4ParticleDefinition * proj
Definition: G4LENDModel.hh:83
void DumpLENDTargetInfo(G4bool force=false)
Definition: G4LENDModel.cc:275
ThreeVector shoot(const G4int Ap, const G4int Af)
const G4LorentzVector & Get4Momentum() const
Hep3Vector unit() const
G4double GetKineticEnergy() const
Hep3Vector getV() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4double GetMass() const
const G4ParticleDefinition * GetDefinition() const
G4int GetN() const
Definition: G4Isotope.hh:94
int G4int
Definition: G4Types.hh:78
double phi() const
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4String default_evaluation
Definition: G4LENDModel.hh:71
static constexpr double k_Boltzmann
G4bool allow_any
Definition: G4LENDModel.hh:73
void SetMomentum(const G4ThreeVector &momentum)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4GLOB_DLL std::ostream G4cout
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:170
Hep3Vector vect() const
G4double GetZ() const
Definition: G4Element.hh:131
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:504
static G4LENDManager * GetInstance()
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void recreate_used_target_map()
Definition: G4LENDModel.cc:77
G4GIDI_target * get_target_from_map(G4int nuclear_code)
Definition: G4LENDModel.cc:267
void setY(double)
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4double GetTemperature() const
Definition: G4Material.hh:183
void create_used_target_map()
Definition: G4LENDModel.cc:93
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
const G4String & GetModelName() const
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:159
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4LENDModel.cc:158