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 "G4LEAntiOmegaMinusInelastic.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "Randomize.hh"
00040
00041 void G4LEAntiOmegaMinusInelastic::ModelDescription(std::ostream& outFile) const
00042 {
00043 outFile << "G4LEAntiOmegaMinusInelastic is one of the Low Energy\n"
00044 << "Parameterized (LEP) models used to implement inelastic\n"
00045 << "antiOmega- scattering from nuclei. It is a re-engineered\n"
00046 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00047 << "initial collision products into backward- and forward-going\n"
00048 << "clusters which are then decayed into final state hadrons. The\n"
00049 << "model does not conserve energy on an event-by-event basis. It\n"
00050 << "may be applied to antiOmega- with initial energies between 0\n"
00051 << "and 25 GeV.\n";
00052 }
00053
00054 G4HadFinalState*
00055 G4LEAntiOmegaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00056 G4Nucleus& targetNucleus)
00057 {
00058 const G4HadProjectile* originalIncident = &aTrack;
00059 if (originalIncident->GetKineticEnergy()<= 0.1*MeV) {
00060 theParticleChange.SetStatusChange(isAlive);
00061 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00062 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00063 return &theParticleChange;
00064 }
00065
00066
00067 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00068
00069 if (verboseLevel > 1) {
00070 const G4Material *targetMaterial = aTrack.GetMaterial();
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
00121 if ((currentParticle.GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
00122 Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00123 incidentHasChanged, targetHasChanged, quasiElastic);
00124
00125 CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00126 modifiedOriginal, targetNucleus, currentParticle,
00127 targetParticle, incidentHasChanged, targetHasChanged,
00128 quasiElastic);
00129
00130 SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00131
00132 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00133
00134 delete originalTarget;
00135 return &theParticleChange;
00136 }
00137
00138 void G4LEAntiOmegaMinusInelastic::Cascade(
00139 G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec,
00140 G4int& vecLen,
00141 const G4HadProjectile* originalIncident,
00142 G4ReactionProduct& currentParticle,
00143 G4ReactionProduct& targetParticle,
00144 G4bool& incidentHasChanged,
00145 G4bool& targetHasChanged,
00146 G4bool& quasiElastic)
00147 {
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00159 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00160 const G4double targetMass = targetParticle.GetMass()/MeV;
00161 G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
00162 targetMass*targetMass +
00163 2.0*targetMass*etOriginal);
00164 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00165 if (availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV) {
00166
00167 quasiElastic = true;
00168 return;
00169 }
00170 static G4bool first = true;
00171 const G4int numMul = 1200;
00172 const G4int numSec = 60;
00173 static G4double protmul[numMul], protnorm[numSec];
00174 static G4double neutmul[numMul], neutnorm[numSec];
00175
00176
00177 G4int counter, nt=0;
00178 G4int npos = 0, nneg = 0, nzero = 0;
00179 G4double test;
00180 const G4double c = 1.25;
00181 const G4double b[] = { 0.7, 0.7 };
00182 if (first) {
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 for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
00190 for (nzero = 0; nzero < numSec/3; ++nzero) {
00191 if (++counter < numMul) {
00192 nt = npos+nneg+nzero;
00193 if (nt > 0 && nt <= numSec) {
00194 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00195 protnorm[nt-1] += protmul[counter];
00196 }
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 for (nneg = npos; nneg <= (npos+2); ++nneg) {
00207 for (nzero = 0; nzero < numSec/3; ++nzero) {
00208 if (++counter < numMul) {
00209 nt = npos+nneg+nzero;
00210 if ( nt>0 && nt<=numSec ) {
00211 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00212 neutnorm[nt-1] += neutmul[counter];
00213 }
00214 }
00215 }
00216 }
00217 }
00218 for (i = 0; i < numSec; ++i) {
00219 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00220 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00221 }
00222 }
00223
00224 const G4double expxu = 82.;
00225 const G4double expxl = -expxu;
00226 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00227 G4ParticleDefinition *aProton = G4Proton::Proton();
00228 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00229 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00230 G4ParticleDefinition *aXiZero = G4XiZero::XiZero();
00231 G4double n, anpn;
00232 GetNormalizationConstant( availableEnergy, n, anpn );
00233 G4double ran = G4UniformRand();
00234 G4double dum, excs = 0.0;
00235 G4int nvefix = 0;
00236 if (targetParticle.GetDefinition() == aProton) {
00237 counter = -1;
00238 for (npos = 0; npos < numSec/3 && ran>=excs; ++npos) {
00239 for (nneg = std::max(0,npos-1); nneg <= (npos+1) && ran>=excs; ++nneg) {
00240 for (nzero = 0; nzero < numSec/3 && ran>=excs; ++nzero) {
00241 if ( ++counter < numMul ) {
00242 nt = npos+nneg+nzero;
00243 if (nt > 0 && nt <= numSec) {
00244 test = std::exp(std::min(expxu,
00245 std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00246 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00247 if (std::fabs(dum) < 1.0) {
00248 if( test >= 1.0e-10 )excs += dum*test;
00249 }
00250 else
00251 excs += dum*test;
00252 }
00253 }
00254 }
00255 }
00256 }
00257 if (ran >= excs) {
00258
00259 quasiElastic = true;
00260 return;
00261 }
00262 npos--; nneg--; nzero--;
00263
00264
00265
00266
00267
00268
00269 if (npos < nneg) {
00270 if (npos+1 == nneg) {
00271 currentParticle.SetDefinitionAndUpdateE(aXiZero);
00272 incidentHasChanged = true;
00273 nvefix = 1;
00274 } else {
00275
00276 currentParticle.SetDefinitionAndUpdateE(aSigmaPlus);
00277 incidentHasChanged = true;
00278 nvefix = 2;
00279 }
00280 } else if (npos > nneg) {
00281 targetParticle.SetDefinitionAndUpdateE(aNeutron);
00282 targetHasChanged = true;
00283 }
00284 } else {
00285
00286 counter = -1;
00287 for (npos = 0; npos < numSec/3 && ran >= excs; ++npos) {
00288 for (nneg = npos; nneg <= (npos+2) && ran>=excs; ++nneg) {
00289 for (nzero = 0; nzero < numSec/3 && ran>=excs; ++nzero) {
00290 if (++counter < numMul) {
00291 nt = npos+nneg+nzero;
00292 if (nt > 0 && nt <= numSec) {
00293 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00294 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00295 if( std::fabs(dum) < 1.0 )
00296 {
00297 if( test >= 1.0e-10 )excs += dum*test;
00298 }
00299 else
00300 excs += dum*test;
00301 }
00302 }
00303 }
00304 }
00305 }
00306 if (ran >= excs) {
00307
00308 quasiElastic = true;
00309 return;
00310 }
00311 npos--; nneg--; nzero--;
00312 if (npos+1 < nneg) {
00313 if( npos+2 == nneg) {
00314 currentParticle.SetDefinitionAndUpdateE( aXiZero );
00315 incidentHasChanged = true;
00316 nvefix = 1;
00317 } else {
00318
00319 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00320 incidentHasChanged = true;
00321 nvefix = 2;
00322 }
00323 targetParticle.SetDefinitionAndUpdateE( aProton );
00324 targetHasChanged = true;
00325 } else if (npos+1 == nneg) {
00326 targetParticle.SetDefinitionAndUpdateE( aProton );
00327 targetHasChanged = true;
00328 }
00329 }
00330
00331 SetUpPions(npos, nneg, nzero, vec, vecLen);
00332 for (G4int i = 0; i < vecLen && nvefix > 0; ++i) {
00333 if (vec[i]->GetDefinition() == G4PionMinus::PionMinus() ) {
00334
00335
00336 if (nvefix >= 1) vec[i]->SetDefinitionAndUpdateE(aKaonMinus);
00337 --nvefix;
00338 }
00339 }
00340 return;
00341 }
00342
00343
00344