44 outFile <<
"G4HEAntiXiZeroInelastic is one of the High Energy\n"
45 <<
"Parameterized (HEP) models used to implement inelastic\n"
46 <<
"anti-Xi0 scattering from nuclei. It is a re-engineered\n"
47 <<
"version of the GHEISHA code of H. Fesefeldt. It divides the\n"
48 <<
"initial collision products into backward- and forward-going\n"
49 <<
"clusters which are then decayed into final state hadrons.\n"
50 <<
"The model does not conserve energy on an event-by-event\n"
51 <<
"basis. It may be applied to anti-Xi0 with initial energies\n"
76 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
78 if (incidentKineticEnergy < 1.)
79 G4cout <<
"GHEAntiXiZeroInelastic: incident energy < 1 GeV" <<
G4endl;
82 G4cout <<
"G4HEAntiXiZeroInelastic::ApplyYourself" <<
G4endl;
84 <<
"mass " << incidentMass
85 <<
"kinetic energy " << incidentKineticEnergy
87 G4cout <<
"target material with (A,Z) = ("
88 << atomicWeight <<
"," << atomicNumber <<
")" <<
G4endl;
92 atomicWeight, atomicNumber);
94 G4cout <<
"nuclear inelasticity = " << inelasticity <<
G4endl;
96 incidentKineticEnergy -= inelasticity;
102 atomicWeight, atomicNumber,
104 excitationEnergyDTA);
106 G4cout <<
"nuclear excitation = " << excitation << excitationEnergyGNP
107 << excitationEnergyDTA <<
G4endl;
109 incidentKineticEnergy -= excitation;
110 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
123 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
124 + targetMass*targetMass
125 + 2.0*targetMass*incidentTotalEnergy);
126 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
132 G4cout <<
"ApplyYourself: CallFirstIntInCascade for particle "
133 << incidentCode <<
G4endl;
135 G4bool successful =
false;
138 incidentParticle, targetParticle, atomicWeight);
141 G4cout <<
"ApplyYourself::StrangeParticlePairProduction" <<
G4endl;
143 if ((
vecLength > 0) && (availableEnergy > 1.))
146 incidentParticle, targetParticle);
149 excitationEnergyGNP, excitationEnergyDTA,
150 incidentParticle, targetParticle,
151 atomicWeight, atomicNumber);
154 excitationEnergyGNP, excitationEnergyDTA,
155 incidentParticle, targetParticle,
156 atomicWeight, atomicNumber);
159 excitationEnergyGNP, excitationEnergyDTA,
160 incidentParticle, targetParticle,
161 atomicWeight, atomicNumber);
165 excitationEnergyGNP, excitationEnergyDTA,
166 incidentParticle, targetParticle,
167 atomicWeight, atomicNumber);
170 excitationEnergyGNP, excitationEnergyDTA,
171 incidentParticle, targetParticle,
172 atomicWeight, atomicNumber);
176 atomicWeight, atomicNumber);
178 G4cout <<
"GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
203 static const G4double expxl = -expxu;
209 static const G4int numMul = 1200;
210 static const G4int numMulAn = 400;
211 static const G4int numSec = 60;
218 static G4bool first =
true;
219 static G4double protmul[numMul], protnorm[numSec];
220 static G4double protmulAn[numMulAn],protnormAn[numSec];
221 static G4double neutmul[numMul], neutnorm[numSec];
222 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
227 G4int i, counter, nt, npos, nneg, nzero;
232 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
233 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
235 for( npos=0; npos<(numSec/3); npos++ )
237 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
239 for( nzero=0; nzero<numSec/3; nzero++ )
241 if( ++counter < numMul )
243 nt = npos+nneg+nzero;
244 if( (nt>0) && (nt<=numSec) )
246 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
247 protnorm[nt-1] += protmul[counter];
253 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
254 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
256 for( npos=0; npos<numSec/3; npos++ )
258 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
260 for( nzero=0; nzero<numSec/3; nzero++ )
262 if( ++counter < numMul )
264 nt = npos+nneg+nzero;
265 if( (nt>0) && (nt<=numSec) )
267 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
268 neutnorm[nt-1] += neutmul[counter];
274 for( i=0; i<numSec; i++ )
276 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
277 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
280 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
281 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
283 for( npos=1; npos<(numSec/3); npos++ )
285 nneg = std::max(0,npos-1);
286 for( nzero=0; nzero<numSec/3; nzero++ )
288 if( ++counter < numMulAn )
290 nt = npos+nneg+nzero;
291 if( (nt>1) && (nt<=numSec) )
293 protmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
294 protnormAn[nt-1] += protmulAn[counter];
299 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
300 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
302 for( npos=0; npos<numSec/3; npos++ )
305 for( nzero=0; nzero<numSec/3; nzero++ )
307 if( ++counter < numMulAn )
309 nt = npos+nneg+nzero;
310 if( (nt>1) && (nt<=numSec) )
312 neutmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
313 neutnormAn[nt-1] += neutmulAn[counter];
318 for( i=0; i<numSec; i++ )
320 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
321 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
328 pv[0] = incidentParticle;
329 pv[1] = targetParticle;
334 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
336 G4int iplab = std::min(9,
G4int( incidentTotalMomentum*2.5 ));
337 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
341 if ( targetCode == protonCode)
403 npos = 0; nneg = 0; nzero = 0;
404 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
405 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
406 0.39, 0.36, 0.33, 0.10, 0.01};
407 G4int iplab =
G4int( incidentTotalMomentum*10.);
408 if ( iplab > 9) iplab = 10 +
G4int( (incidentTotalMomentum -1.)*5. );
409 if ( iplab > 14) iplab = 15 +
G4int( incidentTotalMomentum -2. );
410 if ( iplab > 22) iplab = 23 +
G4int( (incidentTotalMomentum -10.)/10.);
411 iplab = std::min(24, iplab);
418 G4double aleab = std::log(availableEnergy);
419 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
420 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
426 for (nt=1; nt<=numSec; nt++) {
427 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
428 dum = pi*nt/(2.0*n*n);
429 if (std::fabs(dum) < 1.0) {
430 if( test >= 1.0e-10 )anpn += dum*test;
438 if( targetCode == protonCode )
441 for( npos=0; npos<numSec/3; npos++ )
443 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
445 for( nzero=0; nzero<numSec/3; nzero++ )
447 if( ++counter < numMul )
449 nt = npos+nneg+nzero;
450 if ( (nt>0) && (nt<=numSec) ) {
451 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
452 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
453 if (std::fabs(dum) < 1.0) {
454 if( test >= 1.0e-10 )excs += dum*test;
459 if (ran < excs)
goto outOfLoop;
473 for( npos=0; npos<numSec/3; npos++ )
475 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
477 for( nzero=0; nzero<numSec/3; nzero++ )
479 if( ++counter < numMul )
481 nt = npos+nneg+nzero;
482 if ( (nt>0) && (nt<=numSec) ) {
483 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
484 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
485 if (std::fabs(dum) < 1.0) {
486 if( test >= 1.0e-10 )excs += dum*test;
491 if (ran < excs)
goto outOfLoop;
506 if( targetCode == protonCode)
523 else if (npos == (nneg+1))
539 else if (npos == (nneg-1))
566 else if ( npos == (nneg-1))
582 else if (npos == (nneg+1))
599 G4double aleab = std::log(availableEnergy);
600 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
601 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
607 for (nt=2; nt<=numSec; nt++) {
608 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
609 dum = pi*nt/(2.0*n*n);
610 if (std::fabs(dum) < 1.0) {
611 if( test >= 1.0e-10 )anpn += dum*test;
619 if( targetCode == protonCode )
622 for( npos=1; npos<numSec/3; npos++ )
625 for( nzero=0; nzero<numSec/3; nzero++ )
627 if( ++counter < numMulAn )
629 nt = npos+nneg+nzero;
630 if ( (nt>1) && (nt<=numSec) ) {
631 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
632 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
633 if (std::fabs(dum) < 1.0) {
634 if( test >= 1.0e-10 )excs += dum*test;
639 if (ran < excs)
goto outOfLoopAn;
651 for( npos=0; npos<numSec/3; npos++ )
654 for( nzero=0; nzero<numSec/3; nzero++ )
656 if( ++counter < numMulAn )
658 nt = npos+nneg+nzero;
659 if ( (nt>1) && (nt<=numSec) ) {
660 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
661 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
662 if (std::fabs(dum) < 1.0) {
663 if( test >= 1.0e-10 )excs += dum*test;
667 if (ran < excs)
goto outOfLoopAn;
680 nt = npos + nneg + nzero;
691 else if ( ran < (
G4double)(npos+nneg)/nt)
707 nt = npos + nneg + nzero;
711 G4cout <<
"Particles produced: " ;
714 for (i=2; i < vecLen; i++)
G4DLLIMPORT std::ostream G4cout
void FirstIntInCasAntiXiZero(G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)
virtual void ModelDescription(std::ostream &) const
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)
G4HEVector AntiSigmaMinus
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