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