ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGXiMinusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGXiMinusInelastic.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 "G4RPGXiMinusInelastic.hh"
29 #include "G4Exp.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "Randomize.hh"
33 
36  G4Nucleus &targetNucleus )
37 {
38  const G4HadProjectile *originalIncident = &aTrack;
39  if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
40  {
44  return &theParticleChange;
45  }
46 
47  // create the target particle
48 
49  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
50 
51  if( verboseLevel > 1 )
52  {
53  const G4Material *targetMaterial = aTrack.GetMaterial();
54  G4cout << "G4RPGXiMinusInelastic::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  if( currentParticle.GetKineticEnergy()/MeV > cutOff )
108  Cascade( vec, vecLen,
109  originalIncident, currentParticle, targetParticle,
110  incidentHasChanged, targetHasChanged, quasiElastic );
111 
112  CalculateMomenta( vec, vecLen,
113  originalIncident, originalTarget, modifiedOriginal,
114  targetNucleus, currentParticle, targetParticle,
115  incidentHasChanged, targetHasChanged, quasiElastic );
116 
117  SetUpChange( vec, vecLen,
118  currentParticle, targetParticle,
119  incidentHasChanged );
120 
121  delete originalTarget;
122  return &theParticleChange;
123 }
124 
125 
126 void
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 CASXM
137  //
138  // XiMinus 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 targetMass = targetParticle.GetMass()/MeV;
149  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
150  targetMass*targetMass +
151  2.0*targetMass*etOriginal );
152  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
153  if (availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV) {
154  quasiElastic = true;
155  return;
156  }
157  static G4ThreadLocal G4bool first = true;
158  const G4int numMul = 1200;
159  const G4int numSec = 60;
160  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
161  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
162 
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) { // Computation of normalization constants will only be done once
169  first = false;
170  G4int i;
171  for (i = 0; i < numMul; ++i) protmul[i] = 0.0;
172  for (i = 0; i < numSec; ++i) protnorm[i] = 0.0;
173  counter = -1;
174  for (np = 0; np < (numSec/3); ++np) {
175  for (nneg = std::max(0,np-1); nneg <= (np+1); ++nneg) {
176  for (nz = 0; nz < numSec/3; ++nz) {
177  if (++counter < numMul) {
178  nt = np + nneg + nz;
179  if (nt > 0 && nt <= numSec) {
180  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
181  protnorm[nt-1] += protmul[counter];
182  }
183  }
184  }
185  }
186  }
187 
188  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
189  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
190  counter = -1;
191  for( np=0; np<numSec/3; ++np )
192  {
193  for( nneg=np; nneg<=(np+2); ++nneg )
194  {
195  for( nz=0; nz<numSec/3; ++nz )
196  {
197  if( ++counter < numMul )
198  {
199  nt = np+nneg+nz;
200  if( nt>0 && nt<=numSec )
201  {
202  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
203  neutnorm[nt-1] += neutmul[counter];
204  }
205  }
206  }
207  }
208  }
209  for( i=0; i<numSec; ++i )
210  {
211  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
212  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
213  }
214  } // end of initialization
215 
216  const G4double expxu = 82.; // upper bound for arg. of exp
217  const G4double expxl = -expxu; // lower bound for arg. of exp
223  //
224  // energetically possible to produce pion(s) --> inelastic scattering
225  //
226  G4double n, anpn;
227  GetNormalizationConstant( availableEnergy, n, anpn );
229  G4double dum, excs = 0.0;
230  if( targetParticle.GetDefinition() == aProton )
231  {
232  counter = -1;
233  for( np=0; np<numSec/3 && ran>=excs; ++np )
234  {
235  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
236  {
237  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
238  {
239  if( ++counter < numMul )
240  {
241  nt = np+nneg+nz;
242  if( nt>0 && nt<=numSec )
243  {
244  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
245  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
246  if( std::fabs(dum) < 1.0 )
247  {
248  if( test >= 1.0e-10 )excs += dum*test;
249  }
250  else
251  excs += dum*test;
252  }
253  }
254  }
255  }
256  }
257  if( ran >= excs ) // 3 previous loops continued to the end
258  {
259  quasiElastic = true;
260  return;
261  }
262  np--; nneg--; nz--;
263  //
264  // number of secondary mesons determined by kno distribution
265  // check for total charge of final state mesons to determine
266  // the kind of baryons to be produced, taking into account
267  // charge and strangeness conservation
268  //
269  if( np < nneg )
270  {
271  if( np+1 == nneg )
272  {
273  currentParticle.SetDefinitionAndUpdateE( aXiZero );
274  incidentHasChanged = true;
275  }
276  else // charge mismatch
277  {
278  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
279  incidentHasChanged = true;
280  //
281  // correct the strangeness by replacing a pi- by a kaon-
282  //
283  vec.Initialize( 1 );
285  p->SetDefinition( aKaonMinus );
286  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
287  vec.SetElement( vecLen++, p );
288  --nneg;
289  }
290  }
291  else if( np == nneg )
292  {
293  if( G4UniformRand() >= 0.5 )
294  {
295  currentParticle.SetDefinitionAndUpdateE( aXiZero );
296  incidentHasChanged = true;
297  targetParticle.SetDefinitionAndUpdateE( aNeutron );
298  targetHasChanged = true;
299  }
300  }
301  else
302  {
303  targetParticle.SetDefinitionAndUpdateE( aNeutron );
304  targetHasChanged = true;
305  }
306  }
307  else // target must be a neutron
308  {
309  counter = -1;
310  for( np=0; np<numSec/3 && ran>=excs; ++np )
311  {
312  for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
313  {
314  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
315  {
316  if( ++counter < numMul )
317  {
318  nt = np+nneg+nz;
319  if( nt>0 && nt<=numSec )
320  {
321  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
322  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
323  if( std::fabs(dum) < 1.0 )
324  {
325  if( test >= 1.0e-10 )excs += dum*test;
326  }
327  else
328  excs += dum*test;
329  }
330  }
331  }
332  }
333  }
334  if( ran >= excs ) // 3 previous loops continued to the end
335  {
336  quasiElastic = true;
337  return;
338  }
339  np--; nneg--; nz--;
340  if( np+1 < nneg )
341  {
342  if( np+2 == nneg )
343  {
344  currentParticle.SetDefinitionAndUpdateE( aXiZero );
345  incidentHasChanged = true;
346  targetParticle.SetDefinitionAndUpdateE( aProton );
347  targetHasChanged = true;
348  }
349  else // charge mismatch
350  {
351  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
352  incidentHasChanged = true;
353  targetParticle.SetDefinitionAndUpdateE( aProton );
354  targetHasChanged = true;
355  //
356  // correct the strangeness by replacing a pi- by a kaon-
357  //
358  vec.Initialize( 1 );
360  p->SetDefinition( aKaonMinus );
361  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
362  vec.SetElement( vecLen++, p );
363  --nneg;
364  }
365  }
366  else if( np+1 == nneg )
367  {
368  if( G4UniformRand() < 0.5 )
369  {
370  currentParticle.SetDefinitionAndUpdateE( aXiZero );
371  incidentHasChanged = true;
372  }
373  else
374  {
375  targetParticle.SetDefinitionAndUpdateE( aProton );
376  targetHasChanged = true;
377  }
378  }
379  }
380  SetUpPions(np, nneg, nz, vec, vecLen);
381  return;
382 }
383 
384  /* end of file */
385