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