ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiSigmaMinusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiSigmaMinusInelastic.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 "G4RPGAntiSigmaMinusInelastic.hh"
29 #include "G4Exp.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "Randomize.hh"
33 
34 G4HadFinalState*
35 G4RPGAntiSigmaMinusInelastic::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 << "G4RPGAntiSigmaMinusInelastic::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 G4RPGAntiSigmaMinusInelastic::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 CASASM
138  // AntiSigmaMinus 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[2] = { 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-2); nneg<=np; ++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=std::max(0,np-1); nneg<=(np+1); ++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=2; np<(numSec/3); ++np )
226  {
227  nneg = np-2;
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=1; 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  const G4double expxu = 82.; // upper bound for arg. of exp
267  const G4double expxl = -expxu; // lower bound for arg. of exp
268  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
269  G4ParticleDefinition *aProton = G4Proton::Proton();
270  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
271  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
272  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
273  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
274  G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
275  G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
276  const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
277  0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
278  0.39,0.36,0.33,0.10,0.01};
279  G4int iplab = G4int( pOriginal/GeV*10.0 );
280  if( iplab > 9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
281  if( iplab > 14 )iplab = G4int( pOriginal/GeV- 2.0 ) + 15;
282  if( iplab > 23 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
283  if( iplab > 24 )iplab = 24;
284  if( G4UniformRand() > anhl[iplab] )
285  {
286  if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
287  {
288  quasiElastic = true;
289  return;
290  }
291  G4double n, anpn;
292  GetNormalizationConstant( availableEnergy, n, anpn );
293  G4double ran = G4UniformRand();
294  G4double dum, excs = 0.0;
295  if( targetParticle.GetDefinition() == aProton )
296  {
297  counter = -1;
298  for( np=0; np<numSec/3 && ran>=excs; ++np )
299  {
300  for( nneg=std::max(0,np-2); nneg<=np && ran>=excs; ++nneg )
301  {
302  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
303  {
304  if( ++counter < numMul )
305  {
306  nt = np+nneg+nz;
307  if( nt>0 && nt<=numSec )
308  {
309  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
310  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
311  if( std::fabs(dum) < 1.0 )
312  {
313  if( test >= 1.0e-10 )excs += dum*test;
314  }
315  else
316  excs += dum*test;
317  }
318  }
319  }
320  }
321  }
322  if( ran >= excs ) // 3 previous loops continued to the end
323  {
324  quasiElastic = true;
325  return;
326  }
327  np--; nneg--; nz--;
328  G4int ncht = std::min( 3, std::max( 1, np-nneg+1 ) );
329  switch( ncht )
330  {
331  case 1:
332  break;
333  case 2:
334  if( G4UniformRand() < 0.5 )
335  {
336  targetParticle.SetDefinitionAndUpdateE( aNeutron );
337  targetHasChanged = true;
338  }
339  else
340  {
341  if( G4UniformRand() < 0.5 )
342  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
343  else
344  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
345  incidentHasChanged = true;
346  }
347  break;
348  case 3:
349  if( G4UniformRand() < 0.5 )
350  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
351  else
352  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
353  incidentHasChanged = true;
354  targetParticle.SetDefinitionAndUpdateE( aNeutron );
355  targetHasChanged = true;
356  break;
357  }
358  }
359  else // target must be a neutron
360  {
361  counter = -1;
362  for( np=0; np<numSec/3 && ran>=excs; ++np )
363  {
364  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
365  {
366  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
367  {
368  if( ++counter < numMul )
369  {
370  nt = np+nneg+nz;
371  if( nt>0 && nt<=numSec )
372  {
373  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
374  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
375  if( std::fabs(dum) < 1.0 )
376  {
377  if( test >= 1.0e-10 )excs += dum*test;
378  }
379  else
380  excs += dum*test;
381  }
382  }
383  }
384  }
385  }
386  if( ran >= excs ) // 3 previous loops continued to the end
387  {
388  quasiElastic = true;
389  return;
390  }
391  np--; nneg--; nz--;
392  G4int ncht = std::min( 3, std::max( 1, np-nneg+2 ) );
393  switch( ncht )
394  {
395  case 1:
396  {
397  targetParticle.SetDefinitionAndUpdateE( aProton );
398  targetHasChanged = true;
399  }
400  break;
401  case 2:
402  if( G4UniformRand() < 0.5 )
403  {
404  if( G4UniformRand() < 0.5 )
405  {
406  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
407  incidentHasChanged = true;
408  targetParticle.SetDefinitionAndUpdateE( aProton );
409  targetHasChanged = true;
410  }
411  }
412  else
413  {
414  if( G4UniformRand() < 0.5 )
415  {
416  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
417  incidentHasChanged = true;
418  targetParticle.SetDefinitionAndUpdateE( aProton );
419  targetHasChanged = true;
420  }
421  }
422  break;
423  case 3:
424  if( G4UniformRand() < 0.5 )
425  currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
426  else
427  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
428  incidentHasChanged = true;
429  break;
430  }
431  }
432  }
433  else // random number <= anhl[iplab]
434  {
435  if( centerofmassEnergy <= aPiPlus->GetPDGMass()/MeV+aKaonPlus->GetPDGMass()/MeV )
436  {
437  quasiElastic = true;
438  return;
439  }
440  G4double n, anpn;
441  GetNormalizationConstant( -centerofmassEnergy, n, anpn );
442  G4double ran = G4UniformRand();
443  G4double dum, excs = 0.0;
444  if( targetParticle.GetDefinition() == aProton )
445  {
446  counter = -1;
447  for( np=2; np<numSec/3 && ran>=excs; ++np )
448  {
449  nneg=np-2;
450  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
451  {
452  if( ++counter < numMulA )
453  {
454  nt = np+nneg+nz;
455  if( nt>1 && nt<=numSec )
456  {
457  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
458  dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
459  if( std::fabs(dum) < 1.0 )
460  {
461  if( test >= 1.0e-10 )excs += dum*test;
462  }
463  else
464  excs += dum*test;
465  }
466  }
467  }
468  }
469  if( ran >= excs ) // 3 previous loops continued to the end
470  {
471  quasiElastic = true;
472  return;
473  }
474  np--; nz--;
475  }
476  else // target must be a neutron
477  {
478  counter = -1;
479  for( np=1; np<numSec/3 && ran>=excs; ++np )
480  {
481  nneg = np-1;
482  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
483  {
484  if( ++counter < numMulA )
485  {
486  nt = np+nneg+nz;
487  if( nt>1 && nt<=numSec )
488  {
489  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
490  dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
491  if( std::fabs(dum) < 1.0 )
492  {
493  if( test >= 1.0e-10 )excs += dum*test;
494  }
495  else
496  excs += dum*test;
497  }
498  }
499  }
500  }
501  if( ran >= excs ) // 3 previous loops continued to the end
502  {
503  quasiElastic = true;
504  return;
505  }
506  np--; nz--;
507  }
508  if( nz > 0 )
509  {
510  if( nneg > 0 )
511  {
512  if( G4UniformRand() < 0.5 )
513  {
514  vec.Initialize( 1 );
515  G4ReactionProduct *p = new G4ReactionProduct;
516  p->SetDefinition( aKaonMinus );
517  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
518  vec.SetElement( vecLen++, p );
519  --nneg;
520  }
521  else // random number >= 0.5
522  {
523  vec.Initialize( 1 );
524  G4ReactionProduct *p = new G4ReactionProduct;
525  p->SetDefinition( aKaonZL );
526  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
527  vec.SetElement( vecLen++, p );
528  --nz;
529  }
530  }
531  else // nneg == 0
532  {
533  vec.Initialize( 1 );
534  G4ReactionProduct *p = new G4ReactionProduct;
535  p->SetDefinition( aKaonZL );
536  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
537  vec.SetElement( vecLen++, p );
538  --nz;
539  }
540  }
541  else // nz == 0
542  {
543  if( nneg > 0 )
544  {
545  vec.Initialize( 1 );
546  G4ReactionProduct *p = new G4ReactionProduct;
547  p->SetDefinition( aKaonMinus );
548  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
549  vec.SetElement( vecLen++, p );
550  --nneg;
551  }
552  }
553  currentParticle.SetMass( 0.0 );
554  targetParticle.SetMass( 0.0 );
555  }
556 
557  SetUpPions( np, nneg, nz, vec, vecLen );
558  return;
559 }
560 
561  /* end of file */
562