45 outFile <<
"G4HEAntiSigmaMinusInelastic is one of the High Energy\n"
46 <<
"Parameterized (HEP) models used to implement inelastic\n"
47 <<
"anti-Sigma- 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-Sigma- with initial\n"
53 <<
"energies above 20 GeV.\n";
74 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
76 if (incidentKineticEnergy < 1.)
77 G4cout <<
"GHEAntiSigmaMinusInelastic: incident energy < 1 GeV" <<
G4endl;
80 G4cout <<
"G4HEAntiSigmaMinusInelastic::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;
107 incidentKineticEnergy -= excitation;
108 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
121 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
122 + targetMass*targetMass
123 + 2.0*targetMass*incidentTotalEnergy);
124 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
130 G4cout <<
"ApplyYourself: CallFirstIntInCascade for particle "
131 << incidentCode <<
G4endl;
133 G4bool successful =
false;
136 incidentParticle, targetParticle, atomicWeight);
139 G4cout <<
"ApplyYourself::StrangeParticlePairProduction" <<
G4endl;
141 if ((
vecLength > 0) && (availableEnergy > 1.))
144 incidentParticle, targetParticle);
146 excitationEnergyGNP, excitationEnergyDTA,
147 incidentParticle, targetParticle,
148 atomicWeight, atomicNumber);
151 excitationEnergyGNP, excitationEnergyDTA,
152 incidentParticle, targetParticle,
153 atomicWeight, atomicNumber);
156 excitationEnergyGNP, excitationEnergyDTA,
157 incidentParticle, targetParticle,
158 atomicWeight, atomicNumber);
162 excitationEnergyGNP, excitationEnergyDTA,
163 incidentParticle, targetParticle,
164 atomicWeight, atomicNumber);
167 excitationEnergyGNP, excitationEnergyDTA,
168 incidentParticle, targetParticle,
169 atomicWeight, atomicNumber);
173 atomicWeight, atomicNumber);
176 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;
219 static G4bool first =
true;
220 static G4double protmul[numMul], protnorm[numSec];
221 static G4double protmulAn[numMulAn],protnormAn[numSec];
222 static G4double neutmul[numMul], neutnorm[numSec];
223 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
228 G4int i, counter, nt, npos, nneg, nzero;
233 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
234 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
236 for( npos=0; npos<(numSec/3); npos++ )
238 for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ )
240 for( nzero=0; nzero<numSec/3; nzero++ )
242 if( ++counter < numMul )
244 nt = npos+nneg+nzero;
245 if( (nt>0) && (nt<=numSec) )
247 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
248 protnorm[nt-1] += protmul[counter];
254 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
255 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
257 for( npos=0; npos<numSec/3; npos++ )
259 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
261 for( nzero=0; nzero<numSec/3; nzero++ )
263 if( ++counter < numMul )
265 nt = npos+nneg+nzero;
266 if( (nt>0) && (nt<=numSec) )
268 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
269 neutnorm[nt-1] += neutmul[counter];
275 for( i=0; i<numSec; i++ )
277 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
278 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
281 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
282 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
284 for( npos=1; npos<(numSec/3); npos++ )
286 nneg = std::max(0,npos-2);
287 for( nzero=0; nzero<numSec/3; nzero++ )
289 if( ++counter < numMulAn )
291 nt = npos+nneg+nzero;
292 if( (nt>1) && (nt<=numSec) )
294 protmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
295 protnormAn[nt-1] += protmulAn[counter];
300 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
301 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
303 for( npos=0; npos<numSec/3; npos++ )
306 for( nzero=0; nzero<numSec/3; nzero++ )
308 if( ++counter < numMulAn )
310 nt = npos+nneg+nzero;
311 if( (nt>1) && (nt<=numSec) )
313 neutmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
314 neutnormAn[nt-1] += neutmulAn[counter];
319 for( i=0; i<numSec; i++ )
321 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
322 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
329 pv[0] = incidentParticle;
330 pv[1] = targetParticle;
335 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
337 G4int iplab = std::min(9,
G4int( incidentTotalMomentum*2.5 ));
338 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
383 npos = 0; nneg = 0; nzero = 0;
384 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
385 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
386 0.39, 0.36, 0.33, 0.10, 0.01};
387 G4int iplab =
G4int( incidentTotalMomentum*10.);
388 if ( iplab > 9) iplab = 10 +
G4int( (incidentTotalMomentum -1.)*5. );
389 if ( iplab > 14) iplab = 15 +
G4int( incidentTotalMomentum -2. );
390 if ( iplab > 22) iplab = 23 +
G4int( (incidentTotalMomentum -10.)/10.);
391 iplab = std::min(24, iplab);
396 G4double aleab = std::log(availableEnergy);
397 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
398 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
404 for (nt=1; nt<=numSec; nt++) {
405 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
406 dum = pi*nt/(2.0*n*n);
407 if (std::fabs(dum) < 1.0) {
408 if( test >= 1.0e-10 )anpn += dum*test;
416 if (targetCode == protonCode) {
418 for (npos = 0; npos < numSec/3; npos++) {
419 for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) {
420 for (nzero = 0; nzero < numSec/3; nzero++) {
421 if (++counter < numMul) {
422 nt = npos+nneg+nzero;
423 if ((nt > 0) && (nt <= numSec) ) {
424 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
425 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
426 if (std::fabs(dum) < 1.0) {
427 if( test >= 1.0e-10 )excs += dum*test;
432 if (ran < excs)
goto outOfLoop;
443 for( npos=0; npos<numSec/3; npos++ )
445 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
447 for( nzero=0; nzero<numSec/3; nzero++ )
449 if( ++counter < numMul )
451 nt = npos+nneg+nzero;
452 if ( (nt>0) && (nt<=numSec) ) {
453 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
454 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
455 if (std::fabs(dum) < 1.0) {
456 if( test >= 1.0e-10 )excs += dum*test;
461 if (ran < excs)
goto outOfLoop;
476 if( targetCode == protonCode)
481 else if (npos == (nneg+1))
528 else if ( npos == (nneg+1))
551 G4double aleab = std::log(availableEnergy);
552 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
553 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
559 for (nt=2; nt<=numSec; nt++) {
560 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
561 dum = pi*nt/(2.0*n*n);
562 if (std::fabs(dum) < 1.0) {
563 if( test >= 1.0e-10 )anpn += dum*test;
571 if( targetCode == protonCode )
574 for( npos=2; npos<numSec/3; npos++ )
577 for( nzero=0; nzero<numSec/3; nzero++ )
579 if( ++counter < numMulAn )
581 nt = npos+nneg+nzero;
582 if ( (nt>1) && (nt<=numSec) ) {
583 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
584 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
585 if (std::fabs(dum) < 1.0) {
586 if( test >= 1.0e-10 )excs += dum*test;
591 if (ran < excs)
goto outOfLoopAn;
603 for( npos=1; npos<numSec/3; npos++ )
606 for( nzero=0; nzero<numSec/3; nzero++ )
608 if (++counter < numMulAn) {
609 nt = npos+nneg+nzero;
610 if ( (nt>1) && (nt<=numSec) ) {
611 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
612 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
613 if (std::fabs(dum) < 1.0) {
614 if( test >= 1.0e-10 )excs += dum*test;
619 if (ran < excs)
goto outOfLoopAn;
632 nt = npos + nneg + nzero;
643 else if ( ran < (
G4double)(npos+nneg)/nt)
659 nt = npos + nneg + nzero;
663 G4cout <<
"Particles produced: " ;
666 for (i=2; i < vecLen; i++)
G4DLLIMPORT std::ostream G4cout
virtual void ModelDescription(std::ostream &) const
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void FirstIntInCasAntiSigmaMinus(G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)
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