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 "G4LEAntiXiZeroInelastic.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "Randomize.hh"
00040
00041 void G4LEAntiXiZeroInelastic::ModelDescription(std::ostream& outFile) const
00042 {
00043 outFile << "G4LEAntiXiZeroInelastic is one of the Low Energy Parameterized\n"
00044 << "(LEP) models used to implement inelastic antiXi0 scattering\n"
00045 << "from nuclei. It is a re-engineered version of the GHEISHA\n"
00046 << "code of H. Fesefeldt. It divides the initial collision\n"
00047 << "products into backward- and forward-going clusters which are\n"
00048 << "then decayed into final state hadrons. The model does not\n"
00049 << "conserve energy on an event-by-event basis. It may be\n"
00050 << "applied to antiXi0 with initial energies between 0 and 25\n"
00051 << "GeV.\n";
00052 }
00053
00054 G4HadFinalState*
00055 G4LEAntiXiZeroInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00056 G4Nucleus& targetNucleus)
00057 {
00058 const G4HadProjectile *originalIncident = &aTrack;
00059
00060
00061 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00062
00063 if (verboseLevel > 1) {
00064 const G4Material *targetMaterial = aTrack.GetMaterial();
00065 G4cout << "G4LEAntiXiZeroInelastic::ApplyYourself called" << G4endl;
00066 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00067 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00068 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00069 << G4endl;
00070 }
00071
00072
00073
00074 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00075 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00076 G4ReactionProduct modifiedOriginal;
00077 modifiedOriginal = *originalIncident;
00078
00079 G4double tkin = targetNucleus.Cinema( ek );
00080 ek += tkin;
00081 modifiedOriginal.SetKineticEnergy( ek*MeV );
00082 G4double et = ek + amas;
00083 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00084 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00085 if (pp > 0.0) {
00086 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00087 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00088 }
00089
00090
00091 tkin = targetNucleus.EvaporationEffects( ek );
00092 ek -= tkin;
00093 modifiedOriginal.SetKineticEnergy( ek*MeV );
00094 et = ek + amas;
00095 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00096 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00097 if (pp > 0.0) {
00098 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00099 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00100 }
00101 G4ReactionProduct currentParticle = modifiedOriginal;
00102 G4ReactionProduct targetParticle;
00103 targetParticle = *originalTarget;
00104 currentParticle.SetSide(1);
00105 targetParticle.SetSide(-1);
00106 G4bool incidentHasChanged = false;
00107 G4bool targetHasChanged = false;
00108 G4bool quasiElastic = false;
00109 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00110 G4int vecLen = 0;
00111 vec.Initialize(0);
00112
00113 const G4double cutOff = 0.1;
00114 const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
00115 if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
00116 (G4UniformRand() > anni) )
00117 Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00118 incidentHasChanged, targetHasChanged, quasiElastic);
00119
00120 CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00121 modifiedOriginal, targetNucleus, currentParticle,
00122 targetParticle, incidentHasChanged, targetHasChanged,
00123 quasiElastic);
00124
00125 SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00126
00127 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00128
00129 delete originalTarget;
00130 return &theParticleChange;
00131 }
00132
00133
00134 void G4LEAntiXiZeroInelastic::Cascade(
00135 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00136 G4int& vecLen,
00137 const G4HadProjectile *originalIncident,
00138 G4ReactionProduct ¤tParticle,
00139 G4ReactionProduct &targetParticle,
00140 G4bool &incidentHasChanged,
00141 G4bool &targetHasChanged,
00142 G4bool &quasiElastic)
00143 {
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
00240 G4double n, anpn;
00241 GetNormalizationConstant( availableEnergy, n, anpn );
00242 G4double ran = G4UniformRand();
00243 G4double dum, excs = 0.0;
00244 if( targetParticle.GetDefinition() == aProton )
00245 {
00246 counter = -1;
00247 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00248 {
00249 for( nneg=std::max(0,npos-2); nneg<=(npos+1) && ran>=excs; ++nneg )
00250 {
00251 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00252 {
00253 if( ++counter < numMul )
00254 {
00255 nt = npos+nneg+nzero;
00256 if( nt>0 && nt<=numSec )
00257 {
00258 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00259 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00260 if( std::fabs(dum) < 1.0 )
00261 {
00262 if( test >= 1.0e-10 )excs += dum*test;
00263 }
00264 else
00265 excs += dum*test;
00266 }
00267 }
00268 }
00269 }
00270 }
00271 if( ran >= excs )
00272 {
00273 quasiElastic = true;
00274 return;
00275 }
00276 npos--; nneg--; nzero--;
00277
00278
00279
00280
00281
00282
00283 if( npos < nneg+1 )
00284 {
00285 if( npos != nneg )
00286 {
00287 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00288 incidentHasChanged = true;
00289
00290
00291
00292 vec.Initialize( 1 );
00293 G4ReactionProduct *p = new G4ReactionProduct;
00294 p->SetDefinition( aKaonMinus );
00295 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00296 vec.SetElement( vecLen++, p );
00297 --nneg;
00298 }
00299 }
00300 else if( npos == nneg+1 )
00301 {
00302 if( G4UniformRand() < 0.5 )
00303 {
00304 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00305 targetHasChanged = true;
00306 }
00307 else
00308 {
00309 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00310 incidentHasChanged = true;
00311 }
00312 }
00313 else
00314 {
00315 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00316 incidentHasChanged = true;
00317 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00318 targetHasChanged = true;
00319 }
00320 }
00321 else
00322 {
00323 counter = -1;
00324 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00325 {
00326 for( nneg=std::max(0,npos-1); nneg<=(npos+2) && ran>=excs; ++nneg )
00327 {
00328 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00329 {
00330 if( ++counter < numMul )
00331 {
00332 nt = npos+nneg+nzero;
00333 if( nt>0 && nt<=numSec )
00334 {
00335 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00336 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00337 if( std::fabs(dum) < 1.0 )
00338 {
00339 if( test >= 1.0e-10 )excs += dum*test;
00340 }
00341 else
00342 excs += dum*test;
00343 }
00344 }
00345 }
00346 }
00347 }
00348 if( ran >= excs )
00349 {
00350 quasiElastic = true;
00351 return;
00352 }
00353 npos--; nneg--; nzero--;
00354 if( npos < nneg )
00355 {
00356 if( npos+1 == nneg )
00357 {
00358 targetParticle.SetDefinitionAndUpdateE( aProton );
00359 targetHasChanged = true;
00360 }
00361 else
00362 {
00363 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00364 incidentHasChanged = true;
00365 targetParticle.SetDefinitionAndUpdateE( aProton );
00366 targetHasChanged = true;
00367
00368
00369
00370 vec.Initialize( 1 );
00371 G4ReactionProduct *p = new G4ReactionProduct;
00372 p->SetDefinition( aKaonMinus );
00373 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00374 vec.SetElement( vecLen++, p );
00375 --nneg;
00376 }
00377 }
00378 else if( npos == nneg )
00379 {
00380 if( G4UniformRand() >= 0.5 )
00381 {
00382 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00383 incidentHasChanged = true;
00384 targetParticle.SetDefinitionAndUpdateE( aProton );
00385 targetHasChanged = true;
00386 }
00387 }
00388 else
00389 {
00390 currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00391 incidentHasChanged = true;
00392 }
00393 }
00394 SetUpPions( npos, nneg, nzero, vec, vecLen );
00395 return;
00396 }
00397
00398
00399