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 // 26 // 27 // This software was developed by Lawrence Liv 27 // This software was developed by Lawrence Livermore National Laboratory. 28 // 28 // 29 // Redistribution and use in source and binary 29 // Redistribution and use in source and binary forms, with or without 30 // modification, are permitted provided that t 30 // modification, are permitted provided that the following conditions are met: 31 // 31 // 32 // 1. Redistributions of source code must reta 32 // 1. Redistributions of source code must retain the above copyright notice, 33 // this list of conditions and the following 33 // this list of conditions and the following disclaimer. 34 // 2. Redistributions in binary form must repr 34 // 2. Redistributions in binary form must reproduce the above copyright notice, 35 // this list of conditions and the following 35 // this list of conditions and the following disclaimer in the documentation 36 // and/or other materials provided with the 36 // and/or other materials provided with the distribution. 37 // 3. The name of the author may not be used t 37 // 3. The name of the author may not be used to endorse or promote products 38 // derived from this software without specif 38 // derived from this software without specific prior written permission. 39 // 39 // 40 // THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``A 40 // THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED 41 // WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 41 // WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 42 // MERCHANTABILITY AND FITNESS FOR A PARTICULA 42 // MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO 43 // EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DI 43 // EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 44 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGE 44 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 45 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES 45 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 46 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AN 46 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 47 // WHETHER IN CONTRACT, STRICT LIABILITY, OR T 47 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 48 // OTHERWISE) ARISING IN ANY WAY OUT OF THE US 48 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF 49 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 49 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 50 // 50 // 51 // Copyright (c) 2006 The Regents of the Unive 51 // Copyright (c) 2006 The Regents of the University of California. 52 // All rights reserved. 52 // All rights reserved. 53 // UCRL-CODE-224807 53 // UCRL-CODE-224807 54 // 54 // 55 // 55 // >> 56 // $Id: G4FissionLibrary.cc 67966 2013-03-13 09:38:38Z gcosmo $ 56 // 57 // 57 // neutron_hp -- source file 58 // neutron_hp -- source file 58 // J.M. Verbeke, Jan-2007 59 // J.M. Verbeke, Jan-2007 59 // A low energy neutron-induced fission model. 60 // A low energy neutron-induced fission model. 60 // 61 // 61 62 62 #include "G4FissionLibrary.hh" 63 #include "G4FissionLibrary.hh" 63 #include "G4ParticleHPManager.hh" << 64 #include "G4SystemOfUnits.hh" 64 #include "G4SystemOfUnits.hh" 65 #include "G4PhysicsModelCatalog.hh" << 66 65 67 G4FissionLibrary::G4FissionLibrary() 66 G4FissionLibrary::G4FissionLibrary() 68 : G4ParticleHPFinalState(), theIsotope(0), t << 67 : G4NeutronHPFinalState(), theIsotope(0), targetMass(0.0) 69 { 68 { 70 hasXsec = false; 69 hasXsec = false; 71 fe=0; << 72 secID = G4PhysicsModelCatalog::GetModelID( " << 73 } 70 } 74 71 75 G4FissionLibrary::~G4FissionLibrary() 72 G4FissionLibrary::~G4FissionLibrary() 76 {} 73 {} 77 74 78 G4ParticleHPFinalState * G4FissionLibrary::New << 75 G4NeutronHPFinalState * G4FissionLibrary::New() 79 { 76 { 80 G4FissionLibrary * theNew = new G4FissionLib 77 G4FissionLibrary * theNew = new G4FissionLibrary; 81 return theNew; 78 return theNew; 82 } 79 } 83 80 84 //void G4FissionLibrary::Init (G4double A, G4d 81 //void G4FissionLibrary::Init (G4double A, G4double Z, G4String & dirName, G4String &) 85 void G4FissionLibrary::Init (G4double A, G4dou << 82 void G4FissionLibrary::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String &) 86 { 83 { 87 G4String tString = "/FS/"; 84 G4String tString = "/FS/"; 88 G4bool dbool; 85 G4bool dbool; 89 theIsotope = static_cast<G4int>(1000*Z+A); 86 theIsotope = static_cast<G4int>(1000*Z+A); 90 G4ParticleHPDataUsed aFile = theNames.GetNam << 87 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool); 91 G4String filename = aFile.GetName(); 88 G4String filename = aFile.GetName(); 92 89 93 if(!dbool) 90 if(!dbool) 94 { 91 { 95 hasAnyData = false; 92 hasAnyData = false; 96 hasFSData = false; 93 hasFSData = false; 97 hasXsec = false; 94 hasXsec = false; 98 return; 95 return; 99 } 96 } 100 //std::ifstream theData(filename, std::ios:: << 97 std::ifstream theData(filename, std::ios::in); 101 std::istringstream theData(std::ios::in); << 102 G4ParticleHPManager::GetInstance()->GetDataS << 103 98 104 // here it comes 99 // here it comes 105 G4int infoType, dataType; 100 G4int infoType, dataType; 106 hasFSData = false; 101 hasFSData = false; 107 while (theData >> infoType) // Loop checking << 102 while (theData >> infoType) 108 { 103 { 109 hasFSData = true; 104 hasFSData = true; 110 theData >> dataType; 105 theData >> dataType; 111 switch(infoType) 106 switch(infoType) 112 { 107 { 113 case 1: 108 case 1: 114 if(dataType==4) theNeutronAngularDis.I 109 if(dataType==4) theNeutronAngularDis.Init(theData); 115 if(dataType==5) thePromptNeutronEnDis. 110 if(dataType==5) thePromptNeutronEnDis.Init(theData); 116 if(dataType==12) theFinalStatePhotons. 111 if(dataType==12) theFinalStatePhotons.InitMean(theData); 117 if(dataType==14) theFinalStatePhotons. 112 if(dataType==14) theFinalStatePhotons.InitAngular(theData); 118 if(dataType==15) theFinalStatePhotons. 113 if(dataType==15) theFinalStatePhotons.InitEnergies(theData); 119 break; 114 break; 120 case 2: 115 case 2: 121 if(dataType==1) theFinalStateNeutrons. 116 if(dataType==1) theFinalStateNeutrons.InitMean(theData); 122 break; 117 break; 123 case 3: 118 case 3: 124 if(dataType==1) theFinalStateNeutrons. 119 if(dataType==1) theFinalStateNeutrons.InitDelayed(theData); 125 if(dataType==5) theDelayedNeutronEnDis 120 if(dataType==5) theDelayedNeutronEnDis.Init(theData); 126 break; 121 break; 127 case 4: 122 case 4: 128 if(dataType==1) theFinalStateNeutrons. 123 if(dataType==1) theFinalStateNeutrons.InitPrompt(theData); 129 break; 124 break; 130 case 5: 125 case 5: 131 if(dataType==1) theEnergyRelease.Init( 126 if(dataType==1) theEnergyRelease.Init(theData); 132 break; 127 break; 133 default: 128 default: 134 G4cout << "G4FissionLibrary::Init: unk 129 G4cout << "G4FissionLibrary::Init: unknown data type"<<dataType<<G4endl; 135 throw G4HadronicException(__FILE__, __ 130 throw G4HadronicException(__FILE__, __LINE__, "G4FissionLibrary::Init: unknown data type"); 136 break; 131 break; 137 } 132 } 138 } 133 } 139 targetMass = theFinalStateNeutrons.GetTarget 134 targetMass = theFinalStateNeutrons.GetTargetMass(); 140 //theData.close(); << 135 theData.close(); 141 } 136 } 142 137 143 G4HadFinalState* G4FissionLibrary::ApplyYourse << 138 G4HadFinalState * G4FissionLibrary::ApplyYourself(const G4HadProjectile & theTrack) 144 { 139 { >> 140 theResult.Clear(); 145 141 146 if ( theResult.Get() == NULL ) theResult.Put << 142 // prepare neutron 147 theResult.Get()->Clear(); << 148 << 149 // prepare neutron << 150 G4double eKinetic = theTrack.GetKineticEnerg 143 G4double eKinetic = theTrack.GetKineticEnergy(); 151 const G4HadProjectile* incidentParticle = &t << 144 const G4HadProjectile *incidentParticle = &theTrack; 152 G4ReactionProduct theNeutron(incidentParticl << 145 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) ); 153 theNeutron.SetMomentum(incidentParticle->Get << 146 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() ); 154 theNeutron.SetKineticEnergy(eKinetic); << 147 theNeutron.SetKineticEnergy( eKinetic ); 155 148 156 // prepare target << 149 // prepare target 157 G4Nucleus aNucleus; 150 G4Nucleus aNucleus; 158 G4ReactionProduct theTarget; 151 G4ReactionProduct theTarget; 159 G4ThreeVector neuVelo = (1./incidentParticle 152 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum(); 160 theTarget = aNucleus.GetBiasedThermalNucleus 153 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature()); 161 154 162 // set neutron and target in the FS classes << 155 // set neutron and target in the FS classes 163 //theNeutronAngularDis.SetNeutron(theNeutron << 156 theNeutronAngularDis.SetNeutron(theNeutron); 164 theNeutronAngularDis.SetProjectileRP(theNeut << 165 theNeutronAngularDis.SetTarget(theTarget); 157 theNeutronAngularDis.SetTarget(theTarget); 166 158 167 // boost to target rest system << 159 // boost to target rest system 168 theNeutron.Lorentz(theNeutron, -1*theTarget) 160 theNeutron.Lorentz(theNeutron, -1*theTarget); 169 161 170 eKinetic = theNeutron.GetKineticEnergy(); 162 eKinetic = theNeutron.GetKineticEnergy(); 171 163 172 // dice neutron and gamma multiplicities, en << 164 // dice neutron and gamma multiplicities, energies and momenta in Lab. @@ 173 // no energy conservation on an event-to-eve << 165 // no energy conservation on an event-to-event basis. we rely on the data to be ok. @@ 174 // also for mean, we rely on the consistency << 166 // also for mean, we rely on the consistency of the data. @@ 175 167 176 G4int nPrompt=0, gPrompt=0; 168 G4int nPrompt=0, gPrompt=0; 177 SampleMult(theTrack, &nPrompt, &gPrompt, eKi 169 SampleMult(theTrack, &nPrompt, &gPrompt, eKinetic); 178 170 179 // Build neutrons and add them to dynamic pa << 171 // Build neutrons and add them to dynamic particle vector 180 G4double momentum; 172 G4double momentum; 181 for(G4int i=0; i<nPrompt; i++) 173 for(G4int i=0; i<nPrompt; i++) 182 { 174 { 183 G4DynamicParticle * it = new G4DynamicPart 175 G4DynamicParticle * it = new G4DynamicParticle; 184 it->SetDefinition(G4Neutron::Neutron()); 176 it->SetDefinition(G4Neutron::Neutron()); 185 it->SetKineticEnergy(fe->getNeutronEnergy( << 177 it->SetKineticEnergy(getneng_(&i)*MeV); 186 momentum = it->GetTotalMomentum(); 178 momentum = it->GetTotalMomentum(); 187 G4ThreeVector temp(momentum*fe->getNeutron << 179 G4ThreeVector temp(momentum*getndircosu_(&i), 188 momentum*fe->getNeutron << 180 momentum*getndircosv_(&i), 189 momentum*fe->getNeutron << 181 momentum*getndircosw_(&i)); 190 it->SetMomentum( temp ); 182 it->SetMomentum( temp ); 191 // it->SetGlobalTime(fe->getNeutronAge(i)*s << 183 // it->SetGlobalTime(getnage_(&i)*second); 192 theResult.Get()->AddSecondary(it, secID); << 184 theResult.AddSecondary(it); 193 // G4cout <<"G4FissionLibrary::ApplyYoursel 185 // G4cout <<"G4FissionLibrary::ApplyYourself: energy of prompt neutron " << i << " = " << it->GetKineticEnergy()<<G4endl; 194 } 186 } 195 187 196 // Build gammas, lorentz transform them, and << 188 // Build gammas, lorentz transform them, and add them to dynamic particle vector 197 for(G4int i=0; i<gPrompt; i++) 189 for(G4int i=0; i<gPrompt; i++) 198 { 190 { 199 G4ReactionProduct * thePhoton = new G4Reac 191 G4ReactionProduct * thePhoton = new G4ReactionProduct; 200 thePhoton->SetDefinition(G4Gamma::Gamma()) 192 thePhoton->SetDefinition(G4Gamma::Gamma()); 201 thePhoton->SetKineticEnergy(fe->getPhotonE << 193 thePhoton->SetKineticEnergy(getpeng_(&i)*MeV); 202 momentum = thePhoton->GetTotalMomentum(); 194 momentum = thePhoton->GetTotalMomentum(); 203 G4ThreeVector temp(momentum*fe->getPhotonD << 195 G4ThreeVector temp(momentum*getpdircosu_(&i), 204 momentum*fe->getPhotonD << 196 momentum*getpdircosv_(&i), 205 momentum*fe->getPhotonD << 197 momentum*getpdircosw_(&i)); 206 thePhoton->SetMomentum( temp ); 198 thePhoton->SetMomentum( temp ); 207 thePhoton->Lorentz(*thePhoton, -1.*theTarg 199 thePhoton->Lorentz(*thePhoton, -1.*theTarget); 208 200 209 G4DynamicParticle * it = new G4DynamicPart 201 G4DynamicParticle * it = new G4DynamicParticle; 210 it->SetDefinition(thePhoton->GetDefinition 202 it->SetDefinition(thePhoton->GetDefinition()); 211 it->SetMomentum(thePhoton->GetMomentum()); 203 it->SetMomentum(thePhoton->GetMomentum()); 212 // it->SetGlobalTime(fe->getPhotonAge(i)*se << 204 // it->SetGlobalTime(getpage_(&i)*second); 213 // G4cout <<"G4FissionLibrary::ApplyYoursel 205 // G4cout <<"G4FissionLibrary::ApplyYourself: energy of prompt photon " << i << " = " << it->GetKineticEnergy()<<G4endl; 214 theResult.Get()->AddSecondary(it, secID); << 206 theResult.AddSecondary(it); 215 delete thePhoton; 207 delete thePhoton; 216 } 208 } 217 // G4cout <<"G4FissionLibrary::ApplyYourself: 209 // G4cout <<"G4FissionLibrary::ApplyYourself: Number of secondaries = "<<theResult.GetNumberOfSecondaries()<< G4endl; 218 // G4cout <<"G4FissionLibrary::ApplyYourself: 210 // G4cout <<"G4FissionLibrary::ApplyYourself: Number of induced prompt neutron = "<<nPrompt<<G4endl; 219 // G4cout <<"G4FissionLibrary::ApplyYourself: 211 // G4cout <<"G4FissionLibrary::ApplyYourself: Number of induced prompt photons = "<<gPrompt<<G4endl; 220 212 221 // finally deal with local energy deposition << 213 // finally deal with local energy depositions. 222 G4double eDepByFragments = theEnergyRelease. 214 G4double eDepByFragments = theEnergyRelease.GetFragmentKinetic(); 223 theResult.Get()->SetLocalEnergyDeposit(eDepB << 215 theResult.SetLocalEnergyDeposit(eDepByFragments); 224 // G4cout << "G4FissionLibrary::local energy 216 // G4cout << "G4FissionLibrary::local energy deposit" << eDepByFragments<<G4endl; 225 // clean up the primary neutron << 217 // clean up the primary neutron 226 theResult.Get()->SetStatusChange(stopAndKill << 218 theResult.SetStatusChange(stopAndKill); 227 return theResult.Get(); << 219 return &theResult; 228 } 220 } 229 221 230 void G4FissionLibrary::SampleMult(const G4HadP 222 void G4FissionLibrary::SampleMult(const G4HadProjectile & theTrack, G4int* nPrompt, 231 G4int* gPro 223 G4int* gPrompt, G4double eKinetic) 232 { 224 { 233 G4double promptNeutronMulti = 0; 225 G4double promptNeutronMulti = 0; 234 promptNeutronMulti = theFinalStateNeutrons. 226 promptNeutronMulti = theFinalStateNeutrons.GetPrompt(eKinetic); // prompt nubar from Geant 235 G4double delayedNeutronMulti = 0; 227 G4double delayedNeutronMulti = 0; 236 delayedNeutronMulti = theFinalStateNeutrons 228 delayedNeutronMulti = theFinalStateNeutrons.GetDelayed(eKinetic); // delayed nubar from Geant 237 229 238 G4double time = theTrack.GetGlobalTime()/se 230 G4double time = theTrack.GetGlobalTime()/second; 239 G4double totalNeutronMulti = theFinalStateN << 240 if(delayedNeutronMulti==0&&promptNeutronMul 231 if(delayedNeutronMulti==0&&promptNeutronMulti==0) { 241 // no data for prompt and delayed neutron 232 // no data for prompt and delayed neutrons in Geant 242 // but there is perhaps data for the tota 233 // but there is perhaps data for the total neutron multiplicity, in which case 243 // we use it for prompt neutron emission 234 // we use it for prompt neutron emission 244 if (fe != 0) delete fe; << 235 G4double totalNeutronMulti = theFinalStateNeutrons.GetMean(eKinetic); 245 fe = new G4fissionEvent(theIsotope, time, << 236 genfissevt_(&theIsotope, &time, &totalNeutronMulti, &eKinetic); 246 } else { 237 } else { 247 // prompt nubar != 0 || delayed nubar != 238 // prompt nubar != 0 || delayed nubar != 0 248 if (fe != 0) delete fe; << 239 genfissevt_(&theIsotope, &time, &promptNeutronMulti, &eKinetic); 249 fe = new G4fissionEvent(theIsotope, time, << 250 } 240 } 251 *nPrompt = fe->getNeutronNu(); << 241 *nPrompt = getnnu_(); 252 if (*nPrompt == -1) *nPrompt = 0; // the fi 242 if (*nPrompt == -1) *nPrompt = 0; // the fission library libFission.a has no data for neutrons 253 *gPrompt = fe->getPhotonNu(); << 243 *gPrompt = getpnu_(); 254 if (*gPrompt == -1) *gPrompt = 0; // the fi 244 if (*gPrompt == -1) *gPrompt = 0; // the fission library libFission.a has no data for gammas 255 } 245 } 256 246 257 247