Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 29 // 30 // 070523 bug fix for G4FPE_DEBUG on by A. How 31 // 071031 bug fix T. Koi on behalf of A. Howar 32 // 081203 bug fix in Register method by T. Koi 33 // 34 // P. Arce, June-2014 Conversion neutron_hp to 35 // 36 // June-2019 - E. Mendoza --> Modification to 37 // data library if the G4NEUTRONHP_SKIP_MISS 38 // flag is defined. The missing XS are set t 39 40 #include "G4ParticleHPChannel.hh" 41 42 #include "G4HadTmpUtil.hh" 43 #include "G4ParticleHPElasticFS.hh" 44 #include "G4ParticleHPFinalState.hh" 45 #include "G4ParticleHPReactionWhiteBoard.hh" 46 #include "G4ParticleHPThermalBoost.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "globals.hh" 49 50 #include <cstdlib> 51 52 G4ParticleHPChannel::G4ParticleHPChannel(G4Par 53 { 54 fManager = G4ParticleHPManager::GetInstance( 55 if (fManager->GetUseWendtFissionModel()) { 56 wendtFissionGenerator = G4WendtFissionFrag 57 // Make sure both fission fragment models 58 fManager->SetProduceFissionFragments(false 59 } 60 theProjectile = (nullptr == p) ? G4Neutron:: 61 theChannelData = new G4ParticleHPVector; 62 } 63 64 G4ParticleHPChannel::~G4ParticleHPChannel() 65 { 66 delete theChannelData; 67 // Following statement disabled to avoid SEG 68 // theBuffer is also deleted as "theChannelD 69 delete[] theIsotopeWiseData; 70 if (theFinalStates != nullptr) { 71 for (G4int i = 0; i < niso; i++) { 72 delete theFinalStates[i]; 73 } 74 delete[] theFinalStates; 75 } 76 delete[] active; 77 } 78 79 G4double G4ParticleHPChannel::GetXsec(G4double 80 { 81 return std::max(0., theChannelData->GetXsec( 82 } 83 84 G4double G4ParticleHPChannel::GetWeightedXsec( 85 G4int isoNumber) const 86 { 87 return theIsotopeWiseData[isoNumber].GetXsec 88 } 89 90 G4double G4ParticleHPChannel::GetFSCrossSectio 91 G4int isoNumber) const 92 { 93 return theFinalStates[isoNumber]->GetXsec(en 94 } 95 96 void G4ParticleHPChannel::Init(G4Element* anEl 97 const G4String& dirName, const G4 98 { 99 theFSType = aFSType; 100 Init(anElement, dirName); 101 } 102 103 void G4ParticleHPChannel::Init(G4Element* anEl 104 { 105 theDir = dirName; 106 theElement = anElement; 107 } 108 109 G4bool G4ParticleHPChannel::Register(G4Particl 110 { 111 ++registerCount; 112 G4int Z = theElement->GetZasInt(); 113 114 niso = (G4int)theElement->GetNumberOfIsotope 115 const std::size_t nsize = niso > 0 ? niso : 116 117 delete[] theIsotopeWiseData; 118 theIsotopeWiseData = new G4ParticleHPIsoData 119 delete[] active; 120 active = new G4bool[nsize]; 121 122 delete[] theFinalStates; 123 theFinalStates = new G4ParticleHPFinalState* 124 delete theChannelData; 125 theChannelData = new G4ParticleHPVector; 126 for (G4int i = 0; i < niso; ++i) { 127 theFinalStates[i] = theFS->New(); 128 theFinalStates[i]->SetProjectile(theProjec 129 } 130 if (niso != 0 && registerCount == 0) { 131 for (G4int i1 = 0; i1 < niso; ++i1) { 132 G4int A = theElement->GetIsotope(i1)->Ge 133 G4int M = theElement->GetIsotope(i1)->Ge 134 //G4cout <<" Init: normal case i=" << i1 135 // << " Z=" << Z << " A=" << A << G4 136 G4double frac = theElement->GetRelativeA 137 theFinalStates[i1]->SetA_Z(A, Z, M); 138 UpdateData(A, Z, M, i1, frac, theProject 139 } 140 } 141 G4bool result = HasDataInAnyFinalState(); 142 143 // To avoid issuing hash by worker threads 144 if (result) theChannelData->Hash(); 145 146 return result; 147 } 148 149 void G4ParticleHPChannel::UpdateData(G4int A, 150 G4double 151 G4Particl 152 { 153 // Initialze the G4FissionFragment generator 154 if (wendtFissionGenerator != nullptr) { 155 wendtFissionGenerator->InitializeANucleus( 156 } 157 158 theFinalStates[index]->Init(A, Z, M, theDir, 159 if (!theFinalStates[index]->HasAnyData()) re 160 // nothing there for exactly this isotope. 161 162 // the above has put the X-sec into the FS 163 theBuffer = nullptr; 164 if (theFinalStates[index]->HasXsec()) { 165 theBuffer = theFinalStates[index]->GetXsec 166 theBuffer->Times(abundance / 100.); 167 theIsotopeWiseData[index].FillChannelData( 168 } 169 else // get data from CrossSection director 170 { 171 const G4String& tString = "/CrossSection"; 172 active[index] = theIsotopeWiseData[index]. 173 174 if (active[index]) theBuffer = theIsotopeW 175 } 176 if (theBuffer != nullptr) Harmonise(theChann 177 } 178 179 void G4ParticleHPChannel::Harmonise(G4Particle 180 G4Particle 181 { 182 G4int s_tmp = 0, n = 0, m_tmp = 0; 183 auto theMerge = new G4ParticleHPVector; 184 G4ParticleHPVector* anActive = theStore; 185 G4ParticleHPVector* aPassive = theNew; 186 G4ParticleHPVector* tmp; 187 G4int a = s_tmp, p = n, t; 188 while (a < anActive->GetVectorLength() && p 189 // Loop checking, 11.05.2015, T. Koi 190 { 191 if (anActive->GetEnergy(a) <= aPassive->Ge 192 G4double xa = anActive->GetEnergy(a); 193 theMerge->SetData(m_tmp, xa, anActive->G 194 m_tmp++; 195 a++; 196 G4double xp = aPassive->GetEnergy(p); 197 if (std::abs(std::abs(xp - xa) / xa) < 0 198 ++p; 199 } 200 } 201 else { 202 tmp = anActive; 203 t = a; 204 anActive = aPassive; 205 a = p; 206 aPassive = tmp; 207 p = t; 208 } 209 } 210 while (a != anActive->GetVectorLength()) // 211 { 212 theMerge->SetData(m_tmp++, anActive->GetEn 213 ++a; 214 } 215 while (p != aPassive->GetVectorLength()) // 216 { 217 if (std::abs(theMerge->GetEnergy(std::max( 218 aPassive->GetEnergy(p)) / aPassive->GetEn 219 theMerge->SetData(m_tmp++, aPassive->Get 220 ++p; 221 } 222 delete theStore; 223 theStore = theMerge; 224 } 225 226 G4WendtFissionFragmentGenerator* G4ParticleHPC 227 if ( wendtFissionGenerator ) return wendtFis 228 else return nullptr; 229 } 230 231 G4HadFinalState* 232 G4ParticleHPChannel::ApplyYourself(const G4Had 233 G4int anIsotope, G4bool isElastic) 234 { 235 //G4cout << "G4ParticleHPChannel::ApplyYours 236 // << " ni=" << anIsotope << " isElastic=" 237 if (anIsotope != -1 && anIsotope != -2) { 238 // Inelastic Case 239 //G4cout << "G4ParticleHPChannel Inelastic 240 //<< " Z= " << GetZ(anIsotope) << " A = " 241 fManager->GetReactionWhiteBoard()->SetTarg 242 fManager->GetReactionWhiteBoard()->SetTarg 243 return theFinalStates[anIsotope]->ApplyYou 244 } 245 G4double sum = 0; 246 G4int it = 0; 247 auto xsec = new G4double[niso]; 248 G4ParticleHPThermalBoost aThermalE; 249 for (G4int i = 0; i < niso; i++) { 250 if (theFinalStates[i]->HasAnyData()) { 251 /* 252 G4cout << "FS: " << i << theTrack.GetDef 253 << " Z=" << theFinalStates[i]->GetZ() 254 << " A=" << theFinalStates[i]->GetN() 255 << G4endl; 256 */ 257 xsec[i] = theIsotopeWiseData[i].GetXsec( 258 aThermalE.GetThermalEnergy(theTrack, t 259 theFinalSta 260 theTrack.Ge 261 sum += xsec[i]; 262 } 263 else { 264 xsec[i] = 0; 265 } 266 } 267 if (sum == 0) { 268 it = G4lrint(niso * G4UniformRand()); 269 } 270 else { 271 G4double random = G4UniformRand(); 272 G4double running = 0; 273 for (G4int ix = 0; ix < niso; ix++) { 274 running += xsec[ix]; 275 if (sum == 0 || random <= running / sum) 276 it = ix; 277 break; 278 } 279 } 280 if (it == niso) it--; 281 } 282 delete[] xsec; 283 G4HadFinalState* theFinalState = nullptr; 284 const auto A = (G4int)this->GetN(it); 285 const auto Z = (G4int)this->GetZ(it); 286 const auto M = (G4int)this->GetM(it); 287 288 //-2:Marker for Fission 289 if ((wendtFissionGenerator != nullptr) && an 290 theFinalState = wendtFissionGenerator->App 291 } 292 293 // Use the standard procedure if the G4Fissi 294 if (theFinalState == nullptr) { 295 G4int icounter = 0; 296 G4int icounter_max = 1024; 297 while (theFinalState == nullptr) // Loop 298 { 299 icounter++; 300 if (icounter > icounter_max) { 301 G4cout << "Loop-counter exceeded the t 302 << __LINE__ << "th line of " << 303 break; 304 } 305 if (isElastic) { 306 // Register 0 K cross-section for DBRC 307 G4ParticleHPVector* xsforFS = theIsoto 308 // Only G4ParticleHPElasticFS has the 309 static_cast<G4ParticleHPElasticFS*>(th 310 } 311 theFinalState = theFinalStates[it]->Appl 312 } 313 } 314 315 // G4cout <<"THE IMPORTANT RETURN"<<G4endl; 316 // G4cout << "TK G4ParticleHPChannel Elastic 317 //<< " Z= " << this->GetZ(it) << " A = " << 318 fManager->GetReactionWhiteBoard()->SetTargA( 319 fManager->GetReactionWhiteBoard()->SetTargZ( 320 fManager->GetReactionWhiteBoard()->SetTargM( 321 322 return theFinalState; 323 } 324 325 void G4ParticleHPChannel::DumpInfo() const 326 { 327 G4cout << " Element: " << theElement->GetNam 328 G4cout << " Directory name: " << theDir << G 329 G4cout << " FS name: " << theFSType << G4end 330 G4cout << " Number of Isotopes: " << niso << 331 G4cout << " Have cross sections: " << G4endl 332 for (int i = 0; i < niso; i++) { 333 G4cout << theFinalStates[i]->HasXsec() << 334 } 335 G4cout << G4endl; 336 if (theChannelData != nullptr) { 337 G4cout << " Cross Section (total for this 338 int np = theChannelData->GetVectorLength() 339 G4cout << np << G4endl; 340 for (int i = 0; i < np; i++) { 341 G4cout << theChannelData->GetEnergy(i) / 342 } 343 } 344 } 345