44 outFile <<
"G4HEAntiXiMinusInelastic is one of the High Energy\n"
45 <<
"Parameterized (HEP) models used to implement inelastic\n"
46 <<
"anti-Xi- 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-Xi- with initial energies\n"
76 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
78 if (incidentKineticEnergy < 1.)
79 G4cout <<
"GHEAntiXiMinusInelastic: incident energy < 1 GeV" <<
G4endl;
82 G4cout <<
"G4HEAntiXiMinusInelastic::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);
148 excitationEnergyGNP, excitationEnergyDTA,
149 incidentParticle, targetParticle,
150 atomicWeight, atomicNumber);
153 excitationEnergyGNP, excitationEnergyDTA,
154 incidentParticle, targetParticle,
155 atomicWeight, atomicNumber);
158 excitationEnergyGNP, excitationEnergyDTA,
159 incidentParticle, targetParticle,
160 atomicWeight, atomicNumber);
164 excitationEnergyGNP, excitationEnergyDTA,
165 incidentParticle, targetParticle,
166 atomicWeight, atomicNumber);
169 excitationEnergyGNP, excitationEnergyDTA,
170 incidentParticle, targetParticle,
171 atomicWeight, atomicNumber);
175 atomicWeight, atomicNumber);
178 G4cout <<
"GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
202 static const G4double expxl = -expxu;
208 static const G4int numMul = 1200;
209 static const G4int numMulAn = 400;
210 static const G4int numSec = 60;
217 static G4bool first =
true;
218 static G4double protmul[numMul], protnorm[numSec];
219 static G4double protmulAn[numMulAn],protnormAn[numSec];
220 static G4double neutmul[numMul], neutnorm[numSec];
221 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
226 G4int i, counter, nt, npos, nneg, nzero;
231 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
232 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
234 for( npos=0; npos<(numSec/3); npos++ )
236 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
238 for( nzero=0; nzero<numSec/3; nzero++ )
240 if( ++counter < numMul )
242 nt = npos+nneg+nzero;
243 if( (nt>0) && (nt<=numSec) )
245 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
246 protnorm[nt-1] += protmul[counter];
252 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
253 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
255 for( npos=0; npos<numSec/3; npos++ )
257 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
259 for( nzero=0; nzero<numSec/3; nzero++ )
261 if( ++counter < numMul )
263 nt = npos+nneg+nzero;
264 if( (nt>0) && (nt<=numSec) )
266 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
267 neutnorm[nt-1] += neutmul[counter];
273 for( i=0; i<numSec; i++ )
275 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
276 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
279 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
280 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
282 for( npos=1; npos<(numSec/3); npos++ )
284 nneg = std::max(0,npos-1);
285 for( nzero=0; nzero<numSec/3; nzero++ )
287 if( ++counter < numMulAn )
289 nt = npos+nneg+nzero;
290 if( (nt>1) && (nt<=numSec) )
292 protmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
293 protnormAn[nt-1] += protmulAn[counter];
298 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
299 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
301 for( npos=0; npos<numSec/3; npos++ )
304 for( nzero=0; nzero<numSec/3; nzero++ )
306 if( ++counter < numMulAn )
308 nt = npos+nneg+nzero;
309 if( (nt>1) && (nt<=numSec) )
311 neutmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
312 neutnormAn[nt-1] += neutmulAn[counter];
317 for( i=0; i<numSec; i++ )
319 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
320 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
327 pv[0] = incidentParticle;
328 pv[1] = targetParticle;
333 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
335 G4int iplab = std::min(9,
G4int( incidentTotalMomentum*2.5));
336 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
340 if ( targetCode == protonCode)
402 npos = 0; nneg = 0; nzero = 0;
403 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
404 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
405 0.39, 0.36, 0.33, 0.10, 0.01};
406 G4int iplab =
G4int( incidentTotalMomentum*10.);
407 if ( iplab > 9) iplab = 10 +
G4int( (incidentTotalMomentum -1.)*5. );
408 if ( iplab > 14) iplab = 15 +
G4int( incidentTotalMomentum -2. );
409 if ( iplab > 22) iplab = 23 +
G4int( (incidentTotalMomentum -10.)/10.);
410 iplab = std::min(24, iplab);
417 G4double aleab = std::log(availableEnergy);
418 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
419 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
425 for (nt=1; nt<=numSec; nt++) {
426 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
427 dum = pi*nt/(2.0*n*n);
428 if (std::fabs(dum) < 1.0) {
429 if( test >= 1.0e-10 )anpn += dum*test;
437 if( targetCode == protonCode )
440 for( npos=0; npos<numSec/3; npos++ )
442 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
444 for( nzero=0; nzero<numSec/3; nzero++ )
446 if( ++counter < numMul )
448 nt = npos+nneg+nzero;
449 if ( (nt>0) && (nt<=numSec) ) {
450 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
451 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;
472 for( npos=0; npos<numSec/3; npos++ )
474 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
476 for( nzero=0; nzero<numSec/3; nzero++ )
478 if( ++counter < numMul )
480 nt = npos+nneg+nzero;
481 if ( (nt>0) && (nt<=numSec) ) {
482 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
483 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
484 if (std::fabs(dum) < 1.0) {
485 if( test >= 1.0e-10 )excs += dum*test;
490 if (ran < excs)
goto outOfLoop;
505 if( targetCode == protonCode)
522 else if (npos == (nneg+1))
538 else if (npos == (nneg-1))
565 else if ( npos == (nneg-1))
581 else if (npos == (nneg+1))
598 G4double aleab = std::log(availableEnergy);
599 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
600 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
606 for (nt=2; nt<=numSec; nt++) {
607 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
608 dum = pi*nt/(2.0*n*n);
609 if (std::fabs(dum) < 1.0) {
610 if( test >= 1.0e-10 )anpn += dum*test;
618 if( targetCode == protonCode )
621 for( npos=1; npos<numSec/3; npos++ )
624 for( nzero=0; nzero<numSec/3; nzero++ )
626 if( ++counter < numMulAn )
628 nt = npos+nneg+nzero;
629 if ( (nt>1) && (nt<=numSec) ) {
630 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
631 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
632 if (std::fabs(dum) < 1.0) {
633 if( test >= 1.0e-10 )excs += dum*test;
638 if (ran < excs)
goto outOfLoopAn;
650 for( npos=0; npos<numSec/3; npos++ )
653 for( nzero=0; nzero<numSec/3; nzero++ )
655 if( ++counter < numMulAn )
657 nt = npos+nneg+nzero;
658 if ( (nt>1) && (nt<=numSec) ) {
659 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
660 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
661 if (std::fabs(dum) < 1.0) {
662 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 FirstIntInCasAntiXiMinus(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