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 #include "G4LESigmaMinusInelastic.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036
00037 void G4LESigmaMinusInelastic::ModelDescription(std::ostream& outFile) const
00038 {
00039 outFile << "G4LESigmaMinusInelastic is one of the Low Energy Parameterized\n"
00040 << "(LEP) models used to implement inelastic Sigma- 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 Sigma- with initial energies between 0 and 25\n"
00047 << "GeV.\n";
00048 }
00049
00050 G4HadFinalState*
00051 G4LESigmaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00052 G4Nucleus& targetNucleus)
00053 {
00054 const G4HadProjectile *originalIncident = &aTrack;
00055 if (originalIncident->GetKineticEnergy()<= 0.1*MeV) {
00056 theParticleChange.SetStatusChange(isAlive);
00057 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00058 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00059 return &theParticleChange;
00060 }
00061
00062
00063 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00064
00065 if (verboseLevel > 1) {
00066 const G4Material *targetMaterial = aTrack.GetMaterial();
00067 G4cout << "G4LESigmaMinusInelastic::ApplyYourself called" << G4endl;
00068 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00069 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00070 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00071 << G4endl;
00072 }
00073
00074
00075
00076 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00077 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00078 G4ReactionProduct modifiedOriginal;
00079 modifiedOriginal = *originalIncident;
00080
00081 G4double tkin = targetNucleus.Cinema(ek);
00082 ek += tkin;
00083 modifiedOriginal.SetKineticEnergy(ek*MeV);
00084 G4double et = ek + amas;
00085 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00086 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00087 if (pp > 0.0) {
00088 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00089 modifiedOriginal.SetMomentum(momentum * (p/pp) );
00090 }
00091
00092
00093 tkin = targetNucleus.EvaporationEffects(ek);
00094 ek -= tkin;
00095 modifiedOriginal.SetKineticEnergy(ek*MeV);
00096 et = ek + amas;
00097 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00098 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00099 if (pp > 0.0) {
00100 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00101 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00102 }
00103 G4ReactionProduct currentParticle = modifiedOriginal;
00104 G4ReactionProduct targetParticle;
00105 targetParticle = *originalTarget;
00106 currentParticle.SetSide(1);
00107 targetParticle.SetSide(-1);
00108 G4bool incidentHasChanged = false;
00109 G4bool targetHasChanged = false;
00110 G4bool quasiElastic = false;
00111 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00112 G4int vecLen = 0;
00113 vec.Initialize(0);
00114
00115 const G4double cutOff = 0.1;
00116 if (originalIncident->GetKineticEnergy()/MeV > cutOff)
00117 Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00118 incidentHasChanged, targetHasChanged, quasiElastic);
00119
00120 CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00121 modifiedOriginal, targetNucleus, currentParticle,
00122 targetParticle, incidentHasChanged, targetHasChanged,
00123 quasiElastic);
00124
00125 SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00126
00127 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00128
00129 delete originalTarget;
00130 return &theParticleChange;
00131 }
00132
00133
00134 void G4LESigmaMinusInelastic::Cascade(
00135 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00136 G4int& vecLen,
00137 const G4HadProjectile *originalIncident,
00138 G4ReactionProduct ¤tParticle,
00139 G4ReactionProduct &targetParticle,
00140 G4bool &incidentHasChanged,
00141 G4bool &targetHasChanged,
00142 G4bool &quasiElastic )
00143 {
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00155 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00156 const G4double targetMass = targetParticle.GetMass()/MeV;
00157 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00158 targetMass*targetMass +
00159 2.0*targetMass*etOriginal );
00160 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00161 if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV )
00162 {
00163 quasiElastic = true;
00164 return;
00165 }
00166 static G4bool first = true;
00167 const G4int numMul = 1200;
00168 const G4int numSec = 60;
00169 static G4double protmul[numMul], protnorm[numSec];
00170 static G4double neutmul[numMul], neutnorm[numSec];
00171
00172 G4int counter, nt=0, npos=0, nneg=0, nzero=0;
00173 G4double test;
00174 const G4double c = 1.25;
00175 const G4double b[] = { 0.70, 0.70 };
00176 if( first )
00177 {
00178 first = false;
00179 G4int i;
00180 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00181 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00182 counter = -1;
00183 for( npos=0; npos<(numSec/3); ++npos )
00184 {
00185 for( nneg=std::max(0,npos-1); nneg<=(npos+1); ++nneg )
00186 {
00187 for( nzero=0; nzero<numSec/3; ++nzero )
00188 {
00189 if( ++counter < numMul )
00190 {
00191 nt = npos+nneg+nzero;
00192 if( nt>0 && nt<=numSec )
00193 {
00194 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00195 protnorm[nt-1] += protmul[counter];
00196 }
00197 }
00198 }
00199 }
00200 }
00201 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00202 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00203 counter = -1;
00204 for( npos=0; npos<numSec/3; ++npos )
00205 {
00206 for( nneg=npos; nneg<=(npos+2); ++nneg )
00207 {
00208 for( nzero=0; nzero<numSec/3; ++nzero )
00209 {
00210 if( ++counter < numMul )
00211 {
00212 nt = npos+nneg+nzero;
00213 if( nt>0 && nt<=numSec )
00214 {
00215 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00216 neutnorm[nt-1] += neutmul[counter];
00217 }
00218 }
00219 }
00220 }
00221 }
00222 for( i=0; i<numSec; ++i )
00223 {
00224 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00225 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00226 }
00227 }
00228
00229 const G4double expxu = 82.;
00230 const G4double expxl = -expxu;
00231 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00232 G4ParticleDefinition *aProton = G4Proton::Proton();
00233 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
00234 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
00235
00236
00237
00238 G4double n, anpn;
00239 GetNormalizationConstant( availableEnergy, n, anpn );
00240 G4double ran = G4UniformRand();
00241 G4double dum, excs = 0.0;
00242 if( targetParticle.GetDefinition() == aProton )
00243 {
00244 counter = -1;
00245 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00246 {
00247 for( nneg=std::max(0,npos-1); nneg<=(npos+1) && ran>=excs; ++nneg )
00248 {
00249 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00250 {
00251 if( ++counter < numMul )
00252 {
00253 nt = npos+nneg+nzero;
00254 if( nt>0 && nt<=numSec )
00255 {
00256 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00257 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00258 if( std::fabs(dum) < 1.0 )
00259 {
00260 if( test >= 1.0e-10 )excs += dum*test;
00261 }
00262 else
00263 excs += dum*test;
00264 }
00265 }
00266 }
00267 }
00268 }
00269 if( ran >= excs )
00270 {
00271 quasiElastic = true;
00272 return;
00273 }
00274 npos--; nneg--; nzero--;
00275 G4int ncht = std::max( 1, npos-nneg+2 );
00276 switch( ncht )
00277 {
00278 case 1:
00279 if( G4UniformRand() < 0.5 )
00280 currentParticle.SetDefinitionAndUpdateE( aLambda );
00281 else
00282 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00283 incidentHasChanged = true;
00284 break;
00285 case 2:
00286 if( G4UniformRand() >= 0.5 )
00287 {
00288 if( G4UniformRand() < 0.5 )
00289 currentParticle.SetDefinitionAndUpdateE( aLambda );
00290 else
00291 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00292 incidentHasChanged = true;
00293 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00294 targetHasChanged = true;
00295 }
00296 break;
00297 default:
00298 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00299 targetHasChanged = true;
00300 break;
00301 }
00302 }
00303 else
00304 {
00305 counter = -1;
00306 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00307 {
00308 for( nneg=npos; nneg<=(npos+2) && ran>=excs; ++nneg )
00309 {
00310 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00311 {
00312 if( ++counter < numMul )
00313 {
00314 nt = npos+nneg+nzero;
00315 if( nt>0 && nt<=numSec )
00316 {
00317 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00318 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00319 if( std::fabs(dum) < 1.0 )
00320 {
00321 if( test >= 1.0e-10 )excs += dum*test;
00322 }
00323 else
00324 excs += dum*test;
00325 }
00326 }
00327 }
00328 }
00329 }
00330 if( ran >= excs )
00331 {
00332 quasiElastic = true;
00333 return;
00334 }
00335 npos--; nneg--; nzero--;
00336 G4int ncht = std::max( 1, npos-nneg+3 );
00337 switch( ncht )
00338 {
00339 case 1:
00340 if( G4UniformRand() < 0.5 )
00341 currentParticle.SetDefinitionAndUpdateE( aLambda );
00342 else
00343 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00344 incidentHasChanged = true;
00345 targetParticle.SetDefinitionAndUpdateE( aProton );
00346 targetHasChanged = true;
00347 break;
00348 case 2:
00349 if( G4UniformRand() < 0.5 )
00350 {
00351 if( G4UniformRand() < 0.5 )
00352 currentParticle.SetDefinitionAndUpdateE( aLambda );
00353 else
00354 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00355 incidentHasChanged = true;
00356 }
00357 else
00358 {
00359 targetParticle.SetDefinitionAndUpdateE( aProton );
00360 targetHasChanged = true;
00361 }
00362 break;
00363 default:
00364 break;
00365 }
00366 }
00367 SetUpPions( npos, nneg, nzero, vec, vecLen );
00368 return;
00369 }
00370
00371
00372