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 "G4LESigmaPlusInelastic.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036
00037 void G4LESigmaPlusInelastic::ModelDescription(std::ostream& outFile) const
00038 {
00039 outFile << "G4LESigmaPlusInelastic 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 G4LESigmaPlusInelastic::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 << "G4LESigmaPlusInelastic::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 (currentParticle.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 G4LESigmaPlusInelastic::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 quasiElastic = true;
00163 return;
00164 }
00165 static G4bool first = true;
00166 const G4int numMul = 1200;
00167 const G4int numSec = 60;
00168 static G4double protmul[numMul], protnorm[numSec];
00169 static G4double neutmul[numMul], neutnorm[numSec];
00170
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.7, 0.7 };
00176 if (first) {
00177 first = false;
00178 G4int i;
00179 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00180 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00181 counter = -1;
00182 for (npos = 0; npos < (numSec/3); ++npos) {
00183 for (nneg = npos; nneg <= (npos+2); ++nneg) {
00184 for (nzero = 0; nzero < numSec/3; ++nzero) {
00185 if (++counter < numMul) {
00186 nt = npos+nneg+nzero;
00187 if (nt > 0 && nt <= numSec) {
00188 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00189 protnorm[nt-1] += protmul[counter];
00190 }
00191 }
00192 }
00193 }
00194 }
00195
00196 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00197 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00198 counter = -1;
00199 for( npos=0; npos<numSec/3; ++npos )
00200 {
00201 for( nneg=std::max(0,npos-1); nneg<=(npos+1); ++nneg )
00202 {
00203 for( nzero=0; nzero<numSec/3; ++nzero )
00204 {
00205 if( ++counter < numMul )
00206 {
00207 nt = npos+nneg+nzero;
00208 if( nt>0 && nt<=numSec )
00209 {
00210 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00211 neutnorm[nt-1] += neutmul[counter];
00212 }
00213 }
00214 }
00215 }
00216 }
00217 for (i = 0; i < numSec; ++i) {
00218 if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
00219 if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
00220 }
00221 }
00222
00223 const G4double expxu = 82.;
00224 const G4double expxl = -expxu;
00225 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00226 G4ParticleDefinition *aProton = G4Proton::Proton();
00227 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
00228 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
00229
00230
00231 G4double n, anpn;
00232 GetNormalizationConstant(availableEnergy, n, anpn);
00233 G4double ran = G4UniformRand();
00234 G4double dum, excs = 0.0;
00235 if (targetParticle.GetDefinition() == aProton) {
00236 counter = -1;
00237 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00238 {
00239 for( nneg=npos; nneg<=(npos+2) && ran>=excs; ++nneg )
00240 {
00241 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00242 {
00243 if( ++counter < numMul )
00244 {
00245 nt = npos+nneg+nzero;
00246 if( nt>0 && nt<=numSec )
00247 {
00248 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00249 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00250 if( std::fabs(dum) < 1.0 )
00251 {
00252 if( test >= 1.0e-10 )excs += dum*test;
00253 }
00254 else
00255 excs += dum*test;
00256 }
00257 }
00258 }
00259 }
00260 }
00261 if( ran >= excs )
00262 {
00263 quasiElastic = true;
00264 return;
00265 }
00266 npos--; nneg--; nzero--;
00267 switch( std::min( 3, std::max( 1, npos-nneg+3 ) ) )
00268 {
00269 case 1:
00270 if( G4UniformRand() < 0.5 )
00271 currentParticle.SetDefinitionAndUpdateE( aLambda );
00272 else
00273 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00274 incidentHasChanged = true;
00275 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00276 targetHasChanged = true;
00277 break;
00278 case 2:
00279 if( G4UniformRand() < 0.5 )
00280 {
00281 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00282 targetHasChanged = true;
00283 }
00284 else
00285 {
00286 if( G4UniformRand() < 0.5 )
00287 currentParticle.SetDefinitionAndUpdateE( aLambda );
00288 else
00289 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00290 incidentHasChanged = true;
00291 }
00292 break;
00293 case 3:
00294 break;
00295 }
00296 }
00297 else
00298 {
00299 counter = -1;
00300 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00301 {
00302 for( nneg=std::max(0,npos-1); nneg<=(npos+1) && ran>=excs; ++nneg )
00303 {
00304 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00305 {
00306 if( ++counter < numMul )
00307 {
00308 nt = npos+nneg+nzero;
00309 if( nt>0 && nt<=numSec )
00310 {
00311 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00312 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00313 if( std::fabs(dum) < 1.0 )
00314 {
00315 if( test >= 1.0e-10 )excs += dum*test;
00316 }
00317 else
00318 excs += dum*test;
00319 }
00320 }
00321 }
00322 }
00323 }
00324 if( ran >= excs )
00325 {
00326 quasiElastic = true;
00327 return;
00328 }
00329 npos--; nneg--; nzero--;
00330 switch( std::min( 3, std::max( 1, npos-nneg+2 ) ) )
00331 {
00332 case 1:
00333 targetParticle.SetDefinitionAndUpdateE( aProton );
00334 targetHasChanged = true;
00335 break;
00336 case 2:
00337 if( G4UniformRand() < 0.5 )
00338 {
00339 if( G4UniformRand() < 0.5 )
00340 {
00341 currentParticle.SetDefinitionAndUpdateE( aLambda );
00342 incidentHasChanged = true;
00343 targetParticle.SetDefinitionAndUpdateE( aProton );
00344 targetHasChanged = true;
00345 }
00346 else
00347 {
00348 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00349 targetHasChanged = true;
00350 }
00351 }
00352 else
00353 {
00354 if( G4UniformRand() < 0.5 )
00355 {
00356 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00357 incidentHasChanged = true;
00358 targetParticle.SetDefinitionAndUpdateE( aProton );
00359 targetHasChanged = true;
00360 }
00361 }
00362 break;
00363 case 3:
00364 if( G4UniformRand() < 0.5 )
00365 currentParticle.SetDefinitionAndUpdateE( aLambda );
00366 else
00367 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00368 incidentHasChanged = true;
00369 break;
00370 }
00371 }
00372 SetUpPions(npos, nneg, nzero, vec, vecLen);
00373 return;
00374 }
00375
00376
00377