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