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 #include "G4LEAntiSigmaMinusInelastic.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "Randomize.hh"
00037
00038 void G4LEAntiSigmaMinusInelastic::ModelDescription(std::ostream& outFile) const
00039 {
00040 outFile << "G4LEAntiSigmaMinusInelastic is one of the Low Energy\n"
00041 << "Parameterized (LEP) models used to implement inelastic\n"
00042 << "antiSigma- scattering from nuclei. It is a re-engineered\n"
00043 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00044 << "initial collision products into backward- and forward-going\n"
00045 << "clusters which are then decayed into final state hadrons. The\n"
00046 << "model does not conserve energy on an event-by-event basis. It\n"
00047 << "may be applied to antiSigma- with initial energies between 0\n"
00048 << "and 25 GeV.\n";
00049 }
00050
00051
00052 G4HadFinalState*
00053 G4LEAntiSigmaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00054 G4Nucleus& targetNucleus)
00055 {
00056 const G4HadProjectile *originalIncident = &aTrack;
00057 if (originalIncident->GetKineticEnergy()<= 0.1*MeV) {
00058 theParticleChange.SetStatusChange(isAlive);
00059 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00060 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00061 return &theParticleChange;
00062 }
00063
00064
00065 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00066
00067 if (verboseLevel > 1) {
00068 const G4Material *targetMaterial = aTrack.GetMaterial();
00069 G4cout << "G4LEAntiSigmaMinusInelastic::ApplyYourself called" << G4endl;
00070 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00071 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00072 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00073 << G4endl;
00074 }
00075
00076
00077
00078 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00079 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00080 G4ReactionProduct modifiedOriginal;
00081 modifiedOriginal = *originalIncident;
00082
00083 G4double tkin = targetNucleus.Cinema( ek );
00084 ek += tkin;
00085 modifiedOriginal.SetKineticEnergy( ek*MeV );
00086 G4double et = ek + amas;
00087 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00088 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00089 if (pp > 0.0) {
00090 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00091 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00092 }
00093
00094
00095 tkin = targetNucleus.EvaporationEffects( ek );
00096 ek -= tkin;
00097 modifiedOriginal.SetKineticEnergy( ek*MeV );
00098 et = ek + amas;
00099 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00100 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00101 if (pp > 0.0) {
00102 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00103 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00104 }
00105 G4ReactionProduct currentParticle = modifiedOriginal;
00106 G4ReactionProduct targetParticle;
00107 targetParticle = *originalTarget;
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;
00118 const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
00119 if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
00120 (G4UniformRand() > anni) )
00121 Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00122 incidentHasChanged, targetHasChanged, quasiElastic);
00123
00124 CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00125 modifiedOriginal, targetNucleus, currentParticle,
00126 targetParticle, incidentHasChanged, targetHasChanged,
00127 quasiElastic);
00128
00129 SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00130
00131 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00132
00133 delete originalTarget;
00134 return &theParticleChange;
00135 }
00136
00137
00138 void G4LEAntiSigmaMinusInelastic::Cascade(
00139 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00140 G4int& vecLen,
00141 const G4HadProjectile *originalIncident,
00142 G4ReactionProduct ¤tParticle,
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 pOriginal = originalIncident->GetTotalMomentum()/MeV;
00161 const G4double targetMass = targetParticle.GetMass()/MeV;
00162 G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
00163 targetMass*targetMass +
00164 2.0*targetMass*etOriginal);
00165 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00166
00167 static G4bool first = true;
00168 const G4int numMul = 1200;
00169 const G4int numMulA = 400;
00170 const G4int numSec = 60;
00171 static G4double protmul[numMul], protnorm[numSec];
00172 static G4double neutmul[numMul], neutnorm[numSec];
00173 static G4double protmulA[numMulA], protnormA[numSec];
00174 static G4double neutmulA[numMulA], neutnormA[numSec];
00175
00176
00177 G4int counter, nt=0, npos=0, nneg=0, nzero=0;
00178 G4double test;
00179 const G4double c = 1.25;
00180 const G4double b[2] = { 0.7, 0.7 };
00181 if( first )
00182 {
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 {
00190 for( nneg=std::max(0,npos-2); nneg<=npos; ++nneg )
00191 {
00192 for( nzero=0; nzero<numSec/3; ++nzero )
00193 {
00194 if( ++counter < numMul )
00195 {
00196 nt = npos+nneg+nzero;
00197 if( nt>0 && nt<=numSec )
00198 {
00199 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00200 protnorm[nt-1] += protmul[counter];
00201 }
00202 }
00203 }
00204 }
00205 }
00206 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00207 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00208 counter = -1;
00209 for( npos=0; npos<numSec/3; ++npos )
00210 {
00211 for( nneg=std::max(0,npos-1); nneg<=(npos+1); ++nneg )
00212 {
00213 for( nzero=0; nzero<numSec/3; ++nzero )
00214 {
00215 if( ++counter < numMul )
00216 {
00217 nt = npos+nneg+nzero;
00218 if( nt>0 && nt<=numSec )
00219 {
00220 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00221 neutnorm[nt-1] += neutmul[counter];
00222 }
00223 }
00224 }
00225 }
00226 }
00227 for( i=0; i<numSec; ++i )
00228 {
00229 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00230 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00231 }
00232
00233
00234
00235 for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
00236 for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
00237 counter = -1;
00238 for( npos=2; npos<(numSec/3); ++npos )
00239 {
00240 nneg = npos-2;
00241 for( nzero=0; nzero<numSec/3; ++nzero )
00242 {
00243 if( ++counter < numMulA )
00244 {
00245 nt = npos+nneg+nzero;
00246 if( nt>1 && nt<=numSec )
00247 {
00248 protmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00249 protnormA[nt-1] += protmulA[counter];
00250 }
00251 }
00252 }
00253 }
00254 for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
00255 for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
00256 counter = -1;
00257 for( npos=1; npos<numSec/3; ++npos )
00258 {
00259 nneg = npos-1;
00260 for( nzero=0; nzero<numSec/3; ++nzero )
00261 {
00262 if( ++counter < numMulA )
00263 {
00264 nt = npos+nneg+nzero;
00265 if( nt>1 && nt<=numSec )
00266 {
00267 neutmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00268 neutnormA[nt-1] += neutmulA[counter];
00269 }
00270 }
00271 }
00272 }
00273 for( i=0; i<numSec; ++i )
00274 {
00275 if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
00276 if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
00277 }
00278 }
00279 const G4double expxu = 82.;
00280 const G4double expxl = -expxu;
00281 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00282 G4ParticleDefinition *aProton = G4Proton::Proton();
00283 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00284 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00285 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
00286 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
00287 G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
00288 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
00289 const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
00290 0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
00291 0.39,0.36,0.33,0.10,0.01};
00292 G4int iplab = G4int( pOriginal/GeV*10.0 );
00293 if( iplab > 9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
00294 if( iplab > 14 )iplab = G4int( pOriginal/GeV- 2.0 ) + 15;
00295 if( iplab > 23 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
00296 if( iplab > 24 )iplab = 24;
00297 if( G4UniformRand() > anhl[iplab] )
00298 {
00299 if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
00300 {
00301 quasiElastic = true;
00302 return;
00303 }
00304 G4double n, anpn;
00305 GetNormalizationConstant( availableEnergy, n, anpn );
00306 G4double ran = G4UniformRand();
00307 G4double dum, excs = 0.0;
00308 if( targetParticle.GetDefinition() == aProton )
00309 {
00310 counter = -1;
00311 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00312 {
00313 for( nneg=std::max(0,npos-2); nneg<=npos && ran>=excs; ++nneg )
00314 {
00315 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00316 {
00317 if( ++counter < numMul )
00318 {
00319 nt = npos+nneg+nzero;
00320 if( nt>0 && nt<=numSec )
00321 {
00322 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00323 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00324 if( std::fabs(dum) < 1.0 )
00325 {
00326 if( test >= 1.0e-10 )excs += dum*test;
00327 }
00328 else
00329 excs += dum*test;
00330 }
00331 }
00332 }
00333 }
00334 }
00335 if( ran >= excs )
00336 {
00337 quasiElastic = true;
00338 return;
00339 }
00340 npos--; nneg--; nzero--;
00341 G4int ncht = std::min( 3, std::max( 1, npos-nneg+1 ) );
00342 switch( ncht )
00343 {
00344 case 1:
00345 break;
00346 case 2:
00347 if( G4UniformRand() < 0.5 )
00348 {
00349 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00350 targetHasChanged = true;
00351 }
00352 else
00353 {
00354 if( G4UniformRand() < 0.5 )
00355 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00356 else
00357 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00358 incidentHasChanged = true;
00359 }
00360 break;
00361 case 3:
00362 if( G4UniformRand() < 0.5 )
00363 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00364 else
00365 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00366 incidentHasChanged = true;
00367 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00368 targetHasChanged = true;
00369 break;
00370 }
00371 }
00372 else
00373 {
00374 counter = -1;
00375 for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00376 {
00377 for( nneg=std::max(0,npos-1); nneg<=(npos+1) && ran>=excs; ++nneg )
00378 {
00379 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00380 {
00381 if( ++counter < numMul )
00382 {
00383 nt = npos+nneg+nzero;
00384 if( nt>0 && nt<=numSec )
00385 {
00386 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00387 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00388 if( std::fabs(dum) < 1.0 )
00389 {
00390 if( test >= 1.0e-10 )excs += dum*test;
00391 }
00392 else
00393 excs += dum*test;
00394 }
00395 }
00396 }
00397 }
00398 }
00399 if( ran >= excs )
00400 {
00401 quasiElastic = true;
00402 return;
00403 }
00404 npos--; nneg--; nzero--;
00405 G4int ncht = std::min( 3, std::max( 1, npos-nneg+2 ) );
00406 switch( ncht )
00407 {
00408 case 1:
00409 {
00410 targetParticle.SetDefinitionAndUpdateE( aProton );
00411 targetHasChanged = true;
00412 }
00413 break;
00414 case 2:
00415 if( G4UniformRand() < 0.5 )
00416 {
00417 if( G4UniformRand() < 0.5 )
00418 {
00419 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00420 incidentHasChanged = true;
00421 targetParticle.SetDefinitionAndUpdateE( aProton );
00422 targetHasChanged = true;
00423 }
00424 }
00425 else
00426 {
00427 if( G4UniformRand() < 0.5 )
00428 {
00429 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00430 incidentHasChanged = true;
00431 targetParticle.SetDefinitionAndUpdateE( aProton );
00432 targetHasChanged = true;
00433 }
00434 }
00435 break;
00436 case 3:
00437 if( G4UniformRand() < 0.5 )
00438 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00439 else
00440 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00441 incidentHasChanged = true;
00442 break;
00443 }
00444 }
00445 }
00446 else
00447 {
00448 if( centerofmassEnergy <= aPiPlus->GetPDGMass()/MeV+aKaonPlus->GetPDGMass()/MeV )
00449 {
00450 quasiElastic = true;
00451 return;
00452 }
00453 G4double n, anpn;
00454 GetNormalizationConstant( -centerofmassEnergy, n, anpn );
00455 G4double ran = G4UniformRand();
00456 G4double dum, excs = 0.0;
00457 if( targetParticle.GetDefinition() == aProton )
00458 {
00459 counter = -1;
00460 for( npos=2; npos<numSec/3 && ran>=excs; ++npos )
00461 {
00462 nneg=npos-2;
00463 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00464 {
00465 if( ++counter < numMulA )
00466 {
00467 nt = npos+nneg+nzero;
00468 if( nt>1 && nt<=numSec )
00469 {
00470 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00471 dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
00472 if( std::fabs(dum) < 1.0 )
00473 {
00474 if( test >= 1.0e-10 )excs += dum*test;
00475 }
00476 else
00477 excs += dum*test;
00478 }
00479 }
00480 }
00481 }
00482 if( ran >= excs )
00483 {
00484 quasiElastic = true;
00485 return;
00486 }
00487 npos--; nzero--;
00488 }
00489 else
00490 {
00491 counter = -1;
00492 for( npos=1; npos<numSec/3 && ran>=excs; ++npos )
00493 {
00494 nneg = npos-1;
00495 for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00496 {
00497 if( ++counter < numMulA )
00498 {
00499 nt = npos+nneg+nzero;
00500 if( nt>1 && nt<=numSec )
00501 {
00502 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00503 dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
00504 if( std::fabs(dum) < 1.0 )
00505 {
00506 if( test >= 1.0e-10 )excs += dum*test;
00507 }
00508 else
00509 excs += dum*test;
00510 }
00511 }
00512 }
00513 }
00514 if( ran >= excs )
00515 {
00516 quasiElastic = true;
00517 return;
00518 }
00519 npos--; nzero--;
00520 }
00521 if( nzero > 0 )
00522 {
00523 if( nneg > 0 )
00524 {
00525 if( G4UniformRand() < 0.5 )
00526 {
00527 vec.Initialize( 1 );
00528 G4ReactionProduct *p = new G4ReactionProduct;
00529 p->SetDefinition( aKaonMinus );
00530 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00531 vec.SetElement( vecLen++, p );
00532 --nneg;
00533 }
00534 else
00535 {
00536 vec.Initialize( 1 );
00537 G4ReactionProduct *p = new G4ReactionProduct;
00538 p->SetDefinition( aKaonZL );
00539 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00540 vec.SetElement( vecLen++, p );
00541 --nzero;
00542 }
00543 }
00544 else
00545 {
00546 vec.Initialize( 1 );
00547 G4ReactionProduct *p = new G4ReactionProduct;
00548 p->SetDefinition( aKaonZL );
00549 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00550 vec.SetElement( vecLen++, p );
00551 --nzero;
00552 }
00553 }
00554 else
00555 {
00556 if( nneg > 0 )
00557 {
00558 vec.Initialize( 1 );
00559 G4ReactionProduct *p = new G4ReactionProduct;
00560 p->SetDefinition( aKaonMinus );
00561 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00562 vec.SetElement( vecLen++, p );
00563 --nneg;
00564 }
00565 }
00566 currentParticle.SetMass( 0.0 );
00567 targetParticle.SetMass( 0.0 );
00568 }
00569 SetUpPions( npos, nneg, nzero, vec, vecLen );
00570 return;
00571 }
00572
00573
00574