Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4HESigmaPlusInelastic.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id$
27
28// G4 Process: Gheisha High Energy Collision model.
29// This includes the high energy cascading model, the two-body-resonance model
30// and the low energy two-body model. Not included are the low energy stuff
31// like nuclear reactions, nuclear fission without any cascading and all
32// processes for particles at rest.
33// First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
34// H. Fesefeldt, RWTH-Aachen, 23-October-1996
35
37#include "globals.hh"
38#include "G4ios.hh"
40#include "G4SystemOfUnits.hh"
41
42
44 : G4HEInelastic(name)
45{
46 vecLength = 0;
47 theMinEnergy = 20*GeV;
48 theMaxEnergy = 10*TeV;
49 MAXPART = 2048;
50 verboseLevel = 0;
51 G4cout << "WARNING: model G4HESigmaPlusInelastic is being deprecated and will\n"
52 << "disappear in Geant4 version 10.0" << G4endl;
53}
54
55void G4HESigmaPlusInelastic::ModelDescription(std::ostream& outFile) const
56{
57 outFile << "G4HESigmaPlusInelastic is one of the High Energy\n"
58 << "Parameterized (HEP) models used to implement inelastic\n"
59 << "Sigma+ scattering from nuclei. It is a re-engineered\n"
60 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
61 << "initial collision products into backward- and forward-going\n"
62 << "clusters which are then decayed into final state hadrons.\n"
63 << "The model does not conserve energy on an event-by-event\n"
64 << "basis. It may be applied to Sigma+ with initial energies\n"
65 << "above 20 GeV.\n";
66}
67
68
71 G4Nucleus& targetNucleus)
72{
73 G4HEVector* pv = new G4HEVector[MAXPART];
74 const G4HadProjectile* aParticle = &aTrack;
75 const G4double A = targetNucleus.GetA_asInt();
76 const G4double Z = targetNucleus.GetZ_asInt();
77 G4HEVector incidentParticle(aParticle);
78
79 G4double atomicNumber = Z;
80 G4double atomicWeight = A;
81
82 G4int incidentCode = incidentParticle.getCode();
83 G4double incidentMass = incidentParticle.getMass();
84 G4double incidentTotalEnergy = incidentParticle.getEnergy();
85
86 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
87 // DHW 19 May 2011: variable set but not used
88
89 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
90
91 if (incidentKineticEnergy < 1.)
92 G4cout << "GHESigmaPlusInelastic: incident energy < 1 GeV" << G4endl;
93
94 if (verboseLevel > 1) {
95 G4cout << "G4HESigmaPlusInelastic::ApplyYourself" << G4endl;
96 G4cout << "incident particle " << incidentParticle.getName()
97 << "mass " << incidentMass
98 << "kinetic energy " << incidentKineticEnergy
99 << G4endl;
100 G4cout << "target material with (A,Z) = ("
101 << atomicWeight << "," << atomicNumber << ")" << G4endl;
102 }
103
104 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
105 atomicWeight, atomicNumber);
106 if (verboseLevel > 1)
107 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
108
109 incidentKineticEnergy -= inelasticity;
110
111 G4double excitationEnergyGNP = 0.;
112 G4double excitationEnergyDTA = 0.;
113
114 G4double excitation = NuclearExcitation(incidentKineticEnergy,
115 atomicWeight, atomicNumber,
116 excitationEnergyGNP,
117 excitationEnergyDTA);
118 if (verboseLevel > 1)
119 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
120 << excitationEnergyDTA << G4endl;
121
122 incidentKineticEnergy -= excitation;
123 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
124 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
125 // *(incidentTotalEnergy+incidentMass));
126 // DHW 19 May 2011: variable set but not used
127
128 G4HEVector targetParticle;
129 if (G4UniformRand() < atomicNumber/atomicWeight) {
130 targetParticle.setDefinition("Proton");
131 } else {
132 targetParticle.setDefinition("Neutron");
133 }
134
135 G4double targetMass = targetParticle.getMass();
136 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
137 + targetMass*targetMass
138 + 2.0*targetMass*incidentTotalEnergy);
139 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
140
141 G4bool inElastic = true;
142 vecLength = 0;
143
144 if (verboseLevel > 1)
145 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
146 << incidentCode << G4endl;
147
148 G4bool successful = false;
149
150 FirstIntInCasSigmaPlus(inElastic, availableEnergy, pv, vecLength,
151 incidentParticle, targetParticle, atomicWeight);
152
153 if (verboseLevel > 1)
154 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
155
156 if ((vecLength > 0) && (availableEnergy > 1.))
157 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
158 pv, vecLength,
159 incidentParticle, targetParticle);
160
161 HighEnergyCascading(successful, pv, vecLength,
162 excitationEnergyGNP, excitationEnergyDTA,
163 incidentParticle, targetParticle,
164 atomicWeight, atomicNumber);
165 if (!successful)
167 excitationEnergyGNP, excitationEnergyDTA,
168 incidentParticle, targetParticle,
169 atomicWeight, atomicNumber);
170 if (!successful)
171 MediumEnergyCascading(successful, pv, vecLength,
172 excitationEnergyGNP, excitationEnergyDTA,
173 incidentParticle, targetParticle,
174 atomicWeight, atomicNumber);
175
176 if (!successful)
178 excitationEnergyGNP, excitationEnergyDTA,
179 incidentParticle, targetParticle,
180 atomicWeight, atomicNumber);
181 if (!successful)
182 QuasiElasticScattering(successful, pv, vecLength,
183 excitationEnergyGNP, excitationEnergyDTA,
184 incidentParticle, targetParticle,
185 atomicWeight, atomicNumber);
186 if (!successful)
187 ElasticScattering(successful, pv, vecLength,
188 incidentParticle,
189 atomicWeight, atomicNumber);
190
191 if (!successful)
192 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
193 << G4endl;
194
196 delete [] pv;
198 return &theParticleChange;
199}
200
201
202void
204 const G4double availableEnergy,
205 G4HEVector pv[],
206 G4int& vecLen,
207 const G4HEVector& incidentParticle,
208 const G4HEVector& targetParticle,
209 const G4double atomicWeight)
210
211// Sigma+ undergoes interaction with nucleon within a nucleus. Check if it is
212// energetically possible to produce pions/kaons. In not, assume nuclear excitation
213// occurs and input particle is degraded in energy. No other particles are produced.
214// If reaction is possible, find the correct number of pions/protons/neutrons
215// produced using an interpolation to multiplicity data. Replace some pions or
216// protons/neutrons by kaons or strange baryons according to the average
217// multiplicity per inelastic reaction.
218{
219 static const G4double expxu = 82.; // upper bound for arg. of exp
220 static const G4double expxl = -expxu; // lower bound for arg. of exp
221
222 static const G4double protb = 0.7;
223 static const G4double neutb = 0.7;
224 static const G4double c = 1.25;
225
226 static const G4int numMul = 1200;
227 static const G4int numSec = 60;
228
229 G4int protonCode = Proton.getCode();
230
231 G4int targetCode = targetParticle.getCode();
232 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
233
234 static G4bool first = true;
235 static G4double protmul[numMul], protnorm[numSec]; // proton constants
236 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
237
238 // misc. local variables
239 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
240
241 G4int i, counter, nt, npos, nneg, nzero;
242
243 if (first) {
244 // compute normalization constants, this will only be done once
245 first = false;
246 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
247 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
248 counter = -1;
249 for( npos=0; npos<(numSec/3); npos++ )
250 {
251 for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ )
252 {
253 for( nzero=0; nzero<numSec/3; nzero++ )
254 {
255 if( ++counter < numMul )
256 {
257 nt = npos+nneg+nzero;
258 if( (nt>0) && (nt<=numSec) )
259 {
260 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
261 protnorm[nt-1] += protmul[counter];
262 }
263 }
264 }
265 }
266 }
267 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
268 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
269 counter = -1;
270 for( npos=0; npos<numSec/3; npos++ )
271 {
272 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
273 {
274 for( nzero=0; nzero<numSec/3; nzero++ )
275 {
276 if( ++counter < numMul )
277 {
278 nt = npos+nneg+nzero;
279 if( (nt>0) && (nt<=numSec) )
280 {
281 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
282 neutnorm[nt-1] += neutmul[counter];
283 }
284 }
285 }
286 }
287 }
288 for( i=0; i<numSec; i++ )
289 {
290 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
291 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
292 }
293 } // end of initialization
294
295
296 // initialize the first two places
297 // the same as beam and target
298 pv[0] = incidentParticle;
299 pv[1] = targetParticle;
300 vecLen = 2;
301
302 if( !inElastic )
303 { // quasi-elastic scattering, no pions produced
304 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
305 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
306 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
307 {
308 G4double ran = G4UniformRand();
309 if( targetCode == protonCode)
310 {
311 pv[0] = Proton;
312 pv[1] = SigmaPlus;
313 }
314 else
315 {
316 if(ran < 0.2)
317 {
318 pv[0] = SigmaZero;
319 pv[1] = Proton;
320 }
321 else if(ran < 0.4)
322 {
323 pv[0] = Lambda;
324 pv[1] = Proton;
325 }
326 else if(ran < 0.6)
327 {
328 pv[0] = Neutron;
329 pv[1] = SigmaPlus;
330 }
331 else if(ran < 0.8)
332 {
333 pv[0] = Proton;
334 pv[1] = SigmaZero;
335 }
336 else
337 {
338 pv[0] = Proton;
339 pv[1] = Lambda;
340 }
341 }
342 }
343 return;
344 }
345 else if (availableEnergy <= PionPlus.getMass())
346 return;
347
348 // inelastic scattering
349
350 npos = 0; nneg = 0; nzero = 0;
351
352 // number of total particles vs. centre of mass Energy - 2*proton mass
353
354 G4double aleab = std::log(availableEnergy);
355 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
356 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
357
358 // normalization constant for kno-distribution.
359 // calculate first the sum of all constants, check for numerical problems.
360 G4double test, dum, anpn = 0.0;
361
362 for (nt=1; nt<=numSec; nt++) {
363 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
364 dum = pi*nt/(2.0*n*n);
365 if (std::fabs(dum) < 1.0) {
366 if (test >= 1.0e-10) anpn += dum*test;
367 } else {
368 anpn += dum*test;
369 }
370 }
371
372 G4double ran = G4UniformRand();
373 G4double excs = 0.0;
374 if( targetCode == protonCode )
375 {
376 counter = -1;
377 for (npos=0; npos<numSec/3; npos++) {
378 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
379 for (nzero=0; nzero<numSec/3; nzero++) {
380 if (++counter < numMul) {
381 nt = npos+nneg+nzero;
382 if ( (nt>0) && (nt<=numSec) ) {
383 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
384 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
385 if (std::fabs(dum) < 1.0) {
386 if (test >= 1.0e-10) excs += dum*test;
387 } else {
388 excs += dum*test;
389 }
390 if (ran < excs) goto outOfLoop; //----------------------->
391 }
392 }
393 }
394 }
395 }
396
397 // 3 previous loops continued to the end
398 inElastic = false; // quasi-elastic scattering
399 return;
400 }
401 else
402 { // target must be a neutron
403 counter = -1;
404 for (npos=0; npos<numSec/3; npos++) {
405 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
406 for (nzero=0; nzero<numSec/3; nzero++) {
407 if (++counter < numMul) {
408 nt = npos+nneg+nzero;
409 if ( (nt>=1) && (nt<=numSec) ) {
410 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
411 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
412 if (std::fabs(dum) < 1.0) {
413 if( test >= 1.0e-10 )excs += dum*test;
414 } else {
415 excs += dum*test;
416 }
417 if (ran < excs) goto outOfLoop; // ---------------------->
418 }
419 }
420 }
421 }
422 }
423 // 3 previous loops continued to the end
424 inElastic = false; // quasi-elastic scattering.
425 return;
426 }
427
428 outOfLoop: // <------------------------------------------------
429
430 ran = G4UniformRand();
431 if (targetCode == protonCode) {
432 if( npos == nneg)
433 {
434 }
435 else if (npos == (nneg+1))
436 {
437 if( ran < 0.25)
438 {
439 pv[0] = SigmaZero;
440 }
441 else if(ran < 0.5)
442 {
443 pv[0] = Lambda;
444 }
445 else
446 {
447 pv[1] = Neutron;
448 }
449 }
450 else
451 {
452 if(ran < 0.5)
453 {
454 pv[0] = SigmaZero;
455 pv[1] = Neutron;
456 }
457 else
458 {
459 pv[0] = Lambda;
460 pv[1] = Neutron;
461 }
462 }
463 } else {
464 if (npos == nneg)
465 {
466 if (ran < 0.5)
467 {
468 }
469 else if (ran < 0.75)
470 {
471 pv[0] = SigmaZero;
472 pv[1] = Proton;
473 }
474 else
475 {
476 pv[0] = Lambda;
477 pv[1] = Proton;
478 }
479 }
480 else if (npos == (nneg-1))
481 {
482 pv[1] = Proton;
483 }
484 else
485 {
486 if (ran < 0.5)
487 {
488 pv[0] = SigmaZero;
489 }
490 else
491 {
492 pv[0] = Lambda;
493 }
494 }
495 }
496
497 nt = npos + nneg + nzero;
498 while (nt > 0) {
499 G4double rnd = G4UniformRand();
500 if (rnd < (G4double)npos/nt) {
501 if (npos > 0) {
502 pv[vecLen++] = PionPlus;
503 npos--;
504 }
505 } else if (rnd < (G4double)(npos+nneg)/nt) {
506 if (nneg > 0) {
507 pv[vecLen++] = PionMinus;
508 nneg--;
509 }
510 } else {
511 if (nzero > 0) {
512 pv[vecLen++] = PionZero;
513 nzero--;
514 }
515 }
516 nt = npos + nneg + nzero;
517 }
518
519 if (verboseLevel > 1) {
520 G4cout << "Particles produced: " ;
521 G4cout << pv[0].getName() << " " ;
522 G4cout << pv[1].getName() << " " ;
523 for (i = 2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
524 G4cout << G4endl;
525 }
526
527 return;
528}
529
@ stopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
G4HEVector PionPlus
G4HEVector SigmaZero
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4HEVector Lambda
G4HEVector SigmaPlus
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)
G4HEVector Neutron
void FillParticleChange(G4HEVector pv[], G4int aVecLength)
G4HEVector PionMinus
void HighEnergyClusterProduction(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
G4HEVector PionZero
G4double NuclearExcitation(G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber, G4double &excitationEnergyCascade, G4double &excitationEnergyEvaporation)
G4HEVector Proton
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)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4HESigmaPlusInelastic(const G4String &name="G4HESigmaPlusInelastic")
void FirstIntInCasSigmaPlus(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 getEnergy() const
Definition: G4HEVector.cc:313
G4double getMass() const
Definition: G4HEVector.cc:361
G4int getCode() const
Definition: G4HEVector.cc:426
G4double getTotalMomentum() const
Definition: G4HEVector.cc:166
G4String getName() const
Definition: G4HEVector.cc:431
void setDefinition(G4String name)
Definition: G4HEVector.cc:812
void SetStatusChange(G4HadFinalStateStatus aS)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115