Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEAntiSigmaPlusInelastic.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 
29 // G4 Process: Gheisha High Energy Collision model.
30 // This includes the high energy cascading model, the two-body-resonance model
31 // and the low energy two-body model. Not included are the low energy stuff
32 // like nuclear reactions, nuclear fission without any cascading and all
33 // processes for particles at rest.
34 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
35 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
36 
38 #include "globals.hh"
39 #include "G4ios.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 
44 {
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";
54 }
55 
56 
59  G4Nucleus& targetNucleus)
60 {
61  G4HEVector* pv = new G4HEVector[MAXPART];
62  const G4HadProjectile* aParticle = &aTrack;
63  const G4double atomicWeight = targetNucleus.GetA_asInt();
64  const G4double atomicNumber = targetNucleus.GetZ_asInt();
65  G4HEVector incidentParticle(aParticle);
66 
67  G4int incidentCode = incidentParticle.getCode();
68  G4double incidentMass = incidentParticle.getMass();
69  G4double incidentTotalEnergy = incidentParticle.getEnergy();
70 
71  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
72  // DHW 19 May 2011: variable set but not used
73 
74  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
75 
76  if (incidentKineticEnergy < 1.)
77  G4cout << "GHEAntiSigmaPlusInelastic: incident energy < 1 GeV" << G4endl;
78 
79  if (verboseLevel > 1) {
80  G4cout << "G4HEAntiSigmaPlusInelastic::ApplyYourself" << G4endl;
81  G4cout << "incident particle " << incidentParticle.getName()
82  << "mass " << incidentMass
83  << "kinetic energy " << incidentKineticEnergy
84  << G4endl;
85  G4cout << "target material with (A,Z) = ("
86  << atomicWeight << "," << atomicNumber << ")" << G4endl;
87  }
88 
89  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
90  atomicWeight, atomicNumber);
91  if (verboseLevel > 1)
92  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
93 
94  incidentKineticEnergy -= inelasticity;
95 
96  G4double excitationEnergyGNP = 0.;
97  G4double excitationEnergyDTA = 0.;
98 
99  G4double excitation = NuclearExcitation(incidentKineticEnergy,
100  atomicWeight, atomicNumber,
101  excitationEnergyGNP,
102  excitationEnergyDTA);
103  if (verboseLevel > 1)
104  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
105  << excitationEnergyDTA << G4endl;
106 
107  incidentKineticEnergy -= excitation;
108  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
109  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
110  // *(incidentTotalEnergy+incidentMass));
111  // DHW 19 May 2011: variable set but not used
112 
113  G4HEVector targetParticle;
114  if (G4UniformRand() < atomicNumber/atomicWeight) {
115  targetParticle.setDefinition("Proton");
116  } else {
117  targetParticle.setDefinition("Neutron");
118  }
119 
120  G4double targetMass = targetParticle.getMass();
121  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
122  + targetMass*targetMass
123  + 2.0*targetMass*incidentTotalEnergy);
124  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
125 
126  G4bool inElastic = true;
127  vecLength = 0;
128 
129  if (verboseLevel > 1)
130  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
131  << incidentCode << G4endl;
132 
133  G4bool successful = false;
134 
135  FirstIntInCasAntiSigmaPlus(inElastic, availableEnergy, pv, vecLength,
136  incidentParticle, targetParticle, atomicWeight);
137 
138  if (verboseLevel > 1)
139  G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
140 
141  if ((vecLength > 0) && (availableEnergy > 1.))
142  StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
143  pv, vecLength,
144  incidentParticle, targetParticle);
145 
146  HighEnergyCascading(successful, pv, vecLength,
147  excitationEnergyGNP, excitationEnergyDTA,
148  incidentParticle, targetParticle,
149  atomicWeight, atomicNumber);
150  if (!successful)
151  HighEnergyClusterProduction(successful, pv, vecLength,
152  excitationEnergyGNP, excitationEnergyDTA,
153  incidentParticle, targetParticle,
154  atomicWeight, atomicNumber);
155  if (!successful)
156  MediumEnergyCascading(successful, pv, vecLength,
157  excitationEnergyGNP, excitationEnergyDTA,
158  incidentParticle, targetParticle,
159  atomicWeight, atomicNumber);
160 
161  if (!successful)
163  excitationEnergyGNP, excitationEnergyDTA,
164  incidentParticle, targetParticle,
165  atomicWeight, atomicNumber);
166  if (!successful)
167  QuasiElasticScattering(successful, pv, vecLength,
168  excitationEnergyGNP, excitationEnergyDTA,
169  incidentParticle, targetParticle,
170  atomicWeight, atomicNumber);
171  if (!successful)
172  ElasticScattering(successful, pv, vecLength,
173  incidentParticle,
174  atomicWeight, atomicNumber);
175 
176  if (!successful)
177  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
178  << G4endl;
179 
181  delete [] pv;
183  return &theParticleChange;
184 }
185 
186 
187 void
189  const G4double availableEnergy,
190  G4HEVector pv[],
191  G4int& vecLen,
192  const G4HEVector& incidentParticle,
193  const G4HEVector& targetParticle,
194  const G4double atomicWeight)
195 
196 // AntiSigma+ undergoes interaction with nucleon within a nucleus. Check if it is
197 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
198 // occurs and input particle is degraded in energy. No other particles are produced.
199 // If reaction is possible, find the correct number of pions/protons/neutrons
200 // produced using an interpolation to multiplicity data. Replace some pions or
201 // protons/neutrons by kaons or strange baryons according to the average
202 // multiplicity per inelastic reaction.
203 {
204  static const G4double expxu = 82.; // upper bound for arg. of exp
205  static const G4double expxl = -expxu; // lower bound for arg. of exp
206 
207  static const G4double protb = 0.7;
208  static const G4double neutb = 0.7;
209  static const G4double c = 1.25;
210 
211  static const G4int numMul = 1200;
212  static const G4int numMulAn = 400;
213  static const G4int numSec = 60;
214 
215  G4int protonCode = Proton.getCode();
216 
217  G4int targetCode = targetParticle.getCode();
218  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
219 
220  static G4bool first = true;
221  static G4double protmul[numMul], protnorm[numSec]; // proton constants
222  static G4double protmulAn[numMulAn],protnormAn[numSec];
223  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
224  static G4double neutmulAn[numMulAn],neutnormAn[numSec];
225 
226  // misc. local variables
227  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
228 
229  G4int i, counter, nt, npos, nneg, nzero;
230 
231  if (first) {
232  // compute normalization constants, this will only be done once
233  first = false;
234  for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
235  for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
236  counter = -1;
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];
245  }
246  }
247  }
248  }
249  }
250 
251  for (i = 0; i < numMul; i++) neutmul[i] = 0.0;
252  for (i = 0; i < numSec; i++) neutnorm[i] = 0.0;
253  counter = -1;
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];
262  }
263  }
264  }
265  }
266  }
267 
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];
271  }
272 
273  // annihilation
274  for (i = 0; i < numMulAn; i++) protmulAn[i] = 0.0;
275  for (i = 0; i < numSec; i++) protnormAn[i] = 0.0;
276  counter = -1;
277  for (npos = 1; npos < (numSec/3); npos++) {
278  nneg = 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];
285  }
286  }
287  }
288  }
289 
290  for (i = 0; i < numMulAn; i++) neutmulAn[i] = 0.0;
291  for (i = 0; i < numSec; i++) neutnormAn[i] = 0.0;
292  counter = -1;
293  for (npos = 0; npos < numSec/3; npos++) {
294  nneg = npos+1;
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];
301  }
302  }
303  }
304  }
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];
308  }
309  } // end of initialization
310 
311 
312  // initialize the first two places the same as beam and target
313  pv[0] = incidentParticle;
314  pv[1] = targetParticle;
315  vecLen = 2;
316 
317  if (!inElastic) { // some two-body reactions
318  G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
319 
320  G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
321  if (G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
322  G4double ran = G4UniformRand();
323 
324  if (targetCode == protonCode) {
325  if (ran < 0.2) {
326  pv[0] = Proton;
327  pv[1] = AntiSigmaPlus;
328  } else if (ran < 0.4) {
329  pv[0] = AntiLambda;
330  pv[1] = Neutron;
331  } else if (ran < 0.6) {
332  pv[0] = Neutron;
333  pv[1] = AntiLambda;
334  } else if (ran < 0.8) {
335  pv[0] = Neutron;
336  pv[1] = AntiSigmaZero;
337  } else {
338  pv[0] = AntiSigmaZero;
339  pv[1] = Neutron;
340  }
341  } else {
342  pv[0] = Neutron;
343  pv[1] = AntiSigmaPlus;
344  }
345  }
346 
347  return;
348  }
349  else if (availableEnergy <= PionPlus.getMass()) return;
350 
351  // inelastic scattering
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);
361 
362  if ( G4UniformRand() > anhl[iplab] ) { // non- annihilation channels
363 
364  // number of total particles vs. centre of mass Energy - 2*proton mass
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;
368 
369  // normalization constant for kno-distribution.
370  // calculate first the sum of all constants, check for numerical problems.
371  G4double test, dum, anpn = 0.0;
372 
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;
378  } else {
379  anpn += dum*test;
380  }
381  }
382 
383  G4double ran = G4UniformRand();
384  G4double excs = 0.0;
385  if( targetCode == protonCode )
386  {
387  counter = -1;
388  for( npos=0; npos<numSec/3; npos++ )
389  {
390  for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
391  {
392  for( nzero=0; nzero<numSec/3; nzero++ )
393  {
394  if( ++counter < numMul )
395  {
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;
402  } else {
403  excs += dum*test;
404  }
405 
406  if (ran < excs) goto outOfLoop; //----------------------->
407  }
408  }
409  }
410  }
411  }
412 
413  // 3 previous loops continued to the end
414  inElastic = false; // quasi-elastic scattering
415  return;
416  }
417  else
418  { // target must be a neutron
419  counter = -1;
420  for( npos=0; npos<numSec/3; npos++ )
421  {
422  for( nneg=npos; nneg<=(npos+2); nneg++ )
423  {
424  for( nzero=0; nzero<numSec/3; nzero++ )
425  {
426  if( ++counter < numMul )
427  {
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;
434  } else {
435  excs += dum*test;
436  }
437 
438  if (ran < excs) goto outOfLoop; // -------------------------->
439  }
440  }
441  }
442  }
443  }
444  // 3 previous loops continued to the end
445  inElastic = false; // quasi-elastic scattering.
446  return;
447  }
448 
449  outOfLoop: // <------------------------------------------------------------------------
450 
451  ran = G4UniformRand();
452 
453  if( targetCode == protonCode)
454  {
455  if( npos == nneg)
456  {
457  if (ran < 0.50)
458  {
459  }
460  else if (ran < 0.75)
461  {
462  pv[0] = AntiSigmaZero;
463  pv[1] = Neutron;
464  }
465  else
466  {
467  pv[0] = AntiLambda;
468  pv[1] = Neutron;
469  }
470  }
471  else if (npos == (nneg-1))
472  {
473  if( ran < 0.50)
474  {
475  pv[0] = AntiLambda;
476  }
477  else
478  {
479  pv[0] = AntiSigmaZero;
480  }
481  }
482  else
483  {
484  pv[1] = Neutron;
485  }
486  }
487  else
488  {
489  if( npos == nneg)
490  {
491  }
492  else if ( npos == (nneg-1))
493  {
494  if (ran < 0.5)
495  {
496  pv[1] = Proton;
497  }
498  else if (ran < 0.75)
499  {
500  pv[0] = AntiLambda;
501  }
502  else
503  {
504  pv[0] = AntiSigmaZero;
505  }
506  }
507  else
508  {
509  if (ran < 0.5)
510  {
511  pv[0] = AntiLambda;
512  pv[1] = Proton;
513  }
514  else
515  {
516  pv[0] = AntiSigmaZero;
517  pv[1] = Proton;
518  }
519  }
520  }
521  }
522  else // annihilation
523  {
524  if ( availableEnergy > 2. * PionPlus.getMass() )
525  {
526 
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;
530 
531  // normalization constant for kno-distribution.
532  // calculate first the sum of all constants, check for numerical problems.
533  G4double test, dum, anpn = 0.0;
534 
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;
540  } else {
541  anpn += dum*test;
542  }
543  }
544 
545  G4double ran = G4UniformRand();
546  G4double excs = 0.0;
547  if( targetCode == protonCode )
548  {
549  counter = -1;
550  for( npos=1; npos<numSec/3; npos++ )
551  {
552  nneg = npos;
553  for( nzero=0; nzero<numSec/3; nzero++ )
554  {
555  if( ++counter < numMulAn )
556  {
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;
563  } else {
564  excs += dum*test;
565  }
566 
567  if (ran < excs) goto outOfLoopAn; //----------------------->
568  }
569  }
570  }
571  }
572  // 3 previous loops continued to the end
573  inElastic = false; // quasi-elastic scattering
574  return;
575  }
576  else
577  { // target must be a neutron
578  counter = -1;
579  for( npos=0; npos<numSec/3; npos++ )
580  {
581  nneg = npos+1;
582  for( nzero=0; nzero<numSec/3; nzero++ )
583  {
584  if( ++counter < numMulAn )
585  {
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;
592  } else {
593  excs += dum*test;
594  }
595 
596  if (ran < excs) goto outOfLoopAn; // -------------------------->
597  }
598  }
599  }
600  }
601  inElastic = false; // quasi-elastic scattering.
602  return;
603  }
604  outOfLoopAn: // <----------------------------------------
605  vecLen = 0;
606  }
607  }
608 
609  nt = npos + nneg + nzero;
610  while ( nt > 0)
611  {
612  G4double ran = G4UniformRand();
613  if ( ran < (G4double)npos/nt)
614  {
615  if( npos > 0 )
616  { pv[vecLen++] = PionPlus;
617  npos--;
618  }
619  }
620  else if ( ran < (G4double)(npos+nneg)/nt)
621  {
622  if( nneg > 0 )
623  {
624  pv[vecLen++] = PionMinus;
625  nneg--;
626  }
627  }
628  else
629  {
630  if( nzero > 0 )
631  {
632  pv[vecLen++] = PionZero;
633  nzero--;
634  }
635  }
636  nt = npos + nneg + nzero;
637  }
638  if (verboseLevel > 1)
639  {
640  G4cout << "Particles produced: " ;
641  G4cout << pv[0].getName() << " " ;
642  G4cout << pv[1].getName() << " " ;
643  for (i=2; i < vecLen; i++)
644  {
645  G4cout << pv[i].getName() << " " ;
646  }
647  G4cout << G4endl;
648  }
649  return;
650  }
651 
652 
653 
654 
655 
656 
657 
658 
659 
660