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