ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiKZeroInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiKZeroInelastic.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 "G4RPGAntiKZeroInelastic.hh"
29 #include "Randomize.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4HadReentrentException.hh"
33 
34  G4HadFinalState *
35  G4RPGAntiKZeroInelastic::ApplyYourself( const G4HadProjectile &aTrack,
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 << "G4RPGAntiKZeroInelastic::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  try
106  {
107  CalculateMomenta( vec, vecLen,
108  originalIncident, originalTarget, modifiedOriginal,
109  targetNucleus, currentParticle, targetParticle,
110  incidentHasChanged, targetHasChanged, quasiElastic );
111  }
112  catch(G4HadReentrentException & aR)
113  {
114  aR.Report(G4cout);
115  throw G4HadReentrentException(__FILE__, __LINE__, "Bailing out");
116  }
117  SetUpChange( vec, vecLen,
118  currentParticle, targetParticle,
119  incidentHasChanged );
120 
121  delete originalTarget;
122  return &theParticleChange;
123  }
124 
125 void G4RPGAntiKZeroInelastic::Cascade(
126  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
127  G4int& vecLen,
128  const G4HadProjectile *originalIncident,
129  G4ReactionProduct &currentParticle,
130  G4ReactionProduct &targetParticle,
131  G4bool &incidentHasChanged,
132  G4bool &targetHasChanged,
133  G4bool &quasiElastic )
134 {
135  // Derived from H. Fesefeldt's original FORTRAN code CASK0B
136  //
137  // K0Long undergoes interaction with nucleon within a nucleus. Check if it is
138  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
139  // occurs and input particle is degraded in energy. No other particles are produced.
140  // If reaction is possible, find the correct number of pions/protons/neutrons
141  // produced using an interpolation to multiplicity data. Replace some pions or
142  // protons/neutrons by kaons or strange baryons according to the average
143  // multiplicity per Inelastic reaction.
144 
145  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
146  const G4double etOriginal = originalIncident->Get4Momentum().e()/MeV;
147  const G4double pOriginal = originalIncident->Get4Momentum().vect().mag()/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 
154  static G4ThreadLocal G4bool first = true;
155  const G4int numMul = 1200;
156  const G4int numSec = 60;
157  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
158  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
159 
160  // np = number of pi+, nneg = number of pi-, nz = number of pi0
161 
162  G4int counter, nt=0, np=0, nneg=0, nz=0;
163  const G4double c = 1.25;
164  const G4double b[] = { 0.7, 0.7 };
165  if( first ) // compute normalization constants, this will only be Done once
166  {
167  first = false;
168  G4int i;
169  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
170  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
171  counter = -1;
172  for( np=0; np<numSec/3; ++np )
173  {
174  for( nneg=std::max(0,np-2); nneg<=np; ++nneg )
175  {
176  for( nz=0; nz<numSec/3; ++nz )
177  {
178  if( ++counter < numMul )
179  {
180  nt = np+nneg+nz;
181  if( nt>0 && nt<=numSec )
182  {
183  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
184  protnorm[nt-1] += protmul[counter];
185  }
186  }
187  }
188  }
189  }
190  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
191  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
192  counter = -1;
193  for( np=0; np<(numSec/3); ++np )
194  {
195  for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
196  {
197  for( nz=0; nz<numSec/3; ++nz )
198  {
199  if( ++counter < numMul )
200  {
201  nt = np+nneg+nz;
202  if( nt>0 && nt<=numSec )
203  {
204  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
205  neutnorm[nt-1] += neutmul[counter];
206  }
207  }
208  }
209  }
210  }
211  for( i=0; i<numSec; ++i )
212  {
213  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
214  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
215  }
216  } // end of initialization
217 
218  const G4double expxu = 82.; // upper bound for arg. of exp
219  const G4double expxl = -expxu; // lower bound for arg. of exp
220  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
221  G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
222  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
223  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
224  G4ParticleDefinition *aProton = G4Proton::Proton();
225  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
226  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
227  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
228  G4ParticleDefinition *aLambda = G4Lambda::Lambda();
229  G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
230  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
231  G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
232  const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
233  G4int iplab = G4int(std::min( 9.0, 5.0*pOriginal*MeV/GeV ));
234 
235  if ((pOriginal*MeV/GeV <= 2.0) && (G4UniformRand() < cech[iplab]) ) {
236  np = nneg = nz = nt = 0;
237  iplab = G4int(std::min( 19.0, pOriginal*MeV/GeV*10.0 ));
238  const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
239  0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
240  if (G4UniformRand() > cnk0[iplab] ) {
241  G4double ran = G4UniformRand();
242  if (ran < 0.25) { // k0Long n --> pi- s+
243  if (targetParticle.GetDefinition() == aNeutron) {
244  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
245  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
246  incidentHasChanged = true;
247  targetHasChanged = true;
248  }
249  } else if( ran < 0.50 ) { // k0Long p --> pi+ s0 or k0Long n --> pi0 s0
250  if( targetParticle.GetDefinition() == aNeutron )
251  currentParticle.SetDefinitionAndUpdateE( aPiZero );
252  else
253  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
254  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
255  incidentHasChanged = true;
256  targetHasChanged = true;
257  } else if( ran < 0.75 ) { // k0Long n --> pi+ s-
258  if( targetParticle.GetDefinition() == aNeutron )
259  {
260  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
261  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
262  incidentHasChanged = true;
263  targetHasChanged = true;
264  }
265  } else { // k0Long p --> pi+ L or k0Long n --> pi0 L
266  if( targetParticle.GetDefinition() == aNeutron )
267  currentParticle.SetDefinitionAndUpdateE( aPiZero );
268  else
269  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
270  targetParticle.SetDefinitionAndUpdateE( aLambda );
271  incidentHasChanged = true;
272  targetHasChanged = true;
273  }
274  } else { // ran <= cnk0
275  quasiElastic = true;
276  if (targetParticle.GetDefinition() == aNeutron) {
277  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
278  targetParticle.SetDefinitionAndUpdateE( aProton );
279  incidentHasChanged = true;
280  targetHasChanged = true;
281  }
282  }
283  } else { // (pOriginal > 2.0*GeV) || (random number >= cech[iplab])
284  if (availableEnergy < aPiPlus->GetPDGMass()/MeV) {
285  quasiElastic = true;
286  return;
287  }
288  G4double n, anpn;
289  GetNormalizationConstant( availableEnergy, n, anpn );
290  G4double ran = G4UniformRand();
291  G4double dum, test, excs = 0.0;
292  if (targetParticle.GetDefinition() == aProton) {
293  counter = -1;
294  for( np=0; (np<numSec/3) && (ran>=excs); ++np )
295  {
296  for( nneg=std::max(0,np-2); nneg<=np && ran>=excs; ++nneg )
297  {
298  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
299  {
300  if( ++counter < numMul )
301  {
302  nt = np+nneg+nz;
303  if( nt>0 && nt<=numSec )
304  {
305  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
306  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
307  if( std::fabs(dum) < 1.0 )
308  {
309  if( test >= 1.0e-10 )excs += dum*test;
310  }
311  else
312  excs += dum*test;
313  }
314  }
315  }
316  }
317  }
318  if( ran >= excs ) // 3 previous loops continued to the end
319  {
320  quasiElastic = true;
321  return;
322  }
323  np--; nneg--; nz--;
324  switch( np-nneg )
325  {
326  case 1:
327  if( G4UniformRand() < 0.5 )
328  {
329  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
330  incidentHasChanged = true;
331  }
332  else
333  {
334  targetParticle.SetDefinitionAndUpdateE( aNeutron );
335  targetHasChanged = true;
336  }
337  case 0:
338  break;
339  default:
340  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
341  targetParticle.SetDefinitionAndUpdateE( aNeutron );
342  incidentHasChanged = true;
343  targetHasChanged = true;
344  break;
345  }
346  }
347  else // target must be a neutron
348  {
349  counter = -1;
350  for( np=0; np<numSec/3 && ran>=excs; ++np )
351  {
352  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
353  {
354  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
355  {
356  if( ++counter < numMul )
357  {
358  nt = np+nneg+nz;
359  if( nt>0 && nt<=numSec )
360  {
361  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
362  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
363  if( std::fabs(dum) < 1.0 )
364  {
365  if( test >= 1.0e-10 )excs += dum*test;
366  }
367  else
368  excs += dum*test;
369  }
370  }
371  }
372  }
373  }
374  if( ran >= excs ) // 3 previous loops continued to the end
375  {
376  quasiElastic = true;
377  return;
378  }
379  np--; nneg--; nz--;
380  switch( np-nneg )
381  {
382  case 0:
383  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
384  targetParticle.SetDefinitionAndUpdateE( aProton );
385  incidentHasChanged = true;
386  targetHasChanged = true;
387  break;
388  case 1:
389  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
390  incidentHasChanged = true;
391  break;
392  default:
393  targetParticle.SetDefinitionAndUpdateE( aProton );
394  targetHasChanged = true;
395  break;
396  }
397  }
398  if( G4UniformRand() >= 0.5 )
399  {
400  if( currentParticle.GetDefinition() == aKaonMinus &&
401  targetParticle.GetDefinition() == aNeutron )
402  {
403  ran = G4UniformRand();
404  if( ran < 0.68 )
405  {
406  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
407  targetParticle.SetDefinitionAndUpdateE( aLambda );
408  }
409  else if( ran < 0.84 )
410  {
411  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
412  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
413  }
414  else
415  {
416  currentParticle.SetDefinitionAndUpdateE( aPiZero );
417  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
418  }
419  }
420  else if( (currentParticle.GetDefinition() == aKaonZS ||
421  currentParticle.GetDefinition() == aKaonZL ) &&
422  targetParticle.GetDefinition() == aProton )
423  {
424  ran = G4UniformRand();
425  if( ran < 0.68 )
426  {
427  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
428  targetParticle.SetDefinitionAndUpdateE( aLambda );
429  }
430  else if( ran < 0.84 )
431  {
432  currentParticle.SetDefinitionAndUpdateE( aPiZero );
433  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
434  }
435  else
436  {
437  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
438  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
439  }
440  }
441  else
442  {
443  ran = G4UniformRand();
444  if( ran < 0.67 )
445  {
446  currentParticle.SetDefinitionAndUpdateE( aPiZero );
447  targetParticle.SetDefinitionAndUpdateE( aLambda );
448  }
449  else if( ran < 0.78 )
450  {
451  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
452  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
453  }
454  else if( ran < 0.89 )
455  {
456  currentParticle.SetDefinitionAndUpdateE( aPiZero );
457  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
458  }
459  else
460  {
461  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
462  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
463  }
464  }
465  incidentHasChanged = true;
466  targetHasChanged = true;
467  }
468  }
469  if( currentParticle.GetDefinition() == aKaonZL )
470  {
471  if( G4UniformRand() >= 0.5 )
472  {
473  currentParticle.SetDefinitionAndUpdateE( aKaonZS );
474  incidentHasChanged = true;
475  }
476  }
477  if( targetParticle.GetDefinition() == aKaonZL )
478  {
479  if( G4UniformRand() >= 0.5 )
480  {
481  targetParticle.SetDefinitionAndUpdateE( aKaonZS );
482  targetHasChanged = true;
483  }
484  }
485  SetUpPions( np, nneg, nz, vec, vecLen );
486  return;
487 }
488 
489  /* end of file */
490