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