00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include "G4NeutronHPInelasticBaseFS.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "G4Nucleus.hh"
00038 #include "G4NucleiProperties.hh"
00039 #include "G4He3.hh"
00040 #include "G4Alpha.hh"
00041 #include "G4Electron.hh"
00042 #include "G4NeutronHPDataUsed.hh"
00043
00044 #include "G4ParticleTable.hh"
00045
00046 void G4NeutronHPInelasticBaseFS::InitGammas(G4double AR, G4double ZR)
00047 {
00048
00049
00050
00051
00052
00053
00054 std::ostringstream ost;
00055 ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
00056 G4String aName = ost.str();
00057 std::ifstream from(aName, std::ios::in);
00058
00059 if(!from) return;
00060
00061 std::ifstream theGammaData(aName, std::ios::in);
00062
00063 G4double eps = 0.001;
00064 theNuclearMassDifference =
00065 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(AR+eps),static_cast<G4int>(ZR+eps)) -
00066 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps));
00067 theGammas.Init(theGammaData);
00068
00069 }
00070
00071 void G4NeutronHPInelasticBaseFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & bit)
00072 {
00073 gammaPath = "/Inelastic/Gammas/";
00074 if(!getenv("G4NEUTRONHPDATA"))
00075 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
00076 G4String tBase = getenv("G4NEUTRONHPDATA");
00077 gammaPath = tBase+gammaPath;
00078 G4String tString = dirName;
00079 G4bool dbool;
00080 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M,tString, bit, dbool);
00081 G4String filename = aFile.GetName();
00082 SetAZMs( A, Z, M, aFile);
00083
00084
00085
00086
00087
00088 if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
00089 {
00090 if(getenv("NeutronHPNamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
00091 hasAnyData = false;
00092 hasFSData = false;
00093 hasXsec = false;
00094 return;
00095 }
00096
00097
00098 std::ifstream theData(filename, std::ios::in);
00099 if(!(theData))
00100 {
00101 hasAnyData = false;
00102 hasFSData = false;
00103 hasXsec = false;
00104 theData.close();
00105 return;
00106 }
00107
00108 G4int infoType, dataType, dummy=INT_MAX;
00109 hasFSData = false;
00110 while (theData >> infoType)
00111 {
00112 theData >> dataType;
00113 if(dummy==INT_MAX) theData >> dummy >> dummy;
00114 if(dataType==3)
00115 {
00116 G4int total;
00117 theData >> total;
00118 theXsection->Init(theData, total, eV);
00119 }
00120 else if(dataType==4)
00121 {
00122 theAngularDistribution = new G4NeutronHPAngular;
00123 theAngularDistribution->Init(theData);
00124 hasFSData = true;
00125 }
00126 else if(dataType==5)
00127 {
00128 theEnergyDistribution = new G4NeutronHPEnergyDistribution;
00129 theEnergyDistribution->Init(theData);
00130 hasFSData = true;
00131 }
00132 else if(dataType==6)
00133 {
00134 theEnergyAngData = new G4NeutronHPEnAngCorrelation;
00135 theEnergyAngData->Init(theData);
00136 hasFSData = true;
00137 }
00138 else if(dataType==12)
00139 {
00140 theFinalStatePhotons = new G4NeutronHPPhotonDist;
00141 theFinalStatePhotons->InitMean(theData);
00142 hasFSData = true;
00143 }
00144 else if(dataType==13)
00145 {
00146 theFinalStatePhotons = new G4NeutronHPPhotonDist;
00147 theFinalStatePhotons->InitPartials(theData);
00148 hasFSData = true;
00149 }
00150 else if(dataType==14)
00151 {
00152 theFinalStatePhotons->InitAngular(theData);
00153 hasFSData = true;
00154 }
00155 else if(dataType==15)
00156 {
00157 theFinalStatePhotons->InitEnergies(theData);
00158 hasFSData = true;
00159 }
00160 else
00161 {
00162 throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4NeutronHPInelasticBaseFS");
00163 }
00164 }
00165 theData.close();
00166 }
00167
00168 void G4NeutronHPInelasticBaseFS::BaseApply(const G4HadProjectile & theTrack,
00169 G4ParticleDefinition ** theDefs,
00170 G4int nDef)
00171 {
00172
00173
00174 theResult.Clear();
00175 G4double eKinetic = theTrack.GetKineticEnergy();
00176 const G4HadProjectile *incidentParticle = &theTrack;
00177 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00178 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00179 theNeutron.SetKineticEnergy( eKinetic );
00180
00181
00182 G4double targetMass;
00183 G4double eps = 0.0001;
00184 targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
00185 G4Neutron::Neutron()->GetPDGMass();
00186
00187 if(theEnergyAngData!=0)
00188 { targetMass = theEnergyAngData->GetTargetMass(); }
00189 if(theAngularDistribution!=0)
00190 { targetMass = theAngularDistribution->GetTargetMass(); }
00191
00192
00193
00194
00195 if ( targetMass == 0 )
00196 {
00197
00198 targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) / G4Neutron::Neutron()->GetPDGMass();
00199 }
00200
00201 G4Nucleus aNucleus;
00202 G4ReactionProduct theTarget;
00203 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
00204 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
00205
00206
00207 G4ReactionProduct boosted;
00208 boosted.Lorentz(theNeutron, theTarget);
00209 eKinetic = boosted.GetKineticEnergy();
00210 G4double orgMomentum = boosted.GetMomentum().mag();
00211
00212
00213 if(!HasFSData())
00214 {
00215 G4NeutronHPNBodyPhaseSpace thePhaseSpaceDistribution;
00216 G4double aPhaseMass=0;
00217 G4int ii;
00218 for(ii=0; ii<nDef; ii++)
00219 {
00220 aPhaseMass+=theDefs[ii]->GetPDGMass();
00221 }
00222 thePhaseSpaceDistribution.Init(aPhaseMass, nDef);
00223 thePhaseSpaceDistribution.SetNeutron(&theNeutron);
00224 thePhaseSpaceDistribution.SetTarget(&theTarget);
00225 for(ii=0; ii<nDef; ii++)
00226 {
00227 G4double massCode = 1000.*std::abs(theDefs[ii]->GetPDGCharge());
00228 massCode += theDefs[ii]->GetBaryonNumber();
00229 G4double dummy = 0;
00230 G4ReactionProduct * aSec = thePhaseSpaceDistribution.Sample(eKinetic, massCode, dummy);
00231 aSec->Lorentz(*aSec, -1.*theTarget);
00232 G4DynamicParticle * aPart = new G4DynamicParticle();
00233 aPart->SetDefinition(aSec->GetDefinition());
00234 aPart->SetMomentum(aSec->GetMomentum());
00235 delete aSec;
00236 theResult.AddSecondary(aPart);
00237 }
00238 theResult.SetStatusChange(stopAndKill);
00239
00240
00241
00242 G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
00243 G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
00244 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
00245 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
00246 adjust_final_state ( init_4p_lab );
00247
00248 return;
00249 }
00250
00251
00252 if(theAngularDistribution!=0)
00253 {
00254 theAngularDistribution->SetTarget(theTarget);
00255 theAngularDistribution->SetNeutron(theNeutron);
00256 }
00257 else if(theEnergyAngData!=0)
00258 {
00259 theEnergyAngData->SetTarget(theTarget);
00260 theEnergyAngData->SetNeutron(theNeutron);
00261 }
00262
00263 G4ReactionProductVector * tmpHadrons = 0;
00264 G4int ii, dummy;
00265 unsigned int i;
00266 if(theEnergyAngData != 0)
00267 {
00268 tmpHadrons = theEnergyAngData->Sample(eKinetic);
00269 }
00270 else if(theAngularDistribution!= 0)
00271 {
00272 G4bool * Done = new G4bool[nDef];
00273 G4int i0;
00274 for(i0=0; i0<nDef; i0++) Done[i0] = false;
00275 if(tmpHadrons == 0)
00276 {
00277 tmpHadrons = new G4ReactionProductVector;
00278 }
00279 else
00280 {
00281 for(i=0; i<tmpHadrons->size(); i++)
00282 {
00283 for(ii=0; ii<nDef; ii++)
00284 if(!Done[ii] && tmpHadrons->operator[](i)->GetDefinition() == theDefs[ii])
00285 Done[ii] = true;
00286 }
00287 }
00288 G4ReactionProduct * aHadron;
00289 G4double localMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps)));
00290 G4ThreeVector bufferedDirection(0,0,0);
00291 for(i0=0; i0<nDef; i0++)
00292 {
00293 if(!Done[i0])
00294 {
00295 aHadron = new G4ReactionProduct;
00296 if(theEnergyDistribution!=0)
00297 {
00298 aHadron->SetDefinition(theDefs[i0]);
00299 aHadron->SetKineticEnergy(theEnergyDistribution->Sample(eKinetic, dummy));
00300 }
00301 else if(nDef == 1)
00302 {
00303 aHadron->SetDefinition(theDefs[i0]);
00304 aHadron->SetKineticEnergy(eKinetic);
00305 }
00306 else if(nDef == 2)
00307 {
00308 aHadron->SetDefinition(theDefs[i0]);
00309 aHadron->SetKineticEnergy(50*MeV);
00310 }
00311 else
00312 {
00313 throw G4HadronicException(__FILE__, __LINE__, "No energy distribution to sample from in InelasticBaseFS::BaseApply");
00314 }
00315 theAngularDistribution->SampleAndUpdate(*aHadron);
00316 if(theEnergyDistribution==0 && nDef == 2)
00317 {
00318 if(i0==0)
00319 {
00320 G4double mass1 = theDefs[0]->GetPDGMass();
00321 G4double mass2 = theDefs[1]->GetPDGMass();
00322 G4double massn = G4Neutron::Neutron()->GetPDGMass();
00323 G4int z1 = static_cast<G4int>(theBaseZ+eps-theDefs[0]->GetPDGCharge()-theDefs[1]->GetPDGCharge());
00324 G4int a1 = static_cast<G4int>(theBaseA+eps)-theDefs[0]->GetBaryonNumber()-theDefs[1]->GetBaryonNumber();
00325 G4double concreteMass = G4NucleiProperties::GetNuclearMass(a1, z1);
00326 G4double availableEnergy = eKinetic+massn+localMass-mass1-mass2-concreteMass;
00327
00328 G4double emin = availableEnergy+mass1+mass2 - std::sqrt((mass1+mass2)*(mass1+mass2)+orgMomentum*orgMomentum);
00329 G4double p1=std::sqrt(2.*mass2*emin);
00330 bufferedDirection = p1*aHadron->GetMomentum().unit();
00331 if(getenv("HTOKEN"))
00332 {
00333 G4cout << "HTOKEN "<<z1<<" "<<theBaseZ<<" "<<a1<<" "<<theBaseA<<" "<<availableEnergy<<" "
00334 << emin<<G4endl;
00335 }
00336 }
00337 else
00338 {
00339 bufferedDirection = -bufferedDirection;
00340 }
00341
00342 if(getenv("HTOKEN"))
00343 {
00344 G4cout << " HTOKEN "<<bufferedDirection.mag2()<<G4endl;
00345 }
00346 aHadron->SetTotalEnergy( std::sqrt(aHadron->GetMass()*aHadron->GetMass()
00347 +bufferedDirection.mag2()) );
00348 aHadron->SetMomentum(bufferedDirection);
00349 aHadron->Lorentz(*aHadron, -1.*(theTarget+theNeutron));
00350 if(getenv("HTOKEN"))
00351 {
00352 G4cout << " HTOKEN "<<aHadron->GetTotalEnergy()<<" "<<aHadron->GetMomentum()<<G4endl;
00353 }
00354 }
00355 tmpHadrons->push_back(aHadron);
00356 }
00357 }
00358 delete [] Done;
00359 }
00360 else
00361 {
00362 throw G4HadronicException(__FILE__, __LINE__, "No data to create the neutrons in NInelasticFS");
00363 }
00364
00365 G4ReactionProductVector * thePhotons = 0;
00366 if(theFinalStatePhotons!=0)
00367 {
00368
00369 G4ReactionProduct boosted_tmp;
00370 boosted_tmp.Lorentz(theNeutron, theTarget);
00371 G4double anEnergy = boosted_tmp.GetKineticEnergy();
00372 thePhotons = theFinalStatePhotons->GetPhotons(anEnergy);
00373 if(thePhotons!=0)
00374 {
00375 for(i=0; i<thePhotons->size(); i++)
00376 {
00377
00378 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1.*theTarget);
00379 }
00380 }
00381 }
00382 else if(theEnergyAngData!=0)
00383 {
00384
00385 G4double theGammaEnergy = theEnergyAngData->GetTotalMeanEnergy();
00386 G4double anEnergy = boosted.GetKineticEnergy();
00387 theGammaEnergy = anEnergy-theGammaEnergy;
00388 theGammaEnergy += theNuclearMassDifference;
00389 G4double eBindProducts = 0;
00390 G4double eBindN = 0;
00391 G4double eBindP = 0;
00392 G4double eBindD = G4NucleiProperties::GetBindingEnergy(2,1);
00393 G4double eBindT = G4NucleiProperties::GetBindingEnergy(3,1);
00394 G4double eBindHe3 = G4NucleiProperties::GetBindingEnergy(3,2);
00395 G4double eBindA = G4NucleiProperties::GetBindingEnergy(4,2);
00396 G4int ia=0;
00397 for(i=0; i<tmpHadrons->size(); i++)
00398 {
00399 if(tmpHadrons->operator[](i)->GetDefinition() == G4Neutron::Neutron())
00400 {
00401 eBindProducts+=eBindN;
00402 }
00403 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Proton::Proton())
00404 {
00405 eBindProducts+=eBindP;
00406 }
00407 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Deuteron::Deuteron())
00408 {
00409 eBindProducts+=eBindD;
00410 }
00411 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Triton::Triton())
00412 {
00413 eBindProducts+=eBindT;
00414 }
00415 else if(tmpHadrons->operator[](i)->GetDefinition() == G4He3::He3())
00416 {
00417 eBindProducts+=eBindHe3;
00418 }
00419 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Alpha::Alpha())
00420 {
00421 eBindProducts+=eBindA;
00422 ia++;
00423 }
00424 }
00425
00426 theGammaEnergy += eBindProducts;
00427
00428
00429
00430 if ( (G4int)(theBaseZ+eps) == 4 && (G4int)(theBaseA+eps) == 9 )
00431 {
00432
00433 if ( std::abs( theNuclearMassDifference -
00434 ( G4NucleiProperties::GetBindingEnergy( 8 , 4 ) -
00435 G4NucleiProperties::GetBindingEnergy( 9 , 4 ) ) ) < 1*keV
00436 && ia == 2 )
00437 {
00438 theGammaEnergy -= (2*eBindA);
00439 }
00440 }
00441
00442 G4ReactionProductVector * theOtherPhotons = 0;
00443 G4int iLevel;
00444 while(theGammaEnergy>=theGammas.GetLevelEnergy(0))
00445 {
00446 for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
00447 {
00448 if(theGammas.GetLevelEnergy(iLevel)<theGammaEnergy) break;
00449 }
00450 if(iLevel==0||iLevel==theGammas.GetNumberOfLevels()-1)
00451 {
00452 theOtherPhotons = theGammas.GetDecayGammas(iLevel);
00453 }
00454 else
00455 {
00456 G4double random = G4UniformRand();
00457 G4double eLow = theGammas.GetLevelEnergy(iLevel);
00458 G4double eHigh = theGammas.GetLevelEnergy(iLevel+1);
00459 if(random > (eHigh-eLow)/(theGammaEnergy-eLow)) iLevel++;
00460 theOtherPhotons = theGammas.GetDecayGammas(iLevel);
00461 }
00462 if(thePhotons==0) thePhotons = new G4ReactionProductVector;
00463 if(theOtherPhotons != 0)
00464 {
00465 for(unsigned int iii=0; iii<theOtherPhotons->size(); iii++)
00466 {
00467 thePhotons->push_back(theOtherPhotons->operator[](iii));
00468 }
00469 delete theOtherPhotons;
00470 }
00471 theGammaEnergy -= theGammas.GetLevelEnergy(iLevel);
00472 if(iLevel == -1) break;
00473 }
00474 }
00475
00476
00477 unsigned int nSecondaries = tmpHadrons->size();
00478 unsigned int nPhotons = 0;
00479 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
00480 nSecondaries += nPhotons;
00481 G4DynamicParticle * theSec;
00482
00483 for(i=0; i<nSecondaries-nPhotons; i++)
00484 {
00485 theSec = new G4DynamicParticle;
00486 theSec->SetDefinition(tmpHadrons->operator[](i)->GetDefinition());
00487 theSec->SetMomentum(tmpHadrons->operator[](i)->GetMomentum());
00488 theResult.AddSecondary(theSec);
00489 delete tmpHadrons->operator[](i);
00490 }
00491 if(thePhotons != 0)
00492 {
00493 for(i=0; i<nPhotons; i++)
00494 {
00495 theSec = new G4DynamicParticle;
00496 theSec->SetDefinition(thePhotons->operator[](i)->GetDefinition());
00497 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
00498 theResult.AddSecondary(theSec);
00499 delete thePhotons->operator[](i);
00500 }
00501 }
00502
00503
00504 delete thePhotons;
00505 delete tmpHadrons;
00506
00507
00508 G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
00509 G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
00510 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
00511 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
00512 adjust_final_state ( init_4p_lab );
00513
00514
00515 theResult.SetStatusChange(stopAndKill);
00516 }