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