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
00038
00039 #include "G4HEPionMinusInelastic.hh"
00040 #include "globals.hh"
00041 #include "G4ios.hh"
00042 #include "G4PhysicalConstants.hh"
00043
00044 void G4HEPionMinusInelastic::ModelDescription(std::ostream& outFile) const
00045 {
00046 outFile << "G4HEPionMinusInelastic is one of the High Energy\n"
00047 << "Parameterized (HEP) models used to implement inelastic\n"
00048 << "pi- scattering from nuclei. It is a re-engineered\n"
00049 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
00050 << "initial collision products into backward- and forward-going\n"
00051 << "clusters which are then decayed into final state hadrons.\n"
00052 << "The model does not conserve energy on an event-by-event\n"
00053 << "basis. It may be applied to pi- with initial energies\n"
00054 << "above 20 GeV.\n";
00055 }
00056
00057
00058 G4HadFinalState*
00059 G4HEPionMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00060 G4Nucleus& targetNucleus)
00061 {
00062 G4HEVector* pv = new G4HEVector[MAXPART];
00063 const G4HadProjectile* aParticle = &aTrack;
00064 const G4double A = targetNucleus.GetA_asInt();
00065 const G4double Z = targetNucleus.GetZ_asInt();
00066 G4HEVector incidentParticle(aParticle);
00067
00068 G4double atomicNumber = Z;
00069 G4double atomicWeight = A;
00070
00071 G4int incidentCode = incidentParticle.getCode();
00072 G4double incidentMass = incidentParticle.getMass();
00073 G4double incidentTotalEnergy = incidentParticle.getEnergy();
00074
00075
00076
00077
00078 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00079
00080 if (incidentKineticEnergy < 1.)
00081 G4cout << "GHEPionMinusInelastic: incident energy < 1 GeV" << G4endl;
00082
00083 if (verboseLevel > 1) {
00084 G4cout << "G4HEPionMinusInelastic::ApplyYourself" << G4endl;
00085 G4cout << "incident particle " << incidentParticle.getName()
00086 << "mass " << incidentMass
00087 << "kinetic energy " << incidentKineticEnergy
00088 << G4endl;
00089 G4cout << "target material with (A,Z) = ("
00090 << atomicWeight << "," << atomicNumber << ")" << G4endl;
00091 }
00092
00093 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
00094 atomicWeight, atomicNumber);
00095 if (verboseLevel > 1)
00096 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00097
00098 incidentKineticEnergy -= inelasticity;
00099
00100 G4double excitationEnergyGNP = 0.;
00101 G4double excitationEnergyDTA = 0.;
00102
00103 G4double excitation = NuclearExcitation(incidentKineticEnergy,
00104 atomicWeight, atomicNumber,
00105 excitationEnergyGNP,
00106 excitationEnergyDTA);
00107 if (verboseLevel > 1)
00108 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
00109 << excitationEnergyDTA << G4endl;
00110
00111 incidentKineticEnergy -= excitation;
00112 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00113
00114
00115
00116
00117 G4HEVector targetParticle;
00118
00119 if (G4UniformRand() < atomicNumber/atomicWeight) {
00120 targetParticle.setDefinition("Proton");
00121 } else {
00122 targetParticle.setDefinition("Neutron");
00123 }
00124
00125 G4double targetMass = targetParticle.getMass();
00126 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00127 + targetMass*targetMass
00128 + 2.0*targetMass*incidentTotalEnergy);
00129 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00130
00131
00132
00133
00134
00135
00136
00137 G4bool inElastic = true;
00138
00139
00140 vecLength = 0;
00141
00142 if (verboseLevel > 1)
00143 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00144 << incidentCode << G4endl;
00145
00146 G4bool successful = false;
00147
00148 FirstIntInCasPionMinus(inElastic, availableEnergy, pv, vecLength,
00149 incidentParticle, targetParticle);
00150
00151 if (verboseLevel > 1)
00152 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00153
00154 if ((vecLength > 0) && (availableEnergy > 1.))
00155 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00156 pv, vecLength,
00157 incidentParticle, targetParticle);
00158
00159 HighEnergyCascading(successful, pv, vecLength,
00160 excitationEnergyGNP, excitationEnergyDTA,
00161 incidentParticle, targetParticle,
00162 atomicWeight, atomicNumber);
00163 if (!successful)
00164 HighEnergyClusterProduction(successful, pv, vecLength,
00165 excitationEnergyGNP, excitationEnergyDTA,
00166 incidentParticle, targetParticle,
00167 atomicWeight, atomicNumber);
00168 if (!successful)
00169 MediumEnergyCascading(successful, pv, vecLength,
00170 excitationEnergyGNP, excitationEnergyDTA,
00171 incidentParticle, targetParticle,
00172 atomicWeight, atomicNumber);
00173
00174 if (!successful)
00175 MediumEnergyClusterProduction(successful, pv, vecLength,
00176 excitationEnergyGNP, excitationEnergyDTA,
00177 incidentParticle, targetParticle,
00178 atomicWeight, atomicNumber);
00179 if (!successful)
00180 QuasiElasticScattering(successful, pv, vecLength,
00181 excitationEnergyGNP, excitationEnergyDTA,
00182 incidentParticle, targetParticle,
00183 atomicWeight, atomicNumber);
00184 if (!successful)
00185 ElasticScattering(successful, pv, vecLength,
00186 incidentParticle,
00187 atomicWeight, atomicNumber);
00188
00189 if (!successful)
00190 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00191 << G4endl;
00192
00193 FillParticleChange(pv, vecLength);
00194 delete [] pv;
00195 theParticleChange.SetStatusChange(stopAndKill);
00196 return &theParticleChange;
00197 }
00198
00199
00200 void
00201 G4HEPionMinusInelastic::FirstIntInCasPionMinus(G4bool& inElastic,
00202 const G4double availableEnergy,
00203 G4HEVector pv[],
00204 G4int& vecLen,
00205 const G4HEVector& incidentParticle,
00206 const G4HEVector& targetParticle)
00207
00208
00209
00210
00211
00212
00213
00214
00215 {
00216 static const G4double expxu = 82.;
00217 static const G4double expxl = -expxu;
00218
00219 static const G4double protb = 0.7;
00220 static const G4double neutb = 0.7;
00221 static const G4double c = 1.25;
00222
00223 static const G4int numMul = 1200;
00224 static const G4int numSec = 60;
00225
00226 G4int protonCode = Proton.getCode();
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 neutmul[numMul], neutnorm[numSec];
00233
00234
00235
00236
00237 G4int i, counter, nt, npos, nneg, nero;
00238
00239 if (first) {
00240
00241 first = false;
00242 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00243 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00244 counter = -1;
00245 for( npos=0; npos<(numSec/3); npos++ )
00246 {
00247 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ )
00248 {
00249 for( nero=0; nero<numSec/3; nero++ )
00250 {
00251 if( ++counter < numMul )
00252 {
00253 nt = npos+nneg+nero;
00254 if( (nt>0) && (nt<=numSec) )
00255 {
00256 protmul[counter] =
00257 pmltpc(npos,nneg,nero,nt,protb,c) ;
00258 protnorm[nt-1] += protmul[counter];
00259 }
00260 }
00261 }
00262 }
00263 }
00264 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00265 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00266 counter = -1;
00267 for( npos=0; npos<numSec/3; npos++ )
00268 {
00269 for( nneg=npos; nneg<=(npos+2); nneg++ )
00270 {
00271 for( nero=0; nero<numSec/3; nero++ )
00272 {
00273 if( ++counter < numMul )
00274 {
00275 nt = npos+nneg+nero;
00276 if( (nt>0) && (nt<=numSec) )
00277 {
00278 neutmul[counter] =
00279 pmltpc(npos,nneg,nero,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
00293
00294
00295
00296 pv[0] = incidentParticle;
00297 pv[1] = targetParticle;
00298 vecLen = 2;
00299
00300 if (!inElastic || (availableEnergy <= PionPlus.getMass()))
00301 return;
00302
00303
00304
00305
00306 npos = 0, nneg = 0, nero = 0;
00307 G4double eab = availableEnergy;
00308 G4int ieab = G4int( eab*5.0 );
00309
00310 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00311 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
00312 {
00313
00314
00315 G4double cech[] = {1., 0.95, 0.79, 0.32, 0.19, 0.16, 0.14, 0.12, 0.10, 0.08};
00316 G4int iplab = Imin(9, G4int( incidentTotalMomentum*5.));
00317 if( G4UniformRand() < cech[iplab] )
00318 {
00319 if( targetCode == protonCode )
00320 {
00321 pv[0] = PionZero;
00322 pv[1] = Neutron;
00323 return;
00324 }
00325 else
00326 {
00327 return;
00328 }
00329 }
00330
00331 G4double w0, wp, wm, wt, ran;
00332 if( targetCode == protonCode )
00333 {
00334 w0 = - sqr(1.+protb)/(2.*c*c);
00335 wp = w0 = std::exp(w0);
00336 wm = - sqr(-1.+protb)/(2.*c*c);
00337 wm = std::exp(wm);
00338 wp *= 10.;
00339 wt = w0+wp+wm;
00340 wp = w0+wp;
00341 ran = G4UniformRand();
00342 if( ran < w0/wt )
00343 { npos = 0; nneg = 0; nero = 1; }
00344 else if ( ran < wp/wt )
00345 { npos = 1; nneg = 0; nero = 0; }
00346 else
00347 { npos = 0; nneg = 1; nero = 0; }
00348 }
00349 else
00350 {
00351 w0 = -sqr(1.+neutb)/(2.*c*c);
00352 wm = -sqr(-1.+neutb)/(2.*c*c);
00353 w0 = std::exp(w0);
00354 wm = std::exp(wm);
00355 if( G4UniformRand() < w0/(w0+wm) )
00356 { npos = 0; nneg = 0; nero = 1; }
00357 else
00358 { npos = 0; nneg = 1; nero = 0; }
00359 }
00360 }
00361 else
00362 {
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( Amin( expxu, Amax( 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 for (nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++) {
00390 for (nero=0; nero<numSec/3; nero++) {
00391 if (++counter < numMul) {
00392 nt = npos+nneg+nero;
00393 if ( (nt>0) && (nt<=numSec) ) {
00394 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00395 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00396 if (std::fabs(dum) < 1.0) {
00397 if( test >= 1.0e-10 )excs += dum*test;
00398 } else {
00399 excs += dum*test;
00400 }
00401 if (ran < excs) goto outOfLoop;
00402 }
00403 }
00404 }
00405 }
00406 }
00407
00408
00409 inElastic = false;
00410 return;
00411 }
00412 else
00413 {
00414 counter = -1;
00415 for (npos=0; npos<numSec/3; npos++) {
00416 for (nneg=npos; nneg<=(npos+2); nneg++) {
00417 for (nero=0; nero<numSec/3; nero++) {
00418 if (++counter < numMul) {
00419 nt = npos+nneg+nero;
00420 if ( (nt>=1) && (nt<=numSec) ) {
00421 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00422 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00423 if (std::fabs(dum) < 1.0) {
00424 if( test >= 1.0e-10 )excs += dum*test;
00425 } else {
00426 excs += dum*test;
00427 }
00428 if (ran < excs) goto outOfLoop;
00429 }
00430 }
00431 }
00432 }
00433 }
00434
00435 inElastic = false;
00436 return;
00437 }
00438 }
00439 outOfLoop:
00440
00441 if( targetCode == protonCode)
00442 {
00443 if( npos == (1+nneg))
00444 {
00445 pv[1] = Neutron;
00446 }
00447 else if (npos == nneg)
00448 {
00449 if( G4UniformRand() < 0.75)
00450 {
00451 }
00452 else
00453 {
00454 pv[0] = PionZero;
00455 pv[1] = Neutron;
00456 }
00457 }
00458 else
00459 {
00460 pv[0] = PionZero;
00461 }
00462 }
00463 else
00464 {
00465 if( npos == (nneg-1))
00466 {
00467 if( G4UniformRand() < 0.5)
00468 {
00469 pv[1] = Proton;
00470 }
00471 else
00472 {
00473 pv[0] = PionZero;
00474 }
00475 }
00476 else if ( npos == nneg)
00477 {
00478 }
00479 else
00480 {
00481 pv[0] = PionZero;
00482 pv[1] = Proton;
00483 }
00484 }
00485
00486
00487 nt = npos + nneg + nero;
00488 while ( nt > 0)
00489 {
00490 G4double ran = G4UniformRand();
00491 if ( ran < (G4double)npos/nt)
00492 {
00493 if( npos > 0 )
00494 { pv[vecLen++] = PionPlus;
00495 npos--;
00496 }
00497 }
00498 else if ( ran < (G4double)(npos+nneg)/nt)
00499 {
00500 if( nneg > 0 )
00501 {
00502 pv[vecLen++] = PionMinus;
00503 nneg--;
00504 }
00505 }
00506 else
00507 {
00508 if( nero > 0 )
00509 {
00510 pv[vecLen++] = PionZero;
00511 nero--;
00512 }
00513 }
00514 nt = npos + nneg + nero;
00515 }
00516 if (verboseLevel > 1)
00517 {
00518 G4cout << "Particles produced: " ;
00519 G4cout << pv[0].getName() << " " ;
00520 G4cout << pv[1].getName() << " " ;
00521 for (i=2; i < vecLen; i++)
00522 {
00523 G4cout << pv[i].getName() << " " ;
00524 }
00525 G4cout << G4endl;
00526 }
00527 return;
00528 }
00529
00530
00531
00532
00533
00534
00535
00536
00537