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