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
00036 #include "G4NeutronHPCaptureFS.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4Gamma.hh"
00040 #include "G4ReactionProduct.hh"
00041 #include "G4Nucleus.hh"
00042 #include "G4PhotonEvaporation.hh"
00043 #include "G4Fragment.hh"
00044 #include "G4ParticleTable.hh"
00045 #include "G4NeutronHPDataUsed.hh"
00046
00047 G4HadFinalState * G4NeutronHPCaptureFS::ApplyYourself(const G4HadProjectile & theTrack)
00048 {
00049
00050 G4int i;
00051 theResult.Clear();
00052
00053 G4double eKinetic = theTrack.GetKineticEnergy();
00054 const G4HadProjectile *incidentParticle = &theTrack;
00055 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00056 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00057 theNeutron.SetKineticEnergy( eKinetic );
00058
00059
00060 G4ReactionProduct theTarget;
00061 G4Nucleus aNucleus;
00062 G4double eps = 0.0001;
00063 if(targetMass<500*MeV)
00064 targetMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theBaseA+eps) , static_cast<G4int>(theBaseZ+eps) )) /
00065 G4Neutron::Neutron()->GetPDGMass();
00066 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
00067 G4double temperature = theTrack.GetMaterial()->GetTemperature();
00068 theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
00069
00070
00071 theNeutron.Lorentz(theNeutron, -1*theTarget);
00072 eKinetic = theNeutron.GetKineticEnergy();
00073
00074
00075
00076 G4ReactionProductVector * thePhotons = 0;
00077 if ( HasFSData() && !getenv ( "G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION" ) )
00078 {
00079
00080 if ( hasExactMF6 )
00081 {
00082 theMF6FinalState.SetTarget(theTarget);
00083 theMF6FinalState.SetNeutron(theNeutron);
00084 thePhotons = theMF6FinalState.Sample( eKinetic );
00085 }
00086 else
00087 thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
00088 }
00089 else
00090 {
00091 G4ThreeVector aCMSMomentum = theNeutron.GetMomentum()+theTarget.GetMomentum();
00092 G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
00093 G4Fragment nucleus(static_cast<G4int>(theBaseA+1), static_cast<G4int>(theBaseZ) ,p4);
00094 G4PhotonEvaporation photonEvaporation;
00095
00096 photonEvaporation.SetICM( TRUE );
00097 G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
00098 G4FragmentVector::iterator it;
00099 thePhotons = new G4ReactionProductVector;
00100 for(it=products->begin(); it!=products->end(); it++)
00101 {
00102 G4ReactionProduct * theOne = new G4ReactionProduct;
00103
00104 if ( (*it)->GetParticleDefinition() != 0 )
00105 theOne->SetDefinition( (*it)->GetParticleDefinition() );
00106 else
00107 theOne->SetDefinition( G4Gamma::Gamma() );
00108
00109
00110
00111 G4ParticleTable* theTable = G4ParticleTable::GetParticleTable();
00112 if( (*it)->GetMomentum().mag() > 10*MeV ) theOne->SetDefinition( theTable->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)) );
00113
00114
00115 if ( (*it)->GetExcitationEnergy() > 1.0e-2*eV )
00116 {
00117 G4double ex = (*it)->GetExcitationEnergy();
00118 G4ReactionProduct* aPhoton = new G4ReactionProduct;
00119 aPhoton->SetDefinition( G4Gamma::Gamma() );
00120 aPhoton->SetMomentum( (*it)->GetMomentum().vect().unit() * ex );
00121
00122 thePhotons->push_back(aPhoton);
00123 }
00124
00125 theOne->SetMomentum( (*it)->GetMomentum().vect() * ( (*it)->GetMomentum().t() - (*it)->GetExcitationEnergy() ) / (*it)->GetMomentum().t() ) ;
00126
00127 thePhotons->push_back(theOne);
00128 delete *it;
00129 }
00130 delete products;
00131 }
00132
00133
00134
00135
00136
00137 G4int nPhotons = 0;
00138 if(thePhotons!=0) nPhotons=thePhotons->size();
00139
00140
00141
00142
00143
00144
00145
00146
00147 if ( nPhotons == 0 )
00148 {
00149 G4ReactionProduct * theOne = new G4ReactionProduct;
00150 theOne->SetDefinition( G4Gamma::Gamma() );
00151 G4double theta = pi*G4UniformRand();
00152 G4double phi = twopi*G4UniformRand();
00153 G4double sinth = std::sin(theta);
00154 G4ThreeVector direction( sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) );
00155 theOne->SetMomentum( direction ) ;
00156 thePhotons->push_back(theOne);
00157 nPhotons++;
00158 }
00159
00160
00161
00162 if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 )
00163 {
00164 G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();
00165 G4double Q = G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ))->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass()
00166 - G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ))->GetPDGMass();
00167 thePhotons->operator[](0)->SetMomentum( Q*direction );
00168 }
00169
00170
00171
00172 for(i=0; i<nPhotons; i++)
00173 {
00174 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), theTarget);
00175 }
00176
00177
00178
00179 if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 )
00180 {
00181 G4DynamicParticle * theOne = new G4DynamicParticle;
00182 G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
00183 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
00184 theOne->SetDefinition(aRecoil);
00185
00186
00187 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect()
00188 +theTarget.GetMomentum()
00189 -thePhotons->operator[](0)->GetMomentum();
00190
00191 G4ThreeVector theMomUnit = aMomentum.unit();
00192 G4double aKinEnergy = theTrack.GetKineticEnergy()
00193 +theTarget.GetKineticEnergy();
00194 G4double theResMass = aRecoil->GetPDGMass();
00195 G4double theResE = aRecoil->GetPDGMass()+aKinEnergy;
00196 G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass);
00197 G4ThreeVector theMomentum = theAbsMom*theMomUnit;
00198 theOne->SetMomentum(theMomentum);
00199 theResult.AddSecondary(theOne);
00200 }
00201
00202
00203 for(i=0; i<nPhotons; i++)
00204 {
00205
00206 G4DynamicParticle * theOne = new G4DynamicParticle;
00207 theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition());
00208 theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
00209 theResult.AddSecondary(theOne);
00210 delete thePhotons->operator[](i);
00211 }
00212 delete thePhotons;
00213
00214
00215 G4bool residual = false;
00216 G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
00217 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
00218 for ( G4int j = 0 ; j != theResult.GetNumberOfSecondaries() ; j++ )
00219 {
00220 if ( theResult.GetSecondary(j)->GetParticle()->GetDefinition() == aRecoil ) residual = true;
00221 }
00222
00223 if ( residual == false )
00224 {
00225
00226
00227 G4int nNonZero = 0;
00228 G4LorentzVector p_photons(0,0,0,0);
00229 for ( G4int j = 0 ; j != theResult.GetNumberOfSecondaries() ; j++ )
00230 {
00231 p_photons += theResult.GetSecondary(j)->GetParticle()->Get4Momentum();
00232
00233 if ( theResult.GetSecondary(j)->GetParticle()->Get4Momentum().e() > 0 ) nNonZero++;
00234 }
00235
00236
00237 G4double deltaE = ( theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy() )
00238 - ( p_photons.e() + aRecoil->GetPDGMass() );
00239
00240
00241 if ( nPhotons - nNonZero > 0 )
00242 {
00243
00244 std::vector<G4double> vRand;
00245 vRand.push_back( 0.0 );
00246 for ( G4int j = 0 ; j != nPhotons - nNonZero - 1 ; j++ )
00247 {
00248 vRand.push_back( G4UniformRand() );
00249 }
00250 vRand.push_back( 1.0 );
00251 std::sort( vRand.begin(), vRand.end() );
00252
00253 std::vector<G4double> vEPhoton;
00254 for ( G4int j = 0 ; j < (G4int)vRand.size() - 1 ; j++ )
00255 {
00256 vEPhoton.push_back( deltaE * ( vRand[j+1] - vRand[j] ) );
00257 }
00258 std::sort( vEPhoton.begin(), vEPhoton.end() );
00259
00260 for ( G4int j = 0 ; j < nPhotons - nNonZero - 1 ; j++ )
00261 {
00262
00263 G4double theta = pi*G4UniformRand();
00264 G4double phi = twopi*G4UniformRand();
00265 G4double sinth = std::sin(theta);
00266 G4double en = vEPhoton[j];
00267 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00268
00269 p_photons += G4LorentzVector ( tempVector, tempVector.mag() );
00270 G4DynamicParticle * theOne = new G4DynamicParticle;
00271 theOne->SetDefinition( G4Gamma::Gamma() );
00272 theOne->SetMomentum( tempVector );
00273 theResult.AddSecondary(theOne);
00274 }
00275
00276
00277 G4DynamicParticle * theOne = new G4DynamicParticle;
00278 theOne->SetDefinition( G4Gamma::Gamma() );
00279
00280 G4ThreeVector lastPhoton = -p_photons.vect().unit()*vEPhoton.back();
00281 p_photons += G4LorentzVector( lastPhoton , lastPhoton.mag() );
00282 theOne->SetMomentum( lastPhoton );
00283 theResult.AddSecondary(theOne);
00284 }
00285
00286
00287 G4DynamicParticle * theOne = new G4DynamicParticle;
00288 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum()
00289 - p_photons.vect();
00290 theOne->SetDefinition(aRecoil);
00291 theOne->SetMomentum( aMomentum );
00292 theResult.AddSecondary(theOne);
00293
00294 }
00295
00296
00297
00298 theResult.SetStatusChange(stopAndKill);
00299 return &theResult;
00300 }
00301
00302 #include <sstream>
00303 void G4NeutronHPCaptureFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & )
00304 {
00305
00306
00307 std::stringstream ss;
00308 ss << static_cast<G4int>(Z);
00309 G4String sZ;
00310 ss >> sZ;
00311 ss.clear();
00312 ss << static_cast<G4int>(A);
00313 G4String sA;
00314 ss >> sA;
00315
00316 ss.clear();
00317 G4String sM;
00318 if ( M > 0 )
00319 {
00320 ss << "m";
00321 ss << M;
00322 ss >> sM;
00323 ss.clear();
00324 }
00325
00326 G4String element_name = theNames.GetName( static_cast<G4int>(Z)-1 );
00327 G4String filenameMF6 = dirName+"/FSMF6/"+sZ+"_"+sA+sM+"_"+element_name;
00328 std::ifstream dummyIFS(filenameMF6, std::ios::in);
00329 if ( dummyIFS.good() == true ) hasExactMF6=true;
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347 if ( hasExactMF6 == true )
00348 {
00349 std::ifstream theData(filenameMF6, std::ios::in);
00350 theMF6FinalState.Init(theData);
00351 theData.close();
00352 return;
00353 }
00354
00355
00356
00357 G4String tString = "/FS";
00358 G4bool dbool;
00359 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
00360
00361 G4String filename = aFile.GetName();
00362 SetAZMs( A, Z, M, aFile );
00363
00364
00365 if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
00366 {
00367 hasAnyData = false;
00368 hasFSData = false;
00369 hasXsec = false;
00370 return;
00371 }
00372 std::ifstream theData(filename, std::ios::in);
00373
00374 hasFSData = theFinalStatePhotons.InitMean(theData);
00375 if(hasFSData)
00376 {
00377 targetMass = theFinalStatePhotons.GetTargetMass();
00378 theFinalStatePhotons.InitAngular(theData);
00379 theFinalStatePhotons.InitEnergies(theData);
00380 }
00381 theData.close();
00382 }