Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEAntiSigmaMinusInelastic.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 << "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";
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 << "GHEAntiSigmaMinusInelastic: incident energy < 1 GeV" << G4endl;
78 
79  if (verboseLevel > 1) {
80  G4cout << "G4HEAntiSigmaMinusInelastic::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  FirstIntInCasAntiSigmaMinus(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  HighEnergyCascading(successful, pv, vecLength,
146  excitationEnergyGNP, excitationEnergyDTA,
147  incidentParticle, targetParticle,
148  atomicWeight, atomicNumber);
149  if (!successful)
150  HighEnergyClusterProduction(successful, pv, vecLength,
151  excitationEnergyGNP, excitationEnergyDTA,
152  incidentParticle, targetParticle,
153  atomicWeight, atomicNumber);
154  if (!successful)
155  MediumEnergyCascading(successful, pv, vecLength,
156  excitationEnergyGNP, excitationEnergyDTA,
157  incidentParticle, targetParticle,
158  atomicWeight, atomicNumber);
159 
160  if (!successful)
162  excitationEnergyGNP, excitationEnergyDTA,
163  incidentParticle, targetParticle,
164  atomicWeight, atomicNumber);
165  if (!successful)
166  QuasiElasticScattering(successful, pv, vecLength,
167  excitationEnergyGNP, excitationEnergyDTA,
168  incidentParticle, targetParticle,
169  atomicWeight, atomicNumber);
170  if (!successful)
171  ElasticScattering(successful, pv, vecLength,
172  incidentParticle,
173  atomicWeight, atomicNumber);
174 
175  if (!successful)
176  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
177  << G4endl;
179  delete [] pv;
181  return &theParticleChange;
182 }
183 
184 
185 void
187  const G4double availableEnergy,
188  G4HEVector pv[],
189  G4int& vecLen,
190  const G4HEVector& incidentParticle,
191  const G4HEVector& targetParticle,
192  const G4double atomicWeight)
193 
194 // AntiSigma- undergoes interaction with nucleon within a nucleus. Check if it is
195 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
196 // occurs and input particle is degraded in energy. No other particles are produced.
197 // If reaction is possible, find the correct number of pions/protons/neutrons
198 // produced using an interpolation to multiplicity data. Replace some pions or
199 // protons/neutrons by kaons or strange baryons according to the average
200 // multiplicity per inelastic reaction.
201 {
202  static const G4double expxu = 82; // upper bound for arg. of exp
203  static const G4double expxl = -expxu; // lower bound for arg. of exp
204 
205  static const G4double protb = 0.7;
206  static const G4double neutb = 0.7;
207  static const G4double c = 1.25;
208 
209  static const G4int numMul = 1200;
210  static const G4int numMulAn = 400;
211  static const G4int numSec = 60;
212 
214  G4int protonCode = Proton.getCode();
215 
216  G4int targetCode = targetParticle.getCode();
217  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
218 
219  static G4bool first = true;
220  static G4double protmul[numMul], protnorm[numSec]; // proton constants
221  static G4double protmulAn[numMulAn],protnormAn[numSec];
222  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
223  static G4double neutmulAn[numMulAn],neutnormAn[numSec];
224 
225  // misc. local variables
226  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
227 
228  G4int i, counter, nt, npos, nneg, nzero;
229 
230  if( first )
231  { // compute normalization constants, this will only be done once
232  first = false;
233  for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
234  for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
235  counter = -1;
236  for( npos=0; npos<(numSec/3); npos++ )
237  {
238  for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ )
239  {
240  for( nzero=0; nzero<numSec/3; nzero++ )
241  {
242  if( ++counter < numMul )
243  {
244  nt = npos+nneg+nzero;
245  if( (nt>0) && (nt<=numSec) )
246  {
247  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
248  protnorm[nt-1] += protmul[counter];
249  }
250  }
251  }
252  }
253  }
254  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
255  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
256  counter = -1;
257  for( npos=0; npos<numSec/3; npos++ )
258  {
259  for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
260  {
261  for( nzero=0; nzero<numSec/3; nzero++ )
262  {
263  if( ++counter < numMul )
264  {
265  nt = npos+nneg+nzero;
266  if( (nt>0) && (nt<=numSec) )
267  {
268  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
269  neutnorm[nt-1] += neutmul[counter];
270  }
271  }
272  }
273  }
274  }
275  for( i=0; i<numSec; i++ )
276  {
277  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
278  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
279  }
280  // annihilation
281  for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
282  for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
283  counter = -1;
284  for( npos=1; npos<(numSec/3); npos++ )
285  {
286  nneg = std::max(0,npos-2);
287  for( nzero=0; nzero<numSec/3; nzero++ )
288  {
289  if( ++counter < numMulAn )
290  {
291  nt = npos+nneg+nzero;
292  if( (nt>1) && (nt<=numSec) )
293  {
294  protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
295  protnormAn[nt-1] += protmulAn[counter];
296  }
297  }
298  }
299  }
300  for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
301  for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
302  counter = -1;
303  for( npos=0; npos<numSec/3; npos++ )
304  {
305  nneg = npos-1;
306  for( nzero=0; nzero<numSec/3; nzero++ )
307  {
308  if( ++counter < numMulAn )
309  {
310  nt = npos+nneg+nzero;
311  if( (nt>1) && (nt<=numSec) )
312  {
313  neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
314  neutnormAn[nt-1] += neutmulAn[counter];
315  }
316  }
317  }
318  }
319  for( i=0; i<numSec; i++ )
320  {
321  if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
322  if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
323  }
324  } // end of initialization
325 
326 
327  // initialize the first two places
328  // the same as beam and target
329  pv[0] = incidentParticle;
330  pv[1] = targetParticle;
331  vecLen = 2;
332 
333  if( !inElastic )
334  { // some two-body reactions
335  G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
336 
337  G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5 ));
338  if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
339  {
340  G4double ran = G4UniformRand();
341 
342  if ( targetCode == neutronCode)
343  {
344  if(ran < 0.2)
345  {
346  pv[0] = AntiSigmaZero;
347  pv[1] = Proton;
348  }
349  else if (ran < 0.4)
350  {
351  pv[0] = AntiLambda;
352  pv[1] = Proton;
353  }
354  else if (ran < 0.6)
355  {
356  pv[0] = Proton;
357  pv[1] = AntiLambda;
358  }
359  else if (ran < 0.8)
360  {
361  pv[0] = Proton;
362  pv[1] = AntiSigmaZero;
363  }
364  else
365  {
366  pv[0] = Neutron;
367  pv[1] = AntiSigmaMinus;
368  }
369  }
370  else
371  {
372  pv[0] = Proton;
373  pv[1] = AntiSigmaMinus;
374  }
375  }
376  return;
377  }
378  else if (availableEnergy <= PionPlus.getMass())
379  return;
380 
381  // inelastic scattering
382 
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);
392 
393  if (G4UniformRand() > anhl[iplab]) { // non- annihilation channels
394 
395  // number of total particles vs. centre of mass Energy - 2*proton mass
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;
399 
400  // normalization constant for kno-distribution.
401  // calculate first the sum of all constants, check for numerical problems.
402  G4double test, dum, anpn = 0.0;
403 
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;
409  } else {
410  anpn += dum*test;
411  }
412  }
413 
414  G4double ran = G4UniformRand();
415  G4double excs = 0.0;
416  if (targetCode == protonCode) {
417  counter = -1;
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;
428  } else {
429  excs += dum*test;
430  }
431 
432  if (ran < excs) goto outOfLoop; //----------------------->
433  }
434  }
435  }
436  }
437  }
438  // 3 previous loops continued to the end
439  inElastic = false; // quasi-elastic scattering
440  return;
441  } else { // target must be a neutron
442  counter = -1;
443  for( npos=0; npos<numSec/3; npos++ )
444  {
445  for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
446  {
447  for( nzero=0; nzero<numSec/3; nzero++ )
448  {
449  if( ++counter < numMul )
450  {
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;
457  } else {
458  excs += dum*test;
459  }
460 
461  if (ran < excs) goto outOfLoop; // -------------------------->
462  }
463  }
464  }
465  }
466  }
467  // 3 previous loops continued to the end
468  inElastic = false; // quasi-elastic scattering.
469  return;
470  }
471 
472  outOfLoop: // <------------------------------------------------------------------------
473 
474  ran = G4UniformRand();
475 
476  if( targetCode == protonCode)
477  {
478  if( npos == nneg)
479  {
480  }
481  else if (npos == (nneg+1))
482  {
483  if( ran < 0.50)
484  {
485  pv[1] = Neutron;
486  }
487  else if (ran < 0.75)
488  {
489  pv[0] = AntiSigmaZero;
490  }
491  else
492  {
493  pv[0] = AntiLambda;
494  }
495  }
496  else
497  {
498  if (ran < 0.5)
499  {
500  pv[0] = AntiSigmaZero;
501  pv[1] = Neutron;
502  }
503  else
504  {
505  pv[0] = AntiLambda;
506  pv[1] = Neutron;
507  }
508  }
509  }
510  else
511  {
512  if( npos == nneg)
513  {
514  if (ran < 0.5)
515  {
516  }
517  else if(ran < 0.75)
518  {
519  pv[0] = AntiSigmaZero;
520  pv[1] = Proton;
521  }
522  else
523  {
524  pv[0] = AntiLambda;
525  pv[1] = Proton;
526  }
527  }
528  else if ( npos == (nneg+1))
529  {
530  if (ran < 0.5)
531  {
532  pv[0] = AntiSigmaZero;
533  }
534  else
535  {
536  pv[0] = AntiLambda;
537  }
538  }
539  else
540  {
541  pv[1] = Proton;
542  }
543  }
544 
545  }
546  else // annihilation
547  {
548  if ( availableEnergy > 2. * PionPlus.getMass() )
549  {
550 
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;
554 
555  // normalization constant for kno-distribution.
556  // calculate first the sum of all constants, check for numerical problems.
557  G4double test, dum, anpn = 0.0;
558 
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;
564  } else {
565  anpn += dum*test;
566  }
567  }
568 
569  G4double ran = G4UniformRand();
570  G4double excs = 0.0;
571  if( targetCode == protonCode )
572  {
573  counter = -1;
574  for( npos=2; npos<numSec/3; npos++ )
575  {
576  nneg = npos-2;
577  for( nzero=0; nzero<numSec/3; nzero++ )
578  {
579  if( ++counter < numMulAn )
580  {
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;
587  } else {
588  excs += dum*test;
589  }
590 
591  if (ran < excs) goto outOfLoopAn; //----------------------->
592  }
593  }
594  }
595  }
596  // 3 previous loops continued to the end
597  inElastic = false; // quasi-elastic scattering
598  return;
599  }
600  else
601  { // target must be a neutron
602  counter = -1;
603  for( npos=1; npos<numSec/3; npos++ )
604  {
605  nneg = npos-1;
606  for( nzero=0; nzero<numSec/3; nzero++ )
607  {
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;
615  } else {
616  excs += dum*test;
617  }
618 
619  if (ran < excs) goto outOfLoopAn; // -------------------------->
620  }
621  }
622  }
623  }
624  inElastic = false; // quasi-elastic scattering.
625  return;
626  }
627  outOfLoopAn: // <------------------------------------------------------------------
628  vecLen = 0;
629  }
630  }
631 
632  nt = npos + nneg + nzero;
633  while ( nt > 0)
634  {
635  G4double ran = G4UniformRand();
636  if ( ran < (G4double)npos/nt)
637  {
638  if( npos > 0 )
639  { pv[vecLen++] = PionPlus;
640  npos--;
641  }
642  }
643  else if ( ran < (G4double)(npos+nneg)/nt)
644  {
645  if( nneg > 0 )
646  {
647  pv[vecLen++] = PionMinus;
648  nneg--;
649  }
650  }
651  else
652  {
653  if( nzero > 0 )
654  {
655  pv[vecLen++] = PionZero;
656  nzero--;
657  }
658  }
659  nt = npos + nneg + nzero;
660  }
661  if (verboseLevel > 1)
662  {
663  G4cout << "Particles produced: " ;
664  G4cout << pv[0].getName() << " " ;
665  G4cout << pv[1].getName() << " " ;
666  for (i=2; i < vecLen; i++)
667  {
668  G4cout << pv[i].getName() << " " ;
669  }
670  G4cout << G4endl;
671  }
672  return;
673  }
674 
675 
676 
677 
678 
679 
680 
681 
682 
683