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 "G4HEAntiSigmaPlusInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042
00043 void G4HEAntiSigmaPlusInelastic::ModelDescription(std::ostream& outFile) const
00044 {
00045 outFile << "G4HEAntiSigmaPlusInelastic 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 G4HEAntiSigmaPlusInelastic::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 << "GHEAntiSigmaPlusInelastic: incident energy < 1 GeV" << G4endl;
00078
00079 if (verboseLevel > 1) {
00080 G4cout << "G4HEAntiSigmaPlusInelastic::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 FirstIntInCasAntiSigmaPlus(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
00146 HighEnergyCascading(successful, pv, vecLength,
00147 excitationEnergyGNP, excitationEnergyDTA,
00148 incidentParticle, targetParticle,
00149 atomicWeight, atomicNumber);
00150 if (!successful)
00151 HighEnergyClusterProduction(successful, pv, vecLength,
00152 excitationEnergyGNP, excitationEnergyDTA,
00153 incidentParticle, targetParticle,
00154 atomicWeight, atomicNumber);
00155 if (!successful)
00156 MediumEnergyCascading(successful, pv, vecLength,
00157 excitationEnergyGNP, excitationEnergyDTA,
00158 incidentParticle, targetParticle,
00159 atomicWeight, atomicNumber);
00160
00161 if (!successful)
00162 MediumEnergyClusterProduction(successful, pv, vecLength,
00163 excitationEnergyGNP, excitationEnergyDTA,
00164 incidentParticle, targetParticle,
00165 atomicWeight, atomicNumber);
00166 if (!successful)
00167 QuasiElasticScattering(successful, pv, vecLength,
00168 excitationEnergyGNP, excitationEnergyDTA,
00169 incidentParticle, targetParticle,
00170 atomicWeight, atomicNumber);
00171 if (!successful)
00172 ElasticScattering(successful, pv, vecLength,
00173 incidentParticle,
00174 atomicWeight, atomicNumber);
00175
00176 if (!successful)
00177 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00178 << G4endl;
00179
00180 FillParticleChange(pv, vecLength);
00181 delete [] pv;
00182 theParticleChange.SetStatusChange(stopAndKill);
00183 return &theParticleChange;
00184 }
00185
00186
00187 void
00188 G4HEAntiSigmaPlusInelastic::FirstIntInCasAntiSigmaPlus(G4bool& inElastic,
00189 const G4double availableEnergy,
00190 G4HEVector pv[],
00191 G4int& vecLen,
00192 const G4HEVector& incidentParticle,
00193 const G4HEVector& targetParticle,
00194 const G4double atomicWeight)
00195
00196
00197
00198
00199
00200
00201
00202
00203 {
00204 static const G4double expxu = 82.;
00205 static const G4double expxl = -expxu;
00206
00207 static const G4double protb = 0.7;
00208 static const G4double neutb = 0.7;
00209 static const G4double c = 1.25;
00210
00211 static const G4int numMul = 1200;
00212 static const G4int numMulAn = 400;
00213 static const G4int numSec = 60;
00214
00215 G4int protonCode = Proton.getCode();
00216
00217 G4int targetCode = targetParticle.getCode();
00218 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00219
00220 static G4bool first = true;
00221 static G4double protmul[numMul], protnorm[numSec];
00222 static G4double protmulAn[numMulAn],protnormAn[numSec];
00223 static G4double neutmul[numMul], neutnorm[numSec];
00224 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00225
00226
00227
00228
00229 G4int i, counter, nt, npos, nneg, nzero;
00230
00231 if (first) {
00232
00233 first = false;
00234 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
00235 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
00236 counter = -1;
00237 for (npos = 0; npos < (numSec/3); npos++) {
00238 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
00239 for (nzero = 0; nzero < numSec/3; nzero++) {
00240 if (++counter < numMul) {
00241 nt = npos+nneg+nzero;
00242 if ((nt>0) && (nt<=numSec) ) {
00243 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00244 protnorm[nt-1] += protmul[counter];
00245 }
00246 }
00247 }
00248 }
00249 }
00250
00251 for (i = 0; i < numMul; i++) neutmul[i] = 0.0;
00252 for (i = 0; i < numSec; i++) neutnorm[i] = 0.0;
00253 counter = -1;
00254 for (npos = 0; npos < numSec/3; npos++) {
00255 for (nneg = npos; nneg <= (npos+2); nneg++) {
00256 for (nzero = 0; nzero < numSec/3; nzero++) {
00257 if (++counter < numMul) {
00258 nt = npos+nneg+nzero;
00259 if ((nt>0) && (nt<=numSec) ) {
00260 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00261 neutnorm[nt-1] += neutmul[counter];
00262 }
00263 }
00264 }
00265 }
00266 }
00267
00268 for (i = 0; i < numSec; i++) {
00269 if (protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00270 if (neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00271 }
00272
00273
00274 for (i = 0; i < numMulAn; i++) protmulAn[i] = 0.0;
00275 for (i = 0; i < numSec; i++) protnormAn[i] = 0.0;
00276 counter = -1;
00277 for (npos = 1; npos < (numSec/3); npos++) {
00278 nneg = npos;
00279 for (nzero = 0; nzero < numSec/3; nzero++) {
00280 if (++counter < numMulAn) {
00281 nt = npos+nneg+nzero;
00282 if ((nt>1) && (nt<=numSec) ) {
00283 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00284 protnormAn[nt-1] += protmulAn[counter];
00285 }
00286 }
00287 }
00288 }
00289
00290 for (i = 0; i < numMulAn; i++) neutmulAn[i] = 0.0;
00291 for (i = 0; i < numSec; i++) neutnormAn[i] = 0.0;
00292 counter = -1;
00293 for (npos = 0; npos < numSec/3; npos++) {
00294 nneg = npos+1;
00295 for (nzero = 0; nzero < numSec/3; nzero++) {
00296 if (++counter < numMulAn) {
00297 nt = npos+nneg+nzero;
00298 if ((nt>1) && (nt<=numSec) ) {
00299 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00300 neutnormAn[nt-1] += neutmulAn[counter];
00301 }
00302 }
00303 }
00304 }
00305 for (i = 0; i < numSec; i++) {
00306 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00307 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00308 }
00309 }
00310
00311
00312
00313 pv[0] = incidentParticle;
00314 pv[1] = targetParticle;
00315 vecLen = 2;
00316
00317 if (!inElastic) {
00318 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00319
00320 G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
00321 if (G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
00322 G4double ran = G4UniformRand();
00323
00324 if (targetCode == protonCode) {
00325 if (ran < 0.2) {
00326 pv[0] = Proton;
00327 pv[1] = AntiSigmaPlus;
00328 } else if (ran < 0.4) {
00329 pv[0] = AntiLambda;
00330 pv[1] = Neutron;
00331 } else if (ran < 0.6) {
00332 pv[0] = Neutron;
00333 pv[1] = AntiLambda;
00334 } else if (ran < 0.8) {
00335 pv[0] = Neutron;
00336 pv[1] = AntiSigmaZero;
00337 } else {
00338 pv[0] = AntiSigmaZero;
00339 pv[1] = Neutron;
00340 }
00341 } else {
00342 pv[0] = Neutron;
00343 pv[1] = AntiSigmaPlus;
00344 }
00345 }
00346
00347 return;
00348 }
00349 else if (availableEnergy <= PionPlus.getMass()) return;
00350
00351
00352 npos = 0; nneg = 0; nzero = 0;
00353 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
00354 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
00355 0.39, 0.36, 0.33, 0.10, 0.01};
00356 G4int iplab = G4int( incidentTotalMomentum*10.);
00357 if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. );
00358 if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. );
00359 if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.);
00360 iplab = std::min(24, iplab);
00361
00362 if ( G4UniformRand() > anhl[iplab] ) {
00363
00364
00365 G4double aleab = std::log(availableEnergy);
00366 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00367 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00368
00369
00370
00371 G4double test, dum, anpn = 0.0;
00372
00373 for (nt=1; nt<=numSec; nt++) {
00374 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00375 dum = pi*nt/(2.0*n*n);
00376 if (std::fabs(dum) < 1.0) {
00377 if( test >= 1.0e-10 )anpn += dum*test;
00378 } else {
00379 anpn += dum*test;
00380 }
00381 }
00382
00383 G4double ran = G4UniformRand();
00384 G4double excs = 0.0;
00385 if( targetCode == protonCode )
00386 {
00387 counter = -1;
00388 for( npos=0; npos<numSec/3; npos++ )
00389 {
00390 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
00391 {
00392 for( nzero=0; nzero<numSec/3; nzero++ )
00393 {
00394 if( ++counter < numMul )
00395 {
00396 nt = npos+nneg+nzero;
00397 if ( (nt>0) && (nt<=numSec) ) {
00398 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00399 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00400 if (std::fabs(dum) < 1.0) {
00401 if( test >= 1.0e-10 )excs += dum*test;
00402 } else {
00403 excs += dum*test;
00404 }
00405
00406 if (ran < excs) goto outOfLoop;
00407 }
00408 }
00409 }
00410 }
00411 }
00412
00413
00414 inElastic = false;
00415 return;
00416 }
00417 else
00418 {
00419 counter = -1;
00420 for( npos=0; npos<numSec/3; npos++ )
00421 {
00422 for( nneg=npos; nneg<=(npos+2); nneg++ )
00423 {
00424 for( nzero=0; nzero<numSec/3; nzero++ )
00425 {
00426 if( ++counter < numMul )
00427 {
00428 nt = npos+nneg+nzero;
00429 if ( (nt>0) && (nt<=numSec) ) {
00430 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00431 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00432 if (std::fabs(dum) < 1.0) {
00433 if( test >= 1.0e-10 )excs += dum*test;
00434 } else {
00435 excs += dum*test;
00436 }
00437
00438 if (ran < excs) goto outOfLoop;
00439 }
00440 }
00441 }
00442 }
00443 }
00444
00445 inElastic = false;
00446 return;
00447 }
00448
00449 outOfLoop:
00450
00451 ran = G4UniformRand();
00452
00453 if( targetCode == protonCode)
00454 {
00455 if( npos == nneg)
00456 {
00457 if (ran < 0.50)
00458 {
00459 }
00460 else if (ran < 0.75)
00461 {
00462 pv[0] = AntiSigmaZero;
00463 pv[1] = Neutron;
00464 }
00465 else
00466 {
00467 pv[0] = AntiLambda;
00468 pv[1] = Neutron;
00469 }
00470 }
00471 else if (npos == (nneg-1))
00472 {
00473 if( ran < 0.50)
00474 {
00475 pv[0] = AntiLambda;
00476 }
00477 else
00478 {
00479 pv[0] = AntiSigmaZero;
00480 }
00481 }
00482 else
00483 {
00484 pv[1] = Neutron;
00485 }
00486 }
00487 else
00488 {
00489 if( npos == nneg)
00490 {
00491 }
00492 else if ( npos == (nneg-1))
00493 {
00494 if (ran < 0.5)
00495 {
00496 pv[1] = Proton;
00497 }
00498 else if (ran < 0.75)
00499 {
00500 pv[0] = AntiLambda;
00501 }
00502 else
00503 {
00504 pv[0] = AntiSigmaZero;
00505 }
00506 }
00507 else
00508 {
00509 if (ran < 0.5)
00510 {
00511 pv[0] = AntiLambda;
00512 pv[1] = Proton;
00513 }
00514 else
00515 {
00516 pv[0] = AntiSigmaZero;
00517 pv[1] = Proton;
00518 }
00519 }
00520 }
00521 }
00522 else
00523 {
00524 if ( availableEnergy > 2. * PionPlus.getMass() )
00525 {
00526
00527 G4double aleab = std::log(availableEnergy);
00528 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00529 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00530
00531
00532
00533 G4double test, dum, anpn = 0.0;
00534
00535 for (nt=2; nt<=numSec; nt++) {
00536 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00537 dum = pi*nt/(2.0*n*n);
00538 if (std::fabs(dum) < 1.0) {
00539 if( test >= 1.0e-10 )anpn += dum*test;
00540 } else {
00541 anpn += dum*test;
00542 }
00543 }
00544
00545 G4double ran = G4UniformRand();
00546 G4double excs = 0.0;
00547 if( targetCode == protonCode )
00548 {
00549 counter = -1;
00550 for( npos=1; npos<numSec/3; npos++ )
00551 {
00552 nneg = npos;
00553 for( nzero=0; nzero<numSec/3; nzero++ )
00554 {
00555 if( ++counter < numMulAn )
00556 {
00557 nt = npos+nneg+nzero;
00558 if ( (nt>1) && (nt<=numSec) ) {
00559 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00560 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00561 if (std::fabs(dum) < 1.0) {
00562 if( test >= 1.0e-10 )excs += dum*test;
00563 } else {
00564 excs += dum*test;
00565 }
00566
00567 if (ran < excs) goto outOfLoopAn;
00568 }
00569 }
00570 }
00571 }
00572
00573 inElastic = false;
00574 return;
00575 }
00576 else
00577 {
00578 counter = -1;
00579 for( npos=0; npos<numSec/3; npos++ )
00580 {
00581 nneg = npos+1;
00582 for( nzero=0; nzero<numSec/3; nzero++ )
00583 {
00584 if( ++counter < numMulAn )
00585 {
00586 nt = npos+nneg+nzero;
00587 if ( (nt>1) && (nt<=numSec) ) {
00588 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00589 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00590 if (std::fabs(dum) < 1.0) {
00591 if( test >= 1.0e-10 )excs += dum*test;
00592 } else {
00593 excs += dum*test;
00594 }
00595
00596 if (ran < excs) goto outOfLoopAn;
00597 }
00598 }
00599 }
00600 }
00601 inElastic = false;
00602 return;
00603 }
00604 outOfLoopAn:
00605 vecLen = 0;
00606 }
00607 }
00608
00609 nt = npos + nneg + nzero;
00610 while ( nt > 0)
00611 {
00612 G4double ran = G4UniformRand();
00613 if ( ran < (G4double)npos/nt)
00614 {
00615 if( npos > 0 )
00616 { pv[vecLen++] = PionPlus;
00617 npos--;
00618 }
00619 }
00620 else if ( ran < (G4double)(npos+nneg)/nt)
00621 {
00622 if( nneg > 0 )
00623 {
00624 pv[vecLen++] = PionMinus;
00625 nneg--;
00626 }
00627 }
00628 else
00629 {
00630 if( nzero > 0 )
00631 {
00632 pv[vecLen++] = PionZero;
00633 nzero--;
00634 }
00635 }
00636 nt = npos + nneg + nzero;
00637 }
00638 if (verboseLevel > 1)
00639 {
00640 G4cout << "Particles produced: " ;
00641 G4cout << pv[0].getName() << " " ;
00642 G4cout << pv[1].getName() << " " ;
00643 for (i=2; i < vecLen; i++)
00644 {
00645 G4cout << pv[i].getName() << " " ;
00646 }
00647 G4cout << G4endl;
00648 }
00649 return;
00650 }
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660