51 G4cout <<
"WARNING: model G4HEAntiNeutronInelastic is being deprecated and will\n"
52 <<
"disappear in Geant4 version 10.0" <<
G4endl;
58 outFile <<
"G4HEAntiNeutronInelastic is one of the High Energy\n"
59 <<
"Parameterized (HEP) models used to implement inelastic\n"
60 <<
"anti-neutron scattering from nuclei. It is a re-engineered\n"
61 <<
"version of the GHEISHA code of H. Fesefeldt. It divides the\n"
62 <<
"initial collision products into backward- and forward-going\n"
63 <<
"clusters which are then decayed into final state hadrons.\n"
64 <<
"The model does not conserve energy on an event-by-event\n"
65 <<
"basis. It may be applied to anti-neutrons with initial energies\n"
87 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
89 if (incidentKineticEnergy < 1.)
90 G4cout <<
"GHEAntiNeutronInelastic: incident energy < 1 GeV" <<
G4endl;;
93 G4cout <<
"G4HEAntiNeutronInelastic::ApplyYourself" <<
G4endl;
95 <<
"mass " << incidentMass
96 <<
"kinetic energy " << incidentKineticEnergy
98 G4cout <<
"target material with (A,Z) = ("
99 << atomicWeight <<
"," << atomicNumber <<
")" <<
G4endl;
103 atomicWeight, atomicNumber);
105 G4cout <<
"nuclear inelasticity = " << inelasticity <<
G4endl;
107 incidentKineticEnergy -= inelasticity;
113 atomicWeight, atomicNumber,
115 excitationEnergyDTA);
117 G4cout <<
"nuclear excitation = " << excitation << excitationEnergyGNP
118 << excitationEnergyDTA <<
G4endl;
120 incidentKineticEnergy -= excitation;
121 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
134 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
135 + targetMass*targetMass
136 + 2.0*targetMass*incidentTotalEnergy);
137 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
143 G4cout <<
"ApplyYourself: CallFirstIntInCascade for particle "
144 << incidentCode <<
G4endl;
146 G4bool successful =
false;
149 incidentParticle, targetParticle, atomicWeight);
152 G4cout <<
"ApplyYourself::StrangeParticlePairProduction" <<
G4endl;
154 if ((
vecLength > 0) && (availableEnergy > 1.))
157 incidentParticle, targetParticle);
159 excitationEnergyGNP, excitationEnergyDTA,
160 incidentParticle, targetParticle,
161 atomicWeight, atomicNumber);
164 excitationEnergyGNP, excitationEnergyDTA,
165 incidentParticle, targetParticle,
166 atomicWeight, atomicNumber);
169 excitationEnergyGNP, excitationEnergyDTA,
170 incidentParticle, targetParticle,
171 atomicWeight, atomicNumber);
175 excitationEnergyGNP, excitationEnergyDTA,
176 incidentParticle, targetParticle,
177 atomicWeight, atomicNumber);
180 excitationEnergyGNP, excitationEnergyDTA,
181 incidentParticle, targetParticle,
182 atomicWeight, atomicNumber);
186 atomicWeight, atomicNumber);
189 G4cout <<
"GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
218 static const G4double expxl = -expxu;
224 static const G4int numMul = 1200;
225 static const G4int numMulAn = 400;
226 static const G4int numSec = 60;
234 static G4bool first =
true;
235 static G4double protmul[numMul], protnorm[numSec];
236 static G4double protmulAn[numMulAn],protnormAn[numSec];
237 static G4double neutmul[numMul], neutnorm[numSec];
238 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
243 G4int i, counter, nt, npos, nneg, nzero;
248 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
249 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
251 for( npos=0; npos<(numSec/3); npos++ )
253 for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ )
255 for( nzero=0; nzero<numSec/3; nzero++ )
257 if( ++counter < numMul )
259 nt = npos+nneg+nzero;
260 if( (nt>0) && (nt<=numSec) )
262 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
263 protnorm[nt-1] += protmul[counter];
269 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
270 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
272 for( npos=0; npos<numSec/3; npos++ )
274 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
276 for( nzero=0; nzero<numSec/3; nzero++ )
278 if( ++counter < numMul )
280 nt = npos+nneg+nzero;
281 if( (nt>0) && (nt<=numSec) )
283 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
284 neutnorm[nt-1] += neutmul[counter];
290 for( i=0; i<numSec; i++ )
292 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
293 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
296 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
297 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
299 for( npos=1; npos<(numSec/3); npos++ )
301 nneg = std::max(0,npos-1);
302 for( nzero=0; nzero<numSec/3; nzero++ )
304 if( ++counter < numMulAn )
306 nt = npos+nneg+nzero;
307 if( (nt>1) && (nt<=numSec) )
309 protmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
310 protnormAn[nt-1] += protmulAn[counter];
315 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
316 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
318 for( npos=0; npos<numSec/3; npos++ )
321 for( nzero=0; nzero<numSec/3; nzero++ )
323 if( ++counter < numMulAn )
325 nt = npos+nneg+nzero;
326 if( (nt>1) && (nt<=numSec) )
328 neutmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
329 neutnormAn[nt-1] += neutmulAn[counter];
334 for( i=0; i<numSec; i++ )
336 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
337 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
344 pv[0] = incidentParticle;
345 pv[1] = targetParticle;
352 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
354 G4int iplab = std::min(9,
G4int( incidentTotalMomentum*2.5));
355 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
368 npos = 0, nneg = 0, nzero = 0;
369 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
370 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
371 0.39, 0.36, 0.33, 0.10, 0.01};
372 G4int iplab =
G4int( incidentTotalMomentum*10.);
373 if ( iplab > 9) iplab = 10 +
G4int( (incidentTotalMomentum -1.)*5. );
374 if ( iplab > 14) iplab = 15 +
G4int( incidentTotalMomentum -2. );
375 if ( iplab > 22) iplab = 23 +
G4int( (incidentTotalMomentum -10.)/10.);
376 iplab = std::min(24, iplab);
384 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
390 if( targetCode == protonCode )
392 w0 = -
sqr(1.+protb)/(2.*c*c);
393 w0 = wp = std::exp(w0);
395 { npos = 0; nneg = 0; nzero = 1; }
397 { npos = 1; nneg = 0; nzero = 0; }
401 w0 = -
sqr(1.+neutb)/(2.*c*c);
402 w0 = wp = std::exp(w0);
403 wm = -
sqr(-1.+neutb)/(2.*c*c);
409 { npos = 0; nneg = 0; nzero = 1; }
410 else if( ran < wp/wt)
411 { npos = 1; nneg = 0; nzero = 0; }
413 { npos = 0; nneg = 1; nzero = 0; }
420 G4double aleab = std::log(availableEnergy);
421 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
422 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
428 for (nt=1; nt<=numSec; nt++) {
429 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
430 dum = pi*nt/(2.0*n*n);
432 if (std::fabs(dum) < 1.0) {
433 if( test >= 1.0e-10 )anpn += dum*test;
441 if (targetCode == protonCode) {
443 for (npos=0; npos<numSec/3; npos++) {
444 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
445 for (nzero=0; nzero<numSec/3; nzero++) {
446 if (++counter < numMul) {
447 nt = npos+nneg+nzero;
448 if ( (nt>0) && (nt<=numSec) ) {
449 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
450 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
452 if (std::fabs(dum) < 1.0) {
453 if( test >= 1.0e-10 )excs += dum*test;
458 if (ran < excs)
goto outOfLoop;
471 for (npos=0; npos<numSec/3; npos++) {
472 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
473 for (nzero=0; nzero<numSec/3; nzero++) {
474 if (++counter < numMul) {
475 nt = npos+nneg+nzero;
476 if ((nt>0) && (nt<=numSec) ) {
477 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
478 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
479 if (std::fabs(dum) < 1.0) {
480 if( test >= 1.0e-10 )excs += dum*test;
485 if (ran < excs)
goto outOfLoop;
498 if( targetCode == protonCode)
503 else if (npos == (nneg+1))
533 else if ( npos == (nneg-1))
549 G4double aleab = std::log(availableEnergy);
550 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
551 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
557 for (nt=2; nt<=numSec; nt++) {
558 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
559 dum = pi*nt/(2.0*n*n);
560 if (std::fabs(dum) < 1.0) {
561 if( test >= 1.0e-10 )anpn += dum*test;
569 if (targetCode == protonCode) {
571 for (npos=1; npos<numSec/3; npos++) {
573 for (nzero=0; nzero<numSec/3; nzero++) {
574 if (++counter < numMulAn) {
575 nt = npos+nneg+nzero;
576 if ( (nt>1) && (nt<=numSec) ) {
577 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
578 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
580 if (std::fabs(dum) < 1.0) {
581 if( test >= 1.0e-10 )excs += dum*test;
586 if (ran < excs)
goto outOfLoopAn;
597 for (npos=0; npos<numSec/3; npos++) {
599 for (nzero=0; nzero<numSec/3; nzero++) {
600 if (++counter < numMulAn) {
601 nt = npos+nneg+nzero;
602 if ( (nt>1) && (nt<=numSec) ) {
603 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
604 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
606 if (std::fabs(dum) < 1.0) {
607 if( test >= 1.0e-10 )excs += dum*test;
612 if (ran < excs)
goto outOfLoopAn;
625 nt = npos + nneg + nzero;
636 else if ( ran < (
G4double)(npos+nneg)/nt)
652 nt = npos + nneg + nzero;
656 G4cout <<
"Particles produced: " ;
659 for (i=2; i < vecLen; i++)
G4DLLIMPORT std::ostream G4cout
virtual void ModelDescription(std::ostream &) const
G4HEAntiNeutronInelastic(const G4String &name="G4HEAntiNeutronInelastic")
void FirstIntInCasAntiNeutron(G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
void MediumEnergyClusterProduction(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
void ElasticScattering(G4bool &successful, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, G4double atomicWeight, G4double atomicNumber)
void QuasiElasticScattering(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
void FillParticleChange(G4HEVector pv[], G4int aVecLength)
void HighEnergyClusterProduction(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
G4double NuclearExcitation(G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber, G4double &excitationEnergyCascade, G4double &excitationEnergyEvaporation)
void MediumEnergyCascading(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
G4double NuclearInelasticity(G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber)
void StrangeParticlePairProduction(const G4double availableEnergy, const G4double centerOfMassEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle)
void HighEnergyCascading(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
G4double getEnergy() const
G4double getTotalMomentum() const
void setDefinition(G4String name)
void SetStatusChange(G4HadFinalStateStatus aS)
G4HadFinalState theParticleChange