ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LENDModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 
47 G4LENDModel::G4LENDModel( G4String name )
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 
63  lend_manager = G4LENDManager::GetInstance();
64 
65 }
66 
67 G4LENDModel::~G4LENDModel()
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 
77 void G4LENDModel::recreate_used_target_map()
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 
87  create_used_target_map();
88 
89 }
90 
91 
92 
93 void G4LENDModel::create_used_target_map()
94 {
95 
96  lend_manager->RequestChangeOfVerboseLevel( verboseLevel );
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
125  G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder();
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  DumpLENDTargetInfo();
152 }
153 
154 
155 
156 #include "G4IonTable.hh"
157 
158 G4HadFinalState * G4LENDModel::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
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 );
190  G4ReactionProduct theTarget( pd );
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 
253 G4HadFinalState* G4LENDModel::returnUnchanged(const G4HadProjectile& aTrack, G4HadFinalState* theResult ) {
254  if ( lend_manager->GetVerboseLevel() >= 1 ) {
255  G4String message;
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 
267 G4GIDI_target* G4LENDModel::get_target_from_map( G4int nuclear_code ) {
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 
275 void G4LENDModel::DumpLENDTargetInfo( G4bool force ) {
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 }