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