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 #include "G4NeutronHPElasticFS.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "G4ReactionProduct.hh"
00038 #include "G4Nucleus.hh"
00039 #include "G4Proton.hh"
00040 #include "G4Deuteron.hh"
00041 #include "G4Triton.hh"
00042 #include "G4Alpha.hh"
00043 #include "G4ThreeVector.hh"
00044 #include "G4LorentzVector.hh"
00045 #include "G4ParticleTable.hh"
00046 #include "G4NeutronHPDataUsed.hh"
00047
00048 void G4NeutronHPElasticFS::Init (G4double A, G4double Z, G4int M, G4String & dirName, G4String & )
00049 {
00050 G4String tString = "/FS";
00051 G4bool dbool;
00052 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
00053 G4String filename = aFile.GetName();
00054 SetAZMs( A, Z, M, aFile );
00055
00056
00057 if(!dbool)
00058 {
00059 hasAnyData = false;
00060 hasFSData = false;
00061 hasXsec = false;
00062 return;
00063 }
00064 std::ifstream theData(filename, std::ios::in);
00065 theData >> repFlag >> targetMass >> frameFlag;
00066 if(repFlag==1)
00067 {
00068 G4int nEnergy;
00069 theData >> nEnergy;
00070 theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
00071 theCoefficients->InitInterpolation(theData);
00072 G4double temp, energy;
00073 G4int tempdep, nLegendre;
00074 G4int i, ii;
00075 for (i=0; i<nEnergy; i++)
00076 {
00077 theData >> temp >> energy >> tempdep >> nLegendre;
00078 energy *=eV;
00079 theCoefficients->Init(i, energy, nLegendre);
00080 theCoefficients->SetTemperature(i, temp);
00081 G4double coeff=0;
00082 for(ii=0; ii<nLegendre; ii++)
00083 {
00084
00085 theData >> coeff;
00086 theCoefficients->SetCoeff(i, ii+1, coeff);
00087 }
00088 }
00089 }
00090 else if (repFlag==2)
00091 {
00092 G4int nEnergy;
00093 theData >> nEnergy;
00094 theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
00095 theProbArray->InitInterpolation(theData);
00096 G4double temp, energy;
00097 G4int tempdep, nPoints;
00098 for(G4int i=0; i<nEnergy; i++)
00099 {
00100 theData >> temp >> energy >> tempdep >> nPoints;
00101 energy *= eV;
00102 theProbArray->InitInterpolation(i, theData);
00103 theProbArray->SetT(i, temp);
00104 theProbArray->SetX(i, energy);
00105 G4double prob, costh;
00106 for(G4int ii=0; ii<nPoints; ii++)
00107 {
00108
00109 theData >> costh >> prob;
00110 theProbArray->SetX(i, ii, costh);
00111 theProbArray->SetY(i, ii, prob);
00112 }
00113 }
00114 }
00115 else if ( repFlag==3 )
00116 {
00117 G4int nEnergy_Legendre;
00118 theData >> nEnergy_Legendre;
00119 theCoefficients = new G4NeutronHPLegendreStore( nEnergy_Legendre );
00120 theCoefficients->InitInterpolation( theData );
00121 G4double temp, energy;
00122 G4int tempdep, nLegendre;
00123
00124 for ( G4int i = 0 ; i < nEnergy_Legendre ; i++ )
00125 {
00126 theData >> temp >> energy >> tempdep >> nLegendre;
00127 energy *=eV;
00128 theCoefficients->Init( i , energy , nLegendre );
00129 theCoefficients->SetTemperature( i , temp );
00130 G4double coeff = 0;
00131 for (G4int ii = 0 ; ii < nLegendre ; ii++ )
00132 {
00133
00134 theData >> coeff;
00135 theCoefficients->SetCoeff(i, ii+1, coeff);
00136 }
00137 }
00138
00139 tE_of_repFlag3 = energy;
00140
00141 G4int nEnergy_Prob;
00142 theData >> nEnergy_Prob;
00143 theProbArray = new G4NeutronHPPartial( nEnergy_Prob , nEnergy_Prob );
00144 theProbArray->InitInterpolation( theData );
00145 G4int nPoints;
00146 for ( G4int i=0 ; i < nEnergy_Prob ; i++ )
00147 {
00148 theData >> temp >> energy >> tempdep >> nPoints;
00149
00150 energy *= eV;
00151
00152
00153 if ( i == 0 )
00154
00155 if ( std::abs( energy - tE_of_repFlag3 ) / tE_of_repFlag3 > 1.0e-15 )
00156 G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl;
00157
00158 theProbArray->InitInterpolation( i , theData );
00159 theProbArray->SetT( i , temp );
00160 theProbArray->SetX( i , energy );
00161 G4double prob, costh;
00162 for( G4int ii = 0 ; ii < nPoints ; ii++ )
00163 {
00164
00165 theData >> costh >> prob;
00166 theProbArray->SetX( i , ii , costh );
00167 theProbArray->SetY( i , ii , prob );
00168 }
00169 }
00170 }
00171 else if (repFlag==0)
00172 {
00173 theData >> frameFlag;
00174 }
00175 else
00176 {
00177 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
00178 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
00179 }
00180 theData.close();
00181 }
00182 G4HadFinalState * G4NeutronHPElasticFS::ApplyYourself(const G4HadProjectile & theTrack)
00183 {
00184
00185 theResult.Clear();
00186 G4double eKinetic = theTrack.GetKineticEnergy();
00187 const G4HadProjectile *incidentParticle = &theTrack;
00188 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
00189 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
00190 theNeutron.SetKineticEnergy( eKinetic );
00191
00192
00193
00194 G4ReactionProduct theTarget;
00195 G4Nucleus aNucleus;
00196 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
00197 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207 G4ThreeVector the3Neutron = theNeutron.GetMomentum();
00208 G4double nEnergy = theNeutron.GetTotalEnergy();
00209 G4ThreeVector the3Target = theTarget.GetMomentum();
00210
00211 G4double tEnergy = theTarget.GetTotalEnergy();
00212 G4ReactionProduct theCMS;
00213 G4double totE = nEnergy+tEnergy;
00214 G4ThreeVector the3CMS = the3Target+the3Neutron;
00215 theCMS.SetMomentum(the3CMS);
00216 G4double cmsMom = std::sqrt(the3CMS*the3CMS);
00217 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
00218 theCMS.SetMass(sqrts);
00219 theCMS.SetTotalEnergy(totE);
00220
00221
00222 G4ReactionProduct boosted;
00223 boosted.Lorentz(theNeutron, theTarget);
00224 eKinetic = boosted.GetKineticEnergy();
00225 G4double cosTh = -2;
00226 if(repFlag == 1)
00227 {
00228 cosTh = theCoefficients->SampleElastic(eKinetic);
00229 }
00230
00231 else if (repFlag==2)
00232 {
00233 cosTh = theProbArray->Sample(eKinetic);
00234 }
00235 else if (repFlag==3)
00236 {
00237 if ( eKinetic <= tE_of_repFlag3 )
00238 {
00239 cosTh = theCoefficients->SampleElastic(eKinetic);
00240 }
00241 else
00242 {
00243 cosTh = theProbArray->Sample(eKinetic);
00244 }
00245 }
00246 else if (repFlag==0)
00247 {
00248 cosTh = 2.*G4UniformRand()-1.;
00249 }
00250 else
00251 {
00252 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
00253 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
00254 }
00255 if(cosTh<-1.1) { return 0; }
00256 G4double phi = twopi*G4UniformRand();
00257 G4double theta = std::acos(cosTh);
00258 G4double sinth = std::sin(theta);
00259 if (frameFlag == 1)
00260 {
00261
00262
00263
00264 theNeutron.Lorentz(theNeutron, theTarget);
00265 G4double e0 = theNeutron.GetTotalEnergy();
00266 G4double p0 = theNeutron.GetTotalMomentum();
00267 G4double mN = theNeutron.GetMass();
00268 G4double mT = theTarget.GetMass();
00269 G4double eE = e0+mT;
00270 G4double ap = (mT+eE)*(mT-eE) + (p0+mN)*(p0-mN);
00271 G4double a = 4*(eE+p0*cosTh)*(eE-p0*cosTh);
00272 G4double b = 4*ap*p0*cosTh;
00273 G4double c = (2.*eE*mN-ap)*(2.*eE*mN+ap);
00274 G4double en = (-b+std::sqrt(b*b - 4*a*c) )/(2*a);
00275 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
00276 theNeutron.SetMomentum(tempVector);
00277 theNeutron.SetTotalEnergy(std::sqrt(en*en+theNeutron.GetMass()*theNeutron.GetMass()));
00278
00279 theNeutron.Lorentz(theNeutron, -1.*theTarget);
00280
00281 theNeutron.Lorentz(theNeutron, theCMS);
00282 theTarget.SetMomentum(-theNeutron.GetMomentum());
00283 theTarget.SetTotalEnergy(theNeutron.GetTotalEnergy());
00284
00285 theNeutron.Lorentz(theNeutron, -1.*theCMS);
00286 theTarget.Lorentz(theTarget, -1.*theCMS);
00287
00288 if ( theNeutron.GetKineticEnergy() <= 0 ) theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00289 if ( theTarget.GetKineticEnergy() <= 0 ) theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00290 }
00291 else if (frameFlag == 2)
00292 {
00293 theNeutron.Lorentz(theNeutron, theCMS);
00294 theTarget.Lorentz(theTarget, theCMS);
00295 G4double en = theNeutron.GetTotalMomentum();
00296 G4ThreeVector cmsMom_tmp=theNeutron.GetMomentum();
00297 G4double cms_theta=cmsMom_tmp.theta();
00298 G4double cms_phi=cmsMom_tmp.phi();
00299 G4ThreeVector tempVector;
00300 tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
00301 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
00302 -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
00303 tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
00304 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
00305 +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
00306 tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
00307 -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
00308 tempVector *= en;
00309 theNeutron.SetMomentum(tempVector);
00310 theTarget.SetMomentum(-tempVector);
00311 G4double tP = theTarget.GetTotalMomentum();
00312 G4double tM = theTarget.GetMass();
00313 theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 theNeutron.Lorentz(theNeutron, -1.*theCMS);
00329
00330 if ( theNeutron.GetKineticEnergy() <= 0 )
00331 {
00332
00333
00334
00335 theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00336 }
00337
00338 theTarget.Lorentz(theTarget, -1.*theCMS);
00339
00340 if ( theTarget.GetKineticEnergy() < 0 )
00341 {
00342
00343
00344
00345 theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + std::pow( 10 , -15.65 ) ) );
00346 }
00347 }
00348 else
00349 {
00350 G4cout <<"Value of frameFlag (1=LAB, 2=CMS): "<<frameFlag;
00351 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::ApplyYourSelf frameflag incorrect");
00352 }
00353
00354
00355 theResult.SetEnergyChange(theNeutron.GetKineticEnergy());
00356 theResult.SetMomentumChange(theNeutron.GetMomentum().unit());
00357 G4DynamicParticle* theRecoil = new G4DynamicParticle;
00358 if(targetMass<4.5)
00359 {
00360 if(targetMass<1)
00361 {
00362
00363 theRecoil->SetDefinition(G4Proton::Proton());
00364 }
00365 else if(targetMass<2 )
00366 {
00367
00368 theRecoil->SetDefinition(G4Deuteron::Deuteron());
00369 }
00370 else if(targetMass<2.999 )
00371 {
00372
00373 theRecoil->SetDefinition(G4He3::He3());
00374 }
00375 else if(targetMass<3 )
00376 {
00377
00378 theRecoil->SetDefinition(G4Triton::Triton());
00379 }
00380 else
00381 {
00382
00383 theRecoil->SetDefinition(G4Alpha::Alpha());
00384 }
00385 }
00386 else
00387 {
00388 theRecoil->SetDefinition(G4ParticleTable::GetParticleTable()
00389 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ)));
00390 }
00391 theRecoil->SetMomentum(theTarget.GetMomentum());
00392 theResult.AddSecondary(theRecoil);
00393
00394
00395 theResult.SetStatusChange(suspend);
00396 return &theResult;
00397 }