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