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 "G4HEPionPlusInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041
00042 void G4HEPionPlusInelastic::ModelDescription(std::ostream& outFile) const
00043 {
00044 outFile << "G4HEPionPlusInelastic is one of the High Energy\n"
00045 << "Parameterized (HEP) models used to implement inelastic\n"
00046 << "pi+ 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 pi+ with initial energies\n"
00052 << "above 20 GeV.\n";
00053 }
00054
00055
00056 G4HadFinalState*
00057 G4HEPionPlusInelastic::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 << "G4HEPionPlusInelastic: incident energy < 1 GeV" << G4endl;
00080
00081 if (verboseLevel > 1) {
00082 G4cout << "G4HEPionPlusInelastic::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
00130 vecLength = 0;
00131
00132 if (verboseLevel > 1)
00133 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00134 << incidentCode << G4endl;
00135
00136 G4bool successful = false;
00137
00138 FirstIntInCasPionPlus(inElastic, availableEnergy, pv, vecLength,
00139 incidentParticle, targetParticle, atomicWeight);
00140
00141 if (verboseLevel > 1)
00142 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00143
00144 if ((vecLength > 0) && (availableEnergy > 1.))
00145 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00146 pv, vecLength,
00147 incidentParticle, targetParticle);
00148
00149 HighEnergyCascading(successful, pv, vecLength,
00150 excitationEnergyGNP, excitationEnergyDTA,
00151 incidentParticle, targetParticle,
00152 atomicWeight, atomicNumber);
00153 if (!successful)
00154 HighEnergyClusterProduction(successful, pv, vecLength,
00155 excitationEnergyGNP, excitationEnergyDTA,
00156 incidentParticle, targetParticle,
00157 atomicWeight, atomicNumber);
00158 if (!successful)
00159 MediumEnergyCascading(successful, pv, vecLength,
00160 excitationEnergyGNP, excitationEnergyDTA,
00161 incidentParticle, targetParticle,
00162 atomicWeight, atomicNumber);
00163
00164 if (!successful)
00165 MediumEnergyClusterProduction(successful, pv, vecLength,
00166 excitationEnergyGNP, excitationEnergyDTA,
00167 incidentParticle, targetParticle,
00168 atomicWeight, atomicNumber);
00169 if (!successful)
00170 QuasiElasticScattering(successful, pv, vecLength,
00171 excitationEnergyGNP, excitationEnergyDTA,
00172 incidentParticle, targetParticle,
00173 atomicWeight, atomicNumber);
00174 if (!successful)
00175 ElasticScattering(successful, pv, vecLength,
00176 incidentParticle,
00177 atomicWeight, atomicNumber);
00178
00179 if (!successful)
00180 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00181 << G4endl;
00182
00183 FillParticleChange(pv, vecLength);
00184 delete [] pv;
00185 theParticleChange.SetStatusChange(stopAndKill);
00186 return &theParticleChange;
00187 }
00188
00189
00190 void
00191 G4HEPionPlusInelastic::FirstIntInCasPionPlus(G4bool& inElastic,
00192 const G4double availableEnergy,
00193 G4HEVector pv[],
00194 G4int& vecLen,
00195 const G4HEVector& incidentParticle,
00196 const G4HEVector& targetParticle,
00197 const G4double atomicWeight)
00198
00199
00200
00201
00202
00203
00204
00205
00206 {
00207 static const G4double expxu = 82.;
00208 static const G4double expxl = -expxu;
00209
00210 static const G4double protb = 0.7;
00211 static const G4double neutb = 0.7;
00212 static const G4double c = 1.25;
00213
00214 static const G4int numMul = 1200;
00215 static const G4int numSec = 60;
00216
00217 G4int neutronCode = Neutron.getCode();
00218 G4int protonCode = Proton.getCode();
00219 G4double pionMass = PionPlus.getMass();
00220
00221 G4int targetCode = targetParticle.getCode();
00222 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00223
00224 static G4bool first = true;
00225 static G4double protmul[numMul], protnorm[numSec];
00226 static G4double neutmul[numMul], neutnorm[numSec];
00227
00228
00229
00230
00231 G4int i, counter, nt, npos, nneg, nzero;
00232
00233 if( first )
00234 {
00235 first = false;
00236 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
00237 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00238 counter = -1;
00239 for( npos=0; npos<(numSec/3); npos++ )
00240 {
00241 for( nneg=Imax(0,npos-2); nneg<=npos; nneg++ )
00242 {
00243 for( nzero=0; nzero<numSec/3; nzero++ )
00244 {
00245 if( ++counter < numMul )
00246 {
00247 nt = npos+nneg+nzero;
00248 if( (nt>0) && (nt<=numSec) )
00249 {
00250 protmul[counter] =
00251 pmltpc(npos,nneg,nzero,nt,protb,c) ;
00252 protnorm[nt-1] += protmul[counter];
00253 }
00254 }
00255 }
00256 }
00257 }
00258 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
00259 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00260 counter = -1;
00261 for( npos=0; npos<numSec/3; npos++ )
00262 {
00263 for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ )
00264 {
00265 for( nzero=0; nzero<numSec/3; nzero++ )
00266 {
00267 if( ++counter < numMul )
00268 {
00269 nt = npos+nneg+nzero;
00270 if( (nt>0) && (nt<=numSec) )
00271 {
00272 neutmul[counter] =
00273 pmltpc(npos,nneg,nzero,nt,neutb,c);
00274 neutnorm[nt-1] += neutmul[counter];
00275 }
00276 }
00277 }
00278 }
00279 }
00280 for( i=0; i<numSec; i++ )
00281 {
00282 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00283 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00284 }
00285 }
00286
00287
00288
00289
00290 pv[0] = incidentParticle;
00291 pv[1] = targetParticle;
00292 vecLen = 2;
00293
00294 if( !inElastic )
00295 {
00296 if( targetCode == neutronCode )
00297 {
00298 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
00299 G4int iplab = G4int( Amin( 9.0, incidentTotalMomentum*5. ) );
00300 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
00301 {
00302 pv[0] = PionZero;
00303 pv[1] = Proton;
00304 }
00305 }
00306 return;
00307 }
00308 else if (availableEnergy <= pionMass)
00309 return;
00310
00311
00312
00313 npos = 0, nneg = 0, nzero = 0;
00314 G4double eab = availableEnergy;
00315 G4int ieab = G4int( eab*5.0 );
00316
00317 G4double supp[] = {0., 0.2, 0.45, 0.55, 0.65, 0.75, 0.85, 0.90, 0.94, 0.98};
00318 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
00319 {
00320
00321
00322 G4double w0, wp, wm, wt, ran;
00323 if( targetCode == protonCode )
00324 {
00325 w0 = - sqr(1.+protb)/(2.*c*c);
00326 wp = w0 = std::exp(w0);
00327 if( G4UniformRand() < w0/(w0+wp) )
00328 { npos = 0; nneg = 0; nzero = 1; }
00329 else
00330 { npos = 1; nneg = 0; nzero = 0; }
00331 }
00332 else
00333 {
00334 w0 = -sqr(1.+neutb)/(2.*c*c);
00335 wp = w0 = std::exp(w0);
00336 wm = -sqr(-1.+neutb)/(2.*c*c);
00337 wm = std::exp(wm);
00338 wt = w0+wp+wm;
00339 wp = w0+wp;
00340 ran = G4UniformRand();
00341 if( ran < w0/wt)
00342 { npos = 0; nneg = 0; nzero = 1; }
00343 else if( ran < wp/wt)
00344 { npos = 1; nneg = 0; nzero = 0; }
00345 else
00346 { npos = 0; nneg = 1; nzero = 0; }
00347 }
00348 }
00349 else
00350 {
00351
00352
00353 G4double aleab = std::log(availableEnergy);
00354 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00355 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00356
00357
00358
00359 G4double test, dum, anpn = 0.0;
00360
00361 for (nt=1; nt<=numSec; nt++) {
00362 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00363 dum = pi*nt/(2.0*n*n);
00364 if (std::fabs(dum) < 1.0) {
00365 if( test >= 1.0e-10 )anpn += dum*test;
00366 } else {
00367 anpn += dum*test;
00368 }
00369 }
00370
00371 G4double ran = G4UniformRand();
00372 G4double excs = 0.0;
00373 if( targetCode == protonCode )
00374 {
00375 counter = -1;
00376 for (npos=0; npos<numSec/3; npos++) {
00377 for (nneg=Imax(0,npos-2); nneg<=npos; nneg++) {
00378 for (nzero=0; nzero<numSec/3; nzero++) {
00379 if (++counter < numMul) {
00380 nt = npos+nneg+nzero;
00381 if ( (nt>0) && (nt<=numSec) ) {
00382 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00383 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00384 if (std::fabs(dum) < 1.0) {
00385 if( test >= 1.0e-10 )excs += dum*test;
00386 } else {
00387 excs += dum*test;
00388 }
00389 if (ran < excs) goto outOfLoop;
00390 }
00391 }
00392 }
00393 }
00394 }
00395
00396
00397 inElastic = false;
00398 return;
00399 }
00400 else
00401 {
00402 counter = -1;
00403 for (npos=0; npos<numSec/3; npos++) {
00404 for (nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++) {
00405 for (nzero=0; nzero<numSec/3; nzero++) {
00406 if (++counter < numMul) {
00407 nt = npos+nneg+nzero;
00408 if ( (nt>=1) && (nt<=numSec) ) {
00409 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00410 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00411 if (std::fabs(dum) < 1.0) {
00412 if( test >= 1.0e-10 )excs += dum*test;
00413 } else {
00414 excs += dum*test;
00415 }
00416 if (ran < excs) goto outOfLoop;
00417 }
00418 }
00419 }
00420 }
00421 }
00422
00423 inElastic = false;
00424 return;
00425 }
00426 }
00427 outOfLoop:
00428
00429 if( targetCode == protonCode)
00430 {
00431 if( npos == nneg)
00432 {
00433 }
00434 else if (npos == (1+nneg))
00435 {
00436 if( G4UniformRand() < 0.5)
00437 {
00438 pv[1] = Neutron;
00439 }
00440 else
00441 {
00442 pv[0] = PionZero;
00443 }
00444 }
00445 else
00446 {
00447 pv[0] = PionZero;
00448 pv[1] = Neutron;
00449 }
00450 }
00451 else
00452 {
00453 if( npos == nneg)
00454 {
00455 if( G4UniformRand() < 0.25)
00456 {
00457 pv[0] = PionZero;
00458 pv[1] = Proton;
00459 }
00460 else
00461 {
00462 }
00463 }
00464 else if ( npos == (1+nneg))
00465 {
00466 pv[0] = PionZero;
00467 }
00468 else
00469 {
00470 pv[1] = Proton;
00471 }
00472 }
00473
00474
00475 nt = npos + nneg + nzero;
00476 while ( nt > 0)
00477 {
00478 G4double ran = G4UniformRand();
00479 if ( ran < (G4double)npos/nt)
00480 {
00481 if( npos > 0 )
00482 { pv[vecLen++] = PionPlus;
00483 npos--;
00484 }
00485 }
00486 else if ( ran < (G4double)(npos+nneg)/nt)
00487 {
00488 if( nneg > 0 )
00489 {
00490 pv[vecLen++] = PionMinus;
00491 nneg--;
00492 }
00493 }
00494 else
00495 {
00496 if( nzero > 0 )
00497 {
00498 pv[vecLen++] = PionZero;
00499 nzero--;
00500 }
00501 }
00502 nt = npos + nneg + nzero;
00503 }
00504 if (verboseLevel > 1)
00505 {
00506 G4cout << "Particles produced: " ;
00507 G4cout << pv[0].getName() << " " ;
00508 G4cout << pv[1].getName() << " " ;
00509 for (i=2; i < vecLen; i++)
00510 {
00511 G4cout << pv[i].getName() << " " ;
00512 }
00513 G4cout << G4endl;
00514 }
00515 return;
00516 }
00517
00518
00519
00520
00521
00522
00523
00524
00525