ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiXiMinusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiXiMinusInelastic.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 // NOTE: The FORTRAN version of the cascade, CASAXM, simply called the
29 // routine for the XiMinus particle. Hence, the ApplyYourself function
30 // below is just a copy of the ApplyYourself from the XiMinus particle.
31 
33 #include "G4Exp.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh"
37 
40  G4Nucleus &targetNucleus )
41 {
42  const G4HadProjectile *originalIncident = &aTrack;
43  if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
44  {
48  return &theParticleChange;
49  }
50 
51  // create the target particle
52 
53  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
54 
55  if( verboseLevel > 1 )
56  {
57  const G4Material *targetMaterial = aTrack.GetMaterial();
58  G4cout << "G4RPGAntiXiMinusInelastic::ApplyYourself called" << G4endl;
59  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
60  G4cout << "target material = " << targetMaterial->GetName() << ", ";
61  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
62  << G4endl;
63  }
64  //
65  // Fermi motion and evaporation
66  // As of Geant3, the Fermi energy calculation had not been Done
67  //
68  G4double ek = originalIncident->GetKineticEnergy()/MeV;
69  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
70  G4ReactionProduct modifiedOriginal;
71  modifiedOriginal = *originalIncident;
72 
73  G4double tkin = targetNucleus.Cinema( ek );
74  ek += tkin;
75  modifiedOriginal.SetKineticEnergy( ek*MeV );
76  G4double et = ek + amas;
77  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
78  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
79  if( pp > 0.0 )
80  {
81  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
82  modifiedOriginal.SetMomentum( momentum * (p/pp) );
83  }
84  //
85  // calculate black track energies
86  //
87  tkin = targetNucleus.EvaporationEffects( ek );
88  ek -= tkin;
89  modifiedOriginal.SetKineticEnergy( ek*MeV );
90  et = ek + amas;
91  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
92  pp = modifiedOriginal.GetMomentum().mag()/MeV;
93  if( pp > 0.0 )
94  {
95  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
96  modifiedOriginal.SetMomentum( momentum * (p/pp) );
97  }
98  G4ReactionProduct currentParticle = modifiedOriginal;
99  G4ReactionProduct targetParticle;
100  targetParticle = *originalTarget;
101  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
102  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
103  G4bool incidentHasChanged = false;
104  G4bool targetHasChanged = false;
105  G4bool quasiElastic = false;
106  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
107  G4int vecLen = 0;
108  vec.Initialize( 0 );
109 
110  const G4double cutOff = 0.1;
111  const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
112  if( (currentParticle.GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
113  Cascade( vec, vecLen,
114  originalIncident, currentParticle, targetParticle,
115  incidentHasChanged, targetHasChanged, quasiElastic );
116 
117  CalculateMomenta( vec, vecLen,
118  originalIncident, originalTarget, modifiedOriginal,
119  targetNucleus, currentParticle, targetParticle,
120  incidentHasChanged, targetHasChanged, quasiElastic );
121 
122  SetUpChange( vec, vecLen,
123  currentParticle, targetParticle,
124  incidentHasChanged );
125 
126  delete originalTarget;
127  return &theParticleChange;
128 }
129 
130 
133  G4int& vecLen,
134  const G4HadProjectile *originalIncident,
135  G4ReactionProduct &currentParticle,
136  G4ReactionProduct &targetParticle,
137  G4bool &incidentHasChanged,
138  G4bool &targetHasChanged,
139  G4bool &quasiElastic )
140 {
141  // Derived from H. Fesefeldt's original FORTRAN code CASAXM
142  // which is just a copy of casxm (cascade for Xi-).
143  // AntiXiMinus 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 targetMass = targetParticle.GetMass()/MeV;
154  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
155  targetMass*targetMass +
156  2.0*targetMass*etOriginal );
157  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
158  if (availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV) {
159  quasiElastic = true;
160  return;
161  }
162  static G4ThreadLocal G4bool first = true;
163  const G4int numMul = 1200;
164  const G4int numSec = 60;
165  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
166  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
167 
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  if (first) { // Computation of normalization constants will only be done once
174  first = false;
175  G4int i;
176  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
177  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
178  counter = -1;
179  for( np=0; np<(numSec/3); ++np )
180  {
181  for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
182  {
183  for( nz=0; nz<numSec/3; ++nz )
184  {
185  if( ++counter < numMul )
186  {
187  nt = np+nneg+nz;
188  if( nt>0 && nt<=numSec )
189  {
190  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
191  protnorm[nt-1] += protmul[counter];
192  }
193  }
194  }
195  }
196  }
197  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
198  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
199  counter = -1;
200  for( np=0; np<numSec/3; ++np )
201  {
202  for( nneg=np; nneg<=(np+2); ++nneg )
203  {
204  for( nz=0; nz<numSec/3; ++nz )
205  {
206  if( ++counter < numMul )
207  {
208  nt = np+nneg+nz;
209  if( nt>0 && nt<=numSec )
210  {
211  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
212  neutnorm[nt-1] += neutmul[counter];
213  }
214  }
215  }
216  }
217  }
218  for( i=0; i<numSec; ++i )
219  {
220  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
221  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
222  }
223  } // end of initialization
224 
225  const G4double expxu = 82.; // upper bound for arg. of exp
226  const G4double expxl = -expxu; // lower bound for arg. of exp
232  //
233  // energetically possible to produce pion(s) --> inelastic scattering
234  //
235  G4double n, anpn;
236  GetNormalizationConstant( availableEnergy, n, anpn );
238  G4double dum, excs = 0.0;
239  if( targetParticle.GetDefinition() == aProton )
240  {
241  counter = -1;
242  for( np=0; np<numSec/3 && ran>=excs; ++np )
243  {
244  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
245  {
246  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
247  {
248  if( ++counter < numMul )
249  {
250  nt = np+nneg+nz;
251  if( nt>0 && nt<=numSec )
252  {
253  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
254  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
255  if( std::fabs(dum) < 1.0 )
256  {
257  if( test >= 1.0e-10 )excs += dum*test;
258  }
259  else
260  excs += dum*test;
261  }
262  }
263  }
264  }
265  }
266  if( ran >= excs ) // 3 previous loops continued to the end
267  {
268  quasiElastic = true;
269  return;
270  }
271  np--; nneg--; nz--;
272  //
273  // number of secondary mesons determined by kno distribution
274  // check for total charge of final state mesons to determine
275  // the kind of baryons to be produced, taking into account
276  // charge and strangeness conservation
277  //
278  if( np < nneg )
279  {
280  if( np+1 == nneg )
281  {
282  currentParticle.SetDefinitionAndUpdateE( aXiZero );
283  incidentHasChanged = true;
284  }
285  else // charge mismatch
286  {
287  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
288  incidentHasChanged = true;
289  //
290  // correct the strangeness by replacing a pi- by a kaon-
291  //
292  vec.Initialize( 1 );
294  p->SetDefinition( aKaonMinus );
295  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
296  vec.SetElement( vecLen++, p );
297  --nneg;
298  }
299  }
300  else if( np == nneg )
301  {
302  if( G4UniformRand() >= 0.5 )
303  {
304  currentParticle.SetDefinitionAndUpdateE( aXiZero );
305  incidentHasChanged = true;
306  targetParticle.SetDefinitionAndUpdateE( aNeutron );
307  targetHasChanged = true;
308  }
309  }
310  else
311  {
312  targetParticle.SetDefinitionAndUpdateE( aNeutron );
313  targetHasChanged = true;
314  }
315  }
316  else // target must be a neutron
317  {
318  counter = -1;
319  for( np=0; np<numSec/3 && ran>=excs; ++np )
320  {
321  for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
322  {
323  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
324  {
325  if( ++counter < numMul )
326  {
327  nt = np+nneg+nz;
328  if( nt>0 && nt<=numSec )
329  {
330  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
331  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
332  if( std::fabs(dum) < 1.0 )
333  {
334  if( test >= 1.0e-10 )excs += dum*test;
335  }
336  else
337  excs += dum*test;
338  }
339  }
340  }
341  }
342  }
343  if( ran >= excs ) // 3 previous loops continued to the end
344  {
345  quasiElastic = true;
346  return;
347  }
348  np--; nneg--; nz--;
349  if( np+1 < nneg )
350  {
351  if( np+2 == nneg )
352  {
353  currentParticle.SetDefinitionAndUpdateE( aXiZero );
354  incidentHasChanged = true;
355  targetParticle.SetDefinitionAndUpdateE( aProton );
356  targetHasChanged = true;
357  }
358  else // charge mismatch
359  {
360  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
361  incidentHasChanged = true;
362  targetParticle.SetDefinitionAndUpdateE( aProton );
363  targetHasChanged = true;
364  //
365  // correct the strangeness by replacing a pi- by a kaon-
366  //
367  vec.Initialize( 1 );
369  p->SetDefinition( aKaonMinus );
370  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
371  vec.SetElement( vecLen++, p );
372  --nneg;
373  }
374  }
375  else if( np+1 == nneg )
376  {
377  if( G4UniformRand() < 0.5 )
378  {
379  currentParticle.SetDefinitionAndUpdateE( aXiZero );
380  incidentHasChanged = true;
381  }
382  else
383  {
384  targetParticle.SetDefinitionAndUpdateE( aProton );
385  targetHasChanged = true;
386  }
387  }
388  }
389 
390  SetUpPions(np, nneg, nz, vec, vecLen);
391  return;
392 }
393 
394  /* end of file */
395