Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEPionMinusInelastic.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 // 11-OCT-2007 F.W. Jones: fixed incorrect Imax (should be Imin) in
29 // sampling of charge exchange.
30 
31 // G4 Process: Gheisha High Energy Collision model.
32 // This includes the high energy cascading model, the two-body-resonance model
33 // and the low energy two-body model. Not included are the low energy stuff
34 // like nuclear reactions, nuclear fission without any cascading and all
35 // processes for particles at rest.
36 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
37 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
38 
39 #include "G4HEPionMinusInelastic.hh"
40 #include "globals.hh"
41 #include "G4ios.hh"
42 #include "G4PhysicalConstants.hh"
43 
44 void G4HEPionMinusInelastic::ModelDescription(std::ostream& outFile) const
45 {
46  outFile << "G4HEPionMinusInelastic is one of the High Energy\n"
47  << "Parameterized (HEP) models used to implement inelastic\n"
48  << "pi- scattering from nuclei. It is a re-engineered\n"
49  << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
50  << "initial collision products into backward- and forward-going\n"
51  << "clusters which are then decayed into final state hadrons.\n"
52  << "The model does not conserve energy on an event-by-event\n"
53  << "basis. It may be applied to pi- with initial energies\n"
54  << "above 20 GeV.\n";
55 }
56 
57 
58 G4HadFinalState*
59 G4HEPionMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
60  G4Nucleus& targetNucleus)
61 {
62  G4HEVector* pv = new G4HEVector[MAXPART];
63  const G4HadProjectile* aParticle = &aTrack;
64  const G4double A = targetNucleus.GetA_asInt();
65  const G4double Z = targetNucleus.GetZ_asInt();
66  G4HEVector incidentParticle(aParticle);
67 
68  G4double atomicNumber = Z;
69  G4double atomicWeight = A;
70 
71  G4int incidentCode = incidentParticle.getCode();
72  G4double incidentMass = incidentParticle.getMass();
73  G4double incidentTotalEnergy = incidentParticle.getEnergy();
74 
75  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
76  // DHW 19 May 2011: variable set but not used
77 
78  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
79 
80  if (incidentKineticEnergy < 1.)
81  G4cout << "GHEPionMinusInelastic: incident energy < 1 GeV" << G4endl;
82 
83  if (verboseLevel > 1) {
84  G4cout << "G4HEPionMinusInelastic::ApplyYourself" << G4endl;
85  G4cout << "incident particle " << incidentParticle.getName()
86  << "mass " << incidentMass
87  << "kinetic energy " << incidentKineticEnergy
88  << G4endl;
89  G4cout << "target material with (A,Z) = ("
90  << atomicWeight << "," << atomicNumber << ")" << G4endl;
91  }
92 
93  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
94  atomicWeight, atomicNumber);
95  if (verboseLevel > 1)
96  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
97 
98  incidentKineticEnergy -= inelasticity;
99 
100  G4double excitationEnergyGNP = 0.;
101  G4double excitationEnergyDTA = 0.;
102 
103  G4double excitation = NuclearExcitation(incidentKineticEnergy,
104  atomicWeight, atomicNumber,
105  excitationEnergyGNP,
106  excitationEnergyDTA);
107  if (verboseLevel > 1)
108  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
109  << excitationEnergyDTA << G4endl;
110 
111  incidentKineticEnergy -= excitation;
112  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
113  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
114  // *(incidentTotalEnergy+incidentMass));
115  // DHW 19 May 2011: variable set but not used
116 
117  G4HEVector targetParticle;
118 
119  if (G4UniformRand() < atomicNumber/atomicWeight) {
120  targetParticle.setDefinition("Proton");
121  } else {
122  targetParticle.setDefinition("Neutron");
123  }
124 
125  G4double targetMass = targetParticle.getMass();
126  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
127  + targetMass*targetMass
128  + 2.0*targetMass*incidentTotalEnergy);
129  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
130 
131  // The value of the inElastic flag was originally defined in the Gheisha
132  // stand-alone code as follows:
133  // G4bool inElastic = InElasticCrossSectionInFirstInt
134  // (availableEnergy, incidentCode, incidentTotalMomentum);
135  // For unknown reasons, it was replaced by the following code in Geant
136 
137  G4bool inElastic = true;
138  // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;
139 
140  vecLength = 0;
141 
142  if (verboseLevel > 1)
143  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
144  << incidentCode << G4endl;
145 
146  G4bool successful = false;
147 
148  FirstIntInCasPionMinus(inElastic, availableEnergy, pv, vecLength,
149  incidentParticle, targetParticle);
150 
151  if (verboseLevel > 1)
152  G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
153 
154  if ((vecLength > 0) && (availableEnergy > 1.))
155  StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
156  pv, vecLength,
157  incidentParticle, targetParticle);
158 
159  HighEnergyCascading(successful, pv, vecLength,
160  excitationEnergyGNP, excitationEnergyDTA,
161  incidentParticle, targetParticle,
162  atomicWeight, atomicNumber);
163  if (!successful)
164  HighEnergyClusterProduction(successful, pv, vecLength,
165  excitationEnergyGNP, excitationEnergyDTA,
166  incidentParticle, targetParticle,
167  atomicWeight, atomicNumber);
168  if (!successful)
169  MediumEnergyCascading(successful, pv, vecLength,
170  excitationEnergyGNP, excitationEnergyDTA,
171  incidentParticle, targetParticle,
172  atomicWeight, atomicNumber);
173 
174  if (!successful)
175  MediumEnergyClusterProduction(successful, pv, vecLength,
176  excitationEnergyGNP, excitationEnergyDTA,
177  incidentParticle, targetParticle,
178  atomicWeight, atomicNumber);
179  if (!successful)
180  QuasiElasticScattering(successful, pv, vecLength,
181  excitationEnergyGNP, excitationEnergyDTA,
182  incidentParticle, targetParticle,
183  atomicWeight, atomicNumber);
184  if (!successful)
185  ElasticScattering(successful, pv, vecLength,
186  incidentParticle,
187  atomicWeight, atomicNumber);
188 
189  if (!successful)
190  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
191  << G4endl;
192 
193  FillParticleChange(pv, vecLength);
194  delete [] pv;
195  theParticleChange.SetStatusChange(stopAndKill);
196  return &theParticleChange;
197 }
198 
199 
200 void
201 G4HEPionMinusInelastic::FirstIntInCasPionMinus(G4bool& inElastic,
202  const G4double availableEnergy,
203  G4HEVector pv[],
204  G4int& vecLen,
205  const G4HEVector& incidentParticle,
206  const G4HEVector& targetParticle)
207 
208 // Pion- undergoes interaction with nucleon within a nucleus. Check if it is
209 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
210 // occurs and input particle is degraded in energy. No other particles are produced.
211 // If reaction is possible, find the correct number of pions/protons/neutrons
212 // produced using an interpolation to multiplicity data. Replace some pions or
213 // protons/neutrons by kaons or strange baryons according to the average
214 // multiplicity per inelastic reaction.
215 {
216  static const G4double expxu = 82.; // upper bound for arg. of exp
217  static const G4double expxl = -expxu; // lower bound for arg. of exp
218 
219  static const G4double protb = 0.7;
220  static const G4double neutb = 0.7;
221  static const G4double c = 1.25;
222 
223  static const G4int numMul = 1200;
224  static const G4int numSec = 60;
225 
226  G4int protonCode = Proton.getCode();
227  G4int targetCode = targetParticle.getCode();
228  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
229 
230  static G4bool first = true;
231  static G4double protmul[numMul], protnorm[numSec]; // proton constants
232  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
233 
234  // misc. local variables
235  // npos = number of pi+, nneg = number of pi-, nero = number of pi0
236 
237  G4int i, counter, nt, npos, nneg, nero;
238 
239  if (first) {
240  // compute normalization constants, this will only be done once
241  first = false;
242  for( i=0; i<numMul; i++ )protmul[i] = 0.0;
243  for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
244  counter = -1;
245  for( npos=0; npos<(numSec/3); npos++ )
246  {
247  for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ )
248  {
249  for( nero=0; nero<numSec/3; nero++ )
250  {
251  if( ++counter < numMul )
252  {
253  nt = npos+nneg+nero;
254  if( (nt>0) && (nt<=numSec) )
255  {
256  protmul[counter] =
257  pmltpc(npos,nneg,nero,nt,protb,c) ;
258  protnorm[nt-1] += protmul[counter];
259  }
260  }
261  }
262  }
263  }
264  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
265  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
266  counter = -1;
267  for( npos=0; npos<numSec/3; npos++ )
268  {
269  for( nneg=npos; nneg<=(npos+2); nneg++ )
270  {
271  for( nero=0; nero<numSec/3; nero++ )
272  {
273  if( ++counter < numMul )
274  {
275  nt = npos+nneg+nero;
276  if( (nt>0) && (nt<=numSec) )
277  {
278  neutmul[counter] =
279  pmltpc(npos,nneg,nero,nt,neutb,c);
280  neutnorm[nt-1] += neutmul[counter];
281  }
282  }
283  }
284  }
285  }
286  for( i=0; i<numSec; i++ )
287  {
288  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
289  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
290  }
291  } // end of initialization
292 
293 
294  // initialize the first two places
295  // the same as beam and target
296  pv[0] = incidentParticle;
297  pv[1] = targetParticle;
298  vecLen = 2;
299 
300  if (!inElastic || (availableEnergy <= PionPlus.getMass()))
301  return;
302 
303 
304  // inelastic scattering
305 
306  npos = 0, nneg = 0, nero = 0;
307  G4double eab = availableEnergy;
308  G4int ieab = G4int( eab*5.0 );
309 
310  G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
311  if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
312  {
313  // suppress high multiplicity events at low momentum
314  // only one additional pion will be produced
315  G4double cech[] = {1., 0.95, 0.79, 0.32, 0.19, 0.16, 0.14, 0.12, 0.10, 0.08};
316  G4int iplab = Imin(9, G4int( incidentTotalMomentum*5.));
317  if( G4UniformRand() < cech[iplab] )
318  {
319  if( targetCode == protonCode )
320  {
321  pv[0] = PionZero;
322  pv[1] = Neutron;
323  return;
324  }
325  else
326  {
327  return;
328  }
329  }
330 
331  G4double w0, wp, wm, wt, ran;
332  if( targetCode == protonCode ) // target is a proton
333  {
334  w0 = - sqr(1.+protb)/(2.*c*c);
335  wp = w0 = std::exp(w0);
336  wm = - sqr(-1.+protb)/(2.*c*c);
337  wm = std::exp(wm);
338  wp *= 10.;
339  wt = w0+wp+wm;
340  wp = w0+wp;
341  ran = G4UniformRand();
342  if( ran < w0/wt )
343  { npos = 0; nneg = 0; nero = 1; }
344  else if ( ran < wp/wt )
345  { npos = 1; nneg = 0; nero = 0; }
346  else
347  { npos = 0; nneg = 1; nero = 0; }
348  }
349  else
350  { // target is a neutron
351  w0 = -sqr(1.+neutb)/(2.*c*c);
352  wm = -sqr(-1.+neutb)/(2.*c*c);
353  w0 = std::exp(w0);
354  wm = std::exp(wm);
355  if( G4UniformRand() < w0/(w0+wm) )
356  { npos = 0; nneg = 0; nero = 1; }
357  else
358  { npos = 0; nneg = 1; nero = 0; }
359  }
360  }
361  else
362  {
363  // number of total particles vs. centre of mass Energy - 2*proton mass
364 
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( Amin( expxu, Amax( 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  for (nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++) {
390  for (nero=0; nero<numSec/3; nero++) {
391  if (++counter < numMul) {
392  nt = npos+nneg+nero;
393  if ( (nt>0) && (nt<=numSec) ) {
394  test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
395  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
396  if (std::fabs(dum) < 1.0) {
397  if( test >= 1.0e-10 )excs += dum*test;
398  } else {
399  excs += dum*test;
400  }
401  if (ran < excs) goto outOfLoop; //----------------------->
402  }
403  }
404  }
405  }
406  }
407 
408  // 3 previous loops continued to the end
409  inElastic = false; // quasi-elastic scattering
410  return;
411  }
412  else
413  { // target must be a neutron
414  counter = -1;
415  for (npos=0; npos<numSec/3; npos++) {
416  for (nneg=npos; nneg<=(npos+2); nneg++) {
417  for (nero=0; nero<numSec/3; nero++) {
418  if (++counter < numMul) {
419  nt = npos+nneg+nero;
420  if ( (nt>=1) && (nt<=numSec) ) {
421  test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
422  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
423  if (std::fabs(dum) < 1.0) {
424  if( test >= 1.0e-10 )excs += dum*test;
425  } else {
426  excs += dum*test;
427  }
428  if (ran < excs) goto outOfLoop; // ----------------->
429  }
430  }
431  }
432  }
433  }
434  // 3 previous loops continued to the end
435  inElastic = false; // quasi-elastic scattering.
436  return;
437  }
438  }
439  outOfLoop: // <-----------------------------------------------
440 
441  if( targetCode == protonCode)
442  {
443  if( npos == (1+nneg))
444  {
445  pv[1] = Neutron;
446  }
447  else if (npos == nneg)
448  {
449  if( G4UniformRand() < 0.75)
450  {
451  }
452  else
453  {
454  pv[0] = PionZero;
455  pv[1] = Neutron;
456  }
457  }
458  else
459  {
460  pv[0] = PionZero;
461  }
462  }
463  else
464  {
465  if( npos == (nneg-1))
466  {
467  if( G4UniformRand() < 0.5)
468  {
469  pv[1] = Proton;
470  }
471  else
472  {
473  pv[0] = PionZero;
474  }
475  }
476  else if ( npos == nneg)
477  {
478  }
479  else
480  {
481  pv[0] = PionZero;
482  pv[1] = Proton;
483  }
484  }
485 
486 
487  nt = npos + nneg + nero;
488  while ( nt > 0)
489  {
490  G4double ran = G4UniformRand();
491  if ( ran < (G4double)npos/nt)
492  {
493  if( npos > 0 )
494  { pv[vecLen++] = PionPlus;
495  npos--;
496  }
497  }
498  else if ( ran < (G4double)(npos+nneg)/nt)
499  {
500  if( nneg > 0 )
501  {
502  pv[vecLen++] = PionMinus;
503  nneg--;
504  }
505  }
506  else
507  {
508  if( nero > 0 )
509  {
510  pv[vecLen++] = PionZero;
511  nero--;
512  }
513  }
514  nt = npos + nneg + nero;
515  }
516  if (verboseLevel > 1)
517  {
518  G4cout << "Particles produced: " ;
519  G4cout << pv[0].getName() << " " ;
520  G4cout << pv[1].getName() << " " ;
521  for (i=2; i < vecLen; i++)
522  {
523  G4cout << pv[i].getName() << " " ;
524  }
525  G4cout << G4endl;
526  }
527  return;
528  }
529 
530 
531 
532 
533 
534 
535 
536 
537 
Generated on Sat May 25 2013 14:33:56 for Geant4 by   doxygen 1.8.4