Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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 
151 
152 
153  G4cout << "Dump UsedTarget for " << GetModelName() << G4endl;
154  //G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) , Pointer of Target" << G4endl;
155  G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) " << G4endl;
156  for ( std::map< G4int , G4LENDUsedTarget* >::iterator
157  it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ )
158  {
159  G4cout
160  << " " << it->second->GetWantedEvaluation()
161  << ", " << it->second->GetWantedZ()
162  << ", " << it->second->GetWantedA()
163  << " -> " << it->second->GetActualEvaluation()
164  << ", " << it->second->GetActualZ()
165  << ", " << it->second->GetActualA()
166  //<< ", " << it->second->GetTarget()
167  << G4endl;
168  }
169 
170 }
171 
172 
173 
174 #include "G4IonTable.hh"
175 
177 {
178 
179  G4double temp = aTrack.GetMaterial()->GetTemperature();
180 
181  //G4int iZ = int ( aTarg.GetZ() );
182  //G4int iA = int ( aTarg.GetN() );
183  //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
184  G4int iZ = aTarg.GetZ_asInt();
185  G4int iA = aTarg.GetA_asInt();
186  G4int iM = 0;
187  if ( aTarg.GetIsotope() != NULL ) {
188  iM = aTarg.GetIsotope()->Getm();
189  }
190 
191  G4double ke = aTrack.GetKineticEnergy();
192 
193  G4HadFinalState* theResult = new G4HadFinalState();
194 
195  G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA , iM ) )->second->GetTarget();
196 
197  G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, NULL, NULL );
198 
199  G4double phi = twopi*G4UniformRand();
200  G4double theta = std::acos( aMu );
201  //G4double sinth = std::sin( theta );
202 
203  G4ReactionProduct theNeutron( aTrack.GetDefinition() );
204  theNeutron.SetMomentum( aTrack.Get4Momentum().vect() );
205  theNeutron.SetKineticEnergy( ke );
206 
207  G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( iZ , iA , iM );
209 
210  G4double mass = pd->GetPDGMass();
211 
212 // add Thermal motion
213  G4double kT = k_Boltzmann*temp;
214  G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass )
215  , G4RandGauss::shoot() * std::sqrt( kT*mass )
216  , G4RandGauss::shoot() * std::sqrt( kT*mass ) );
217 
218  theTarget.SetMomentum( v );
219 
220 
221  G4ThreeVector the3Neutron = theNeutron.GetMomentum();
222  G4double nEnergy = theNeutron.GetTotalEnergy();
223  G4ThreeVector the3Target = theTarget.GetMomentum();
224  G4double tEnergy = theTarget.GetTotalEnergy();
225  G4ReactionProduct theCMS;
226  G4double totE = nEnergy+tEnergy;
227  G4ThreeVector the3CMS = the3Target+the3Neutron;
228  theCMS.SetMomentum(the3CMS);
229  G4double cmsMom = std::sqrt(the3CMS*the3CMS);
230  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
231  theCMS.SetMass(sqrts);
232  theCMS.SetTotalEnergy(totE);
233 
234  theNeutron.Lorentz(theNeutron, theCMS);
235  theTarget.Lorentz(theTarget, theCMS);
236  G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
237  G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS
238  G4double cms_theta=cms3Mom.theta();
239  G4double cms_phi=cms3Mom.phi();
240  G4ThreeVector tempVector;
241  tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
242  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
243  -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
244  tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
245  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
246  +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
247  tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
248  -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
249  tempVector *= en;
250  theNeutron.SetMomentum(tempVector);
251  theTarget.SetMomentum(-tempVector);
252  G4double tP = theTarget.GetTotalMomentum();
253  G4double tM = theTarget.GetMass();
254  theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
255  theNeutron.Lorentz(theNeutron, -1.*theCMS);
256  theTarget.Lorentz(theTarget, -1.*theCMS);
257 
258  theResult->SetEnergyChange(theNeutron.GetKineticEnergy());
259  theResult->SetMomentumChange(theNeutron.GetMomentum().unit());
260  G4DynamicParticle* theRecoil = new G4DynamicParticle;
261 
262  theRecoil->SetDefinition( G4IonTable::GetIonTable()->GetIon( iZ , iA , iM , iZ ) );
263  theRecoil->SetMomentum( theTarget.GetMomentum() );
264 
265  theResult->AddSecondary( theRecoil );
266 
267  return theResult;
268 
269 }
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:159
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
const XML_Char * name
Definition: expat.h:151
ThreeVector shoot(const G4int Ap, const G4int Af)
void SetMomentum(const G4ThreeVector &momentum)
G4double GetTotalMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
double MyRNG(void *)
Definition: G4LENDModel.cc:45
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4NistElementBuilder * GetNistElementBuilder()
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
G4double GetZ() const
Definition: G4Element.hh:131
void recreate_used_target_map()
Definition: G4LENDModel.cc:77
G4double GetIsotopeAbundance(G4int Z, G4int N) const
G4int GetNumberOfNistIsotopes(G4int Z) const
static constexpr double second
Definition: G4SIunits.hh:157
G4ParticleDefinition * proj
Definition: G4LENDModel.hh:79
const G4String & GetModelName() const
int G4int
Definition: G4Types.hh:78
void setY(double)
G4LENDModel(G4String name="LENDModel")
Definition: G4LENDModel.cc:47
void setZ(double)
void setX(double)
static constexpr double twopi
Definition: G4SIunits.hh:76
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:97
G4GLOB_DLL std::ostream G4cout
float k_Boltzmann
Definition: hepunit.py:299
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
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:81
G4double GetKineticEnergy() const
G4int GetNistFirstIsotopeN(G4int Z) const
void SetTotalEnergy(const G4double en)
G4ErrorTarget * theTarget
Definition: errprop.cc:59
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
static constexpr double eV
Definition: G4SIunits.hh:215
double phi() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4LENDModel.cc:176
const G4LorentzVector & Get4Momentum() const
tuple v
Definition: test.py:18
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
double theta() const
G4int GetZ() const
Definition: G4Isotope.hh:91
void SetEnergyChange(G4double anEnergy)
G4double GetTotalEnergy() const
static G4LENDManager * GetInstance()
G4double GetPDGMass() const
double mass
Definition: G4GIDI_mass.cc:44
void create_used_target_map()
Definition: G4LENDModel.cc:93
const G4Isotope * GetIsotope()
Definition: G4Nucleus.hh:119
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:170
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4Material * GetMaterial() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4bool RequestChangeOfVerboseLevel(G4int)
G4LENDManager * lend_manager
Definition: G4LENDModel.hh:80
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
void SetMomentumChange(const G4ThreeVector &aV)
G4double GetMass() const