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 "G4LEXiZeroInelastic.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036
00037 void G4LEXiZeroInelastic::ModelDescription(std::ostream& outFile) const
00038 {
00039 outFile << "G4LEXiZeroInelastic is one of the Low Energy Parameterized\n"
00040 << "(LEP) models used to implement inelastic X0 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 X0 with initial energies between 0 and 25\n"
00047 << "GeV.\n";
00048 }
00049
00050
00051 G4HadFinalState*
00052 G4LEXiZeroInelastic::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 << "G4LEXiZeroInelastic::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 G4LEXiZeroInelastic::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-2); 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=std::max(0,npos-1); 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 *aXiMinus = G4XiMinus::XiMinus();
00237
00238
00239 G4double n, anpn;
00240 GetNormalizationConstant(availableEnergy, n, anpn);
00241 G4double ran = G4UniformRand();
00242 G4double dum, excs = 0.0;
00243 if( targetParticle.GetDefinition() == aProton )
00244 {
00245 counter = -1;
00246 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00247 {
00248 for( nneg=std::max(0,npos-2); nneg<=(npos+1) && ran>=excs; ++nneg )
00249 {
00250 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00251 {
00252 if( ++counter < numMul )
00253 {
00254 nt = npos+nneg+nzero;
00255 if( nt>0 && nt<=numSec )
00256 {
00257 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00258 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00259 if( std::fabs(dum) < 1.0 )
00260 {
00261 if( test >= 1.0e-10 )excs += dum*test;
00262 }
00263 else
00264 excs += dum*test;
00265 }
00266 }
00267 }
00268 }
00269 }
00270 if( ran >= excs )
00271 {
00272 quasiElastic = true;
00273 return;
00274 }
00275 npos--; nneg--; nzero--;
00276
00277
00278
00279
00280
00281
00282 if( npos < nneg+1 )
00283 {
00284 if( npos != nneg )
00285 {
00286 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00287 incidentHasChanged = true;
00288
00289
00290
00291 vec.Initialize( 1 );
00292 G4ReactionProduct *p = new G4ReactionProduct;
00293 p->SetDefinition( aKaonMinus );
00294 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00295 vec.SetElement( vecLen++, p );
00296 --nneg;
00297 }
00298 }
00299 else if( npos == nneg+1 )
00300 {
00301 if( G4UniformRand() < 0.5 )
00302 {
00303 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00304 targetHasChanged = true;
00305 }
00306 else
00307 {
00308 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00309 incidentHasChanged = true;
00310 }
00311 }
00312 else
00313 {
00314 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00315 incidentHasChanged = true;
00316 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00317 targetHasChanged = true;
00318 }
00319 }
00320 else
00321 {
00322 counter = -1;
00323 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00324 {
00325 for( nneg=std::max(0,npos-1); nneg<=(npos+2) && ran>=excs; ++nneg )
00326 {
00327 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00328 {
00329 if( ++counter < numMul )
00330 {
00331 nt = npos+nneg+nzero;
00332 if( nt>0 && nt<=numSec )
00333 {
00334 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00335 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00336 if( std::fabs(dum) < 1.0 )
00337 {
00338 if( test >= 1.0e-10 )excs += dum*test;
00339 }
00340 else
00341 excs += dum*test;
00342 }
00343 }
00344 }
00345 }
00346 }
00347 if( ran >= excs )
00348 {
00349 quasiElastic = true;
00350 return;
00351 }
00352 npos--; nneg--; nzero--;
00353 if( npos < nneg )
00354 {
00355 if( npos+1 == nneg )
00356 {
00357 targetParticle.SetDefinitionAndUpdateE( aProton );
00358 targetHasChanged = true;
00359 }
00360 else
00361 {
00362 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00363 incidentHasChanged = true;
00364 targetParticle.SetDefinitionAndUpdateE( aProton );
00365 targetHasChanged = true;
00366
00367
00368
00369 vec.Initialize( 1 );
00370 G4ReactionProduct *p = new G4ReactionProduct;
00371 p->SetDefinition( aKaonMinus );
00372 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00373 vec.SetElement( vecLen++, p );
00374 --nneg;
00375 }
00376 }
00377 else if( npos == nneg )
00378 {
00379 if( G4UniformRand() >= 0.5 )
00380 {
00381 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00382 incidentHasChanged = true;
00383 targetParticle.SetDefinitionAndUpdateE( aProton );
00384 targetHasChanged = true;
00385 }
00386 }
00387 else
00388 {
00389 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00390 incidentHasChanged = true;
00391 }
00392 }
00393 SetUpPions(npos, nneg, nzero, vec, vecLen);
00394 return;
00395 }
00396
00397
00398