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