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