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 "G4HEAntiOmegaMinusInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042
00043 void G4HEAntiOmegaMinusInelastic::ModelDescription(std::ostream& outFile) const
00044 {
00045 outFile << "G4HEAntiOmegaMinusInelastic is one of the High Energy\n"
00046 << "Parameterized (HEP) models used to implement inelastic\n"
00047 << "anti-Omega- 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-Omega- with initial\n"
00053 << "energies above 20 GeV.\n";
00054 }
00055
00056
00057 G4HadFinalState*
00058 G4HEAntiOmegaMinusInelastic::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 << "GHEAntiOmegaMinusInelastic: incident energy < 1 GeV" << G4endl;
00078
00079 if (verboseLevel > 1) {
00080 G4cout << "G4HEAntiOmegaMinusInelastic::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
00108 incidentKineticEnergy -= excitation;
00109 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00110
00111
00112
00113
00114 G4HEVector targetParticle;
00115 if (G4UniformRand() < atomicNumber/atomicWeight) {
00116 targetParticle.setDefinition("Proton");
00117 } else {
00118 targetParticle.setDefinition("Neutron");
00119 }
00120
00121 G4double targetMass = targetParticle.getMass();
00122 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00123 + targetMass*targetMass
00124 + 2.0*targetMass*incidentTotalEnergy);
00125 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00126
00127 G4bool inElastic = true;
00128 vecLength = 0;
00129
00130 if (verboseLevel > 1)
00131 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00132 << incidentCode << G4endl;
00133
00134 G4bool successful = false;
00135
00136 FirstIntInCasAntiOmegaMinus(inElastic, availableEnergy, pv, vecLength,
00137 incidentParticle, targetParticle, atomicWeight);
00138
00139 if (verboseLevel > 1)
00140 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00141
00142 if ((vecLength > 0) && (availableEnergy > 1.))
00143 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00144 pv, vecLength,
00145 incidentParticle, targetParticle);
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 G4HEAntiOmegaMinusInelastic::FirstIntInCasAntiOmegaMinus(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 static const G4double expxu = 82.;
00201 static const G4double expxl = -expxu;
00202
00203 static const G4double protb = 0.7;
00204 static const G4double neutb = 0.7;
00205 static const G4double c = 1.25;
00206
00207 static const G4int numMul = 1200;
00208 static const G4int numMulAn = 400;
00209 static const G4int numSec = 60;
00210
00211 G4int protonCode = Proton.getCode();
00212
00213 G4int targetCode = targetParticle.getCode();
00214 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00215
00216 static G4bool first = true;
00217 static G4double protmul[numMul], protnorm[numSec];
00218 static G4double protmulAn[numMulAn],protnormAn[numSec];
00219 static G4double neutmul[numMul], neutnorm[numSec];
00220 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00221
00222
00223
00224
00225 G4int i, counter, nt, npos, nneg, nzero;
00226
00227 if( first )
00228 {
00229 first = false;
00230 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
00231 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
00232 counter = -1;
00233 for( npos=0; npos<(numSec/3); npos++ )
00234 {
00235 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
00236 {
00237 for( nzero=0; nzero<numSec/3; nzero++ )
00238 {
00239 if( ++counter < numMul )
00240 {
00241 nt = npos+nneg+nzero;
00242 if( (nt>0) && (nt<=numSec) )
00243 {
00244 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00245 protnorm[nt-1] += protmul[counter];
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 {
00256 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
00257 {
00258 for( nzero=0; nzero<numSec/3; nzero++ )
00259 {
00260 if( ++counter < numMul )
00261 {
00262 nt = npos+nneg+nzero;
00263 if( (nt>0) && (nt<=numSec) )
00264 {
00265 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00266 neutnorm[nt-1] += neutmul[counter];
00267 }
00268 }
00269 }
00270 }
00271 }
00272 for( i=0; i<numSec; i++ )
00273 {
00274 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00275 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00276 }
00277
00278 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
00279 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
00280 counter = -1;
00281 for( npos=1; npos<(numSec/3); npos++ )
00282 {
00283 nneg = std::max(0,npos-1);
00284 for( nzero=0; nzero<numSec/3; nzero++ )
00285 {
00286 if( ++counter < numMulAn )
00287 {
00288 nt = npos+nneg+nzero;
00289 if( (nt>1) && (nt<=numSec) )
00290 {
00291 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00292 protnormAn[nt-1] += protmulAn[counter];
00293 }
00294 }
00295 }
00296 }
00297 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
00298 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
00299 counter = -1;
00300 for( npos=0; npos<numSec/3; npos++ )
00301 {
00302 nneg = npos;
00303 for( nzero=0; nzero<numSec/3; nzero++ )
00304 {
00305 if( ++counter < numMulAn )
00306 {
00307 nt = npos+nneg+nzero;
00308 if( (nt>1) && (nt<=numSec) )
00309 {
00310 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00311 neutnormAn[nt-1] += neutmulAn[counter];
00312 }
00313 }
00314 }
00315 }
00316 for( i=0; i<numSec; i++ )
00317 {
00318 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00319 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00320 }
00321 }
00322
00323
00324
00325
00326 pv[0] = incidentParticle;
00327 pv[1] = targetParticle;
00328 vecLen = 2;
00329
00330 if( !inElastic )
00331 {
00332 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00333
00334 G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5 ));
00335 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00336 {
00337 G4double ran = G4UniformRand();
00338
00339 if ( targetCode == protonCode)
00340 {
00341 if(ran < 0.2)
00342 {
00343 pv[0] = AntiSigmaZero;
00344 }
00345 else if (ran < 0.4)
00346 {
00347 pv[0] = AntiSigmaMinus;
00348 pv[1] = Neutron;
00349 }
00350 else if (ran < 0.6)
00351 {
00352 pv[0] = Proton;
00353 pv[1] = AntiLambda;
00354 }
00355 else if (ran < 0.8)
00356 {
00357 pv[0] = Proton;
00358 pv[1] = AntiSigmaZero;
00359 }
00360 else
00361 {
00362 pv[0] = Neutron;
00363 pv[1] = AntiSigmaMinus;
00364 }
00365 }
00366 else
00367 {
00368 if (ran < 0.2)
00369 {
00370 pv[0] = AntiSigmaZero;
00371 }
00372 else if (ran < 0.4)
00373 {
00374 pv[0] = AntiSigmaPlus;
00375 pv[1] = Proton;
00376 }
00377 else if (ran < 0.6)
00378 {
00379 pv[0] = Neutron;
00380 pv[1] = AntiLambda;
00381 }
00382 else if (ran < 0.8)
00383 {
00384 pv[0] = Neutron;
00385 pv[1] = AntiSigmaZero;
00386 }
00387 else
00388 {
00389 pv[0] = Proton;
00390 pv[1] = AntiSigmaPlus;
00391 }
00392 }
00393 }
00394 return;
00395 }
00396 else if (availableEnergy <= PionPlus.getMass())
00397 return;
00398
00399
00400
00401 npos = 0; nneg = 0; nzero = 0;
00402 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
00403 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
00404 0.39, 0.36, 0.33, 0.10, 0.01};
00405 G4int iplab = G4int( incidentTotalMomentum*10.);
00406 if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. );
00407 if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. );
00408 if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.);
00409 iplab = std::min(24, iplab);
00410
00411 if (G4UniformRand() > anhl[iplab]) {
00412
00413
00414 G4double aleab = std::log(availableEnergy);
00415 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00416 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00417
00418
00419
00420 G4double test, dum, anpn = 0.0;
00421
00422 for (nt=1; nt<=numSec; nt++) {
00423 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00424 dum = pi*nt/(2.0*n*n);
00425 if (std::fabs(dum) < 1.0) {
00426 if( test >= 1.0e-10 )anpn += dum*test;
00427 } else {
00428 anpn += dum*test;
00429 }
00430 }
00431
00432 G4double ran = G4UniformRand();
00433 G4double excs = 0.0;
00434 if( targetCode == protonCode )
00435 {
00436 counter = -1;
00437 for (npos=0; npos<numSec/3; npos++) {
00438 for (nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++) {
00439 for (nzero=0; nzero<numSec/3; nzero++) {
00440 if (++counter < numMul) {
00441 nt = npos+nneg+nzero;
00442 if ( (nt>0) && (nt<=numSec) ) {
00443 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00444 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00445 if (std::fabs(dum) < 1.0) {
00446 if( test >= 1.0e-10 )excs += dum*test;
00447 } else {
00448 excs += dum*test;
00449 }
00450
00451 if (ran < excs) goto outOfLoop;
00452 }
00453 }
00454 }
00455 }
00456 }
00457
00458
00459 inElastic = false;
00460 return;
00461 }
00462 else
00463 {
00464 counter = -1;
00465 for (npos=0; npos<numSec/3; npos++) {
00466 for (nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++) {
00467 for (nzero=0; nzero<numSec/3; nzero++) {
00468 if (++counter < numMul) {
00469 nt = npos+nneg+nzero;
00470 if ( (nt>0) && (nt<=numSec) ) {
00471 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00472 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00473 if (std::fabs(dum) < 1.0) {
00474 if( test >= 1.0e-10 )excs += dum*test;
00475 } else {
00476 excs += dum*test;
00477 }
00478
00479 if (ran < excs) goto outOfLoop;
00480 }
00481 }
00482 }
00483 }
00484 }
00485
00486 inElastic = false;
00487 return;
00488 }
00489
00490 outOfLoop:
00491
00492 ran = G4UniformRand();
00493
00494 if( targetCode == protonCode)
00495 {
00496 if( npos == nneg)
00497 {
00498 if (ran < 0.40)
00499 {
00500 }
00501 else if (ran < 0.8)
00502 {
00503 pv[0] = AntiSigmaZero;
00504 }
00505 else
00506 {
00507 pv[0] = AntiSigmaMinus;
00508 pv[1] = Neutron;
00509 }
00510 }
00511 else if (npos == (nneg+1))
00512 {
00513 if( ran < 0.25)
00514 {
00515 pv[1] = Neutron;
00516 }
00517 else if (ran < 0.5)
00518 {
00519 pv[0] = AntiSigmaZero;
00520 pv[1] = Neutron;
00521 }
00522 else
00523 {
00524 pv[0] = AntiSigmaPlus;
00525 }
00526 }
00527 else if (npos == (nneg-1))
00528 {
00529 pv[0] = AntiSigmaMinus;
00530 }
00531 else
00532 {
00533 pv[0] = AntiSigmaPlus;
00534 pv[1] = Neutron;
00535 }
00536 }
00537 else
00538 {
00539 if( npos == nneg)
00540 {
00541 if (ran < 0.4)
00542 {
00543 }
00544 else if(ran < 0.8)
00545 {
00546 pv[0] = AntiSigmaZero;
00547 }
00548 else
00549 {
00550 pv[0] = AntiSigmaPlus;
00551 pv[1] = Proton;
00552 }
00553 }
00554 else if ( npos == (nneg-1))
00555 {
00556 if (ran < 0.5)
00557 {
00558 pv[0] = AntiSigmaMinus;
00559 }
00560 else if (ran < 0.75)
00561 {
00562 pv[1] = Proton;
00563 }
00564 else
00565 {
00566 pv[0] = AntiSigmaZero;
00567 pv[1] = Proton;
00568 }
00569 }
00570 else if (npos == (nneg+1))
00571 {
00572 pv[0] = AntiSigmaPlus;
00573 }
00574 else
00575 {
00576 pv[0] = AntiSigmaMinus;
00577 pv[1] = Proton;
00578 }
00579 }
00580
00581 }
00582 else
00583 {
00584 if ( availableEnergy > 2. * PionPlus.getMass() )
00585 {
00586
00587 G4double aleab = std::log(availableEnergy);
00588 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00589 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00590
00591
00592
00593 G4double test, dum, anpn = 0.0;
00594
00595 for (nt=2; nt<=numSec; nt++) {
00596 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00597 dum = pi*nt/(2.0*n*n);
00598 if (std::fabs(dum) < 1.0) {
00599 if( test >= 1.0e-10 )anpn += dum*test;
00600 } else {
00601 anpn += dum*test;
00602 }
00603 }
00604
00605 G4double ran = G4UniformRand();
00606 G4double excs = 0.0;
00607 if( targetCode == protonCode )
00608 {
00609 counter = -1;
00610 for( npos=1; npos<numSec/3; npos++ )
00611 {
00612 nneg = npos-1;
00613 for( nzero=0; nzero<numSec/3; nzero++ )
00614 {
00615 if( ++counter < numMulAn )
00616 {
00617 nt = npos+nneg+nzero;
00618 if ( (nt>1) && (nt<=numSec) ) {
00619 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00620 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00621 if (std::fabs(dum) < 1.0) {
00622 if( test >= 1.0e-10 )excs += dum*test;
00623 } else {
00624 excs += dum*test;
00625 }
00626
00627 if (ran < excs) goto outOfLoopAn;
00628 }
00629 }
00630 }
00631 }
00632
00633 inElastic = false;
00634 return;
00635 }
00636 else
00637 {
00638 counter = -1;
00639 for( npos=0; npos<numSec/3; npos++ )
00640 {
00641 nneg = npos;
00642 for( nzero=0; nzero<numSec/3; nzero++ )
00643 {
00644 if( ++counter < numMulAn )
00645 {
00646 nt = npos+nneg+nzero;
00647 if ( (nt>1) && (nt<=numSec) ) {
00648 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00649 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00650 if (std::fabs(dum) < 1.0) {
00651 if( test >= 1.0e-10 )excs += dum*test;
00652 } else {
00653 excs += dum*test;
00654 }
00655
00656 if (ran < excs) goto outOfLoopAn;
00657 }
00658 }
00659 }
00660 }
00661 inElastic = false;
00662 return;
00663 }
00664 outOfLoopAn:
00665 vecLen = 0;
00666 }
00667 }
00668
00669 nt = npos + nneg + nzero;
00670 while ( nt > 0)
00671 {
00672 G4double ran = G4UniformRand();
00673 if ( ran < (G4double)npos/nt)
00674 {
00675 if( npos > 0 )
00676 { pv[vecLen++] = PionPlus;
00677 npos--;
00678 }
00679 }
00680 else if ( ran < (G4double)(npos+nneg)/nt)
00681 {
00682 if( nneg > 0 )
00683 {
00684 pv[vecLen++] = PionMinus;
00685 nneg--;
00686 }
00687 }
00688 else
00689 {
00690 if( nzero > 0 )
00691 {
00692 pv[vecLen++] = PionZero;
00693 nzero--;
00694 }
00695 }
00696 nt = npos + nneg + nzero;
00697 }
00698 if (verboseLevel > 1)
00699 {
00700 G4cout << "Particles produced: " ;
00701 G4cout << pv[0].getName() << " " ;
00702 G4cout << pv[1].getName() << " " ;
00703 for (i=2; i < vecLen; i++)
00704 {
00705 G4cout << pv[i].getName() << " " ;
00706 }
00707 G4cout << G4endl;
00708 }
00709 return;
00710 }
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720