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
00037 #include "G4HEAntiSigmaMinusInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042
00043 void G4HEAntiSigmaMinusInelastic::ModelDescription(std::ostream& outFile) const
00044 {
00045 outFile << "G4HEAntiSigmaMinusInelastic is one of the High Energy\n"
00046 << "Parameterized (HEP) models used to implement inelastic\n"
00047 << "anti-Sigma- scattering from nuclei. It is a re-engineered\n"
00048 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00049 << "initial collision products into backward- and forward-going\n"
00050 << "clusters which are then decayed into final state hadrons.\n"
00051 << "The model does not conserve energy on an event-by-event\n"
00052 << "basis. It may be applied to anti-Sigma- with initial\n"
00053 << "energies above 20 GeV.\n";
00054 }
00055
00056
00057 G4HadFinalState*
00058 G4HEAntiSigmaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00059 G4Nucleus& targetNucleus)
00060 {
00061 G4HEVector* pv = new G4HEVector[MAXPART];
00062 const G4HadProjectile* aParticle = &aTrack;
00063 const G4double atomicWeight = targetNucleus.GetA_asInt();
00064 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00065 G4HEVector incidentParticle(aParticle);
00066
00067 G4int incidentCode = incidentParticle.getCode();
00068 G4double incidentMass = incidentParticle.getMass();
00069 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00070
00071
00072
00073
00074 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00075
00076 if (incidentKineticEnergy < 1.)
00077 G4cout << "GHEAntiSigmaMinusInelastic: incident energy < 1 GeV" << G4endl;
00078
00079 if (verboseLevel > 1) {
00080 G4cout << "G4HEAntiSigmaMinusInelastic::ApplyYourself" << G4endl;
00081 G4cout << "incident particle " << incidentParticle.getName()
00082 << "mass " << incidentMass
00083 << "kinetic energy " << incidentKineticEnergy
00084 << G4endl;
00085 G4cout << "target material with (A,Z) = ("
00086 << atomicWeight << "," << atomicNumber << ")" << G4endl;
00087 }
00088
00089 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
00090 atomicWeight, atomicNumber);
00091 if (verboseLevel > 1)
00092 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00093
00094 incidentKineticEnergy -= inelasticity;
00095
00096 G4double excitationEnergyGNP = 0.;
00097 G4double excitationEnergyDTA = 0.;
00098
00099 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00100 atomicWeight, atomicNumber,
00101 excitationEnergyGNP,
00102 excitationEnergyDTA);
00103 if (verboseLevel > 1)
00104 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
00105 << excitationEnergyDTA << G4endl;
00106
00107 incidentKineticEnergy -= excitation;
00108 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00109
00110
00111
00112
00113 G4HEVector targetParticle;
00114 if (G4UniformRand() < atomicNumber/atomicWeight) {
00115 targetParticle.setDefinition("Proton");
00116 } else {
00117 targetParticle.setDefinition("Neutron");
00118 }
00119
00120 G4double targetMass = targetParticle.getMass();
00121 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00122 + targetMass*targetMass
00123 + 2.0*targetMass*incidentTotalEnergy);
00124 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00125
00126 G4bool inElastic = true;
00127 vecLength = 0;
00128
00129 if (verboseLevel > 1)
00130 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00131 << incidentCode << G4endl;
00132
00133 G4bool successful = false;
00134
00135 FirstIntInCasAntiSigmaMinus(inElastic, availableEnergy, pv, vecLength,
00136 incidentParticle, targetParticle, atomicWeight);
00137
00138 if (verboseLevel > 1)
00139 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00140
00141 if ((vecLength > 0) && (availableEnergy > 1.))
00142 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00143 pv, vecLength,
00144 incidentParticle, targetParticle);
00145 HighEnergyCascading(successful, pv, vecLength,
00146 excitationEnergyGNP, excitationEnergyDTA,
00147 incidentParticle, targetParticle,
00148 atomicWeight, atomicNumber);
00149 if (!successful)
00150 HighEnergyClusterProduction(successful, pv, vecLength,
00151 excitationEnergyGNP, excitationEnergyDTA,
00152 incidentParticle, targetParticle,
00153 atomicWeight, atomicNumber);
00154 if (!successful)
00155 MediumEnergyCascading(successful, pv, vecLength,
00156 excitationEnergyGNP, excitationEnergyDTA,
00157 incidentParticle, targetParticle,
00158 atomicWeight, atomicNumber);
00159
00160 if (!successful)
00161 MediumEnergyClusterProduction(successful, pv, vecLength,
00162 excitationEnergyGNP, excitationEnergyDTA,
00163 incidentParticle, targetParticle,
00164 atomicWeight, atomicNumber);
00165 if (!successful)
00166 QuasiElasticScattering(successful, pv, vecLength,
00167 excitationEnergyGNP, excitationEnergyDTA,
00168 incidentParticle, targetParticle,
00169 atomicWeight, atomicNumber);
00170 if (!successful)
00171 ElasticScattering(successful, pv, vecLength,
00172 incidentParticle,
00173 atomicWeight, atomicNumber);
00174
00175 if (!successful)
00176 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00177 << G4endl;
00178 FillParticleChange(pv, vecLength);
00179 delete [] pv;
00180 theParticleChange.SetStatusChange(stopAndKill);
00181 return &theParticleChange;
00182 }
00183
00184
00185 void
00186 G4HEAntiSigmaMinusInelastic::FirstIntInCasAntiSigmaMinus(G4bool& inElastic,
00187 const G4double availableEnergy,
00188 G4HEVector pv[],
00189 G4int& vecLen,
00190 const G4HEVector& incidentParticle,
00191 const G4HEVector& targetParticle,
00192 const G4double atomicWeight)
00193
00194
00195
00196
00197
00198
00199
00200
00201 {
00202 static const G4double expxu = 82;
00203 static const G4double expxl = -expxu;
00204
00205 static const G4double protb = 0.7;
00206 static const G4double neutb = 0.7;
00207 static const G4double c = 1.25;
00208
00209 static const G4int numMul = 1200;
00210 static const G4int numMulAn = 400;
00211 static const G4int numSec = 60;
00212
00213 G4int neutronCode = Neutron.getCode();
00214 G4int protonCode = Proton.getCode();
00215
00216 G4int targetCode = targetParticle.getCode();
00217 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00218
00219 static G4bool first = true;
00220 static G4double protmul[numMul], protnorm[numSec];
00221 static G4double protmulAn[numMulAn],protnormAn[numSec];
00222 static G4double neutmul[numMul], neutnorm[numSec];
00223 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00224
00225
00226
00227
00228 G4int i, counter, nt, npos, nneg, nzero;
00229
00230 if( first )
00231 {
00232 first = false;
00233 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
00234 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
00235 counter = -1;
00236 for( npos=0; npos<(numSec/3); npos++ )
00237 {
00238 for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ )
00239 {
00240 for( nzero=0; nzero<numSec/3; nzero++ )
00241 {
00242 if( ++counter < numMul )
00243 {
00244 nt = npos+nneg+nzero;
00245 if( (nt>0) && (nt<=numSec) )
00246 {
00247 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00248 protnorm[nt-1] += protmul[counter];
00249 }
00250 }
00251 }
00252 }
00253 }
00254 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00255 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00256 counter = -1;
00257 for( npos=0; npos<numSec/3; npos++ )
00258 {
00259 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
00260 {
00261 for( nzero=0; nzero<numSec/3; nzero++ )
00262 {
00263 if( ++counter < numMul )
00264 {
00265 nt = npos+nneg+nzero;
00266 if( (nt>0) && (nt<=numSec) )
00267 {
00268 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00269 neutnorm[nt-1] += neutmul[counter];
00270 }
00271 }
00272 }
00273 }
00274 }
00275 for( i=0; i<numSec; i++ )
00276 {
00277 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00278 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00279 }
00280
00281 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
00282 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
00283 counter = -1;
00284 for( npos=1; npos<(numSec/3); npos++ )
00285 {
00286 nneg = std::max(0,npos-2);
00287 for( nzero=0; nzero<numSec/3; nzero++ )
00288 {
00289 if( ++counter < numMulAn )
00290 {
00291 nt = npos+nneg+nzero;
00292 if( (nt>1) && (nt<=numSec) )
00293 {
00294 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00295 protnormAn[nt-1] += protmulAn[counter];
00296 }
00297 }
00298 }
00299 }
00300 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
00301 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
00302 counter = -1;
00303 for( npos=0; npos<numSec/3; npos++ )
00304 {
00305 nneg = npos-1;
00306 for( nzero=0; nzero<numSec/3; nzero++ )
00307 {
00308 if( ++counter < numMulAn )
00309 {
00310 nt = npos+nneg+nzero;
00311 if( (nt>1) && (nt<=numSec) )
00312 {
00313 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00314 neutnormAn[nt-1] += neutmulAn[counter];
00315 }
00316 }
00317 }
00318 }
00319 for( i=0; i<numSec; i++ )
00320 {
00321 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00322 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00323 }
00324 }
00325
00326
00327
00328
00329 pv[0] = incidentParticle;
00330 pv[1] = targetParticle;
00331 vecLen = 2;
00332
00333 if( !inElastic )
00334 {
00335 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00336
00337 G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5 ));
00338 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00339 {
00340 G4double ran = G4UniformRand();
00341
00342 if ( targetCode == neutronCode)
00343 {
00344 if(ran < 0.2)
00345 {
00346 pv[0] = AntiSigmaZero;
00347 pv[1] = Proton;
00348 }
00349 else if (ran < 0.4)
00350 {
00351 pv[0] = AntiLambda;
00352 pv[1] = Proton;
00353 }
00354 else if (ran < 0.6)
00355 {
00356 pv[0] = Proton;
00357 pv[1] = AntiLambda;
00358 }
00359 else if (ran < 0.8)
00360 {
00361 pv[0] = Proton;
00362 pv[1] = AntiSigmaZero;
00363 }
00364 else
00365 {
00366 pv[0] = Neutron;
00367 pv[1] = AntiSigmaMinus;
00368 }
00369 }
00370 else
00371 {
00372 pv[0] = Proton;
00373 pv[1] = AntiSigmaMinus;
00374 }
00375 }
00376 return;
00377 }
00378 else if (availableEnergy <= PionPlus.getMass())
00379 return;
00380
00381
00382
00383 npos = 0; nneg = 0; nzero = 0;
00384 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
00385 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
00386 0.39, 0.36, 0.33, 0.10, 0.01};
00387 G4int iplab = G4int( incidentTotalMomentum*10.);
00388 if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. );
00389 if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. );
00390 if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.);
00391 iplab = std::min(24, iplab);
00392
00393 if (G4UniformRand() > anhl[iplab]) {
00394
00395
00396 G4double aleab = std::log(availableEnergy);
00397 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00398 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00399
00400
00401
00402 G4double test, dum, anpn = 0.0;
00403
00404 for (nt=1; nt<=numSec; nt++) {
00405 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00406 dum = pi*nt/(2.0*n*n);
00407 if (std::fabs(dum) < 1.0) {
00408 if( test >= 1.0e-10 )anpn += dum*test;
00409 } else {
00410 anpn += dum*test;
00411 }
00412 }
00413
00414 G4double ran = G4UniformRand();
00415 G4double excs = 0.0;
00416 if (targetCode == protonCode) {
00417 counter = -1;
00418 for (npos = 0; npos < numSec/3; npos++) {
00419 for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) {
00420 for (nzero = 0; nzero < numSec/3; nzero++) {
00421 if (++counter < numMul) {
00422 nt = npos+nneg+nzero;
00423 if ((nt > 0) && (nt <= numSec) ) {
00424 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00425 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00426 if (std::fabs(dum) < 1.0) {
00427 if( test >= 1.0e-10 )excs += dum*test;
00428 } else {
00429 excs += dum*test;
00430 }
00431
00432 if (ran < excs) goto outOfLoop;
00433 }
00434 }
00435 }
00436 }
00437 }
00438
00439 inElastic = false;
00440 return;
00441 } else {
00442 counter = -1;
00443 for( npos=0; npos<numSec/3; npos++ )
00444 {
00445 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
00446 {
00447 for( nzero=0; nzero<numSec/3; nzero++ )
00448 {
00449 if( ++counter < numMul )
00450 {
00451 nt = npos+nneg+nzero;
00452 if ( (nt>0) && (nt<=numSec) ) {
00453 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00454 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00455 if (std::fabs(dum) < 1.0) {
00456 if( test >= 1.0e-10 )excs += dum*test;
00457 } else {
00458 excs += dum*test;
00459 }
00460
00461 if (ran < excs) goto outOfLoop;
00462 }
00463 }
00464 }
00465 }
00466 }
00467
00468 inElastic = false;
00469 return;
00470 }
00471
00472 outOfLoop:
00473
00474 ran = G4UniformRand();
00475
00476 if( targetCode == protonCode)
00477 {
00478 if( npos == nneg)
00479 {
00480 }
00481 else if (npos == (nneg+1))
00482 {
00483 if( ran < 0.50)
00484 {
00485 pv[1] = Neutron;
00486 }
00487 else if (ran < 0.75)
00488 {
00489 pv[0] = AntiSigmaZero;
00490 }
00491 else
00492 {
00493 pv[0] = AntiLambda;
00494 }
00495 }
00496 else
00497 {
00498 if (ran < 0.5)
00499 {
00500 pv[0] = AntiSigmaZero;
00501 pv[1] = Neutron;
00502 }
00503 else
00504 {
00505 pv[0] = AntiLambda;
00506 pv[1] = Neutron;
00507 }
00508 }
00509 }
00510 else
00511 {
00512 if( npos == nneg)
00513 {
00514 if (ran < 0.5)
00515 {
00516 }
00517 else if(ran < 0.75)
00518 {
00519 pv[0] = AntiSigmaZero;
00520 pv[1] = Proton;
00521 }
00522 else
00523 {
00524 pv[0] = AntiLambda;
00525 pv[1] = Proton;
00526 }
00527 }
00528 else if ( npos == (nneg+1))
00529 {
00530 if (ran < 0.5)
00531 {
00532 pv[0] = AntiSigmaZero;
00533 }
00534 else
00535 {
00536 pv[0] = AntiLambda;
00537 }
00538 }
00539 else
00540 {
00541 pv[1] = Proton;
00542 }
00543 }
00544
00545 }
00546 else
00547 {
00548 if ( availableEnergy > 2. * PionPlus.getMass() )
00549 {
00550
00551 G4double aleab = std::log(availableEnergy);
00552 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00553 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00554
00555
00556
00557 G4double test, dum, anpn = 0.0;
00558
00559 for (nt=2; nt<=numSec; nt++) {
00560 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00561 dum = pi*nt/(2.0*n*n);
00562 if (std::fabs(dum) < 1.0) {
00563 if( test >= 1.0e-10 )anpn += dum*test;
00564 } else {
00565 anpn += dum*test;
00566 }
00567 }
00568
00569 G4double ran = G4UniformRand();
00570 G4double excs = 0.0;
00571 if( targetCode == protonCode )
00572 {
00573 counter = -1;
00574 for( npos=2; npos<numSec/3; npos++ )
00575 {
00576 nneg = npos-2;
00577 for( nzero=0; nzero<numSec/3; nzero++ )
00578 {
00579 if( ++counter < numMulAn )
00580 {
00581 nt = npos+nneg+nzero;
00582 if ( (nt>1) && (nt<=numSec) ) {
00583 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00584 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00585 if (std::fabs(dum) < 1.0) {
00586 if( test >= 1.0e-10 )excs += dum*test;
00587 } else {
00588 excs += dum*test;
00589 }
00590
00591 if (ran < excs) goto outOfLoopAn;
00592 }
00593 }
00594 }
00595 }
00596
00597 inElastic = false;
00598 return;
00599 }
00600 else
00601 {
00602 counter = -1;
00603 for( npos=1; npos<numSec/3; npos++ )
00604 {
00605 nneg = npos-1;
00606 for( nzero=0; nzero<numSec/3; nzero++ )
00607 {
00608 if (++counter < numMulAn) {
00609 nt = npos+nneg+nzero;
00610 if ( (nt>1) && (nt<=numSec) ) {
00611 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00612 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00613 if (std::fabs(dum) < 1.0) {
00614 if( test >= 1.0e-10 )excs += dum*test;
00615 } else {
00616 excs += dum*test;
00617 }
00618
00619 if (ran < excs) goto outOfLoopAn;
00620 }
00621 }
00622 }
00623 }
00624 inElastic = false;
00625 return;
00626 }
00627 outOfLoopAn:
00628 vecLen = 0;
00629 }
00630 }
00631
00632 nt = npos + nneg + nzero;
00633 while ( nt > 0)
00634 {
00635 G4double ran = G4UniformRand();
00636 if ( ran < (G4double)npos/nt)
00637 {
00638 if( npos > 0 )
00639 { pv[vecLen++] = PionPlus;
00640 npos--;
00641 }
00642 }
00643 else if ( ran < (G4double)(npos+nneg)/nt)
00644 {
00645 if( nneg > 0 )
00646 {
00647 pv[vecLen++] = PionMinus;
00648 nneg--;
00649 }
00650 }
00651 else
00652 {
00653 if( nzero > 0 )
00654 {
00655 pv[vecLen++] = PionZero;
00656 nzero--;
00657 }
00658 }
00659 nt = npos + nneg + nzero;
00660 }
00661 if (verboseLevel > 1)
00662 {
00663 G4cout << "Particles produced: " ;
00664 G4cout << pv[0].getName() << " " ;
00665 G4cout << pv[1].getName() << " " ;
00666 for (i=2; i < vecLen; i++)
00667 {
00668 G4cout << pv[i].getName() << " " ;
00669 }
00670 G4cout << G4endl;
00671 }
00672 return;
00673 }
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683