Geant4_10
G4LENDModel.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 // 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 
46 :G4HadronicInteraction( name )
47 {
48 
49  proj = NULL; //will be set in an inherited class
50 
51  SetMinEnergy( 0.*eV );
52  SetMaxEnergy( 20.*MeV );
53 
54  //default_evaluation = "endl99";
55  default_evaluation = "ENDF.B-VII.0";
56 
57  allow_nat = false;
58  allow_any = false;
59 
61 
62 }
63 
65 {
66  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
67  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
68  {
69  delete it->second;
70  }
71 }
72 
73 
75 {
76 
77  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
78  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
79  {
80  delete it->second;
81  }
82  usedTarget_map.clear();
83 
85 
86 }
87 
88 
89 
91 {
92 
94 
95  size_t numberOfElements = G4Element::GetNumberOfElements();
96  static const G4ElementTable* theElementTable = G4Element::GetElementTable();
97 
98  for ( size_t i = 0 ; i < numberOfElements ; ++i )
99  {
100 
101  const G4Element* anElement = (*theElementTable)[i];
102  G4int numberOfIsotope = anElement->GetNumberOfIsotopes();
103 
104  if ( numberOfIsotope > 0 )
105  {
106  // User Defined Abundances
107  for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ )
108  {
109  G4int iZ = anElement->GetIsotope( i_iso )->GetZ();
110  G4int iA = anElement->GetIsotope( i_iso )->GetN();
111  G4int iIsomer = anElement->GetIsotope( i_iso )->Getm();
112 
113  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iA );
114  if ( allow_nat == true ) aTarget->AllowNat();
115  if ( allow_any == true ) aTarget->AllowAny();
116  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iA , iIsomer ) , aTarget ) );
117  }
118  }
119  else
120  {
121  // Natural Abundances
123  G4int iZ = int ( anElement->GetZ() );
124  //G4cout << nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ) << G4endl;
125  G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) );
126 
127  for ( G4int ii = 0 ; ii < numberOfNistIso ; ii++ )
128  {
129  //G4cout << nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + i ) << G4endl;
130  if ( nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii ) > 0 )
131  {
132  G4int iMass = nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii;
133  //G4cout << iZ << " " << nistElementBuild->GetNistFirstIsotopeN( iZ ) + i << " " << nistElementBuild->GetIsotopeAbundance ( iZ , iMass ) << G4endl;
134  G4int iIsomer = 0;
135 
136  G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iMass );
137  if ( allow_nat == true ) aTarget->AllowNat();
138  if ( allow_any == true ) aTarget->AllowAny();
139  usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iMass , iIsomer ) , aTarget ) );
140 
141  }
142 
143  }
144 
145  }
146  }
147 
148 
149 
150  G4cout << "Dump UsedTarget for " << GetModelName() << G4endl;
151  G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) , Pointer of Target" << G4endl;
152  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
153  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
154  {
155  G4cout
156  << " " << it->second->GetWantedEvaluation()
157  << ", " << it->second->GetWantedZ()
158  << ", " << it->second->GetWantedA()
159  << " -> " << it->second->GetActualEvaluation()
160  << ", " << it->second->GetActualZ()
161  << ", " << it->second->GetActualA()
162  << ", " << it->second->GetTarget()
163  << G4endl;
164  }
165 
166 }
167 
168 
169 
170 #include "G4ParticleTable.hh"
171 
173 {
174 
175  G4double temp = aTrack.GetMaterial()->GetTemperature();
176 
177  //G4int iZ = int ( aTarg.GetZ() );
178  //G4int iA = int ( aTarg.GetN() );
179  //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
180  G4int iZ = aTarg.GetZ_asInt();
181  G4int iA = aTarg.GetA_asInt();
182  G4int iM = 0;
183  if ( aTarg.GetIsotope() != NULL ) {
184  iM = aTarg.GetIsotope()->Getm();
185  }
186 
187  G4double ke = aTrack.GetKineticEnergy();
188 
189  G4HadFinalState* theResult = new G4HadFinalState();
190 
191  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA , iM ) )->second->GetTarget();
192 
193  G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, NULL, NULL );
194 
195  G4double phi = twopi*G4UniformRand();
196  G4double theta = std::acos( aMu );
197  //G4double sinth = std::sin( theta );
198 
199  G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>( aTrack.GetDefinition() ) );
200  theNeutron.SetMomentum( aTrack.Get4Momentum().vect() );
201  theNeutron.SetKineticEnergy( ke );
202 
203  G4ReactionProduct theTarget( G4ParticleTable::GetParticleTable()->FindIon( iZ , iA , 0 , iZ ) );
204 
206 
207 // add Thermal motion
208  G4double kT = k_Boltzmann*temp;
209  G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass )
210  , G4RandGauss::shoot() * std::sqrt( kT*mass )
211  , G4RandGauss::shoot() * std::sqrt( kT*mass ) );
212 
213  theTarget.SetMomentum( v );
214 
215 
216  G4ThreeVector the3Neutron = theNeutron.GetMomentum();
217  G4double nEnergy = theNeutron.GetTotalEnergy();
218  G4ThreeVector the3Target = theTarget.GetMomentum();
219  G4double tEnergy = theTarget.GetTotalEnergy();
220  G4ReactionProduct theCMS;
221  G4double totE = nEnergy+tEnergy;
222  G4ThreeVector the3CMS = the3Target+the3Neutron;
223  theCMS.SetMomentum(the3CMS);
224  G4double cmsMom = std::sqrt(the3CMS*the3CMS);
225  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
226  theCMS.SetMass(sqrts);
227  theCMS.SetTotalEnergy(totE);
228 
229  theNeutron.Lorentz(theNeutron, theCMS);
230  theTarget.Lorentz(theTarget, theCMS);
231  G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
232  G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS
233  G4double cms_theta=cms3Mom.theta();
234  G4double cms_phi=cms3Mom.phi();
235  G4ThreeVector tempVector;
236  tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
237  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
238  -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
239  tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
240  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
241  +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
242  tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
243  -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
244  tempVector *= en;
245  theNeutron.SetMomentum(tempVector);
246  theTarget.SetMomentum(-tempVector);
247  G4double tP = theTarget.GetTotalMomentum();
248  G4double tM = theTarget.GetMass();
249  theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
250  theNeutron.Lorentz(theNeutron, -1.*theCMS);
251  theTarget.Lorentz(theTarget, -1.*theCMS);
252 
253  theResult->SetEnergyChange(theNeutron.GetKineticEnergy());
254  theResult->SetMomentumChange(theNeutron.GetMomentum().unit());
255  G4DynamicParticle* theRecoil = new G4DynamicParticle;
256 
257  theRecoil->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( iZ, iA , 0, iZ ) );
258  theRecoil->SetMomentum( theTarget.GetMomentum() );
259 
260  theResult->AddSecondary( theRecoil );
261 
262  return theResult;
263 
264 }
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:158
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
ThreeVector shoot(const G4int Ap, const G4int Af)
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4NistElementBuilder * GetNistElementBuilder()
G4double GetZ() const
Definition: G4Element.hh:131
void recreate_used_target_map()
Definition: G4LENDModel.cc:74
G4double GetIsotopeAbundance(G4int Z, G4int N) const
const XML_Char * name
Definition: expat.h:151
G4int GetNumberOfNistIsotopes(G4int Z) const
G4ParticleDefinition * proj
Definition: G4LENDModel.hh:77
const G4String & GetModelName() const
int G4int
Definition: G4Types.hh:78
void setY(double)
G4LENDModel(G4String name="LENDModel")
Definition: G4LENDModel.cc:45
void setZ(double)
void setX(double)
G4int GetNucleusEncoding(G4int iZ, G4int iA, G4int iM)
void SetMinEnergy(G4double anEnergy)
void SetMass(const G4double mas)
Hep3Vector vect() const
G4int GetN() const
Definition: G4Isotope.hh:94
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
float k_Boltzmann
Definition: hepunit.py:299
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
const G4ParticleDefinition * GetDefinition() const
double getElasticFinalState(double e_in, double temperature, double(*rng)(void *), void *rngState)
G4int Getm() const
Definition: G4Isotope.hh:100
std::map< G4int, G4LENDUsedTarget * > usedTarget_map
Definition: G4LENDModel.hh:79
G4double GetKineticEnergy() const
G4int GetNistFirstIsotopeN(G4int Z) const
void SetTotalEnergy(const G4double en)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
double phi() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4LENDModel.cc:172
const G4LorentzVector & Get4Momentum() const
tuple v
Definition: test.py:18
double theta() const
G4double GetKineticEnergy() const
G4int GetZ() const
Definition: G4Isotope.hh:91
void SetEnergyChange(G4double anEnergy)
G4double GetTotalEnergy() const
static G4LENDManager * GetInstance()
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
double mass
Definition: G4GIDI_mass.cc:72
void create_used_target_map()
Definition: G4LENDModel.cc:90
const G4Isotope * GetIsotope()
Definition: G4Nucleus.hh:119
Hep3Vector unit() const
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void SetMaxEnergy(const G4double anEnergy)
G4ThreeVector GetMomentum() const
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:169
G4double GetTemperature() const
Definition: G4Material.hh:180
#define G4endl
Definition: G4ios.hh:61
const G4Material * GetMaterial() const
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4bool RequestChangeOfVerboseLevel(G4int)
G4LENDManager * lend_manager
Definition: G4LENDModel.hh:78
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
void SetMomentumChange(const G4ThreeVector &aV)
G4double GetMass() const
void AddSecondary(G4DynamicParticle *aP)