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