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