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" 44 #include "G4PhysicsModelCatalog.hh" 45 45 46 double MyRNG(void*) { return G4Random::getThe 46 double MyRNG(void*) { return G4Random::getTheEngine()->flat(); } 47 47 48 G4LENDModel::G4LENDModel( G4String name ) 48 G4LENDModel::G4LENDModel( G4String name ) 49 :G4HadronicInteraction( name ), secID(-1) 49 :G4HadronicInteraction( name ), secID(-1) 50 { 50 { 51 51 52 proj = NULL; //will be set in an inherited 52 proj = NULL; //will be set in an inherited class 53 53 54 SetMinEnergy( 0.*eV ); 54 SetMinEnergy( 0.*eV ); 55 SetMaxEnergy( 20.*MeV ); 55 SetMaxEnergy( 20.*MeV ); 56 56 57 //default_evaluation = "endl99"; 57 //default_evaluation = "endl99"; 58 //default_evaluation = "ENDF.B-VII.0"; 58 //default_evaluation = "ENDF.B-VII.0"; 59 default_evaluation = "ENDF/BVII.1"; 59 default_evaluation = "ENDF/BVII.1"; 60 60 61 allow_nat = false; 61 allow_nat = false; 62 allow_any = false; 62 allow_any = false; 63 63 64 lend_manager = G4LENDManager::GetInstance() 64 lend_manager = G4LENDManager::GetInstance(); 65 65 66 secID = G4PhysicsModelCatalog::GetModelID( 66 secID = G4PhysicsModelCatalog::GetModelID( "model_" + GetModelName() ); 67 } 67 } 68 68 69 G4LENDModel::~G4LENDModel() 69 G4LENDModel::~G4LENDModel() 70 { 70 { 71 for ( std::map< G4int , G4LENDUsedTarget* > 71 for ( std::map< G4int , G4LENDUsedTarget* >::iterator 72 it = usedTarget_map.begin() ; it != u 72 it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ ) 73 { 73 { 74 delete it->second; 74 delete it->second; 75 } 75 } 76 } 76 } 77 77 78 78 79 void G4LENDModel::recreate_used_target_map() 79 void G4LENDModel::recreate_used_target_map() 80 { 80 { 81 81 82 for ( std::map< G4int , G4LENDUsedTarget* > 82 for ( std::map< G4int , G4LENDUsedTarget* >::iterator 83 it = usedTarget_map.begin() ; it != u 83 it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ ) 84 { 84 { 85 delete it->second; 85 delete it->second; 86 } 86 } 87 usedTarget_map.clear(); 87 usedTarget_map.clear(); 88 88 89 create_used_target_map(); 89 create_used_target_map(); 90 90 91 } 91 } 92 92 93 93 94 94 95 void G4LENDModel::create_used_target_map() 95 void G4LENDModel::create_used_target_map() 96 { 96 { 97 97 98 lend_manager->RequestChangeOfVerboseLevel( 98 lend_manager->RequestChangeOfVerboseLevel( verboseLevel ); 99 99 100 std::size_t numberOfElements = G4Element::G << 100 size_t numberOfElements = G4Element::GetNumberOfElements(); 101 static const G4ElementTable* theElementTabl 101 static const G4ElementTable* theElementTable = G4Element::GetElementTable(); 102 102 103 for ( std::size_t i = 0 ; i < numberOfEleme << 103 for ( size_t i = 0 ; i < numberOfElements ; ++i ) 104 { 104 { 105 105 106 const G4Element* anElement = (*theElemen 106 const G4Element* anElement = (*theElementTable)[i]; 107 G4int numberOfIsotope = (G4int)anElement << 107 G4int numberOfIsotope = anElement->GetNumberOfIsotopes(); 108 108 109 if ( numberOfIsotope > 0 ) 109 if ( numberOfIsotope > 0 ) 110 { 110 { 111 // User Defined Abundances 111 // User Defined Abundances 112 for ( G4int i_iso = 0 ; i_iso < numbe << 112 for ( G4int i_iso = 0 ; i_iso < numberOfIsotope ; i_iso++ ) 113 { 113 { 114 G4int iZ = anElement->GetIsotope( 114 G4int iZ = anElement->GetIsotope( i_iso )->GetZ(); 115 G4int iA = anElement->GetIsotope( 115 G4int iA = anElement->GetIsotope( i_iso )->GetN(); 116 G4int iIsomer = anElement->GetIsot 116 G4int iIsomer = anElement->GetIsotope( i_iso )->Getm(); 117 117 118 G4LENDUsedTarget* aTarget = new G4 118 G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iA , iIsomer ); 119 if ( allow_nat == true ) aTarget-> 119 if ( allow_nat == true ) aTarget->AllowNat(); 120 if ( allow_any == true ) aTarget-> 120 if ( allow_any == true ) aTarget->AllowAny(); 121 usedTarget_map.insert( std::pair< 121 usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iA , iIsomer ) , aTarget ) ); 122 } 122 } 123 } 123 } 124 else 124 else 125 { 125 { 126 // Natural Abundances 126 // Natural Abundances 127 G4NistElementBuilder* nistElementBuil 127 G4NistElementBuilder* nistElementBuild = lend_manager->GetNistElementBuilder(); 128 G4int iZ = int ( anElement->GetZ() ); 128 G4int iZ = int ( anElement->GetZ() ); 129 //G4cout << nistElementBuild->GetNumb 129 //G4cout << nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ) << G4endl; 130 G4int numberOfNistIso = nistElementBu 130 G4int numberOfNistIso = nistElementBuild->GetNumberOfNistIsotopes( int ( anElement->GetZ() ) ); 131 131 132 for ( G4int ii = 0 ; ii < numberOfNis 132 for ( G4int ii = 0 ; ii < numberOfNistIso ; ii++ ) 133 { 133 { 134 //G4cout << nistElementBuild->GetI 134 //G4cout << nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + i ) << G4endl; 135 if ( nistElementBuild->GetIsotopeA 135 if ( nistElementBuild->GetIsotopeAbundance( iZ , nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii ) > 0 ) 136 { 136 { 137 G4int iMass = nistElementBuild- 137 G4int iMass = nistElementBuild->GetNistFirstIsotopeN( iZ ) + ii; 138 //G4cout << iZ << " " << nistEl 138 //G4cout << iZ << " " << nistElementBuild->GetNistFirstIsotopeN( iZ ) + i << " " << nistElementBuild->GetIsotopeAbundance ( iZ , iMass ) << G4endl; 139 G4int iIsomer = 0; 139 G4int iIsomer = 0; 140 140 141 G4LENDUsedTarget* aTarget = new 141 G4LENDUsedTarget* aTarget = new G4LENDUsedTarget ( proj , default_evaluation , iZ , iMass ); 142 if ( allow_nat == true ) aTarge 142 if ( allow_nat == true ) aTarget->AllowNat(); 143 if ( allow_any == true ) aTarge 143 if ( allow_any == true ) aTarget->AllowAny(); 144 usedTarget_map.insert( std::pai 144 usedTarget_map.insert( std::pair< G4int , G4LENDUsedTarget* > ( lend_manager->GetNucleusEncoding( iZ , iMass , iIsomer ) , aTarget ) ); 145 145 146 } 146 } 147 147 148 } 148 } 149 149 150 } 150 } 151 } 151 } 152 152 153 DumpLENDTargetInfo(); 153 DumpLENDTargetInfo(); 154 } 154 } 155 155 156 156 157 157 158 #include "G4IonTable.hh" 158 #include "G4IonTable.hh" 159 159 160 G4HadFinalState * G4LENDModel::ApplyYourself(c 160 G4HadFinalState * G4LENDModel::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg ) 161 { 161 { 162 162 163 G4double temp = aTrack.GetMaterial()->GetTe 163 G4double temp = aTrack.GetMaterial()->GetTemperature(); 164 164 165 //G4int iZ = int ( aTarg.GetZ() ); 165 //G4int iZ = int ( aTarg.GetZ() ); 166 //G4int iA = int ( aTarg.GetN() ); 166 //G4int iA = int ( aTarg.GetN() ); 167 //migrate to integer A and Z (GetN_asInt re 167 //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 168 G4int iZ = aTarg.GetZ_asInt(); 168 G4int iZ = aTarg.GetZ_asInt(); 169 G4int iA = aTarg.GetA_asInt(); 169 G4int iA = aTarg.GetA_asInt(); 170 G4int iM = 0; 170 G4int iM = 0; 171 if ( aTarg.GetIsotope() != NULL ) { 171 if ( aTarg.GetIsotope() != NULL ) { 172 iM = aTarg.GetIsotope()->Getm(); 172 iM = aTarg.GetIsotope()->Getm(); 173 } 173 } 174 174 175 G4double ke = aTrack.GetKineticEnergy(); 175 G4double ke = aTrack.GetKineticEnergy(); 176 176 177 G4HadFinalState* theResult = new G4HadFinal 177 G4HadFinalState* theResult = new G4HadFinalState(); 178 178 179 G4GIDI_target* aTarget = usedTarget_map.fin 179 G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA , iM ) )->second->GetTarget(); 180 180 181 G4double aMu = aTarget->getElasticFinalStat 181 G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, NULL, NULL ); 182 182 183 G4double phi = twopi*G4UniformRand(); 183 G4double phi = twopi*G4UniformRand(); 184 G4double theta = std::acos( aMu ); 184 G4double theta = std::acos( aMu ); 185 //G4double sinth = std::sin( theta ); 185 //G4double sinth = std::sin( theta ); 186 186 187 G4ReactionProduct theNeutron( aTrack.GetDef 187 G4ReactionProduct theNeutron( aTrack.GetDefinition() ); 188 theNeutron.SetMomentum( aTrack.Get4Momentum 188 theNeutron.SetMomentum( aTrack.Get4Momentum().vect() ); 189 theNeutron.SetKineticEnergy( ke ); 189 theNeutron.SetKineticEnergy( ke ); 190 190 191 G4ParticleDefinition* pd = G4IonTable::GetI 191 G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( iZ , iA , iM ); 192 G4ReactionProduct theTarget( pd ); 192 G4ReactionProduct theTarget( pd ); 193 193 194 G4double mass = pd->GetPDGMass(); 194 G4double mass = pd->GetPDGMass(); 195 195 196 // add Thermal motion 196 // add Thermal motion 197 G4double kT = k_Boltzmann*temp; 197 G4double kT = k_Boltzmann*temp; 198 G4ThreeVector v ( G4RandGauss::shoot() * st 198 G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass ) 199 , G4RandGauss::shoot() * st 199 , G4RandGauss::shoot() * std::sqrt( kT*mass ) 200 , G4RandGauss::shoot() * st 200 , G4RandGauss::shoot() * std::sqrt( kT*mass ) ); 201 201 202 theTarget.SetMomentum( v ); 202 theTarget.SetMomentum( v ); 203 203 204 204 205 G4ThreeVector the3Neutron = theNeutron.Ge 205 G4ThreeVector the3Neutron = theNeutron.GetMomentum(); 206 G4double nEnergy = theNeutron.GetTotalEne 206 G4double nEnergy = theNeutron.GetTotalEnergy(); 207 G4ThreeVector the3Target = theTarget.GetM 207 G4ThreeVector the3Target = theTarget.GetMomentum(); 208 G4double tEnergy = theTarget.GetTotalEner 208 G4double tEnergy = theTarget.GetTotalEnergy(); 209 G4ReactionProduct theCMS; 209 G4ReactionProduct theCMS; 210 G4double totE = nEnergy+tEnergy; 210 G4double totE = nEnergy+tEnergy; 211 G4ThreeVector the3CMS = the3Target+the3Ne 211 G4ThreeVector the3CMS = the3Target+the3Neutron; 212 theCMS.SetMomentum(the3CMS); 212 theCMS.SetMomentum(the3CMS); 213 G4double cmsMom = std::sqrt(the3CMS*the3C 213 G4double cmsMom = std::sqrt(the3CMS*the3CMS); 214 G4double sqrts = std::sqrt((totE-cmsMom)* 214 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom)); 215 theCMS.SetMass(sqrts); 215 theCMS.SetMass(sqrts); 216 theCMS.SetTotalEnergy(totE); 216 theCMS.SetTotalEnergy(totE); 217 217 218 theNeutron.Lorentz(theNeutron, theCMS); 218 theNeutron.Lorentz(theNeutron, theCMS); 219 theTarget.Lorentz(theTarget, theCMS); 219 theTarget.Lorentz(theTarget, theCMS); 220 G4double en = theNeutron.GetTotalMoment 220 G4double en = theNeutron.GetTotalMomentum(); // already in CMS. 221 G4ThreeVector cms3Mom=theNeutron.GetMom 221 G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS 222 G4double cms_theta=cms3Mom.theta(); 222 G4double cms_theta=cms3Mom.theta(); 223 G4double cms_phi=cms3Mom.phi(); 223 G4double cms_phi=cms3Mom.phi(); 224 G4ThreeVector tempVector; 224 G4ThreeVector tempVector; 225 tempVector.setX(std::cos(theta)*std::si 225 tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi) 226 +std::sin(theta)*std::c 226 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi) 227 -std::sin(theta)*std::s 227 -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) ); 228 tempVector.setY(std::cos(theta)*std::si 228 tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi) 229 +std::sin(theta)*std::c 229 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi) 230 +std::sin(theta)*std::s 230 +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) ); 231 tempVector.setZ(std::cos(theta)*std::co 231 tempVector.setZ(std::cos(theta)*std::cos(cms_theta) 232 -std::sin(theta)*std::c 232 -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) ); 233 tempVector *= en; 233 tempVector *= en; 234 theNeutron.SetMomentum(tempVector); 234 theNeutron.SetMomentum(tempVector); 235 theTarget.SetMomentum(-tempVector); 235 theTarget.SetMomentum(-tempVector); 236 G4double tP = theTarget.GetTotalMomentu 236 G4double tP = theTarget.GetTotalMomentum(); 237 G4double tM = theTarget.GetMass(); 237 G4double tM = theTarget.GetMass(); 238 theTarget.SetTotalEnergy(std::sqrt((tP+ 238 theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM)); 239 theNeutron.Lorentz(theNeutron, -1.*theC 239 theNeutron.Lorentz(theNeutron, -1.*theCMS); 240 theTarget.Lorentz(theTarget, -1.*theCMS 240 theTarget.Lorentz(theTarget, -1.*theCMS); 241 241 242 theResult->SetEnergyChange(theNeutron.Get 242 theResult->SetEnergyChange(theNeutron.GetKineticEnergy()); 243 theResult->SetMomentumChange(theNeutron.G 243 theResult->SetMomentumChange(theNeutron.GetMomentum().unit()); 244 G4DynamicParticle* theRecoil = new G4Dyna 244 G4DynamicParticle* theRecoil = new G4DynamicParticle; 245 245 246 theRecoil->SetDefinition( G4IonTable::Get 246 theRecoil->SetDefinition( G4IonTable::GetIonTable()->GetIon( iZ , iA , iM , iZ ) ); 247 theRecoil->SetMomentum( theTarget.GetMome 247 theRecoil->SetMomentum( theTarget.GetMomentum() ); 248 248 249 theResult->AddSecondary( theRecoil ); 249 theResult->AddSecondary( theRecoil ); 250 250 251 return theResult; 251 return theResult; 252 252 253 } 253 } 254 254 255 G4HadFinalState* G4LENDModel::returnUnchanged( 255 G4HadFinalState* G4LENDModel::returnUnchanged(const G4HadProjectile& aTrack, G4HadFinalState* theResult ) { 256 if ( lend_manager->GetVerboseLevel() >= 1 ) 256 if ( lend_manager->GetVerboseLevel() >= 1 ) { 257 G4String message; 257 G4String message; 258 message = "Produce unchanged final state 258 message = "Produce unchanged final state is requested in "; 259 message += this->GetModelName(); 259 message += this->GetModelName(); 260 message += ". Cross section and model li 260 message += ". Cross section and model likely have an inconsistency."; 261 G4Exception( "G4LENDModel::returnUnchang 261 G4Exception( "G4LENDModel::returnUnchanged(,)" , "LENDModel-01" , JustWarning , 262 message ); 262 message ); 263 } 263 } 264 theResult->SetEnergyChange( aTrack.GetKinet 264 theResult->SetEnergyChange( aTrack.GetKineticEnergy() ); 265 theResult->SetMomentumChange( aTrack.Get4Mo 265 theResult->SetMomentumChange( aTrack.Get4Momentum().getV().unit() ); 266 return theResult; 266 return theResult; 267 } 267 } 268 268 269 G4GIDI_target* G4LENDModel::get_target_from_ma 269 G4GIDI_target* G4LENDModel::get_target_from_map( G4int nuclear_code ) { 270 G4GIDI_target* target = NULL; 270 G4GIDI_target* target = NULL; 271 if ( usedTarget_map.find( nuclear_code ) != 271 if ( usedTarget_map.find( nuclear_code ) != usedTarget_map.end() ) { 272 target = usedTarget_map.find( nuclear_co 272 target = usedTarget_map.find( nuclear_code )->second->GetTarget(); 273 } 273 } 274 return target; 274 return target; 275 } 275 } 276 276 277 void G4LENDModel::DumpLENDTargetInfo( G4bool f 277 void G4LENDModel::DumpLENDTargetInfo( G4bool force ) { 278 278 279 if ( lend_manager->GetVerboseLevel() >= 1 | 279 if ( lend_manager->GetVerboseLevel() >= 1 || force ) { 280 if ( usedTarget_map.size() == 0 ) create 280 if ( usedTarget_map.size() == 0 ) create_used_target_map(); 281 G4cout << "Dumping UsedTarget of " << Ge 281 G4cout << "Dumping UsedTarget of " << GetModelName() << " for " << proj->GetParticleName() << G4endl; 282 G4cout << "Requested Evaluation, Z , A - 282 G4cout << "Requested Evaluation, Z , A -> Actual Evaluation, Z , A(0=Nat) " << G4endl; 283 for ( std::map< G4int , G4LENDUsedTarget 283 for ( std::map< G4int , G4LENDUsedTarget* >::iterator 284 it = usedTarget_map.begin() ; it != u 284 it = usedTarget_map.begin() ; it != usedTarget_map.end() ; it ++ ) { 285 G4cout 285 G4cout 286 << " " << it->second->GetWantedEvalua 286 << " " << it->second->GetWantedEvaluation() 287 << ", " << it->second->GetWantedZ() 287 << ", " << it->second->GetWantedZ() 288 << ", " << it->second->GetWantedA() 288 << ", " << it->second->GetWantedA() 289 << " -> " << it->second->GetActualEva 289 << " -> " << it->second->GetActualEvaluation() 290 << ", " << it->second->GetActualZ() 290 << ", " << it->second->GetActualZ() 291 << ", " << it->second->GetActualA() 291 << ", " << it->second->GetActualA() 292 << G4endl; 292 << G4endl; 293 } 293 } 294 } 294 } 295 } 295 } 296 296