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