45 outFile <<
"G4HEAntiOmegaMinusInelastic is one of the High Energy\n"
46 <<
"Parameterized (HEP) models used to implement inelastic\n"
47 <<
"anti-Omega- scattering from nuclei. It is a re-engineered\n"
48 <<
"version of the GHEISHA code of H. Fesefeldt. It divides the\n"
49 <<
"initial collision products into backward- and forward-going\n"
50 <<
"clusters which are then decayed into final state hadrons.\n"
51 <<
"The model does not conserve energy on an event-by-event\n"
52 <<
"basis. It may be applied to anti-Omega- with initial\n"
53 <<
"energies above 20 GeV.\n";
74 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
76 if (incidentKineticEnergy < 1.)
77 G4cout <<
"GHEAntiOmegaMinusInelastic: incident energy < 1 GeV" <<
G4endl;
80 G4cout <<
"G4HEAntiOmegaMinusInelastic::ApplyYourself" <<
G4endl;
82 <<
"mass " << incidentMass
83 <<
"kinetic energy " << incidentKineticEnergy
85 G4cout <<
"target material with (A,Z) = ("
86 << atomicWeight <<
"," << atomicNumber <<
")" <<
G4endl;
90 atomicWeight, atomicNumber);
92 G4cout <<
"nuclear inelasticity = " << inelasticity <<
G4endl;
94 incidentKineticEnergy -= inelasticity;
100 atomicWeight, atomicNumber,
102 excitationEnergyDTA);
104 G4cout <<
"nuclear excitation = " << excitation << excitationEnergyGNP
105 << excitationEnergyDTA <<
G4endl;
108 incidentKineticEnergy -= excitation;
109 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
122 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
123 + targetMass*targetMass
124 + 2.0*targetMass*incidentTotalEnergy);
125 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
131 G4cout <<
"ApplyYourself: CallFirstIntInCascade for particle "
132 << incidentCode <<
G4endl;
134 G4bool successful =
false;
137 incidentParticle, targetParticle, atomicWeight);
140 G4cout <<
"ApplyYourself::StrangeParticlePairProduction" <<
G4endl;
142 if ((
vecLength > 0) && (availableEnergy > 1.))
145 incidentParticle, targetParticle);
147 excitationEnergyGNP, excitationEnergyDTA,
148 incidentParticle, targetParticle,
149 atomicWeight, atomicNumber);
152 excitationEnergyGNP, excitationEnergyDTA,
153 incidentParticle, targetParticle,
154 atomicWeight, atomicNumber);
157 excitationEnergyGNP, excitationEnergyDTA,
158 incidentParticle, targetParticle,
159 atomicWeight, atomicNumber);
163 excitationEnergyGNP, excitationEnergyDTA,
164 incidentParticle, targetParticle,
165 atomicWeight, atomicNumber);
168 excitationEnergyGNP, excitationEnergyDTA,
169 incidentParticle, targetParticle,
170 atomicWeight, atomicNumber);
174 atomicWeight, atomicNumber);
177 G4cout <<
"GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
201 static const G4double expxl = -expxu;
207 static const G4int numMul = 1200;
208 static const G4int numMulAn = 400;
209 static const G4int numSec = 60;
216 static G4bool first =
true;
217 static G4double protmul[numMul], protnorm[numSec];
218 static G4double protmulAn[numMulAn],protnormAn[numSec];
219 static G4double neutmul[numMul], neutnorm[numSec];
220 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
225 G4int i, counter, nt, npos, nneg, nzero;
230 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
231 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
233 for( npos=0; npos<(numSec/3); npos++ )
235 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
237 for( nzero=0; nzero<numSec/3; nzero++ )
239 if( ++counter < numMul )
241 nt = npos+nneg+nzero;
242 if( (nt>0) && (nt<=numSec) )
244 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
245 protnorm[nt-1] += protmul[counter];
251 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
252 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
254 for( npos=0; npos<numSec/3; npos++ )
256 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
258 for( nzero=0; nzero<numSec/3; nzero++ )
260 if( ++counter < numMul )
262 nt = npos+nneg+nzero;
263 if( (nt>0) && (nt<=numSec) )
265 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
266 neutnorm[nt-1] += neutmul[counter];
272 for( i=0; i<numSec; i++ )
274 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
275 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
278 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
279 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
281 for( npos=1; npos<(numSec/3); npos++ )
283 nneg = std::max(0,npos-1);
284 for( nzero=0; nzero<numSec/3; nzero++ )
286 if( ++counter < numMulAn )
288 nt = npos+nneg+nzero;
289 if( (nt>1) && (nt<=numSec) )
291 protmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
292 protnormAn[nt-1] += protmulAn[counter];
297 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
298 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
300 for( npos=0; npos<numSec/3; npos++ )
303 for( nzero=0; nzero<numSec/3; nzero++ )
305 if( ++counter < numMulAn )
307 nt = npos+nneg+nzero;
308 if( (nt>1) && (nt<=numSec) )
310 neutmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
311 neutnormAn[nt-1] += neutmulAn[counter];
316 for( i=0; i<numSec; i++ )
318 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
319 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
326 pv[0] = incidentParticle;
327 pv[1] = targetParticle;
332 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
334 G4int iplab = std::min(9,
G4int( incidentTotalMomentum*2.5 ));
335 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
339 if ( targetCode == protonCode)
401 npos = 0; nneg = 0; nzero = 0;
402 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
403 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
404 0.39, 0.36, 0.33, 0.10, 0.01};
405 G4int iplab =
G4int( incidentTotalMomentum*10.);
406 if ( iplab > 9) iplab = 10 +
G4int( (incidentTotalMomentum -1.)*5. );
407 if ( iplab > 14) iplab = 15 +
G4int( incidentTotalMomentum -2. );
408 if ( iplab > 22) iplab = 23 +
G4int( (incidentTotalMomentum -10.)/10.);
409 iplab = std::min(24, iplab);
414 G4double aleab = std::log(availableEnergy);
415 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
416 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
422 for (nt=1; nt<=numSec; nt++) {
423 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
424 dum = pi*nt/(2.0*n*n);
425 if (std::fabs(dum) < 1.0) {
426 if( test >= 1.0e-10 )anpn += dum*test;
434 if( targetCode == protonCode )
437 for (npos=0; npos<numSec/3; npos++) {
438 for (nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++) {
439 for (nzero=0; nzero<numSec/3; nzero++) {
440 if (++counter < numMul) {
441 nt = npos+nneg+nzero;
442 if ( (nt>0) && (nt<=numSec) ) {
443 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
444 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
445 if (std::fabs(dum) < 1.0) {
446 if( test >= 1.0e-10 )excs += dum*test;
451 if (ran < excs)
goto outOfLoop;
465 for (npos=0; npos<numSec/3; npos++) {
466 for (nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++) {
467 for (nzero=0; nzero<numSec/3; nzero++) {
468 if (++counter < numMul) {
469 nt = npos+nneg+nzero;
470 if ( (nt>0) && (nt<=numSec) ) {
471 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
472 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
473 if (std::fabs(dum) < 1.0) {
474 if( test >= 1.0e-10 )excs += dum*test;
479 if (ran < excs)
goto outOfLoop;
494 if( targetCode == protonCode)
511 else if (npos == (nneg+1))
527 else if (npos == (nneg-1))
554 else if ( npos == (nneg-1))
570 else if (npos == (nneg+1))
587 G4double aleab = std::log(availableEnergy);
588 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
589 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
595 for (nt=2; nt<=numSec; nt++) {
596 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
597 dum = pi*nt/(2.0*n*n);
598 if (std::fabs(dum) < 1.0) {
599 if( test >= 1.0e-10 )anpn += dum*test;
607 if( targetCode == protonCode )
610 for( npos=1; npos<numSec/3; npos++ )
613 for( nzero=0; nzero<numSec/3; nzero++ )
615 if( ++counter < numMulAn )
617 nt = npos+nneg+nzero;
618 if ( (nt>1) && (nt<=numSec) ) {
619 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
620 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
621 if (std::fabs(dum) < 1.0) {
622 if( test >= 1.0e-10 )excs += dum*test;
627 if (ran < excs)
goto outOfLoopAn;
639 for( npos=0; npos<numSec/3; npos++ )
642 for( nzero=0; nzero<numSec/3; nzero++ )
644 if( ++counter < numMulAn )
646 nt = npos+nneg+nzero;
647 if ( (nt>1) && (nt<=numSec) ) {
648 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
649 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
650 if (std::fabs(dum) < 1.0) {
651 if( test >= 1.0e-10 )excs += dum*test;
656 if (ran < excs)
goto outOfLoopAn;
669 nt = npos + nneg + nzero;
680 else if ( ran < (
G4double)(npos+nneg)/nt)
696 nt = npos + nneg + nzero;
700 G4cout <<
"Particles produced: " ;
703 for (i=2; i < vecLen; i++)
G4DLLIMPORT std::ostream G4cout
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void FirstIntInCasAntiOmegaMinus(G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)
virtual void ModelDescription(std::ostream &) const
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