ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiSigmaPlusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiSigmaPlusInelastic.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 "G4RPGAntiSigmaPlusInelastic.hh"
29 #include "G4Exp.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "Randomize.hh"
33 
34 G4HadFinalState*
35 G4RPGAntiSigmaPlusInelastic::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  // Choose the target particle
48 
49  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
50 
51  if( verboseLevel > 1 )
52  {
53  const G4Material *targetMaterial = aTrack.GetMaterial();
54  G4cout << "G4RPGAntiSigmaPlusInelastic::ApplyYourself called" << G4endl;
55  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
56  G4cout << "target material = " << targetMaterial->GetName() << ", ";
57  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
58  << G4endl;
59  }
60 
61  // Fermi motion and evaporation
62  // As of Geant3, the Fermi energy calculation had not been Done
63 
64  G4double ek = originalIncident->GetKineticEnergy()/MeV;
65  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
66  G4ReactionProduct modifiedOriginal;
67  modifiedOriginal = *originalIncident;
68 
69  G4double tkin = targetNucleus.Cinema( ek );
70  ek += tkin;
71  modifiedOriginal.SetKineticEnergy( ek*MeV );
72  G4double et = ek + amas;
73  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
74  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
75  if( pp > 0.0 )
76  {
77  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
78  modifiedOriginal.SetMomentum( momentum * (p/pp) );
79  }
80  //
81  // calculate black track energies
82  //
83  tkin = targetNucleus.EvaporationEffects( ek );
84  ek -= tkin;
85  modifiedOriginal.SetKineticEnergy( ek*MeV );
86  et = ek + amas;
87  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
88  pp = modifiedOriginal.GetMomentum().mag()/MeV;
89  if( pp > 0.0 )
90  {
91  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
92  modifiedOriginal.SetMomentum( momentum * (p/pp) );
93  }
94  G4ReactionProduct currentParticle = modifiedOriginal;
95  G4ReactionProduct targetParticle;
96  targetParticle = *originalTarget;
97  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
98  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
99  G4bool incidentHasChanged = false;
100  G4bool targetHasChanged = false;
101  G4bool quasiElastic = false;
102  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
103  G4int vecLen = 0;
104  vec.Initialize( 0 );
105 
106  const G4double cutOff = 0.1;
107  const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
108  if( (currentParticle.GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
109  Cascade( vec, vecLen,
110  originalIncident, currentParticle, targetParticle,
111  incidentHasChanged, targetHasChanged, quasiElastic );
112 
113  CalculateMomenta( vec, vecLen,
114  originalIncident, originalTarget, modifiedOriginal,
115  targetNucleus, currentParticle, targetParticle,
116  incidentHasChanged, targetHasChanged, quasiElastic );
117 
118  SetUpChange( vec, vecLen,
119  currentParticle, targetParticle,
120  incidentHasChanged );
121 
122  delete originalTarget;
123  return &theParticleChange;
124 }
125 
126 
127 void G4RPGAntiSigmaPlusInelastic::Cascade(
128  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
129  G4int& vecLen,
130  const G4HadProjectile *originalIncident,
131  G4ReactionProduct &currentParticle,
132  G4ReactionProduct &targetParticle,
133  G4bool &incidentHasChanged,
134  G4bool &targetHasChanged,
135  G4bool &quasiElastic )
136 {
137  // Derived from H. Fesefeldt's original FORTRAN code CASASP
138  // AntiSigmaPlus 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->GetTotalEnergy()/MeV;
148  const G4double pOriginal = originalIncident->GetTotalMomentum()/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 numMulA = 400;
158  const G4int numSec = 60;
159  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
160  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
161  static G4ThreadLocal G4double protmulA[numMulA], protnormA[numSec]; // proton constants
162  static G4ThreadLocal G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
163  // np = number of pi+, nneg = number of pi-, nz = number of pi0
164  G4int counter, nt=0, np=0, nneg=0, nz=0;
165  G4double test;
166  const G4double c = 1.25;
167  const G4double b[] = { 0.7, 0.7 };
168  if( first ) // compute normalization constants, this will only be Done once
169  {
170  first = false;
171  G4int i;
172  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
173  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
174  counter = -1;
175  for( np=0; np<(numSec/3); ++np )
176  {
177  for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
178  {
179  for( nz=0; nz<numSec/3; ++nz )
180  {
181  if( ++counter < numMul )
182  {
183  nt = np+nneg+nz;
184  if( nt>0 && nt<=numSec )
185  {
186  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
187  protnorm[nt-1] += protmul[counter];
188  }
189  }
190  }
191  }
192  }
193  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
194  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
195  counter = -1;
196  for( np=0; np<numSec/3; ++np )
197  {
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  for( i=0; i<numSec; ++i )
215  {
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  for( np=1; np<(numSec/3); ++np )
226  {
227  nneg = np;
228  for( nz=0; nz<numSec/3; ++nz )
229  {
230  if( ++counter < numMulA )
231  {
232  nt = np+nneg+nz;
233  if( nt>1 && nt<=numSec )
234  {
235  protmulA[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
236  protnormA[nt-1] += protmulA[counter];
237  }
238  }
239  }
240  }
241  for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
242  for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
243  counter = -1;
244  for( np=0; np<numSec/3; ++np )
245  {
246  nneg = np+1;
247  for( nz=0; nz<numSec/3; ++nz )
248  {
249  if( ++counter < numMulA )
250  {
251  nt = np+nneg+nz;
252  if( nt>1 && nt<=numSec )
253  {
254  neutmulA[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
255  neutnormA[nt-1] += neutmulA[counter];
256  }
257  }
258  }
259  }
260  for( i=0; i<numSec; ++i )
261  {
262  if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
263  if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
264  }
265  } // end of initialization
266 
267  const G4double expxu = 82.; // upper bound for arg. of exp
268  const G4double expxl = -expxu; // lower bound for arg. of exp
269  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
270  G4ParticleDefinition *aProton = G4Proton::Proton();
271  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
272  G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
273  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
274  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
275  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
276  G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
277  const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
278  0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
279  0.39,0.36,0.33,0.10,0.01};
280  G4int iplab = G4int( pOriginal/GeV*10.0 );
281  if( iplab > 9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
282  if( iplab > 14 )iplab = G4int( pOriginal/GeV- 2.0 ) + 15;
283  if( iplab > 22 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
284  if( iplab > 24 )iplab = 24;
285  if( G4UniformRand() > anhl[iplab] )
286  {
287  if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
288  {
289  quasiElastic = true;
290  return;
291  }
292  G4double n, anpn;
293  GetNormalizationConstant( availableEnergy, n, anpn );
294  G4double ran = G4UniformRand();
295  G4double dum, excs = 0.0;
296  if( targetParticle.GetDefinition() == aProton )
297  {
298  counter = -1;
299  for( np=0; np<numSec/3 && ran>=excs; ++np )
300  {
301  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
302  {
303  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
304  {
305  if( ++counter < numMul )
306  {
307  nt = np+nneg+nz;
308  if( (nt>0) && (nt<=numSec) )
309  {
310  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
311  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
312  if( std::fabs(dum) < 1.0 )
313  {
314  if( test >= 1.0e-10 )excs += dum*test;
315  }
316  else
317  excs += dum*test;
318  }
319  }
320  }
321  }
322  }
323  if( ran >= excs ) // 3 previous loops continued to the end
324  {
325  quasiElastic = true;
326  return;
327  }
328  np--; nneg--; nz--;
329  G4int ncht = std::min( 3, std::max( 1, np-nneg+2 ) );
330  switch( ncht )
331  {
332  case 1:
333  if( G4UniformRand() < 0.5 )
334  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
335  else
336  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
337  incidentHasChanged = true;
338  break;
339  case 2:
340  if( G4UniformRand() >= 0.5 )
341  {
342  if( G4UniformRand() < 0.5 )
343  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
344  else
345  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
346  incidentHasChanged = true;
347  }
348  targetParticle.SetDefinitionAndUpdateE( aNeutron );
349  targetHasChanged = true;
350  break;
351  case 3:
352  targetParticle.SetDefinitionAndUpdateE( aNeutron );
353  targetHasChanged = true;
354  break;
355  }
356  }
357  else // target must be a neutron
358  {
359  counter = -1;
360  for( np=0; np<numSec/3 && ran>=excs; ++np )
361  {
362  for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
363  {
364  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
365  {
366  if( ++counter < numMul )
367  {
368  nt = np+nneg+nz;
369  if( (nt>0) && (nt<=numSec) )
370  {
371  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
372  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
373  if( std::fabs(dum) < 1.0 )
374  {
375  if( test >= 1.0e-10 )excs += dum*test;
376  }
377  else
378  excs += dum*test;
379  }
380  }
381  }
382  }
383  }
384  if( ran >= excs ) // 3 previous loops continued to the end
385  {
386  quasiElastic = true;
387  return;
388  }
389  np--; nneg--; nz--;
390  G4int ncht = std::min( 3, std::max( 1, np-nneg+3 ) );
391  switch( ncht )
392  {
393  case 1:
394  if( G4UniformRand() < 0.5 )
395  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
396  else
397  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
398  incidentHasChanged = true;
399  targetParticle.SetDefinitionAndUpdateE( aProton );
400  targetHasChanged = true;
401  break;
402  case 2:
403  if( G4UniformRand() < 0.5 )
404  {
405  if( G4UniformRand() < 0.5 )
406  {
407  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
408  incidentHasChanged = true;
409  }
410  else
411  {
412  targetParticle.SetDefinitionAndUpdateE( aProton );
413  targetHasChanged = true;
414  }
415  }
416  else
417  {
418  if( G4UniformRand() < 0.5 )
419  {
420  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
421  incidentHasChanged = true;
422  }
423  else
424  {
425  targetParticle.SetDefinitionAndUpdateE( aProton );
426  targetHasChanged = true;
427  }
428  }
429  break;
430  case 3:
431  break;
432  }
433  }
434  }
435  else // random number <= anhl[iplab]
436  {
437  if( centerofmassEnergy <= aPiPlus->GetPDGMass()/MeV+aKaonPlus->GetPDGMass()/MeV )
438  {
439  quasiElastic = true;
440  return;
441  }
442  G4double n, anpn;
443  GetNormalizationConstant( -centerofmassEnergy, n, anpn );
444  G4double ran = G4UniformRand();
445  G4double dum, excs = 0.0;
446  if( targetParticle.GetDefinition() == aProton )
447  {
448  counter = -1;
449  for( np=1; np<numSec/3 && ran>=excs; ++np )
450  {
451  nneg = np;
452  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
453  {
454  if( ++counter < numMulA )
455  {
456  nt = np+nneg+nz;
457  if( nt>1 && nt<=numSec )
458  {
459  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
460  dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
461  if( std::fabs(dum) < 1.0 )
462  {
463  if( test >= 1.0e-10 )excs += dum*test;
464  }
465  else
466  excs += dum*test;
467  }
468  }
469  }
470  }
471  if( ran >= excs ) // 3 previous loops continued to the end
472  {
473  quasiElastic = true;
474  return;
475  }
476  np--; nz--;
477  }
478  else // target must be a neutron
479  {
480  counter = -1;
481  for( np=0; np<numSec/3 && ran>=excs; ++np )
482  {
483  nneg = np+1;
484  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
485  {
486  if( ++counter < numMulA )
487  {
488  nt = np+nneg+nz;
489  if( nt>1 && nt<=numSec )
490  {
491  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
492  dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
493  if( std::fabs(dum) < 1.0 )
494  {
495  if( test >= 1.0e-10 )excs += dum*test;
496  }
497  else
498  excs += dum*test;
499  }
500  }
501  }
502  }
503  if( ran >= excs ) // 3 previous loops continued to the end
504  {
505  quasiElastic = true;
506  return;
507  }
508  np--; nz--;
509  }
510  if( nz > 0 )
511  {
512  if( nneg > 0 )
513  {
514  if( G4UniformRand() < 0.5 )
515  {
516  vec.Initialize( 1 );
517  G4ReactionProduct *p= new G4ReactionProduct;
518  p->SetDefinition( aKaonMinus );
519  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
520  vec.SetElement( vecLen++, p );
521  --nneg;
522  }
523  else
524  {
525  vec.Initialize( 1 );
526  G4ReactionProduct *p= new G4ReactionProduct ;
527  p->SetDefinition( aKaonZL );
528  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
529  vec.SetElement( vecLen++, p );
530  --nz;
531  }
532  }
533  else // nneg == 0
534  {
535  vec.Initialize( 1 );
536  G4ReactionProduct *p = new G4ReactionProduct;
537  p->SetDefinition( aKaonZL );
538  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
539  vec.SetElement( vecLen++, p );
540  --nz;
541  }
542  }
543  else // nz == 0
544  {
545  if( nneg > 0 )
546  {
547  vec.Initialize( 1 );
548  G4ReactionProduct *p = new G4ReactionProduct;
549  p->SetDefinition( aKaonMinus );
550  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
551  vec.SetElement( vecLen++, p );
552  --nneg;
553  }
554  }
555  currentParticle.SetMass( 0.0 );
556  targetParticle.SetMass( 0.0 );
557  }
558 
559  SetUpPions( np, nneg, nz, vec, vecLen );
560  return;
561 }
562 
563  /* end of file */
564