ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGKPlusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGKPlusInelastic.cc
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 //
27 //
28 
29 #include "G4RPGKPlusInelastic.hh"
30 #include "G4Exp.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "Randomize.hh"
34 
37  G4Nucleus &targetNucleus )
38 {
39  const G4HadProjectile *originalIncident = &aTrack;
40  if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
41  {
45  return &theParticleChange;
46  }
47 
48  // create the target particle
49 
50  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
51  G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
52 
53  if( verboseLevel > 1 )
54  {
55  const G4Material *targetMaterial = aTrack.GetMaterial();
56  G4cout << "G4RPGKPlusInelastic::ApplyYourself called" << G4endl;
57  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
58  G4cout << "target material = " << targetMaterial->GetName() << ", ";
59  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
60  << G4endl;
61  }
62 
63  G4ReactionProduct currentParticle(originalIncident->GetDefinition() );
64  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
65  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
66 
67  // Fermi motion and evaporation
68  // As of Geant3, the Fermi energy calculation had not been Done
69 
70  G4double ek = originalIncident->GetKineticEnergy();
71  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
72 
73  G4double tkin = targetNucleus.Cinema( ek );
74  ek += tkin;
75  currentParticle.SetKineticEnergy( ek );
76  G4double et = ek + amas;
77  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
78  G4double pp = currentParticle.GetMomentum().mag();
79  if( pp > 0.0 )
80  {
81  G4ThreeVector momentum = currentParticle.GetMomentum();
82  currentParticle.SetMomentum( momentum * (p/pp) );
83  }
84 
85  // calculate black track energies
86 
87  tkin = targetNucleus.EvaporationEffects( ek );
88  ek -= tkin;
89  currentParticle.SetKineticEnergy( ek );
90  et = ek + amas;
91  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
92  pp = currentParticle.GetMomentum().mag();
93  if( pp > 0.0 )
94  {
95  G4ThreeVector momentum = currentParticle.GetMomentum();
96  currentParticle.SetMomentum( momentum * (p/pp) );
97  }
98 
99  G4ReactionProduct modifiedOriginal = currentParticle;
100 
101  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
102  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
103  G4bool incidentHasChanged = false;
104  G4bool targetHasChanged = false;
105  G4bool quasiElastic = false;
106  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
107  G4int vecLen = 0;
108  vec.Initialize( 0 );
109 
110  const G4double cutOff = 0.1*MeV;
111  if( currentParticle.GetKineticEnergy() > cutOff )
112  Cascade( vec, vecLen,
113  originalIncident, currentParticle, targetParticle,
114  incidentHasChanged, targetHasChanged, quasiElastic );
115 
116  CalculateMomenta( vec, vecLen,
117  originalIncident, originalTarget, modifiedOriginal,
118  targetNucleus, currentParticle, targetParticle,
119  incidentHasChanged, targetHasChanged, quasiElastic );
120 
121  SetUpChange( vec, vecLen,
122  currentParticle, targetParticle,
123  incidentHasChanged );
124 
125  delete originalTarget;
126 
127  return &theParticleChange;
128 }
129 
130 
133  G4int &vecLen,
134  const G4HadProjectile *originalIncident,
135  G4ReactionProduct &currentParticle,
136  G4ReactionProduct &targetParticle,
137  G4bool &incidentHasChanged,
138  G4bool &targetHasChanged,
139  G4bool &quasiElastic )
140 {
141  // Derived from H. Fesefeldt's original FORTRAN code CASKP
142  //
143  // K+ undergoes interaction with nucleon within a nucleus. Check if it is
144  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
145  // occurs and input particle is degraded in energy. No other particles are produced.
146  // If reaction is possible, find the correct number of pions/protons/neutrons
147  // produced using an interpolation to multiplicity data. Replace some pions or
148  // protons/neutrons by kaons or strange baryons according to the average
149  // multiplicity per Inelastic reaction.
150  //
151  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
152  const G4double etOriginal = originalIncident->GetTotalEnergy();
153  const G4double targetMass = targetParticle.GetMass();
154  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
155  targetMass*targetMass +
156  2.0*targetMass*etOriginal );
157  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
158  if( availableEnergy < G4PionPlus::PionPlus()->GetPDGMass() )
159  {
160  quasiElastic = true;
161  return;
162  }
163  static G4ThreadLocal G4bool first = true;
164  const G4int numMul = 1200;
165  const G4int numSec = 60;
166  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
167  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
168 
169  // np = number of pi+, nneg = number of pi-, nz = number of pi0
170 
171  G4int nt=0, np=0, nneg=0, nz=0;
172  const G4double c = 1.25;
173  const G4double b[] = { 0.70, 0.70 };
174  if( first ) // compute normalization constants, this will only be Done once
175  {
176  first = false;
177  G4int i;
178  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
179  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
180  G4int counter = -1;
181  for( np=0; np<(numSec/3); ++np )
182  {
183  for( nneg=std::max(0,np-2); nneg<=np; ++nneg )
184  {
185  for( nz=0; nz<numSec/3; ++nz )
186  {
187  if( ++counter < numMul )
188  {
189  nt = np+nneg+nz;
190  if( nt > 0 )
191  {
192  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
193  protnorm[nt-1] += protmul[counter];
194  }
195  }
196  }
197  }
198  }
199  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
200  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
201  counter = -1;
202  for( np=0; np<numSec/3; ++np )
203  {
204  for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
205  {
206  for( nz=0; nz<numSec/3; ++nz )
207  {
208  if( ++counter < numMul )
209  {
210  nt = np+nneg+nz;
211  if( (nt>0) && (nt<=numSec) )
212  {
213  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
214  neutnorm[nt-1] += neutmul[counter];
215  }
216  }
217  }
218  }
219  }
220  for( i=0; i<numSec; ++i )
221  {
222  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
223  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
224  }
225  } // end of initialization
226 
227  const G4double expxu = 82.; // upper bound for arg. of exp
228  const G4double expxl = -expxu; // lower bound for arg. of exp
233  G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
234  const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
235  G4double test, w0, wp, wt, wm;
236  if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
237  {
238  // suppress high multiplicity events at low momentum
239  // only one pion will be produced
240 
241  nneg = np = nz = 0;
242  if( targetParticle.GetDefinition() == aProton )
243  {
244  test = G4Exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[0])/(2.0*c*c) ) ) );
245  w0 = test;
246  wp = test*2.0;
247  if( G4UniformRand() < w0/(w0+wp) )
248  nz = 1;
249  else
250  np = 1;
251  }
252  else // target is a neutron
253  {
254  test = G4Exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[1])/(2.0*c*c) ) ) );
255  w0 = test;
256  wp = test;
257  test = G4Exp( std::min( expxu, std::max( expxl, -sqr(-1.0+b[1])/(2.0*c*c) ) ) );
258  wm = test;
259  wt = w0+wp+wm;
260  wp += w0;
262  if( ran < w0/wt )
263  nz = 1;
264  else if( ran < wp/wt )
265  np = 1;
266  else
267  nneg = 1;
268  }
269  }
270  else
271  {
272  G4double n, anpn;
273  GetNormalizationConstant( availableEnergy, n, anpn );
275  G4double dum, excs = 0.0;
276  if( targetParticle.GetDefinition() == aProton )
277  {
278  G4int counter = -1;
279  for( np=0; (np<numSec/3) && (ran>=excs); ++np )
280  {
281  for( nneg=std::max(0,np-2); (nneg<=np) && (ran>=excs); ++nneg )
282  {
283  for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz )
284  {
285  if( ++counter < numMul )
286  {
287  nt = np+nneg+nz;
288  if( nt > 0 )
289  {
290  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
291  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
292  if( std::fabs(dum) < 1.0 )
293  {
294  if( test >= 1.0e-10 )excs += dum*test;
295  }
296  else
297  excs += dum*test;
298  }
299  }
300  }
301  }
302  }
303  if( ran >= excs )return; // 3 previous loops continued to the end
304  np--; nneg--; nz--;
305  }
306  else // target must be a neutron
307  {
308  G4int counter = -1;
309  for( np=0; (np<numSec/3) && (ran>=excs); ++np )
310  {
311  for( nneg=std::max(0,np-1); (nneg<=(np+1)) && (ran>=excs); ++nneg )
312  {
313  for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz )
314  {
315  if( ++counter < numMul )
316  {
317  nt = np+nneg+nz;
318  if( (nt>=1) && (nt<=numSec) )
319  {
320  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
321  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
322  if( std::fabs(dum) < 1.0 )
323  {
324  if( test >= 1.0e-10 )excs += dum*test;
325  }
326  else
327  excs += dum*test;
328  }
329  }
330  }
331  }
332  }
333  if( ran >= excs )return; // 3 previous loops continued to the end
334  np--; nneg--; nz--;
335  }
336  }
337 
338  if( targetParticle.GetDefinition() == aProton )
339  {
340  switch( np-nneg )
341  {
342  case 1:
343  if( G4UniformRand() < 0.5 )
344  {
345  if( G4UniformRand() < 0.5 )
346  currentParticle.SetDefinitionAndUpdateE( aKaonZS );
347  else
348  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
349  incidentHasChanged = true;
350  }
351  else
352  {
353  targetParticle.SetDefinitionAndUpdateE( aNeutron );
354  targetHasChanged = true;
355  }
356  break;
357  case 2:
358  if( G4UniformRand() < 0.5 )
359  currentParticle.SetDefinitionAndUpdateE( aKaonZS );
360  else
361  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
362  incidentHasChanged = true;
363  targetParticle.SetDefinitionAndUpdateE( aNeutron );
364  incidentHasChanged = true;
365  targetHasChanged = true;
366  break;
367  default:
368  break;
369  }
370  }
371  else // target is a neutron
372  {
373  switch( np-nneg )
374  {
375  case 0:
376  if( G4UniformRand() < 0.25 )
377  {
378  if( G4UniformRand() < 0.5 )
379  currentParticle.SetDefinitionAndUpdateE( aKaonZS );
380  else
381  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
382  targetParticle.SetDefinitionAndUpdateE( aProton );
383  incidentHasChanged = true;
384  targetHasChanged = true;
385  }
386  break;
387  case 1:
388  if( G4UniformRand() < 0.5 )
389  currentParticle.SetDefinitionAndUpdateE( aKaonZS );
390  else
391  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
392  incidentHasChanged = true;
393  break;
394  default: // assumes nneg = np+1 so charge is conserved
395  targetParticle.SetDefinitionAndUpdateE( aProton );
396  targetHasChanged = true;
397  break;
398  }
399  }
400 
401  SetUpPions(np, nneg, nz, vec, vecLen);
402  return;
403 }
404 
405  /* end of file */
406