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