Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LEAntiProtonInelastic.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 // Hadronic Process: AntiProton Inelastic Process
29 // J.L. Chuma, TRIUMF, 13-Feb-1997
30 // J.P. Wellisch: 23-Apr-97: Bug hunting
31 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
32 //
33 
34 #include <iostream>
35 
36 #include "G4LEAntiProtonInelastic.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "Randomize.hh"
40 
41 
42 G4LEAntiProtonInelastic::G4LEAntiProtonInelastic(const G4String& name)
43  :G4InelasticInteraction(name)
44 {
45  SetMinEnergy(0.0);
46  SetMaxEnergy(25.*GeV);
47  G4cout << "WARNING: model G4LEAntiProtonInelastic is being deprecated and will\n"
48  << "disappear in Geant4 version 10.0" << G4endl;
49 }
50 
51 
52 void G4LEAntiProtonInelastic::ModelDescription(std::ostream& outFile) const
53 {
54  outFile << "G4LEAntiProtonInelastic is one of the Low Energy Parameterized\n"
55  << "(LEP) models used to implement inelastic anti-proton scattering\n"
56  << "from nuclei. It is a re-engineered version of the GHEISHA\n"
57  << "code of H. Fesefeldt. It divides the initial collision\n"
58  << "products into backward- and forward-going clusters which are\n"
59  << "then decayed into final state hadrons. The model does not\n"
60  << "conserve energy on an event-by-event basis. It may be\n"
61  << "applied to anti-protons with initial energies between 0 and 25\n"
62  << "GeV.\n";
63 }
64 
65 
66 G4HadFinalState*
67 G4LEAntiProtonInelastic::ApplyYourself(const G4HadProjectile& aTrack,
68  G4Nucleus& targetNucleus)
69 {
70  const G4HadProjectile* originalIncident = &aTrack;
71 
72  if (originalIncident->GetKineticEnergy() <= 0.1*MeV) {
73  theParticleChange.SetStatusChange(isAlive);
74  theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
75  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
76  return &theParticleChange;
77  }
78 
79  // create the target particle
80 
81  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
82 
83  if (verboseLevel > 1) {
84  const G4Material *targetMaterial = aTrack.GetMaterial();
85  G4cout << "G4LEAntiProtonInelastic::ApplyYourself called" << G4endl;
86  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
87  G4cout << "target material = " << targetMaterial->GetName() << ", ";
88  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
89  << G4endl;
90  }
91 
92  // Fermi motion and evaporation
93  // As of Geant3, the Fermi energy calculation had not been Done
94 
95  G4double ek = originalIncident->GetKineticEnergy()/MeV;
96  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
97  G4ReactionProduct modifiedOriginal;
98  modifiedOriginal = *originalIncident;
99 
100  G4double tkin = targetNucleus.Cinema( ek );
101  ek += tkin;
102  modifiedOriginal.SetKineticEnergy( ek*MeV );
103  G4double et = ek + amas;
104  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
105  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
106  if (pp > 0.0) {
107  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
108  modifiedOriginal.SetMomentum( momentum * (p/pp) );
109  }
110 
111  // calculate black track energies
112  tkin = targetNucleus.EvaporationEffects( ek );
113  ek -= tkin;
114  modifiedOriginal.SetKineticEnergy( ek*MeV );
115  et = ek + amas;
116  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
117  pp = modifiedOriginal.GetMomentum().mag()/MeV;
118  if (pp > 0.0) {
119  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
120  modifiedOriginal.SetMomentum( momentum * (p/pp) );
121  }
122  G4ReactionProduct currentParticle = modifiedOriginal;
123  G4ReactionProduct targetParticle;
124  targetParticle = *originalTarget;
125  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
126  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
127  G4bool incidentHasChanged = false;
128  G4bool targetHasChanged = false;
129  G4bool quasiElastic = false;
130  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
131  G4int vecLen = 0;
132  vec.Initialize( 0 );
133 
134  const G4double cutOff = 0.1;
135  const G4double anni = std::min( 1.3*originalIncident->GetTotalMomentum()/GeV, 0.4 );
136 
137  if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
138  (G4UniformRand() > anni) )
139  Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
140  incidentHasChanged, targetHasChanged, quasiElastic);
141  else
142  quasiElastic = true;
143 
144  CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
145  modifiedOriginal, targetNucleus, currentParticle,
146  targetParticle, incidentHasChanged, targetHasChanged,
147  quasiElastic);
148 
149  SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
150 
151  if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
152 
153  delete originalTarget;
154  return &theParticleChange;
155 }
156 
157 
158 void G4LEAntiProtonInelastic::Cascade(
159  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
160  G4int& vecLen,
161  const G4HadProjectile *originalIncident,
162  G4ReactionProduct &currentParticle,
163  G4ReactionProduct &targetParticle,
164  G4bool &incidentHasChanged,
165  G4bool &targetHasChanged,
166  G4bool &quasiElastic )
167 {
168  // derived from original FORTRAN code CASPB by H. Fesefeldt (13-Sep-1987)
169  //
170  // AntiProton undergoes interaction with nucleon within a nucleus. Check if
171  // it is energetically possible to produce pions/kaons. In not, assume
172  // nuclear excitation occurs and input particle is degraded in energy. No
173  // other particles are produced. If reaction is possible, find the correct
174  // number of pions/protons/neutrons produced using an interpolation to
175  // multiplicity data. Replace some pions or protons/neutrons by kaons or
176  // strange baryons according to the average multiplicity per inelastic
177  // reaction.
178 
179  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
180  const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
181  const G4double pOriginal = originalIncident->GetTotalMomentum()/MeV;
182  const G4double targetMass = targetParticle.GetMass()/MeV;
183  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
184  targetMass*targetMass +
185  2.0*targetMass*etOriginal );
186  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
187 
188  static G4bool first = true;
189  const G4int numMul = 1200;
190  const G4int numMulA = 400;
191  const G4int numSec = 60;
192  static G4double protmul[numMul], protnorm[numSec]; // proton constants
193  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
194  static G4double protmulA[numMulA], protnormA[numSec]; // proton constants
195  static G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
196 
197  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
198  G4int counter, nt=0;
199  G4int npos=0, nneg = 0, nzero=0;
200  G4double test;
201  const G4double c = 1.25;
202  const G4double b[] = { 0.7, 0.7 };
203  if (first) { // Computation of normalization constants will only be done once
204  first = false;
205  G4int i;
206  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
207  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
208  counter = -1;
209  for (npos = 0; npos < (numSec/3); ++npos) {
210  for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
211  for (nzero = 0; nzero<numSec/3; ++nzero) {
212  if ( ++counter < numMul ) {
213  nt = npos+nneg+nzero;
214  if (nt > 0 && nt <= numSec) {
215  protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
216  protnorm[nt-1] += protmul[counter];
217  }
218  }
219  }
220  }
221  }
222 
223  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
224  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
225  counter = -1;
226  for (npos = 0; npos < numSec/3; ++npos) {
227  for (nneg = npos; nneg <= (npos+2); ++nneg) {
228  for (nzero = 0; nzero < numSec/3; ++nzero) {
229  if (++counter < numMul) {
230  nt = npos+nneg+nzero;
231  if ((nt > 0) && (nt <= numSec) ) {
232  neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
233  neutnorm[nt-1] += neutmul[counter];
234  }
235  }
236  }
237  }
238  }
239  for (i = 0; i < numSec; ++i) {
240  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
241  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
242  }
243 
244  // do the same for annihilation channels
245  for (i = 0; i < numMulA; ++i) protmulA[i] = 0.0;
246  for (i = 0; i < numSec; ++i) protnormA[i] = 0.0;
247  counter = -1;
248  for (npos = 0; npos < (numSec/3); ++npos) {
249  nneg = npos;
250  for (nzero=0; nzero<numSec/3; ++nzero) {
251  if (++counter < numMulA) {
252  nt = npos+nneg+nzero;
253  if (nt > 1 && nt <= numSec) {
254  protmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
255  protnormA[nt-1] += protmulA[counter];
256  }
257  }
258  }
259  }
260  for (i = 0; i < numMulA; ++i) neutmulA[i] = 0.0;
261  for (i = 0; i < numSec; ++i) neutnormA[i] = 0.0;
262  counter = -1;
263  for (npos=0; npos < numSec/3; ++npos) {
264  nneg = npos+1;
265  for( nzero=0; nzero<numSec/3; ++nzero ) {
266  if( ++counter < numMulA ) {
267  nt = npos+nneg+nzero;
268  if( (nt>1) && (nt<=numSec) ) {
269  neutmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
270  neutnormA[nt-1] += neutmulA[counter];
271  }
272  }
273  }
274  }
275  for (i=0; i<numSec; ++i ) {
276  if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
277  if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
278  }
279  } // end of initialization
280 
281  const G4double expxu = 82.; // upper bound for arg. of exp
282  const G4double expxl = -expxu; // lower bound for arg. of exp
283 
284  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
285  G4ParticleDefinition *aProton = G4Proton::Proton();
286  G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
287  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
288 
289  const G4double anhl[] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.90,
290  0.6,0.52,0.47,0.44,0.41,0.39,0.37,0.35,0.34,0.24,
291  0.19,0.15,0.12,0.10,0.09,0.07,0.06,0.05,0.0};
292 
293  G4int iplab = G4int( pOriginal/GeV*10.0 );
294  if( iplab > 9 )iplab = G4int( pOriginal/GeV ) + 9;
295  if( iplab > 18 )iplab = G4int( pOriginal/GeV/10.0 ) + 18;
296  if( iplab > 27 )iplab = 28;
297  if (G4UniformRand() > anhl[iplab]) {
298  if (availableEnergy <= aPiPlus->GetPDGMass()/MeV) {
299  quasiElastic = true;
300  return;
301  }
302  G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
303  const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
304  G4double w0, wp, wt, wm;
305  if ((availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) ) {
306  // suppress high multiplicity events at low momentum
307  // only one pion will be produced
308  npos = nneg = nzero = 0;
309  if (targetParticle.GetDefinition() == aProton) {
310  test = std::exp(std::min(expxu,
311  std::max(expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
312  w0 = test;
313  wp = test;
314  test = std::exp(std::min(expxu,
315  std::max(expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
316  wm = test;
317  wt = w0+wp+wm;
318  wp += w0;
319  G4double ran = G4UniformRand();
320  if( ran < w0/wt )
321  nzero = 1;
322  else if( ran < wp/wt )
323  npos = 1;
324  else
325  nneg = 1;
326  } else {
327  // target is a neutron
328  test = std::exp(std::min(expxu,
329  std::max(expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
330  w0 = test;
331  test = std::exp(std::min(expxu,
332  std::max(expxl, -(-1.0+b[0])*(-1.0+b[0])/(2.0*c*c) ) ) );
333  wm = test;
334  G4double ran = G4UniformRand();
335  if (ran < w0/(w0+wm) )
336  nzero = 1;
337  else
338  nneg = 1;
339  }
340  } else {
341  // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
342  G4double n, anpn;
343  GetNormalizationConstant( availableEnergy, n, anpn );
344  G4double ran = G4UniformRand();
345  G4double dum, excs = 0.0;
346  if (targetParticle.GetDefinition() == aProton) {
347  counter = -1;
348  for( npos=0; npos<numSec/3 && ran>=excs; ++npos ) {
349  for( nneg=std::max(0,npos-1); nneg<=(npos+1) && ran>=excs; ++nneg) {
350  for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
351  {
352  if( ++counter < numMul )
353  {
354  nt = npos+nneg+nzero;
355  if( (nt>0) && (nt<=numSec) )
356  {
357  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
358  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
359  if( std::fabs(dum) < 1.0 )
360  {
361  if( test >= 1.0e-10 )excs += dum*test;
362  }
363  else
364  excs += dum*test;
365  }
366  }
367  }
368  }
369  }
370  if( ran >= excs ) // 3 previous loops continued to the end
371  {
372  quasiElastic = true;
373  return;
374  }
375  }
376  else // target must be a neutron
377  {
378  counter = -1;
379  for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
380  {
381  for (nneg = npos; nneg <= (npos+2) && ran >= excs; ++nneg) {
382  for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
383  {
384  if( ++counter < numMul )
385  {
386  nt = npos+nneg+nzero;
387  if( (nt>0) && (nt<=numSec) )
388  {
389  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
390  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
391  if( std::fabs(dum) < 1.0 )
392  {
393  if( test >= 1.0e-10 )excs += dum*test;
394  }
395  else
396  excs += dum*test;
397  }
398  }
399  }
400  }
401  }
402  if( ran >= excs ) // 3 previous loops continued to the end
403  {
404  quasiElastic = true;
405  return;
406  }
407  }
408  npos--; nneg--; nzero--;
409  }
410  if( targetParticle.GetDefinition() == aProton )
411  {
412  switch (npos-nneg)
413  {
414  case 0:
415  if( G4UniformRand() < 0.33 )
416  {
417  currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
418  targetParticle.SetDefinitionAndUpdateE( aNeutron );
419  incidentHasChanged = true;
420  targetHasChanged = true;
421  }
422  break;
423  case 1:
424  targetParticle.SetDefinitionAndUpdateE( aNeutron );
425  targetHasChanged = true;
426  break;
427  default:
428  currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
429  incidentHasChanged = true;
430  break;
431  }
432  }
433  else // target must be a neutron
434  {
435  switch (npos-nneg)
436  {
437  case -1:
438  if( G4UniformRand() < 0.5 )
439  {
440  targetParticle.SetDefinitionAndUpdateE( aProton );
441  targetHasChanged = true;
442  }
443  else
444  {
445  currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
446  incidentHasChanged = true;
447  }
448  break;
449  case 0:
450  break;
451  default:
452  currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
453  targetParticle.SetDefinitionAndUpdateE( aProton );
454  incidentHasChanged = true;
455  targetHasChanged = true;
456  break;
457  }
458  }
459  }
460  else // random number <= anhl[iplab]
461  {
462  if( centerofmassEnergy <= 2*aPiPlus->GetPDGMass()/MeV )
463  {
464  quasiElastic = true;
465  return;
466  }
467  //
468  // annihilation channels
469  //
470  G4double n, anpn;
471  GetNormalizationConstant( -centerofmassEnergy, n, anpn );
472  G4double ran = G4UniformRand();
473  G4double dum, excs = 0.0;
474  if( targetParticle.GetDefinition() == aProton )
475  {
476  counter = -1;
477  for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos )
478  {
479  nneg = npos;
480  for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero )
481  {
482  if( ++counter < numMulA )
483  {
484  nt = npos+nneg+nzero;
485  if( (nt>1) && (nt<=numSec) )
486  {
487  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
488  dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
489  if( std::abs(dum) < 1.0 )
490  {
491  if( test >= 1.0e-10 )excs += dum*test;
492  }
493  else
494  excs += dum*test;
495  }
496  }
497  }
498  }
499  }
500  else // target must be a neutron
501  {
502  counter = -1;
503  for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos )
504  {
505  nneg = npos+1;
506  for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero )
507  {
508  if( ++counter < numMulA )
509  {
510  nt = npos+nneg+nzero;
511  if( (nt>1) && (nt<=numSec) )
512  {
513  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
514  dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
515  if( std::fabs(dum) < 1.0 )
516  {
517  if( test >= 1.0e-10 )excs += dum*test;
518  }
519  else
520  excs += dum*test;
521  }
522  }
523  }
524  }
525  }
526  if (ran >= excs) {
527  // 3 previous loops continued to the end
528  quasiElastic = true;
529  return;
530  }
531  npos--; nzero--;
532  currentParticle.SetMass( 0.0 );
533  targetParticle.SetMass( 0.0 );
534  }
535 
536  while(npos+nneg+nzero<3) nzero++;
537  SetUpPions( npos, nneg, nzero, vec, vecLen );
538  return;
539 }
Generated on Sat May 25 2013 14:34:01 for Geant4 by   doxygen 1.8.4