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 "G4RPGLambdaInelastic.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "Randomize.hh"
00033
00034 G4HadFinalState*
00035 G4RPGLambdaInelastic::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 << "G4RPGLambdaInelastic::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
00088 G4ReactionProduct currentParticle = modifiedOriginal;
00089 G4ReactionProduct targetParticle;
00090 targetParticle = *originalTarget;
00091 currentParticle.SetSide( 1 );
00092 targetParticle.SetSide( -1 );
00093 G4bool incidentHasChanged = false;
00094 G4bool targetHasChanged = false;
00095 G4bool quasiElastic = false;
00096 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00097 G4int vecLen = 0;
00098 vec.Initialize( 0 );
00099
00100 const G4double cutOff = 0.1;
00101 if( currentParticle.GetKineticEnergy()/MeV > cutOff )
00102 Cascade( vec, vecLen,
00103 originalIncident, currentParticle, targetParticle,
00104 incidentHasChanged, targetHasChanged, quasiElastic );
00105
00106 CalculateMomenta( vec, vecLen,
00107 originalIncident, originalTarget, modifiedOriginal,
00108 targetNucleus, currentParticle, targetParticle,
00109 incidentHasChanged, targetHasChanged, quasiElastic );
00110
00111 SetUpChange( vec, vecLen,
00112 currentParticle, targetParticle,
00113 incidentHasChanged );
00114
00115 delete originalTarget;
00116 return &theParticleChange;
00117 }
00118
00119
00120 void G4RPGLambdaInelastic::Cascade(
00121 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00122 G4int& vecLen,
00123 const G4HadProjectile *originalIncident,
00124 G4ReactionProduct ¤tParticle,
00125 G4ReactionProduct &targetParticle,
00126 G4bool &incidentHasChanged,
00127 G4bool &targetHasChanged,
00128 G4bool &quasiElastic )
00129 {
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00141 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00142 const G4double targetMass = targetParticle.GetMass()/MeV;
00143 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00144 targetMass*targetMass +
00145 2.0*targetMass*etOriginal );
00146 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00147 if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV )
00148 {
00149 quasiElastic = true;
00150 return;
00151 }
00152 static G4bool first = true;
00153 const G4int numMul = 1200;
00154 const G4int numSec = 60;
00155 static G4double protmul[numMul], protnorm[numSec];
00156 static G4double neutmul[numMul], neutnorm[numSec];
00157
00158
00159
00160 G4int counter, nt=0, np=0, nneg=0, nz=0;
00161 G4double test;
00162 const G4double c = 1.25;
00163 const G4double b[] = { 0.70, 0.35 };
00164 if( first ) {
00165 first = false;
00166 G4int i;
00167 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00168 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00169 counter = -1;
00170 for( np=0; np<(numSec/3); ++np ) {
00171 for( nneg=std::max(0,np-2); nneg<=(np+1); ++nneg ) {
00172 for( nz=0; nz<numSec/3; ++nz ) {
00173 if( ++counter < numMul ) {
00174 nt = np+nneg+nz;
00175 if( nt>0 && nt<=numSec ) {
00176 protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
00177 protnorm[nt-1] += protmul[counter];
00178 }
00179 }
00180 }
00181 }
00182 }
00183 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00184 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00185 counter = -1;
00186 for( np=0; np<numSec/3; ++np ) {
00187 for( nneg=std::max(0,np-1); nneg<=(np+2); ++nneg ) {
00188 for( nz=0; nz<numSec/3; ++nz ) {
00189 if( ++counter < numMul ) {
00190 nt = np+nneg+nz;
00191 if( nt>0 && nt<=numSec ) {
00192 neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
00193 neutnorm[nt-1] += neutmul[counter];
00194 }
00195 }
00196 }
00197 }
00198 }
00199 for( i=0; i<numSec; ++i ) {
00200 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00201 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00202 }
00203 }
00204
00205 const G4double expxu = 82.;
00206 const G4double expxl = -expxu;
00207 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00208 G4ParticleDefinition *aProton = G4Proton::Proton();
00209 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00210 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
00211 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
00212
00213
00214
00215
00216 G4double n, anpn;
00217 GetNormalizationConstant( availableEnergy, n, anpn );
00218 G4double ran = G4UniformRand();
00219 G4double dum, excs = 0.0;
00220 if( targetParticle.GetDefinition() == aProton ) {
00221 counter = -1;
00222 for( np=0; np<numSec/3 && ran>=excs; ++np ) {
00223 for( nneg=std::max(0,np-2); nneg<=(np+1) && ran>=excs; ++nneg ) {
00224 for( nz=0; nz<numSec/3 && ran>=excs; ++nz ) {
00225 if( ++counter < numMul ) {
00226 nt = np+nneg+nz;
00227 if( nt>0 && nt<=numSec ) {
00228 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00229 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00230 if( std::fabs(dum) < 1.0 ) {
00231 if( test >= 1.0e-10 )excs += dum*test;
00232 } else {
00233 excs += dum*test;
00234 }
00235 }
00236 }
00237 }
00238 }
00239 }
00240 if( ran >= excs )
00241 {
00242 quasiElastic = true;
00243 return;
00244 }
00245 np--; nneg--; nz--;
00246 G4int ncht = std::max( 1, np-nneg );
00247 switch( ncht ) {
00248 case 1:
00249 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00250 incidentHasChanged = true;
00251 break;
00252 case 2:
00253 if( G4UniformRand() < 0.5 ) {
00254 if( G4UniformRand() < 0.5 ) {
00255 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00256 incidentHasChanged = true;
00257 }
00258 } else {
00259 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00260 incidentHasChanged = true;
00261 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00262 targetHasChanged = true;
00263 }
00264 break;
00265 case 3:
00266 if( G4UniformRand() < 0.5 ) {
00267 if( G4UniformRand() < 0.5 ) {
00268 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00269 incidentHasChanged = true;
00270 }
00271 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00272 targetHasChanged = true;
00273 } else {
00274 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00275 incidentHasChanged = true;
00276 }
00277 break;
00278 default:
00279 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00280 incidentHasChanged = true;
00281 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00282 targetHasChanged = true;
00283 break;
00284 }
00285 }
00286 else
00287 {
00288 counter = -1;
00289 for( np=0; np<numSec/3 && ran>=excs; ++np ) {
00290 for( nneg=std::max(0,np-1); nneg<=(np+2) && ran>=excs; ++nneg ) {
00291 for( nz=0; nz<numSec/3 && ran>=excs; ++nz ) {
00292 if( ++counter < numMul ) {
00293 nt = np+nneg+nz;
00294 if( nt>0 && nt<=numSec ) {
00295 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00296 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00297 if( std::fabs(dum) < 1.0 ) {
00298 if( test >= 1.0e-10 )excs += dum*test;
00299 } else {
00300 excs += dum*test;
00301 }
00302 }
00303 }
00304 }
00305 }
00306 }
00307 if( ran >= excs )
00308 {
00309 quasiElastic = true;
00310 return;
00311 }
00312 np--; nneg--; nz--;
00313 G4int ncht = std::max( 1, np-nneg+3 );
00314 switch( ncht ) {
00315 case 1:
00316 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00317 incidentHasChanged = true;
00318 targetParticle.SetDefinitionAndUpdateE( aProton );
00319 targetHasChanged = true;
00320 break;
00321 case 2:
00322 if( G4UniformRand() < 0.5 ) {
00323 if( G4UniformRand() < 0.5 ) {
00324 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00325 incidentHasChanged = true;
00326 }
00327 targetParticle.SetDefinitionAndUpdateE( aProton );
00328 targetHasChanged = true;
00329 } else {
00330 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00331 incidentHasChanged = true;
00332 }
00333 break;
00334 case 3:
00335 if( G4UniformRand() < 0.5 ) {
00336 if( G4UniformRand() < 0.5 ) {
00337 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
00338 incidentHasChanged = true;
00339 }
00340 } else {
00341 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00342 incidentHasChanged = true;
00343 targetParticle.SetDefinitionAndUpdateE( aProton );
00344 targetHasChanged = true;
00345 }
00346 break;
00347 default:
00348 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00349 incidentHasChanged = true;
00350 break;
00351 }
00352 }
00353
00354 SetUpPions(np, nneg, nz, vec, vecLen);
00355 return;
00356 }
00357
00358
00359