ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGKMinusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGKMinusInelastic.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 "G4RPGKMinusInelastic.hh"
29 #include "G4Exp.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "Randomize.hh"
33 
34 G4HadFinalState*
35 G4RPGKMinusInelastic::ApplyYourself( const G4HadProjectile &aTrack,
36  G4Nucleus &targetNucleus )
37 {
38  const G4HadProjectile *originalIncident = &aTrack;
39  if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
40  {
41  theParticleChange.SetStatusChange(isAlive);
42  theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
43  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
44  return &theParticleChange;
45  }
46 
47  // create the target particle
48 
49  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
50  G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
51 
52  if( verboseLevel > 1 )
53  {
54  const G4Material *targetMaterial = aTrack.GetMaterial();
55  G4cout << "G4RPGKMinusInelastic::ApplyYourself called" << G4endl;
56  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
57  G4cout << "target material = " << targetMaterial->GetName() << ", ";
58  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
59  << G4endl;
60  }
61 
62  G4ReactionProduct currentParticle(originalIncident->GetDefinition() );
63  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
64  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
65 
66  // Fermi motion and evaporation
67  // As of Geant3, the Fermi energy calculation had not been Done
68 
69  G4double ek = originalIncident->GetKineticEnergy();
70  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
71 
72  G4double tkin = targetNucleus.Cinema( ek );
73  ek += tkin;
74  currentParticle.SetKineticEnergy( ek );
75  G4double et = ek + amas;
76  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
77  G4double pp = currentParticle.GetMomentum().mag();
78  if( pp > 0.0 )
79  {
80  G4ThreeVector momentum = currentParticle.GetMomentum();
81  currentParticle.SetMomentum( momentum * (p/pp) );
82  }
83 
84  // calculate black track energies
85 
86  tkin = targetNucleus.EvaporationEffects( ek );
87  ek -= tkin;
88  currentParticle.SetKineticEnergy( ek );
89  et = ek + amas;
90  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
91  pp = currentParticle.GetMomentum().mag();
92  if( pp > 0.0 )
93  {
94  G4ThreeVector momentum = currentParticle.GetMomentum();
95  currentParticle.SetMomentum( momentum * (p/pp) );
96  }
97 
98  G4ReactionProduct modifiedOriginal = currentParticle;
99 
100  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
101  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
102  G4bool incidentHasChanged = false;
103  G4bool targetHasChanged = false;
104  G4bool quasiElastic = false;
105  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
106  G4int vecLen = 0;
107  vec.Initialize( 0 );
108 
109  const G4double cutOff = 0.1*MeV;
110  if( currentParticle.GetKineticEnergy() > cutOff )
111  Cascade( vec, vecLen,
112  originalIncident, currentParticle, targetParticle,
113  incidentHasChanged, targetHasChanged, quasiElastic );
114 
115  CalculateMomenta( vec, vecLen,
116  originalIncident, originalTarget, modifiedOriginal,
117  targetNucleus, currentParticle, targetParticle,
118  incidentHasChanged, targetHasChanged, quasiElastic );
119 
120  SetUpChange( vec, vecLen,
121  currentParticle, targetParticle,
122  incidentHasChanged );
123 
124  delete originalTarget;
125  return &theParticleChange;
126 }
127 
128 
129 void G4RPGKMinusInelastic::Cascade(
130  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
131  G4int& vecLen,
132  const G4HadProjectile *originalIncident,
133  G4ReactionProduct &currentParticle,
134  G4ReactionProduct &targetParticle,
135  G4bool &incidentHasChanged,
136  G4bool &targetHasChanged,
137  G4bool &quasiElastic )
138 {
139  // Derived from H. Fesefeldt's original FORTRAN code CASKM
140  //
141  // K- undergoes interaction with nucleon within a nucleus. Check if it is
142  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
143  // occurs and input particle is degraded in energy. No other particles are produced.
144  // If reaction is possible, find the correct number of pions/protons/neutrons
145  // produced using an interpolation to multiplicity data. Replace some pions or
146  // protons/neutrons by kaons or strange baryons according to the average
147  // multiplicity per Inelastic reaction.
148  //
149  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
150  const G4double etOriginal = originalIncident->GetTotalEnergy();
151  const G4double pOriginal = originalIncident->GetTotalMomentum();
152  const G4double targetMass = targetParticle.GetMass();
153  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
154  targetMass*targetMass +
155  2.0*targetMass*etOriginal );
156  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
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 nt(0), np(0), nneg(0), nz(0);
165  const G4double c = 1.25;
166  const G4double b[] = { 0.70, 0.70 };
167  if( first ) // compute normalization constants, this will only be Done once
168  {
169  first = false;
170  G4int i;
171  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
172  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
173  G4int counter = -1;
174  for( np=0; np<(numSec/3); ++np )
175  {
176  for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
177  {
178  for( nz=0; nz<numSec/3; ++nz )
179  {
180  if( ++counter < numMul )
181  {
182  nt = np+nneg+nz;
183  if( (nt>0) && (nt<=numSec) )
184  {
185  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
186  protnorm[nt-1] += protmul[counter];
187  }
188  }
189  }
190  }
191  }
192  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
193  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
194  counter = -1;
195  for( np=0; np<numSec/3; ++np )
196  {
197  for( nneg=np; nneg<=(np+2); ++nneg )
198  {
199  for( nz=0; nz<numSec/3; ++nz )
200  {
201  if( ++counter < numMul )
202  {
203  nt = np+nneg+nz;
204  if( (nt>0) && (nt<=numSec) )
205  {
206  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
207  neutnorm[nt-1] += neutmul[counter];
208  }
209  }
210  }
211  }
212  }
213  for( i=0; i<numSec; ++i )
214  {
215  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
216  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
217  }
218  } // end of initialization
219 
220  const G4double expxu = 82.; // upper bound for arg. of exp
221  const G4double expxl = -expxu; // lower bound for arg. of exp
222  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
223  G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
224  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
225  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
226  G4ParticleDefinition *aProton = G4Proton::Proton();
227  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
228  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
229  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
230  G4ParticleDefinition *aLambda = G4Lambda::Lambda();
231  G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
232  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
233  G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
234  const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
235  G4int iplab = G4int(std::min( 9.0, pOriginal/GeV*5.0 ));
236  if( (pOriginal <= 2.0*GeV) && (G4UniformRand() < cech[iplab]) )
237  {
238  np = nneg = nz = nt = 0;
239  iplab = G4int(std::min( 19.0, pOriginal/GeV*10.0 ));
240  const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
241  0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
242  if( G4UniformRand() <= cnk0[iplab] )
243  {
244  quasiElastic = true;
245  if( targetParticle.GetDefinition() == aProton )
246  {
247  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
248  incidentHasChanged = true;
249  targetParticle.SetDefinitionAndUpdateE( aNeutron );
250  targetHasChanged = true;
251  }
252  }
253  else // random number > cnk0[iplab]
254  {
255  G4double ran = G4UniformRand();
256  if( ran < 0.25 ) // k- p --> pi- s+
257  {
258  if( targetParticle.GetDefinition() == aProton )
259  {
260  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
261  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
262  incidentHasChanged = true;
263  targetHasChanged = true;
264  }
265  }
266  else if( ran < 0.50 ) // k- p --> pi0 s0 or k- n --> pi- s0
267  {
268  if( targetParticle.GetDefinition() == aNeutron )
269  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
270  else
271  currentParticle.SetDefinitionAndUpdateE( aPiZero );
272  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
273  incidentHasChanged = true;
274  targetHasChanged = true;
275  }
276  else if( ran < 0.75 ) // k- p --> pi+ s- or k- n --> pi0 s-
277  {
278  if( targetParticle.GetDefinition() == aNeutron )
279  currentParticle.SetDefinitionAndUpdateE( aPiZero );
280  else
281  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
282  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
283  incidentHasChanged = true;
284  targetHasChanged = true;
285  }
286  else // k- p --> pi0 L or k- n --> pi- L
287  {
288  if( targetParticle.GetDefinition() == aNeutron )
289  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
290  else
291  currentParticle.SetDefinitionAndUpdateE( aPiZero );
292  targetParticle.SetDefinitionAndUpdateE( aLambda );
293  incidentHasChanged = true;
294  targetHasChanged = true;
295  }
296  }
297  }
298  else // (pOriginal > 2.0*GeV) || (random number >= cech[iplab])
299  {
300  if( availableEnergy < aPiPlus->GetPDGMass() )
301  { // not energetically possible to produce pion(s)
302  quasiElastic = true;
303  return;
304  }
305  G4double n, anpn;
306  GetNormalizationConstant( availableEnergy, n, anpn );
307  G4double ran = G4UniformRand();
308  G4double dum, test, excs = 0.0;
309  if( targetParticle.GetDefinition() == aProton )
310  {
311  G4int counter = -1;
312  for( np=0; np<numSec/3 && ran>=excs; ++np )
313  {
314  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
315  {
316  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
317  {
318  if( ++counter < numMul )
319  {
320  nt = np+nneg+nz;
321  if( nt > 0 )
322  {
323  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
324  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
325  if( std::fabs(dum) < 1.0 )
326  {
327  if( test >= 1.0e-10 )excs += dum*test;
328  }
329  else
330  excs += dum*test;
331  }
332  }
333  }
334  }
335  }
336  if( ran >= excs ) // 3 previous loops continued to the end
337  {
338  quasiElastic = true;
339  return;
340  }
341  np--; nneg--; nz--;
342  if( np == nneg )
343  {
344  if( G4UniformRand() >= 0.75 )
345  {
346  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
347  targetParticle.SetDefinitionAndUpdateE( aNeutron );
348  incidentHasChanged = true;
349  targetHasChanged = true;
350  }
351  }
352  else if( np == nneg+1 )
353  {
354  targetParticle.SetDefinitionAndUpdateE( aNeutron );
355  targetHasChanged = true;
356  }
357  else
358  {
359  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
360  incidentHasChanged = true;
361  }
362  }
363  else // target must be a neutron
364  {
365  G4int counter = -1;
366  for( np=0; np<numSec/3 && ran>=excs; ++np )
367  {
368  for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
369  {
370  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
371  {
372  if( ++counter < numMul )
373  {
374  nt = np+nneg+nz;
375  if( (nt>=1) && (nt<=numSec) )
376  {
377  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
378  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
379  if( std::fabs(dum) < 1.0 )
380  {
381  if( test >= 1.0e-10 )excs += dum*test;
382  }
383  else
384  excs += dum*test;
385  }
386  }
387  }
388  }
389  }
390  if( ran >= excs ) // 3 previous loops continued to the end
391  {
392  quasiElastic = true;
393  return;
394  }
395  np--; nneg--; nz--;
396  if( np == nneg-1 )
397  {
398  if( G4UniformRand() < 0.5 )
399  {
400  targetParticle.SetDefinitionAndUpdateE( aProton );
401  targetHasChanged = true;
402  }
403  else
404  {
405  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
406  incidentHasChanged = true;
407  }
408  }
409  else if( np != nneg )
410  {
411  currentParticle.SetDefinitionAndUpdateE( aKaonZL );
412  incidentHasChanged = true;
413  }
414  }
415  if( G4UniformRand() >= 0.5 )
416  {
417  if( (currentParticle.GetDefinition() == aKaonMinus &&
418  targetParticle.GetDefinition() == aNeutron ) ||
419  (currentParticle.GetDefinition() == aKaonZL &&
420  targetParticle.GetDefinition() == aProton ) )
421  {
422  ran = G4UniformRand();
423  if( ran < 0.68 )
424  {
425  if( targetParticle.GetDefinition() == aProton )
426  {
427  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
428  targetParticle.SetDefinitionAndUpdateE( aLambda );
429  incidentHasChanged = true;
430  targetHasChanged = true;
431  }
432  else
433  {
434  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
435  targetParticle.SetDefinitionAndUpdateE( aLambda );
436  incidentHasChanged = true;
437  targetHasChanged = true;
438  }
439  }
440  else if( ran < 0.84 )
441  {
442  if( targetParticle.GetDefinition() == aProton )
443  {
444  currentParticle.SetDefinitionAndUpdateE( aPiZero );
445  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
446  incidentHasChanged = true;
447  targetHasChanged = true;
448  }
449  else
450  {
451  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
452  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
453  incidentHasChanged = true;
454  targetHasChanged = true;
455  }
456  }
457  else
458  {
459  if( targetParticle.GetDefinition() == aProton )
460  {
461  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
462  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
463  incidentHasChanged = true;
464  targetHasChanged = true;
465  }
466  else
467  {
468  currentParticle.SetDefinitionAndUpdateE( aPiZero );
469  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
470  incidentHasChanged = true;
471  targetHasChanged = true;
472  }
473  }
474  }
475  else // ( current != aKaonMinus || target != aNeutron ) &&
476  // ( current != aKaonZL || target != aProton )
477  {
478  ran = G4UniformRand();
479  if( ran < 0.67 )
480  {
481  currentParticle.SetDefinitionAndUpdateE( aPiZero );
482  targetParticle.SetDefinitionAndUpdateE( aLambda );
483  incidentHasChanged = true;
484  targetHasChanged = true;
485  }
486  else if( ran < 0.78 )
487  {
488  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
489  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
490  incidentHasChanged = true;
491  targetHasChanged = true;
492  }
493  else if( ran < 0.89 )
494  {
495  currentParticle.SetDefinitionAndUpdateE( aPiZero );
496  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
497  incidentHasChanged = true;
498  targetHasChanged = true;
499  }
500  else
501  {
502  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
503  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
504  incidentHasChanged = true;
505  targetHasChanged = true;
506  }
507  }
508  }
509  }
510 
511  if (currentParticle.GetDefinition() == aKaonZL) {
512  if (G4UniformRand() >= 0.5) {
513  currentParticle.SetDefinitionAndUpdateE(aKaonZS);
514  incidentHasChanged = true;
515  }
516  }
517 
518  if (targetParticle.GetDefinition() == aKaonZL) {
519  if (G4UniformRand() >= 0.5) {
520  targetParticle.SetDefinitionAndUpdateE(aKaonZS);
521  targetHasChanged = true;
522  }
523  }
524 
525  SetUpPions(np, nneg, nz, vec, vecLen);
526  return;
527 }
528 
529  /* end of file */
530 
531 
532