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 // 27 // // 28 // File: G4LENDInelastic.cc 28 // File: G4LENDInelastic.cc // 29 // Date: 24 March 2020 29 // Date: 24 March 2020 // 30 // Author: Dennis Wright 30 // Author: Dennis Wright // 31 // 31 // // 32 // Description: model for inelastic scatterin 32 // Description: model for inelastic scattering of neutrons, light ions and // 33 // gammas at energies of order 2 33 // gammas at energies of order 20 MeV and lower. // 34 // This model uses GIDI particle 34 // This model uses GIDI particle data which are stored mostly // 35 // in spectrum mode. In this mo 35 // in spectrum mode. In this mode, spectra are reproduced // 36 // for each possible particle ty 36 // for each possible particle type which can result from a // 37 // given interaction. Unlike Ge 37 // given interaction. Unlike Geant4, this is done without // 38 // consideration of event-by-eve 38 // consideration of event-by-event conservation rules. // 39 // Indeed, forcing such conserva 39 // Indeed, forcing such conservation on GIDI output products // 40 // introduces correlations and d 40 // introduces correlations and distortions in the resulting // 41 // spectra which are not present 41 // spectra which are not present in the data. // 42 // 42 // // 43 // In order to use GIDI data wit 43 // In order to use GIDI data within the Geant4 framework, a // 44 // minimal event-by-event baryon 44 // minimal event-by-event baryon number conservation is // 45 // enforced which allows deviati 45 // enforced which allows deviations of up to 1 GeV without // 46 // giving warnings. Neither cha 46 // giving warnings. Neither charge, nor energy, nor momentum // 47 // conservation is enforced. Un 47 // conservation is enforced. Under this scheme, light // 48 // fragment (n, p, d, t, alpha) 48 // fragment (n, p, d, t, alpha) spectra are well reproduced // 49 // after a large number of event 49 // after a large number of events. Charge and energy // 50 // conservation also approach th 50 // conservation also approach their event-by-event values in // 51 // this limit. The mass, charge 51 // this limit. The mass, charge and energy distributions of // 52 // large fragments, however, are 52 // large fragments, however, are not expected to reproduce the // 53 // data very well. This is a re 53 // data very well. This is a result of forcing the crude // 54 // baryon number conservation an 54 // baryon number conservation and ensuring that the light // 55 // fragment spectra are correct. 55 // fragment spectra are correct. // 56 // 56 // // 57 ////////////////////////////////////////////// 57 //////////////////////////////////////////////////////////////////////////////// 58 58 59 #include "G4LENDInelastic.hh" 59 #include "G4LENDInelastic.hh" 60 #include "G4SystemOfUnits.hh" 60 #include "G4SystemOfUnits.hh" 61 #include "G4Nucleus.hh" 61 #include "G4Nucleus.hh" 62 #include "G4IonTable.hh" 62 #include "G4IonTable.hh" 63 #include <algorithm> 63 #include <algorithm> 64 #include <random> 64 #include <random> 65 65 66 G4HadFinalState* G4LENDInelastic::ApplyYoursel 66 G4HadFinalState* G4LENDInelastic::ApplyYourself(const G4HadProjectile& aTrack, 67 67 G4Nucleus& aTarg) 68 { 68 { 69 G4ThreeVector projMom = aTrack.Get4Momentum( 69 G4ThreeVector projMom = aTrack.Get4Momentum().vect(); 70 G4double temp = aTrack.GetMaterial()->GetTem 70 G4double temp = aTrack.GetMaterial()->GetTemperature(); 71 71 72 G4int iZ = aTarg.GetZ_asInt(); 72 G4int iZ = aTarg.GetZ_asInt(); 73 G4int iA = aTarg.GetA_asInt(); 73 G4int iA = aTarg.GetA_asInt(); 74 G4int iM = 0; 74 G4int iM = 0; 75 if (aTarg.GetIsotope() != nullptr) iM = aTar 75 if (aTarg.GetIsotope() != nullptr) iM = aTarg.GetIsotope()->Getm(); 76 76 77 G4double ke = aTrack.GetKineticEnergy(); 77 G4double ke = aTrack.GetKineticEnergy(); 78 78 79 G4HadFinalState* theResult = &theParticleCha 79 G4HadFinalState* theResult = &theParticleChange; 80 theResult->Clear(); 80 theResult->Clear(); 81 81 82 G4GIDI_target* aGIDITarget = 82 G4GIDI_target* aGIDITarget = 83 get_target_from_map(lend_manager->GetNucle 83 get_target_from_map(lend_manager->GetNucleusEncoding(iZ, iA, iM) ); 84 if (aGIDITarget == nullptr) { 84 if (aGIDITarget == nullptr) { 85 // G4cout << " No target found " << G4endl; 85 // G4cout << " No target found " << G4endl; 86 theParticleChange.Clear(); 86 theParticleChange.Clear(); 87 theParticleChange.SetStatusChange(isAlive) 87 theParticleChange.SetStatusChange(isAlive); 88 theParticleChange.SetEnergyChange(aTrack.G 88 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy()); 89 theParticleChange.SetMomentumChange(aTrack 89 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 90 return &theParticleChange; 90 return &theParticleChange; 91 } 91 } 92 // return returnUnchanged(aTrack, theResult); 92 // return returnUnchanged(aTrack, theResult); 93 93 94 // Get GIDI final state products for givern 94 // Get GIDI final state products for givern target and projectile 95 G4int loop(0); 95 G4int loop(0); 96 G4int loopMax = 1000; 96 G4int loopMax = 1000; 97 std::vector<G4GIDI_Product>* products; 97 std::vector<G4GIDI_Product>* products; 98 do { 98 do { 99 products = aGIDITarget->getOthersFinalStat 99 products = aGIDITarget->getOthersFinalState(ke*MeV, temp, MyRNG, NULL); 100 loop++; 100 loop++; 101 } while (products == nullptr && loop < loopM 101 } while (products == nullptr && loop < loopMax); 102 102 103 // G4LENDInelastic accepts all light fragmen 103 // G4LENDInelastic accepts all light fragments and gammas from GIDI (A < 5) 104 // and removes any heavy fragments which cau 104 // and removes any heavy fragments which cause large baryon number violation. 105 // Charge and energy non-conservation still 105 // Charge and energy non-conservation still occur, but over a large number 106 // of events, this improves on average. 106 // of events, this improves on average. 107 107 108 if (loop > loopMax - 1) { 108 if (loop > loopMax - 1) { 109 // G4cout << " too many loops, return initi 109 // G4cout << " too many loops, return initial state " << G4endl; 110 110 111 theParticleChange.Clear(); 111 theParticleChange.Clear(); 112 theParticleChange.SetStatusChange(isAlive) 112 theParticleChange.SetStatusChange(isAlive); 113 theParticleChange.SetEnergyChange(aTrack.G 113 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy()); 114 theParticleChange.SetMomentumChange(aTrack 114 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 115 return &theParticleChange; 115 return &theParticleChange; 116 116 117 // if (aTrack.GetDefinition() == G4Proton:: 117 // if (aTrack.GetDefinition() == G4Proton::Proton() || 118 // aTrack.GetDefinition() == G4Neutron: 118 // aTrack.GetDefinition() == G4Neutron::Neutron() ) { 119 // theResult = preco->ApplyYourself(aTra 119 // theResult = preco->ApplyYourself(aTrack, aTarg); 120 // } else { 120 // } else { 121 // theResult = returnUnchanged(aTrack, t 121 // theResult = returnUnchanged(aTrack, theResult); 122 // } 122 // } 123 123 124 } else { 124 } else { 125 G4int iTotZ = iZ + aTrack.GetDefinition()- 125 G4int iTotZ = iZ + aTrack.GetDefinition()->GetAtomicNumber(); 126 G4int iTotA = iA + aTrack.GetDefinition()- 126 G4int iTotA = iA + aTrack.GetDefinition()->GetAtomicMass(); 127 127 128 // Loop over GIDI products and separate li 128 // Loop over GIDI products and separate light from heavy fragments 129 G4int GZtot(0); 129 G4int GZtot(0); 130 G4int GAtot(0); 130 G4int GAtot(0); 131 G4int productA(0); 131 G4int productA(0); 132 G4int productZ(0); 132 G4int productZ(0); 133 std::vector<G4int> lightProductIndex; 133 std::vector<G4int> lightProductIndex; 134 std::vector<G4int> heavyProductIndex; 134 std::vector<G4int> heavyProductIndex; 135 for (G4int i = 0; i < int( products->size( 135 for (G4int i = 0; i < int( products->size() ); i++ ) { 136 productA = (*products)[i].A; 136 productA = (*products)[i].A; 137 if (productA < 5) { 137 if (productA < 5) { 138 lightProductIndex.push_back(i); 138 lightProductIndex.push_back(i); 139 GZtot += (*products)[i].Z; 139 GZtot += (*products)[i].Z; 140 GAtot += productA; 140 GAtot += productA; 141 } else { 141 } else { 142 heavyProductIndex.push_back(i); 142 heavyProductIndex.push_back(i); 143 } 143 } 144 } 144 } 145 145 146 // Randomize order of heavies to correct s 146 // Randomize order of heavies to correct somewhat for sampling bias 147 // std::random_shuffle(heavyProductIndex.b 147 // std::random_shuffle(heavyProductIndex.begin(), heavyProductIndex.end() ); 148 // std::cout << " Heavy product index befo 148 // std::cout << " Heavy product index before shuffle : " ; 149 // for (G4int i = 0; i < int(heavyProductI 149 // for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) std::cout << heavyProductIndex[i] << ", " ; 150 // std::cout << std::endl; 150 // std::cout << std::endl; 151 151 152 auto rng = std::default_random_engine {}; 152 auto rng = std::default_random_engine {}; 153 std::shuffle(heavyProductIndex.begin(), he 153 std::shuffle(heavyProductIndex.begin(), heavyProductIndex.end(), rng); 154 154 155 // std::cout << " Heavy product index afte 155 // std::cout << " Heavy product index after shuffle : " ; 156 // for (G4int i = 0; i < int(heavyProductI 156 // for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) std::cout << heavyProductIndex[i] << ", " ; 157 // std::cout << std::endl; 157 // std::cout << std::endl; 158 158 159 std::vector<G4int> savedHeavyIndex; 159 std::vector<G4int> savedHeavyIndex; 160 G4int itest(0); 160 G4int itest(0); 161 for (G4int i = 0; i < int(heavyProductInde 161 for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) { 162 itest = heavyProductIndex[i]; 162 itest = heavyProductIndex[i]; 163 productA = (*products)[itest].A; 163 productA = (*products)[itest].A; 164 productZ = (*products)[itest].Z; 164 productZ = (*products)[itest].Z; 165 if ((GAtot + productA <= iTotA) && (GZto 165 if ((GAtot + productA <= iTotA) && (GZtot + productZ <= iTotZ) ) { 166 savedHeavyIndex.push_back(itest); 166 savedHeavyIndex.push_back(itest); 167 GZtot += productZ; 167 GZtot += productZ; 168 GAtot += productA; 168 GAtot += productA; 169 } 169 } 170 } 170 } 171 171 172 /* 172 /* 173 G4cout << " saved light products = "; 173 G4cout << " saved light products = "; 174 for (G4int k = 0; k < int(lightProductInde 174 for (G4int k = 0; k < int(lightProductIndex.size() ); k++ ) { 175 itest = lightProductIndex[k]; 175 itest = lightProductIndex[k]; 176 G4cout << "(" << (*products)[itest].Z << 176 G4cout << "(" << (*products)[itest].Z << ", " << (*products)[itest].A << "), "; 177 } 177 } 178 G4cout << G4endl; 178 G4cout << G4endl; 179 179 180 G4cout << " saved heavy products = "; 180 G4cout << " saved heavy products = "; 181 for (G4int k = 0; k < int(savedHeavyIndex. 181 for (G4int k = 0; k < int(savedHeavyIndex.size() ); k++ ) { 182 itest = savedHeavyIndex[k]; 182 itest = savedHeavyIndex[k]; 183 G4cout << "(" << (*products)[itest].Z << 183 G4cout << "(" << (*products)[itest].Z << ", " << (*products)[itest].A << "), "; 184 } 184 } 185 G4cout << G4endl; 185 G4cout << G4endl; 186 */ 186 */ 187 // Now convert saved products to Geant4 pa 187 // Now convert saved products to Geant4 particles 188 // Note that, at least for heavy fragments 188 // Note that, at least for heavy fragments, GIDI masses and Geant4 masses 189 // have slightly different values. 189 // have slightly different values. 190 190 191 G4DynamicParticle* theSec = nullptr; 191 G4DynamicParticle* theSec = nullptr; 192 G4ThreeVector Psum; 192 G4ThreeVector Psum; 193 for (G4int i = 0; i < int(lightProductInde 193 for (G4int i = 0; i < int(lightProductIndex.size() ); i++) { 194 itest = lightProductIndex[i]; 194 itest = lightProductIndex[i]; 195 productZ = (*products)[itest].Z; 195 productZ = (*products)[itest].Z; 196 productA = (*products)[itest].A; 196 productA = (*products)[itest].A; 197 theSec = new G4DynamicParticle(); 197 theSec = new G4DynamicParticle(); 198 if (productA == 1 && productZ == 0) { 198 if (productA == 1 && productZ == 0) { 199 theSec->SetDefinition(G4Neutron::Neutr 199 theSec->SetDefinition(G4Neutron::Neutron() ); 200 } else if (productA == 1 && productZ == 200 } else if (productA == 1 && productZ == 1) { 201 theSec->SetDefinition(G4Proton::Proton 201 theSec->SetDefinition(G4Proton::Proton() ); 202 } else if (productA == 2 && productZ == 202 } else if (productA == 2 && productZ == 1) { 203 theSec->SetDefinition(G4Deuteron::Deut 203 theSec->SetDefinition(G4Deuteron::Deuteron() ); 204 } else if (productA == 3 && productZ == 204 } else if (productA == 3 && productZ == 1) { 205 theSec->SetDefinition(G4Triton::Triton 205 theSec->SetDefinition(G4Triton::Triton() ); 206 } else if (productA == 4 && productZ == 206 } else if (productA == 4 && productZ == 2) { 207 theSec->SetDefinition(G4Alpha::Alpha() 207 theSec->SetDefinition(G4Alpha::Alpha() ); 208 } else { 208 } else { 209 theSec->SetDefinition(G4Gamma::Gamma() 209 theSec->SetDefinition(G4Gamma::Gamma() ); 210 } 210 } 211 211 212 G4ThreeVector momentum((*products)[itest 212 G4ThreeVector momentum((*products)[itest].px*MeV, 213 (*products)[itest 213 (*products)[itest].py*MeV, 214 (*products)[itest 214 (*products)[itest].pz*MeV ); 215 Psum += momentum; 215 Psum += momentum; 216 theSec->SetMomentum(momentum); 216 theSec->SetMomentum(momentum); 217 // theResult->AddSecondary(theSec); 217 // theResult->AddSecondary(theSec); 218 theParticleChange.AddSecondary(theSec, s << 218 theParticleChange.AddSecondary(theSec); 219 } 219 } 220 220 221 G4int productM(0); 221 G4int productM(0); 222 for (G4int i = 0; i < int(savedHeavyIndex. 222 for (G4int i = 0; i < int(savedHeavyIndex.size() ); i++) { 223 itest = savedHeavyIndex[i]; 223 itest = savedHeavyIndex[i]; 224 productZ = (*products)[itest].Z; 224 productZ = (*products)[itest].Z; 225 productA = (*products)[itest].A; 225 productA = (*products)[itest].A; 226 productM = (*products)[itest].m; 226 productM = (*products)[itest].m; 227 theSec = new G4DynamicParticle(); 227 theSec = new G4DynamicParticle(); 228 theSec->SetDefinition(G4IonTable::GetIon 228 theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(productZ, 229 229 productA, 230 230 productM) ); 231 G4ThreeVector momentum((*products)[itest 231 G4ThreeVector momentum((*products)[itest].px*MeV, 232 (*products)[itest 232 (*products)[itest].py*MeV, 233 (*products)[itest 233 (*products)[itest].pz*MeV ); 234 Psum += momentum; 234 Psum += momentum; 235 theSec->SetMomentum(momentum); 235 theSec->SetMomentum(momentum); 236 // theResult->AddSecondary(theSec); 236 // theResult->AddSecondary(theSec); 237 theParticleChange.AddSecondary(theSec, s << 237 theParticleChange.AddSecondary(theSec); 238 } 238 } 239 239 240 // Create heavy fragment if necessary to t 240 // Create heavy fragment if necessary to try to balance A, Z 241 // Note: this step is only required to pre 241 // Note: this step is only required to prevent warnings at the process level 242 // where "catastrophic" non-conserva 242 // where "catastrophic" non-conservation tolerances are set to 1 GeV. 243 // The residual generated will not n 243 // The residual generated will not necessarily be the one that would 244 // occur in the actual reaction. 244 // occur in the actual reaction. 245 if (iTotA - GAtot > 1) { 245 if (iTotA - GAtot > 1) { 246 theSec = new G4DynamicParticle(); 246 theSec = new G4DynamicParticle(); 247 if (iTotZ == GZtot) { 247 if (iTotZ == GZtot) { 248 // Special case when a nucleus of only 248 // Special case when a nucleus of only neutrons is requested 249 // Violate charge conservation and set 249 // Violate charge conservation and set Z = 1 250 // G4cout << " Z = 1, A = "<< iTotA - 250 // G4cout << " Z = 1, A = "<< iTotA - GAtot << " created " << G4endl; 251 theSec->SetDefinition(G4IonTable::GetI 251 theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(1, iTotA-GAtot, 0) ); 252 } else { 252 } else { 253 theSec->SetDefinition(G4IonTable::GetI 253 theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(iTotZ-GZtot, iTotA-GAtot, 0) ); 254 } 254 } 255 theSec->SetMomentum(projMom - Psum); 255 theSec->SetMomentum(projMom - Psum); 256 // theResult->AddSecondary(theSec); 256 // theResult->AddSecondary(theSec); 257 theParticleChange.AddSecondary(theSec, s << 257 theParticleChange.AddSecondary(theSec); 258 } 258 } 259 } // loop OK 259 } // loop OK 260 260 261 delete products; 261 delete products; 262 // theResult->SetStatusChange( stopAndKill ); 262 // theResult->SetStatusChange( stopAndKill ); 263 theParticleChange.SetStatusChange( stopAndKi 263 theParticleChange.SetStatusChange( stopAndKill ); 264 // return theResult; 264 // return theResult; 265 return &theParticleChange; 265 return &theParticleChange; 266 } 266 } 267 267