Geant4 Cross Reference |
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 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron transport model. 29 // 30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi) 31 // 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 33 // V. Ivanchenko, July-2023 Basic revision of particle HP classes 34 // 35 36 #include "G4ParticleHPChannelList.hh" 37 38 #include "G4Element.hh" 39 #include "G4HadFinalState.hh" 40 #include "G4HadProjectile.hh" 41 #include "G4ParticleHPFinalState.hh" 42 #include "G4ParticleHPManager.hh" 43 #include "G4ParticleHPThermalBoost.hh" 44 #include "G4Neutron.hh" 45 46 G4ParticleHPChannelList::G4ParticleHPChannelList(G4int n, G4ParticleDefinition* p) 47 { 48 nChannels = n; 49 theChannels = new G4ParticleHPChannel*[n]; 50 theProjectile = (nullptr == p) ? G4Neutron::Neutron() : p; 51 } 52 53 G4ParticleHPChannelList::~G4ParticleHPChannelList() 54 { 55 if (theChannels != nullptr) { 56 for (G4int i = 0; i < nChannels; ++i) { 57 delete theChannels[i]; 58 } 59 delete[] theChannels; 60 } 61 } 62 63 G4HadFinalState* G4ParticleHPChannelList::ApplyYourself(const G4Element*, 64 const G4HadProjectile& aTrack) 65 { 66 G4ParticleHPThermalBoost aThermalE; 67 G4int i, ii; 68 69 // decide on the isotope 70 G4int numberOfIsos(0); 71 for (ii = 0; ii < nChannels; ++ii) { 72 numberOfIsos = theChannels[ii]->GetNiso(); 73 if (numberOfIsos != 0) break; 74 } 75 auto running = new G4double[numberOfIsos]; 76 running[0] = 0; 77 for (i = 0; i < numberOfIsos; i++) { 78 if (i != 0) running[i] = running[i - 1]; 79 for (ii = 0; ii < nChannels; ii++) { 80 if (theChannels[ii]->HasAnyData(i)) { 81 running[i] += theChannels[ii]->GetWeightedXsec( 82 aThermalE.GetThermalEnergy(aTrack, theChannels[ii]->GetN(i), theChannels[ii]->GetZ(i), 83 aTrack.GetMaterial()->GetTemperature()), i); 84 } 85 } 86 } 87 G4int isotope = nChannels - 1; 88 G4double random = G4UniformRand(); 89 for (i = 0; i < numberOfIsos; ++i) { 90 isotope = i; 91 if (running[numberOfIsos - 1] != 0) 92 if (random < running[i] / running[numberOfIsos - 1]) break; 93 } 94 delete[] running; 95 96 // decide on the channel 97 running = new G4double[nChannels]; 98 running[0] = 0; 99 G4int targA = -1; // For production of unChanged 100 G4int targZ = -1; 101 for (i = 0; i < nChannels; i++) { 102 if (i != 0) running[i] = running[i - 1]; 103 if (theChannels[i]->HasAnyData(isotope)) { 104 targA = (G4int)theChannels[i]->GetN(isotope); // Will be simply used the last valid value 105 targZ = (G4int)theChannels[i]->GetZ(isotope); 106 running[i] += theChannels[i]->GetFSCrossSection( 107 aThermalE.GetThermalEnergy(aTrack, targA, targZ, aTrack.GetMaterial()->GetTemperature()), 108 isotope); 109 targA = (G4int)theChannels[i]->GetN(isotope); // Will be simply used the last valid value 110 targZ = (G4int)theChannels[i]->GetZ(isotope); 111 // G4cout << " G4ParticleHPChannelList " << i << " targA " << targA << " targZ " << targZ << 112 //" running " << running[i] << G4endl; 113 } 114 } 115 116 // TK120607 117 if (running[nChannels - 1] == 0) { 118 // This happened usually by the miss match between the cross section data and model 119 if (targA == -1 && targZ == -1) { 120 throw G4HadronicException( 121 __FILE__, __LINE__, 122 "ParticleHP model encounter lethal discrepancy with cross section data"); 123 } 124 125 // TK121106 126 G4cout << "Warning from NeutronHP: could not find proper reaction channel. This may cause by " 127 "inconsistency between cross section and model. Unchanged final states are returned." 128 << G4endl; 129 unChanged.Clear(); 130 131 // For Ep Check create unchanged final state including rest target 132 G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon(targZ, targA, 0.0); 133 auto targ_dp = new G4DynamicParticle(targ_pd, G4ThreeVector(1, 0, 0), 0.0); 134 unChanged.SetEnergyChange(aTrack.GetKineticEnergy()); 135 unChanged.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 136 unChanged.AddSecondary(targ_dp); 137 // TK121106 138 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargA(targA); 139 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargZ(targZ); 140 delete[] running; 141 return &unChanged; 142 } 143 // TK120607 144 145 G4int lChan = 0; 146 random = G4UniformRand(); 147 for (i = 0; i < nChannels; i++) { 148 lChan = i; 149 if (running[nChannels - 1] != 0) 150 if (random < running[i] / running[nChannels - 1]) break; 151 } 152 delete[] running; 153 #ifdef G4PHPDEBUG 154 if (G4ParticleHPManager::GetInstance()->GetDEBUG()) 155 G4cout << " G4ParticleHPChannelList SELECTED ISOTOPE " << isotope << " SELECTED CHANNEL " 156 << lChan << G4endl; 157 #endif 158 return theChannels[lChan]->ApplyYourself(aTrack, isotope); 159 } 160 161 G4HadFinalState* G4ParticleHPChannelList::ApplyYourself( G4int isotope, G4int anZ, G4int anA, const G4HadProjectile & aTrack ) { 162 // do not choose Z and A again and go to selection of channel 163 G4ParticleHPThermalBoost aThermalE; 164 G4int i; 165 G4double random; 166 // decide on the channel 167 G4double* running = new G4double[nChannels]; 168 running[0] = 0.0; 169 // targA and targZ does not set to -1 170 G4int targA = anA; 171 G4int targZ = anZ; 172 G4double energy = aThermalE.GetThermalEnergy( aTrack, targA, targZ, aTrack.GetMaterial()->GetTemperature() ); 173 for ( i = 0; i < nChannels; i++ ) { 174 if ( i != 0 ) running[i] = running[i-1]; 175 if ( theChannels[i]->HasAnyData( isotope ) ) { 176 targA = (G4int) theChannels[i]->GetN( isotope ); // Will be simply used the last valid value 177 targZ = (G4int) theChannels[i]->GetZ( isotope ); 178 running[i] += theChannels[i]->GetFSCrossSection( energy, isotope ); 179 } 180 } 181 if ( running[nChannels-1] == 0 ) { 182 // This happened usually by the miss match between the cross section data and model 183 if ( targA == -1 && targZ == -1 ) { 184 throw G4HadronicException( __FILE__, __LINE__, "ParticleHP model encounter lethal discrepancy with cross section data" ); 185 } 186 G4cout << "Warning from NeutronHP: could not find proper reaction channel. " 187 << "This may be caused by inconsistency between cross section and model. " 188 << "Unchanged final states are returned." << G4endl; 189 unChanged.Clear(); 190 // For Ep Check create unchanged final state including rest target 191 G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon( targZ, targA, 0.0 ); 192 G4DynamicParticle* targ_dp = new G4DynamicParticle( targ_pd, G4ThreeVector(1.0, 0.0, 0.0), 0.0 ); 193 unChanged.SetEnergyChange( aTrack.GetKineticEnergy() ); 194 unChanged.SetMomentumChange( aTrack.Get4Momentum().vect().unit() ); 195 unChanged.AddSecondary( targ_dp ); 196 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargA( targA ); 197 G4ParticleHPManager::GetInstance()->GetReactionWhiteBoard()->SetTargZ( targZ ); 198 delete [] running; 199 return &unChanged; 200 } 201 G4int lChan = 0; 202 random = G4UniformRand(); 203 for ( i = 0; i < nChannels; i++ ) { 204 lChan = i; 205 if ( running[nChannels-1] != 0 ) if ( random < running[i]/running[nChannels-1] ) break; 206 } 207 delete [] running; 208 #ifdef G4PHPDEBUG 209 if ( G4FindDataDir( "G4ParticleHPDebug" ) != nullptr ) G4cout << " G4ParticleHPChannelList SELECTED ISOTOPE " << isotope 210 << " SELECTED CHANNEL " << lChan << G4endl; 211 #endif 212 return theChannels[lChan]->ApplyYourself( aTrack, isotope ); 213 } 214 215 void G4ParticleHPChannelList::Init(G4Element* anElement, const G4String& dirName, 216 G4ParticleDefinition* projectile) 217 { 218 theDir = dirName; 219 // G4cout << theDir << G4endl; 220 theElement = anElement; 221 // G4cout << theElement->GetName() << G4endl; 222 theProjectile = projectile; 223 } 224 225 void G4ParticleHPChannelList::Register(G4ParticleHPFinalState* theFS, const G4String& aName) 226 { 227 theChannels[idx] = new G4ParticleHPChannel(theProjectile); 228 theChannels[idx]->Init(theElement, theDir, aName); 229 theChannels[idx]->Register(theFS); 230 ++idx; 231 } 232 233 void G4ParticleHPChannelList::DumpInfo() 234 { 235 G4cout << "================================================================" << G4endl; 236 G4cout << " Element: " << theElement->GetName() << G4endl; 237 G4cout << " Number of channels: " << nChannels << G4endl; 238 G4cout << " Projectile: " << theProjectile->GetParticleName() << G4endl; 239 G4cout << " Directory name: " << theDir << G4endl; 240 for (G4int i = 0; i < nChannels; ++i) { 241 if (theChannels[i]->HasDataInAnyFinalState()) { 242 G4cout << "----------------------------------------------------------------" << G4endl; 243 theChannels[i]->DumpInfo(); 244 G4cout << "----------------------------------------------------------------" << G4endl; 245 } 246 } 247 G4cout << "================================================================" << G4endl; 248 } 249