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