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