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