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 <iostream>
00033
00034 #include "G4LEKaonPlusInelastic.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "Randomize.hh"
00038
00039 G4LEKaonPlusInelastic::G4LEKaonPlusInelastic(const G4String& name)
00040 :G4InelasticInteraction(name)
00041 {
00042 SetMinEnergy(0.0);
00043 SetMaxEnergy(25.*GeV);
00044 G4cout << "WARNING: model G4LEKaonPlusInelastic is being deprecated and will\n"
00045 << "disappear in Geant4 version 10.0" << G4endl;
00046 }
00047
00048
00049 void G4LEKaonPlusInelastic::ModelDescription(std::ostream& outFile) const
00050 {
00051 outFile << "G4LEKaonPlusInelastic is one of the Low Energy Parameterized\n"
00052 << "(LEP) models used to implement inelastic K+ scattering\n"
00053 << "from nuclei. It is a re-engineered version of the GHEISHA\n"
00054 << "code of H. Fesefeldt. It divides the initial collision\n"
00055 << "products into backward- and forward-going clusters which are\n"
00056 << "then decayed into final state hadrons. The model does not\n"
00057 << "conserve energy on an event-by-event basis. It may be\n"
00058 << "applied to kaons with initial energies between 0 and 25\n"
00059 << "GeV.\n";
00060 }
00061
00062
00063 G4HadFinalState*
00064 G4LEKaonPlusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00065 G4Nucleus& targetNucleus)
00066 {
00067 const G4HadProjectile *originalIncident = &aTrack;
00068 if (originalIncident->GetKineticEnergy()<= 0.1*MeV) {
00069 theParticleChange.SetStatusChange(isAlive);
00070 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00071 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00072 return &theParticleChange;
00073 }
00074
00075
00076 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00077 G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
00078
00079 if (verboseLevel > 1) {
00080 const G4Material *targetMaterial = aTrack.GetMaterial();
00081 G4cout << "G4LEKaonPlusInelastic::ApplyYourself called" << G4endl;
00082 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
00083 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00084 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00085 << G4endl;
00086 }
00087 G4ReactionProduct currentParticle( const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition()));
00088 currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
00089 currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
00090
00091
00092
00093 G4double ek = originalIncident->GetKineticEnergy();
00094 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00095
00096 G4double tkin = targetNucleus.Cinema(ek);
00097 ek += tkin;
00098 currentParticle.SetKineticEnergy( ek );
00099 G4double et = ek + amas;
00100 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00101 G4double pp = currentParticle.GetMomentum().mag();
00102 if (pp > 0.0) {
00103 G4ThreeVector momentum = currentParticle.GetMomentum();
00104 currentParticle.SetMomentum( momentum * (p/pp) );
00105 }
00106
00107
00108 tkin = targetNucleus.EvaporationEffects( ek );
00109 ek -= tkin;
00110 currentParticle.SetKineticEnergy( ek );
00111 et = ek + amas;
00112 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00113 pp = currentParticle.GetMomentum().mag();
00114 if (pp > 0.0) {
00115 G4ThreeVector momentum = currentParticle.GetMomentum();
00116 currentParticle.SetMomentum( momentum * (p/pp) );
00117 }
00118
00119 G4ReactionProduct modifiedOriginal = currentParticle;
00120
00121 currentParticle.SetSide(1);
00122 targetParticle.SetSide(-1);
00123 G4bool incidentHasChanged = false;
00124 G4bool targetHasChanged = false;
00125 G4bool quasiElastic = false;
00126 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00127 G4int vecLen = 0;
00128 vec.Initialize( 0 );
00129
00130 const G4double cutOff = 0.1*MeV;
00131 if (currentParticle.GetKineticEnergy() > cutOff)
00132 Cascade(vec, vecLen, originalIncident, currentParticle,
00133 targetParticle, incidentHasChanged, targetHasChanged,
00134 quasiElastic);
00135
00136 CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00137 modifiedOriginal, targetNucleus, currentParticle,
00138 targetParticle, incidentHasChanged, targetHasChanged,
00139 quasiElastic);
00140
00141 SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00142
00143 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00144
00145 delete originalTarget;
00146 return &theParticleChange;
00147 }
00148
00149
00150 void G4LEKaonPlusInelastic::Cascade(
00151 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00152 G4int &vecLen,
00153 const G4HadProjectile *originalIncident,
00154 G4ReactionProduct ¤tParticle,
00155 G4ReactionProduct &targetParticle,
00156 G4bool &incidentHasChanged,
00157 G4bool &targetHasChanged,
00158 G4bool &quasiElastic)
00159 {
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
00171 const G4double etOriginal = originalIncident->GetTotalEnergy();
00172 const G4double targetMass = targetParticle.GetMass();
00173 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00174 targetMass*targetMass +
00175 2.0*targetMass*etOriginal );
00176 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00177 if( availableEnergy < G4PionPlus::PionPlus()->GetPDGMass() )
00178 {
00179 quasiElastic = true;
00180 return;
00181 }
00182 static G4bool first = true;
00183 const G4int numMul = 1200;
00184 const G4int numSec = 60;
00185 static G4double protmul[numMul], protnorm[numSec];
00186 static G4double neutmul[numMul], neutnorm[numSec];
00187
00188 G4int nt=0, npos=0, nneg=0, nzero=0;
00189 const G4double c = 1.25;
00190 const G4double b[] = { 0.70, 0.70 };
00191 if( first )
00192 {
00193 first = false;
00194 G4int i;
00195 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00196 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00197 G4int counter = -1;
00198 for( npos=0; npos<(numSec/3); ++npos )
00199 {
00200 for( nneg=std::max(0,npos-2); nneg<=npos; ++nneg )
00201 {
00202 for( nzero=0; nzero<numSec/3; ++nzero )
00203 {
00204 if( ++counter < numMul )
00205 {
00206 nt = npos+nneg+nzero;
00207 if( nt > 0 )
00208 {
00209 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00210 protnorm[nt-1] += protmul[counter];
00211 }
00212 }
00213 }
00214 }
00215 }
00216 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00217 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00218 counter = -1;
00219 for( npos=0; npos<numSec/3; ++npos )
00220 {
00221 for( nneg=std::max(0,npos-1); nneg<=(npos+1); ++nneg )
00222 {
00223 for( nzero=0; nzero<numSec/3; ++nzero )
00224 {
00225 if( ++counter < numMul )
00226 {
00227 nt = npos+nneg+nzero;
00228 if( (nt>0) && (nt<=numSec) )
00229 {
00230 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00231 neutnorm[nt-1] += neutmul[counter];
00232 }
00233 }
00234 }
00235 }
00236 }
00237 for( i=0; i<numSec; ++i )
00238 {
00239 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00240 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00241 }
00242 }
00243
00244 const G4double expxu = 82.;
00245 const G4double expxl = -expxu;
00246 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
00247 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
00248 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00249 G4ParticleDefinition *aProton = G4Proton::Proton();
00250 G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
00251 const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
00252 G4double test, w0, wp, wt, wm;
00253 if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
00254 {
00255
00256
00257
00258 nneg = npos = nzero = 0;
00259 if( targetParticle.GetDefinition() == aProton )
00260 {
00261 test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[0])/(2.0*c*c) ) ) );
00262 w0 = test;
00263 wp = test*2.0;
00264 if( G4UniformRand() < w0/(w0+wp) )
00265 nzero = 1;
00266 else
00267 npos = 1;
00268 }
00269 else
00270 {
00271 test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[1])/(2.0*c*c) ) ) );
00272 w0 = test;
00273 wp = test;
00274 test = std::exp( std::min( expxu, std::max( expxl, -sqr(-1.0+b[1])/(2.0*c*c) ) ) );
00275 wm = test;
00276 wt = w0+wp+wm;
00277 wp += w0;
00278 G4double ran = G4UniformRand();
00279 if( ran < w0/wt )
00280 nzero = 1;
00281 else if( ran < wp/wt )
00282 npos = 1;
00283 else
00284 nneg = 1;
00285 }
00286 }
00287 else
00288 {
00289 G4double n, anpn;
00290 GetNormalizationConstant( availableEnergy, n, anpn );
00291 G4double ran = G4UniformRand();
00292 G4double dum, excs = 0.0;
00293 if( targetParticle.GetDefinition() == aProton )
00294 {
00295 G4int counter = -1;
00296 for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos )
00297 {
00298 for( nneg=std::max(0,npos-2); (nneg<=npos) && (ran>=excs); ++nneg )
00299 {
00300 for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero )
00301 {
00302 if( ++counter < numMul )
00303 {
00304 nt = npos+nneg+nzero;
00305 if( nt > 0 )
00306 {
00307 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00308 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00309 if( std::fabs(dum) < 1.0 )
00310 {
00311 if( test >= 1.0e-10 )excs += dum*test;
00312 }
00313 else
00314 excs += dum*test;
00315 }
00316 }
00317 }
00318 }
00319 }
00320 if( ran >= excs )return;
00321 npos--; nneg--; nzero--;
00322 }
00323 else
00324 {
00325 G4int counter = -1;
00326 for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos )
00327 {
00328 for( nneg=std::max(0,npos-1); (nneg<=(npos+1)) && (ran>=excs); ++nneg )
00329 {
00330 for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero )
00331 {
00332 if( ++counter < numMul )
00333 {
00334 nt = npos+nneg+nzero;
00335 if( (nt>=1) && (nt<=numSec) )
00336 {
00337 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00338 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00339 if( std::fabs(dum) < 1.0 )
00340 {
00341 if( test >= 1.0e-10 )excs += dum*test;
00342 }
00343 else
00344 excs += dum*test;
00345 }
00346 }
00347 }
00348 }
00349 }
00350 if( ran >= excs )return;
00351 npos--; nneg--; nzero--;
00352 }
00353 }
00354 if( targetParticle.GetDefinition() == aProton )
00355 {
00356 switch( npos-nneg )
00357 {
00358 case 1:
00359 if( G4UniformRand() < 0.5 )
00360 {
00361 if( G4UniformRand() < 0.5 )
00362 currentParticle.SetDefinitionAndUpdateE( aKaonZS );
00363 else
00364 currentParticle.SetDefinitionAndUpdateE( aKaonZL );
00365 incidentHasChanged = true;
00366 }
00367 else
00368 {
00369 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00370 targetHasChanged = true;
00371 }
00372 break;
00373 case 2:
00374 if( G4UniformRand() < 0.5 )
00375 currentParticle.SetDefinitionAndUpdateE( aKaonZS );
00376 else
00377 currentParticle.SetDefinitionAndUpdateE( aKaonZL );
00378 incidentHasChanged = true;
00379 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00380 incidentHasChanged = true;
00381 targetHasChanged = true;
00382 break;
00383 default:
00384 break;
00385 }
00386 }
00387 else
00388 {
00389 switch( npos-nneg )
00390 {
00391 case 0:
00392 if( G4UniformRand() < 0.25 )
00393 {
00394 if( G4UniformRand() < 0.5 )
00395 currentParticle.SetDefinitionAndUpdateE( aKaonZS );
00396 else
00397 currentParticle.SetDefinitionAndUpdateE( aKaonZL );
00398 targetParticle.SetDefinitionAndUpdateE( aProton );
00399 incidentHasChanged = true;
00400 targetHasChanged = true;
00401 }
00402 break;
00403 case 1:
00404 if( G4UniformRand() < 0.5 )
00405 currentParticle.SetDefinitionAndUpdateE( aKaonZS );
00406 else
00407 currentParticle.SetDefinitionAndUpdateE( aKaonZL );
00408 incidentHasChanged = true;
00409 break;
00410 default:
00411 targetParticle.SetDefinitionAndUpdateE( aProton );
00412 targetHasChanged = true;
00413 break;
00414 }
00415 }
00416 SetUpPions( npos, nneg, nzero, vec, vecLen );
00417 return;
00418 }
00419
00420
00421