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