45 outFile <<
"G4HEAntiSigmaPlusInelastic 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 <<
"GHEAntiSigmaPlusInelastic: incident energy < 1 GeV" <<
G4endl;
80 G4cout <<
"G4HEAntiSigmaPlusInelastic::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);
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"
205 static const G4double expxl = -expxu;
211 static const G4int numMul = 1200;
212 static const G4int numMulAn = 400;
213 static const G4int numSec = 60;
220 static G4bool first =
true;
221 static G4double protmul[numMul], protnorm[numSec];
222 static G4double protmulAn[numMulAn],protnormAn[numSec];
223 static G4double neutmul[numMul], neutnorm[numSec];
224 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
229 G4int i, counter, nt, npos, nneg, nzero;
234 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
235 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
237 for (npos = 0; npos < (numSec/3); npos++) {
238 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
239 for (nzero = 0; nzero < numSec/3; nzero++) {
240 if (++counter < numMul) {
241 nt = npos+nneg+nzero;
242 if ((nt>0) && (nt<=numSec) ) {
243 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
244 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++) {
255 for (nneg = npos; nneg <= (npos+2); nneg++) {
256 for (nzero = 0; nzero < numSec/3; nzero++) {
257 if (++counter < numMul) {
258 nt = npos+nneg+nzero;
259 if ((nt>0) && (nt<=numSec) ) {
260 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
261 neutnorm[nt-1] += neutmul[counter];
268 for (i = 0; i < numSec; i++) {
269 if (protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
270 if (neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
274 for (i = 0; i < numMulAn; i++) protmulAn[i] = 0.0;
275 for (i = 0; i < numSec; i++) protnormAn[i] = 0.0;
277 for (npos = 1; npos < (numSec/3); npos++) {
279 for (nzero = 0; nzero < numSec/3; nzero++) {
280 if (++counter < numMulAn) {
281 nt = npos+nneg+nzero;
282 if ((nt>1) && (nt<=numSec) ) {
283 protmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
284 protnormAn[nt-1] += protmulAn[counter];
290 for (i = 0; i < numMulAn; i++) neutmulAn[i] = 0.0;
291 for (i = 0; i < numSec; i++) neutnormAn[i] = 0.0;
293 for (npos = 0; npos < numSec/3; npos++) {
295 for (nzero = 0; nzero < numSec/3; nzero++) {
296 if (++counter < numMulAn) {
297 nt = npos+nneg+nzero;
298 if ((nt>1) && (nt<=numSec) ) {
299 neutmulAn[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
300 neutnormAn[nt-1] += neutmulAn[counter];
305 for (i = 0; i < numSec; i++) {
306 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
307 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
313 pv[0] = incidentParticle;
314 pv[1] = targetParticle;
318 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
320 G4int iplab = std::min(9,
G4int( incidentTotalMomentum*2.5));
321 if (
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
324 if (targetCode == protonCode) {
328 }
else if (ran < 0.4) {
331 }
else if (ran < 0.6) {
334 }
else if (ran < 0.8) {
352 npos = 0; nneg = 0; nzero = 0;
353 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
354 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
355 0.39, 0.36, 0.33, 0.10, 0.01};
356 G4int iplab =
G4int( incidentTotalMomentum*10.);
357 if ( iplab > 9) iplab = 10 +
G4int( (incidentTotalMomentum -1.)*5. );
358 if ( iplab > 14) iplab = 15 +
G4int( incidentTotalMomentum -2. );
359 if ( iplab > 22) iplab = 23 +
G4int( (incidentTotalMomentum -10.)/10.);
360 iplab = std::min(24, iplab);
365 G4double aleab = std::log(availableEnergy);
366 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
367 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
373 for (nt=1; nt<=numSec; nt++) {
374 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
375 dum = pi*nt/(2.0*n*n);
376 if (std::fabs(dum) < 1.0) {
377 if( test >= 1.0e-10 )anpn += dum*test;
385 if( targetCode == protonCode )
388 for( npos=0; npos<numSec/3; npos++ )
390 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
392 for( nzero=0; nzero<numSec/3; nzero++ )
394 if( ++counter < numMul )
396 nt = npos+nneg+nzero;
397 if ( (nt>0) && (nt<=numSec) ) {
398 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
399 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
400 if (std::fabs(dum) < 1.0) {
401 if( test >= 1.0e-10 )excs += dum*test;
406 if (ran < excs)
goto outOfLoop;
420 for( npos=0; npos<numSec/3; npos++ )
422 for( nneg=npos; nneg<=(npos+2); nneg++ )
424 for( nzero=0; nzero<numSec/3; nzero++ )
426 if( ++counter < numMul )
428 nt = npos+nneg+nzero;
429 if ( (nt>0) && (nt<=numSec) ) {
430 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
431 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
432 if (std::fabs(dum) < 1.0) {
433 if( test >= 1.0e-10 )excs += dum*test;
438 if (ran < excs)
goto outOfLoop;
453 if( targetCode == protonCode)
471 else if (npos == (nneg-1))
492 else if ( npos == (nneg-1))
527 G4double aleab = std::log(availableEnergy);
528 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
529 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
535 for (nt=2; nt<=numSec; nt++) {
536 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
537 dum = pi*nt/(2.0*n*n);
538 if (std::fabs(dum) < 1.0) {
539 if( test >= 1.0e-10 )anpn += dum*test;
547 if( targetCode == protonCode )
550 for( npos=1; npos<numSec/3; npos++ )
553 for( nzero=0; nzero<numSec/3; nzero++ )
555 if( ++counter < numMulAn )
557 nt = npos+nneg+nzero;
558 if ( (nt>1) && (nt<=numSec) ) {
559 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
560 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
561 if (std::fabs(dum) < 1.0) {
562 if( test >= 1.0e-10 )excs += dum*test;
567 if (ran < excs)
goto outOfLoopAn;
579 for( npos=0; npos<numSec/3; npos++ )
582 for( nzero=0; nzero<numSec/3; nzero++ )
584 if( ++counter < numMulAn )
586 nt = npos+nneg+nzero;
587 if ( (nt>1) && (nt<=numSec) ) {
588 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
589 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
590 if (std::fabs(dum) < 1.0) {
591 if( test >= 1.0e-10 )excs += dum*test;
596 if (ran < excs)
goto outOfLoopAn;
609 nt = npos + nneg + nzero;
620 else if ( ran < (
G4double)(npos+nneg)/nt)
636 nt = npos + nneg + nzero;
640 G4cout <<
"Particles produced: " ;
643 for (i=2; i < vecLen; i++)
G4DLLIMPORT std::ostream G4cout
virtual void ModelDescription(std::ostream &) const
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void FirstIntInCasAntiSigmaPlus(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)
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