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 #include <iostream>
00030
00031 #include "G4LENeutronInelastic.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "Randomize.hh"
00035 #include "G4Electron.hh"
00036
00037 void G4LENeutronInelastic::ModelDescription(std::ostream& outFile) const
00038 {
00039 outFile << "G4LENeutronInelastic is one of the Low Energy Parameterized\n"
00040 << "(LEP) models used to implement inelastic neutron scattering\n"
00041 << "from nuclei. It is a re-engineered version of the GHEISHA\n"
00042 << "code of H. Fesefeldt. It divides the initial collision\n"
00043 << "products into backward- and forward-going clusters which are\n"
00044 << "then decayed into final state hadrons. The model does not\n"
00045 << "conserve energy on an event-by-event basis. It may be\n"
00046 << "applied to neutrons with initial energies between 0 and 25\n"
00047 << "GeV.\n";
00048 }
00049
00050
00051 G4HadFinalState*
00052 G4LENeutronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00053 G4Nucleus& targetNucleus)
00054 {
00055 theParticleChange.Clear();
00056 const G4HadProjectile *originalIncident = &aTrack;
00057
00058
00059 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00060
00061 if (verboseLevel > 1) {
00062 const G4Material* targetMaterial = aTrack.GetMaterial();
00063 G4cout << "G4LENeutronInelastic::ApplyYourself called" << G4endl;
00064 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00065 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00066 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00067 << G4endl;
00068 }
00069
00070 G4ReactionProduct modifiedOriginal;
00071 modifiedOriginal = *originalIncident;
00072 G4ReactionProduct targetParticle;
00073 targetParticle = *originalTarget;
00074 if (originalIncident->GetKineticEnergy()/GeV < 0.01 + 2.*G4UniformRand()/9.) {
00075 SlowNeutron(originalIncident, modifiedOriginal, targetParticle, targetNucleus);
00076 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00077 delete originalTarget;
00078 return &theParticleChange;
00079 }
00080
00081
00082
00083 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00084 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00085
00086 G4double tkin = targetNucleus.Cinema(ek);
00087 ek += tkin;
00088 modifiedOriginal.SetKineticEnergy( ek*MeV );
00089 G4double et = ek + amas;
00090 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00091 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00092 if (pp > 0.0) {
00093 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00094 modifiedOriginal.SetMomentum(momentum * (p/pp) );
00095 }
00096
00097
00098 tkin = targetNucleus.EvaporationEffects( ek );
00099 ek -= tkin;
00100 modifiedOriginal.SetKineticEnergy( ek*MeV );
00101 et = ek + amas;
00102 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00103 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00104 if (pp > 0.0) {
00105 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00106 modifiedOriginal.SetMomentum(momentum * (p/pp) );
00107 }
00108 const G4double cutOff = 0.1;
00109 if (modifiedOriginal.GetKineticEnergy()/MeV <= cutOff) {
00110 SlowNeutron(originalIncident, modifiedOriginal, targetParticle, targetNucleus);
00111 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00112 delete originalTarget;
00113 return &theParticleChange;
00114 }
00115
00116 G4ReactionProduct currentParticle = modifiedOriginal;
00117 currentParticle.SetSide(1);
00118 targetParticle.SetSide(-1);
00119 G4bool incidentHasChanged = false;
00120 G4bool targetHasChanged = false;
00121 G4bool quasiElastic = false;
00122 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00123 G4int vecLen = 0;
00124 vec.Initialize(0);
00125
00126 Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00127 incidentHasChanged, targetHasChanged, quasiElastic);
00128
00129 CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00130 modifiedOriginal, targetNucleus, currentParticle,
00131 targetParticle, incidentHasChanged, targetHasChanged,
00132 quasiElastic);
00133
00134 SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00135
00136 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00137 delete originalTarget;
00138 return &theParticleChange;
00139 }
00140
00141
00142 void G4LENeutronInelastic::SlowNeutron(const G4HadProjectile* originalIncident,
00143 G4ReactionProduct& modifiedOriginal,
00144 G4ReactionProduct& targetParticle,
00145 G4Nucleus& targetNucleus)
00146 {
00147 const G4double A = targetNucleus.GetA_asInt();
00148 const G4double Z = targetNucleus.GetZ_asInt();
00149
00150 G4double currentKinetic = modifiedOriginal.GetKineticEnergy()/MeV;
00151 G4double currentMass = modifiedOriginal.GetMass()/MeV;
00152 if( A < 1.5 )
00153 {
00154
00155
00156
00157
00158
00159 G4double cost1, eka = 0.0;
00160 while (eka <= 0.0)
00161 {
00162 cost1 = -1.0 + 2.0*G4UniformRand();
00163 eka = 1.0 + 2.0*cost1*A + A*A;
00164 }
00165 G4double cost = std::min( 1.0, std::max( -1.0, (A*cost1+1.0)/std::sqrt(eka) ) );
00166 eka /= (1.0+A)*(1.0+A);
00167 G4double ek = currentKinetic*MeV/GeV;
00168 G4double amas = currentMass*MeV/GeV;
00169 ek *= eka;
00170 G4double en = ek + amas;
00171 G4double p = std::sqrt(std::abs(en*en-amas*amas));
00172 G4double sint = std::sqrt(std::abs(1.0-cost*cost));
00173 G4double phi = G4UniformRand()*twopi;
00174 G4double px = sint*std::sin(phi);
00175 G4double py = sint*std::cos(phi);
00176 G4double pz = cost;
00177 targetParticle.SetMomentum( px*GeV, py*GeV, pz*GeV );
00178 G4double pxO = originalIncident->Get4Momentum().x()/GeV;
00179 G4double pyO = originalIncident->Get4Momentum().y()/GeV;
00180 G4double pzO = originalIncident->Get4Momentum().z()/GeV;
00181 G4double ptO = pxO*pxO + pyO+pyO;
00182 if( ptO > 0.0 )
00183 {
00184 G4double pO = std::sqrt(pxO*pxO+pyO*pyO+pzO*pzO);
00185 cost = pzO/pO;
00186 sint = 0.5*(std::sqrt(std::abs((1.0-cost)*(1.0+cost)))+std::sqrt(ptO)/pO);
00187 G4double ph = pi/2.0;
00188 if( pyO < 0.0 )ph = ph*1.5;
00189 if( std::abs(pxO) > 0.000001 )ph = std::atan2(pyO,pxO);
00190 G4double cosp = std::cos(ph);
00191 G4double sinp = std::sin(ph);
00192 px = cost*cosp*px - sinp*py+sint*cosp*pz;
00193 py = cost*sinp*px + cosp*py+sint*sinp*pz;
00194 pz = -sint*px + cost*pz;
00195 }
00196 else
00197 {
00198 if( pz < 0.0 )pz *= -1.0;
00199 }
00200 G4double pu = std::sqrt(px*px+py*py+pz*pz);
00201 modifiedOriginal.SetMomentum( targetParticle.GetMomentum() * (p/pu) );
00202 modifiedOriginal.SetKineticEnergy( ek*GeV );
00203
00204 targetParticle.SetMomentum(
00205 originalIncident->Get4Momentum().vect() - modifiedOriginal.GetMomentum() );
00206 G4double pp = targetParticle.GetMomentum().mag();
00207 G4double tarmas = targetParticle.GetMass();
00208 targetParticle.SetTotalEnergy( std::sqrt( pp*pp + tarmas*tarmas ) );
00209
00210 theParticleChange.SetEnergyChange( modifiedOriginal.GetKineticEnergy() );
00211 G4DynamicParticle *pd = new G4DynamicParticle;
00212 pd->SetDefinition( targetParticle.GetDefinition() );
00213 pd->SetMomentum( targetParticle.GetMomentum() );
00214 theParticleChange.AddSecondary( pd );
00215 return;
00216 }
00217 G4FastVector<G4ReactionProduct,4> vec;
00218 G4int vecLen = 0;
00219 vec.Initialize( 0 );
00220
00221 G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
00222 G4double massVec[9];
00223 massVec[0] = targetNucleus.AtomicMass( A+1.0, Z );
00224 massVec[1] = theAtomicMass;
00225 massVec[2] = 0.;
00226 if (Z > 1.0)
00227 massVec[2] = targetNucleus.AtomicMass( A , Z-1.0 );
00228 massVec[3] = 0.;
00229 if (Z > 1.0 && A > 1.0)
00230 massVec[3] = targetNucleus.AtomicMass( A-1.0, Z-1.0 );
00231 massVec[4] = 0.;
00232 if (Z > 1.0 && A > 2.0 && A-2.0 > Z-1.0)
00233 massVec[4] = targetNucleus.AtomicMass( A-2.0, Z-1.0 );
00234 massVec[5] = 0.;
00235 if (Z > 2.0 && A > 3.0 && A-3.0 > Z-2.0)
00236 massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-2.0 );
00237 massVec[6] = 0.;
00238 if (A > 1.0 && A-1.0 > Z)
00239 massVec[6] = targetNucleus.AtomicMass( A-1.0, Z );
00240 massVec[7] = massVec[3];
00241 massVec[8] = 0.;
00242 if (Z > 2.0 && A > 1.0)
00243 massVec[8] = targetNucleus.AtomicMass( A-1.0, Z-2.0 );
00244
00245 theReactionDynamics.NuclearReaction( vec, vecLen, originalIncident,
00246 targetNucleus, theAtomicMass, massVec );
00247
00248 theParticleChange.SetStatusChange( stopAndKill );
00249 theParticleChange.SetEnergyChange( 0.0 );
00250
00251 G4DynamicParticle * pd;
00252 for( G4int i=0; i<vecLen; ++i )
00253 {
00254 pd = new G4DynamicParticle();
00255 pd->SetDefinition( vec[i]->GetDefinition() );
00256 pd->SetMomentum( vec[i]->GetMomentum() );
00257 theParticleChange.AddSecondary( pd );
00258 delete vec[i];
00259 }
00260 }
00261
00262
00263 void G4LENeutronInelastic::Cascade(
00264 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00265 G4int& vecLen,
00266 const G4HadProjectile *originalIncident,
00267 G4ReactionProduct ¤tParticle,
00268 G4ReactionProduct &targetParticle,
00269 G4bool &incidentHasChanged,
00270 G4bool &targetHasChanged,
00271 G4bool &quasiElastic )
00272 {
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00284 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00285 const G4double targetMass = targetParticle.GetMass()/MeV;
00286 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00287 targetMass*targetMass +
00288 2.0*targetMass*etOriginal );
00289 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00290 if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV )
00291 {
00292 quasiElastic = true;
00293 return;
00294 }
00295 static G4bool first = true;
00296 const G4int numMul = 1200;
00297 const G4int numSec = 60;
00298 static G4double protmul[numMul], protnorm[numSec];
00299 static G4double neutmul[numMul], neutnorm[numSec];
00300
00301 G4int counter, nt=0, npos=0, nneg=0, nzero=0;
00302 const G4double c = 1.25;
00303 const G4double b[] = { 0.35, 0.0 };
00304 if( first )
00305 {
00306 first = false;
00307 G4int i;
00308 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00309 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00310 counter = -1;
00311 for( npos=0; npos<numSec/3; ++npos )
00312 {
00313 for( nneg=std::max(0,npos-1); nneg<=(npos+1); ++nneg )
00314 {
00315 for( nzero=0; nzero<numSec/3; ++nzero )
00316 {
00317 if( ++counter < numMul )
00318 {
00319 nt = npos+nneg+nzero;
00320 if( nt > 0 )
00321 {
00322 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c) /
00323 ( theReactionDynamics.Factorial(1-npos+nneg)*
00324 theReactionDynamics.Factorial(1+npos-nneg) );
00325 protnorm[nt-1] += protmul[counter];
00326 }
00327 }
00328 }
00329 }
00330 }
00331 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00332 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00333 counter = -1;
00334 for( npos=0; npos<(numSec/3); ++npos )
00335 {
00336 for( nneg=npos; nneg<=(npos+2); ++nneg )
00337 {
00338 for( nzero=0; nzero<numSec/3; ++nzero )
00339 {
00340 if( ++counter < numMul )
00341 {
00342 nt = npos+nneg+nzero;
00343 if( (nt>0) && (nt<=numSec) )
00344 {
00345 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c) /
00346 ( theReactionDynamics.Factorial(nneg-npos)*
00347 theReactionDynamics.Factorial(2-nneg+npos) );
00348 neutnorm[nt-1] += neutmul[counter];
00349 }
00350 }
00351 }
00352 }
00353 }
00354 for( i=0; i<numSec; ++i )
00355 {
00356 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00357 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00358 }
00359 }
00360
00361 const G4double expxu = 82.;
00362 const G4double expxl = -expxu;
00363 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00364 G4ParticleDefinition *aProton = G4Proton::Proton();
00365 G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
00366 const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
00367 G4double test, w0, wp, wt, wm;
00368 if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
00369 {
00370
00371
00372
00373 nneg = npos = nzero = 0;
00374 if( targetParticle.GetDefinition() == aNeutron )
00375 {
00376 test = std::exp( std::min( expxu, std::max( expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
00377 w0 = test/2.0;
00378 wm = test;
00379 if( G4UniformRand() < w0/(w0+wm) )
00380 nzero = 1;
00381 else
00382 nneg = 1;
00383 } else {
00384 test = std::exp( std::min( expxu, std::max( expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
00385 w0 = test;
00386 wp = test/2.0;
00387 test = std::exp( std::min( expxu, std::max( expxl, -(-1.0+b[0])*(-1.0+b[0])/(2.0*c*c) ) ) );
00388 wm = test/2.0;
00389 wt = w0+wp+wm;
00390 wp += w0;
00391 G4double ran = G4UniformRand();
00392 if( ran < w0/wt )
00393 nzero = 1;
00394 else if( ran < wp/wt )
00395 npos = 1;
00396 else
00397 nneg = 1;
00398 }
00399 } else {
00400 G4double n, anpn;
00401 GetNormalizationConstant( availableEnergy, n, anpn );
00402 G4double ran = G4UniformRand();
00403 G4double dum, excs = 0.0;
00404 if( targetParticle.GetDefinition() == aProton )
00405 {
00406 counter = -1;
00407 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00408 {
00409 for( nneg=std::max(0,npos-1); nneg<=(npos+1) && ran>=excs; ++nneg )
00410 {
00411 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00412 {
00413 if( ++counter < numMul )
00414 {
00415 nt = npos+nneg+nzero;
00416 if( nt > 0 )
00417 {
00418 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00419 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00420 if( std::fabs(dum) < 1.0 ) {
00421 if( test >= 1.0e-10 )excs += dum*test;
00422 } else {
00423 excs += dum*test;
00424 }
00425 }
00426 }
00427 }
00428 }
00429 }
00430 if( ran >= excs )
00431 {
00432 quasiElastic = true;
00433 return;
00434 }
00435 npos--; nneg--; nzero--;
00436 } else {
00437 counter = -1;
00438 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00439 {
00440 for( nneg=npos; nneg<=(npos+2) && ran>=excs; ++nneg )
00441 {
00442 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00443 {
00444 if( ++counter < numMul )
00445 {
00446 nt = npos+nneg+nzero;
00447 if( (nt>=1) && (nt<=numSec) )
00448 {
00449 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00450 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00451 if( std::fabs(dum) < 1.0 ) {
00452 if( test >= 1.0e-10 )excs += dum*test;
00453 } else {
00454 excs += dum*test;
00455 }
00456 }
00457 }
00458 }
00459 }
00460 }
00461 if( ran >= excs )
00462 {
00463 quasiElastic = true;
00464 return;
00465 }
00466 npos--; nneg--; nzero--;
00467 }
00468 }
00469 if( targetParticle.GetDefinition() == aProton )
00470 {
00471 switch( npos-nneg )
00472 {
00473 case 0:
00474 if( G4UniformRand() < 0.33 )
00475 {
00476 currentParticle.SetDefinitionAndUpdateE( aProton );
00477 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00478 incidentHasChanged = true;
00479 targetHasChanged = true;
00480 }
00481 break;
00482 case 1:
00483 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00484 targetHasChanged = true;
00485 break;
00486 default:
00487 currentParticle.SetDefinitionAndUpdateE( aProton );
00488 incidentHasChanged = true;
00489 break;
00490 }
00491 } else {
00492 switch( npos-nneg )
00493 {
00494 case -1:
00495 if( G4UniformRand() < 0.5 )
00496 {
00497 currentParticle.SetDefinitionAndUpdateE( aProton );
00498 incidentHasChanged = true;
00499 } else {
00500 targetParticle.SetDefinitionAndUpdateE( aProton );
00501 targetHasChanged = true;
00502 }
00503 break;
00504 case 0:
00505 break;
00506 default:
00507 currentParticle.SetDefinitionAndUpdateE( aProton );
00508 targetParticle.SetDefinitionAndUpdateE( aProton );
00509 incidentHasChanged = true;
00510 targetHasChanged = true;
00511 break;
00512 }
00513 }
00514 SetUpPions( npos, nneg, nzero, vec, vecLen );
00515
00516 return;
00517 }
00518
00519
00520