Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiKZeroInelastic.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 
30 #include "Randomize.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
34 
37  G4Nucleus &targetNucleus )
38  {
39  const G4HadProjectile *originalIncident = &aTrack;
40  //
41  // create the target particle
42  //
43  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
44 
45  if( verboseLevel > 1 )
46  {
47  const G4Material *targetMaterial = aTrack.GetMaterial();
48  G4cout << "G4RPGAntiKZeroInelastic::ApplyYourself called" << G4endl;
49  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
50  G4cout << "target material = " << targetMaterial->GetName() << ", ";
51  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
52  << G4endl;
53  }
54  //
55  // Fermi motion and evaporation
56  // As of Geant3, the Fermi energy calculation had not been Done
57  //
58  G4double ek = originalIncident->GetKineticEnergy()/MeV;
59  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
60  G4ReactionProduct modifiedOriginal;
61  modifiedOriginal = *originalIncident;
62 
63  G4double tkin = targetNucleus.Cinema( ek );
64  ek += tkin;
65  modifiedOriginal.SetKineticEnergy( ek*MeV );
66  G4double et = ek + amas;
67  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
68  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
69  if( pp > 0.0 )
70  {
71  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
72  modifiedOriginal.SetMomentum( momentum * (p/pp) );
73  }
74  //
75  // calculate black track energies
76  //
77  tkin = targetNucleus.EvaporationEffects( ek );
78  ek -= tkin;
79  modifiedOriginal.SetKineticEnergy( ek*MeV );
80  et = ek + amas;
81  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
82  pp = modifiedOriginal.GetMomentum().mag()/MeV;
83  if( pp > 0.0 )
84  {
85  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
86  modifiedOriginal.SetMomentum( momentum * (p/pp) );
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  try
107  {
108  CalculateMomenta( vec, vecLen,
109  originalIncident, originalTarget, modifiedOriginal,
110  targetNucleus, currentParticle, targetParticle,
111  incidentHasChanged, targetHasChanged, quasiElastic );
112  }
113  catch(G4HadReentrentException aR)
114  {
115  aR.Report(G4cout);
116  throw G4HadReentrentException(__FILE__, __LINE__, "Bailing out");
117  }
118  SetUpChange( vec, vecLen,
119  currentParticle, targetParticle,
120  incidentHasChanged );
121 
122  delete originalTarget;
123  return &theParticleChange;
124  }
125 
126 void G4RPGAntiKZeroInelastic::Cascade(
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 CASK0B
137  //
138  // K0Long 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->Get4Momentum().e()/MeV;
148  const G4double pOriginal = originalIncident->Get4Momentum().vect().mag()/MeV;
149  const G4double targetMass = targetParticle.GetMass()/MeV;
150  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
151  targetMass*targetMass +
152  2.0*targetMass*etOriginal );
153  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
154 
155  static G4ThreadLocal G4bool first = true;
156  const G4int numMul = 1200;
157  const G4int numSec = 60;
158  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
159  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
160 
161  // np = number of pi+, nneg = number of pi-, nz = number of pi0
162 
163  G4int counter, nt=0, np=0, nneg=0, nz=0;
164  const G4double c = 1.25;
165  const G4double b[] = { 0.7, 0.7 };
166  if( first ) // compute normalization constants, this will only be Done once
167  {
168  first = false;
169  G4int i;
170  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
171  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
172  counter = -1;
173  for( np=0; np<numSec/3; ++np )
174  {
175  for( nneg=std::max(0,np-2); nneg<=np; ++nneg )
176  {
177  for( nz=0; nz<numSec/3; ++nz )
178  {
179  if( ++counter < numMul )
180  {
181  nt = np+nneg+nz;
182  if( nt>0 && nt<=numSec )
183  {
184  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
185  protnorm[nt-1] += protmul[counter];
186  }
187  }
188  }
189  }
190  }
191  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
192  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
193  counter = -1;
194  for( np=0; np<(numSec/3); ++np )
195  {
196  for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
197  {
198  for( nz=0; nz<numSec/3; ++nz )
199  {
200  if( ++counter < numMul )
201  {
202  nt = np+nneg+nz;
203  if( nt>0 && nt<=numSec )
204  {
205  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
206  neutnorm[nt-1] += neutmul[counter];
207  }
208  }
209  }
210  }
211  }
212  for( i=0; i<numSec; ++i )
213  {
214  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
215  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
216  }
217  } // end of initialization
218 
219  const G4double expxu = 82.; // upper bound for arg. of exp
220  const G4double expxl = -expxu; // lower bound for arg. of exp
233  const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
234  G4int iplab = G4int(std::min( 9.0, 5.0*pOriginal*MeV/GeV ));
235 
236  if ((pOriginal*MeV/GeV <= 2.0) && (G4UniformRand() < cech[iplab]) ) {
237  np = nneg = nz = nt = 0;
238  iplab = G4int(std::min( 19.0, pOriginal*MeV/GeV*10.0 ));
239  const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
240  0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
241  if (G4UniformRand() > cnk0[iplab] ) {
242  G4double ran = G4UniformRand();
243  if (ran < 0.25) { // k0Long n --> pi- s+
244  if (targetParticle.GetDefinition() == aNeutron) {
245  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
246  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
247  incidentHasChanged = true;
248  targetHasChanged = true;
249  }
250  } else if( ran < 0.50 ) { // k0Long p --> pi+ s0 or k0Long n --> pi0 s0
251  if( targetParticle.GetDefinition() == aNeutron )
252  currentParticle.SetDefinitionAndUpdateE( aPiZero );
253  else
254  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
255  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
256  incidentHasChanged = true;
257  targetHasChanged = true;
258  } else if( ran < 0.75 ) { // k0Long n --> pi+ s-
259  if( targetParticle.GetDefinition() == aNeutron )
260  {
261  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
262  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
263  incidentHasChanged = true;
264  targetHasChanged = true;
265  }
266  } else { // k0Long p --> pi+ L or k0Long n --> pi0 L
267  if( targetParticle.GetDefinition() == aNeutron )
268  currentParticle.SetDefinitionAndUpdateE( aPiZero );
269  else
270  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
271  targetParticle.SetDefinitionAndUpdateE( aLambda );
272  incidentHasChanged = true;
273  targetHasChanged = true;
274  }
275  } else { // ran <= cnk0
276  quasiElastic = true;
277  if (targetParticle.GetDefinition() == aNeutron) {
278  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
279  targetParticle.SetDefinitionAndUpdateE( aProton );
280  incidentHasChanged = true;
281  targetHasChanged = true;
282  }
283  }
284  } else { // (pOriginal > 2.0*GeV) || (random number >= cech[iplab])
285  if (availableEnergy < aPiPlus->GetPDGMass()/MeV) {
286  quasiElastic = true;
287  return;
288  }
289  G4double n, anpn;
290  GetNormalizationConstant( availableEnergy, n, anpn );
291  G4double ran = G4UniformRand();
292  G4double dum, test, excs = 0.0;
293  if (targetParticle.GetDefinition() == aProton) {
294  counter = -1;
295  for( np=0; (np<numSec/3) && (ran>=excs); ++np )
296  {
297  for( nneg=std::max(0,np-2); nneg<=np && ran>=excs; ++nneg )
298  {
299  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
300  {
301  if( ++counter < numMul )
302  {
303  nt = np+nneg+nz;
304  if( nt>0 && nt<=numSec )
305  {
306  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
307  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
308  if( std::fabs(dum) < 1.0 )
309  {
310  if( test >= 1.0e-10 )excs += dum*test;
311  }
312  else
313  excs += dum*test;
314  }
315  }
316  }
317  }
318  }
319  if( ran >= excs ) // 3 previous loops continued to the end
320  {
321  quasiElastic = true;
322  return;
323  }
324  np--; nneg--; nz--;
325  switch( np-nneg )
326  {
327  case 1:
328  if( G4UniformRand() < 0.5 )
329  {
330  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
331  incidentHasChanged = true;
332  }
333  else
334  {
335  targetParticle.SetDefinitionAndUpdateE( aNeutron );
336  targetHasChanged = true;
337  }
338  case 0:
339  break;
340  default:
341  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
342  targetParticle.SetDefinitionAndUpdateE( aNeutron );
343  incidentHasChanged = true;
344  targetHasChanged = true;
345  break;
346  }
347  }
348  else // target must be a neutron
349  {
350  counter = -1;
351  for( np=0; np<numSec/3 && ran>=excs; ++np )
352  {
353  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
354  {
355  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
356  {
357  if( ++counter < numMul )
358  {
359  nt = np+nneg+nz;
360  if( nt>0 && nt<=numSec )
361  {
362  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
363  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
364  if( std::fabs(dum) < 1.0 )
365  {
366  if( test >= 1.0e-10 )excs += dum*test;
367  }
368  else
369  excs += dum*test;
370  }
371  }
372  }
373  }
374  }
375  if( ran >= excs ) // 3 previous loops continued to the end
376  {
377  quasiElastic = true;
378  return;
379  }
380  np--; nneg--; nz--;
381  switch( np-nneg )
382  {
383  case 0:
384  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
385  targetParticle.SetDefinitionAndUpdateE( aProton );
386  incidentHasChanged = true;
387  targetHasChanged = true;
388  break;
389  case 1:
390  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
391  incidentHasChanged = true;
392  break;
393  default:
394  targetParticle.SetDefinitionAndUpdateE( aProton );
395  targetHasChanged = true;
396  break;
397  }
398  }
399  if( G4UniformRand() >= 0.5 )
400  {
401  if( currentParticle.GetDefinition() == aKaonMinus &&
402  targetParticle.GetDefinition() == aNeutron )
403  {
404  ran = G4UniformRand();
405  if( ran < 0.68 )
406  {
407  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
408  targetParticle.SetDefinitionAndUpdateE( aLambda );
409  }
410  else if( ran < 0.84 )
411  {
412  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
413  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
414  }
415  else
416  {
417  currentParticle.SetDefinitionAndUpdateE( aPiZero );
418  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
419  }
420  }
421  else if( (currentParticle.GetDefinition() == aKaonZS ||
422  currentParticle.GetDefinition() == aKaonZL ) &&
423  targetParticle.GetDefinition() == aProton )
424  {
425  ran = G4UniformRand();
426  if( ran < 0.68 )
427  {
428  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
429  targetParticle.SetDefinitionAndUpdateE( aLambda );
430  }
431  else if( ran < 0.84 )
432  {
433  currentParticle.SetDefinitionAndUpdateE( aPiZero );
434  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
435  }
436  else
437  {
438  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
439  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
440  }
441  }
442  else
443  {
444  ran = G4UniformRand();
445  if( ran < 0.67 )
446  {
447  currentParticle.SetDefinitionAndUpdateE( aPiZero );
448  targetParticle.SetDefinitionAndUpdateE( aLambda );
449  }
450  else if( ran < 0.78 )
451  {
452  currentParticle.SetDefinitionAndUpdateE( aPiMinus );
453  targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
454  }
455  else if( ran < 0.89 )
456  {
457  currentParticle.SetDefinitionAndUpdateE( aPiZero );
458  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
459  }
460  else
461  {
462  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
463  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
464  }
465  }
466  incidentHasChanged = true;
467  targetHasChanged = true;
468  }
469  }
470  if( currentParticle.GetDefinition() == aKaonZL )
471  {
472  if( G4UniformRand() >= 0.5 )
473  {
474  currentParticle.SetDefinitionAndUpdateE( aKaonZS );
475  incidentHasChanged = true;
476  }
477  }
478  if( targetParticle.GetDefinition() == aKaonZL )
479  {
480  if( G4UniformRand() >= 0.5 )
481  {
482  targetParticle.SetDefinitionAndUpdateE( aKaonZS );
483  targetHasChanged = true;
484  }
485  }
486  SetUpPions( np, nneg, nz, vec, vecLen );
487  return;
488 }
489 
490  /* end of file */
491 
G4double EvaporationEffects(G4double kineticEnergy)
Definition: G4Nucleus.cc:264
void SetUpChange(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:176
static G4KaonZeroLong * KaonZeroLong()
void SetSide(const G4int sid)
void CalculateMomenta(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4DynamicParticle *originalTarget, G4ReactionProduct &modifiedOriginal, G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool quasiElastic)
G4ParticleDefinition * GetDefinition() const
#define G4ThreadLocal
Definition: tls.hh:52
void Initialize(G4int items)
Definition: G4FastVector.hh:63
int G4int
Definition: G4Types.hh:78
G4DynamicParticle * ReturnTargetParticle() const
Definition: G4Nucleus.cc:227
void SetDefinitionAndUpdateE(G4ParticleDefinition *aParticleDefinition)
const G4String & GetParticleName() const
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4ParticleDefinition * GetDefinition() const
G4double Pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
static G4KaonZeroShort * KaonZeroShort()
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4int n
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
static G4SigmaMinus * SigmaMinus()
G4double GetKineticEnergy() const
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double Cinema(G4double kineticEnergy)
Definition: G4Nucleus.cc:368
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
const G4Material * GetMaterial() const
def test
Definition: mcscore.py:117
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
void GetNormalizationConstant(const G4double availableEnergy, G4double &n, G4double &anpn)
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
double G4double
Definition: G4Types.hh:76
void SetUpPions(const G4int np, const G4int nm, const G4int nz, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
double mag() const
void Report(std::ostream &aS)
G4double GetMass() const