ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGOmegaMinusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGOmegaMinusInelastic.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 
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 // G4double targetMass = originalTarget->GetDefinition()->GetPDGMass();
51  G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
52 
53  if( verboseLevel > 1 )
54  {
55  const G4Material *targetMaterial = aTrack.GetMaterial();
56  G4cout << "G4RPGOmegaMinusInelastic::ApplyYourself called" << G4endl;
57  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
58  G4cout << "target material = " << targetMaterial->GetName() << ", ";
59  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
60  << G4endl;
61  }
62 
63  G4ReactionProduct currentParticle(originalIncident->GetDefinition() );
64  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
65  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
66 
67  // Fermi motion and evaporation
68  // As of Geant3, the Fermi energy calculation had not been Done
69 
70  G4double ek = originalIncident->GetKineticEnergy();
71  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
72 
73  G4double tkin = targetNucleus.Cinema( ek );
74  ek += tkin;
75  currentParticle.SetKineticEnergy( ek );
76  G4double et = ek + amas;
77  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
78  G4double pp = currentParticle.GetMomentum().mag();
79  if( pp > 0.0 )
80  {
81  G4ThreeVector momentum = currentParticle.GetMomentum();
82  currentParticle.SetMomentum( momentum * (p/pp) );
83  }
84 
85  // calculate black track energies
86 
87  tkin = targetNucleus.EvaporationEffects( ek );
88  ek -= tkin;
89  currentParticle.SetKineticEnergy( ek );
90  et = ek + amas;
91  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
92  pp = currentParticle.GetMomentum().mag();
93  if( pp > 0.0 )
94  {
95  G4ThreeVector momentum = currentParticle.GetMomentum();
96  currentParticle.SetMomentum( momentum * (p/pp) );
97  }
98 
99  G4ReactionProduct modifiedOriginal = currentParticle;
100 
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*MeV;
111  if( currentParticle.GetKineticEnergy() > cutOff )
112  Cascade( vec, vecLen,
113  originalIncident, currentParticle, targetParticle,
114  incidentHasChanged, targetHasChanged, quasiElastic );
115 
116  CalculateMomenta( vec, vecLen,
117  originalIncident, originalTarget, modifiedOriginal,
118  targetNucleus, currentParticle, targetParticle,
119  incidentHasChanged, targetHasChanged, quasiElastic );
120 
121  SetUpChange( vec, vecLen,
122  currentParticle, targetParticle,
123  incidentHasChanged );
124 
125  delete originalTarget;
126  return &theParticleChange;
127 }
128 
129 
132  G4int& vecLen,
133  const G4HadProjectile *originalIncident,
134  G4ReactionProduct &currentParticle,
135  G4ReactionProduct &targetParticle,
136  G4bool &incidentHasChanged,
137  G4bool &targetHasChanged,
138  G4bool &quasiElastic )
139 {
140  // Derived from H. Fesefeldt's original FORTRAN code CASOM
141  // OmegaMinus undergoes interaction with nucleon within a nucleus. Check if it is
142  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
143  // occurs and input particle is degraded in energy. No other particles are produced.
144  // If reaction is possible, find the correct number of pions/protons/neutrons
145  // produced using an interpolation to multiplicity data. Replace some pions or
146  // protons/neutrons by kaons or strange baryons according to the average
147  // multiplicity per Inelastic reaction.
148 
149  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
150  const G4double etOriginal = originalIncident->GetTotalEnergy();
151 // const G4double pOriginal = originalIncident->GetTotalMomentum();
152  const G4double targetMass = targetParticle.GetMass();
153  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
154  targetMass*targetMass +
155  2.0*targetMass*etOriginal );
156  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
157  if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass() )
158  {
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  // np = number of pi+, nneg = number of pi-, nz = number of pi0
168  G4int counter, nt=0, np=0, nneg=0, nz=0;
169  G4double test;
170  const G4double c = 1.25;
171  const G4double b[] = { 0.70, 0.70 };
172  if( first ) // compute normalization constants, this will only be Done once
173  {
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 )
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
233 
234  // energetically possible to produce pion(s) --> inelastic scattering
235 
236  G4double n, anpn;
237  GetNormalizationConstant( availableEnergy, n, anpn );
239  G4double dum, excs = 0.0;
240  if( targetParticle.GetDefinition() == aProton )
241  {
242  counter = -1;
243  for( np=0; np<numSec/3 && ran>=excs; ++np )
244  {
245  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
246  {
247  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
248  {
249  if( ++counter < numMul )
250  {
251  nt = np+nneg+nz;
252  if( nt > 0 )
253  {
254  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
255  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
256  if( std::fabs(dum) < 1.0 )
257  {
258  if( test >= 1.0e-10 )excs += dum*test;
259  }
260  else
261  excs += dum*test;
262  }
263  }
264  }
265  }
266  }
267  if( ran >= excs ) // 3 previous loops continued to the end
268  {
269  quasiElastic = true;
270  return;
271  }
272  np--; nneg--; nz--;
273  }
274  else // target must be a neutron
275  {
276  counter = -1;
277  for( np=0; np<numSec/3 && ran>=excs; ++np )
278  {
279  for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
280  {
281  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
282  {
283  if( ++counter < numMul )
284  {
285  nt = np+nneg+nz;
286  if( (nt>=1) && (nt<=numSec) )
287  {
288  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
289  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
290  if( std::fabs(dum) < 1.0 )
291  {
292  if( test >= 1.0e-10 )excs += dum*test;
293  }
294  else
295  excs += dum*test;
296  }
297  }
298  }
299  }
300  }
301  if( ran >= excs ) // 3 previous loops continued to the end
302  {
303  quasiElastic = true;
304  return;
305  }
306  np--; nneg--; nz--;
307  }
308  // number of secondary mesons determined by kno distribution
309  // check for total charge of final state mesons to determine
310  // the kind of baryons to be produced, taking into account
311  // charge and strangeness conservation
312  //
313  G4int nvefix = 0;
314  if( targetParticle.GetDefinition() == aProton )
315  {
316  if( nneg > np )
317  {
318  if( nneg == np+1 )
319  {
320  currentParticle.SetDefinitionAndUpdateE( aXiZero );
321  nvefix = 1;
322  }
323  else
324  {
325  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
326  nvefix = 2;
327  }
328  incidentHasChanged = true;
329  }
330  else if( nneg < np )
331  {
332  targetParticle.SetDefinitionAndUpdateE( aNeutron );
333  targetHasChanged = true;
334  }
335  }
336  else // target is a neutron
337  {
338  if( np+1 < nneg )
339  {
340  if( nneg == np+2 )
341  {
342  currentParticle.SetDefinitionAndUpdateE( aXiZero );
343  incidentHasChanged = true;
344  nvefix = 1;
345  }
346  else // charge mismatch
347  {
348  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
349  incidentHasChanged = true;
350  nvefix = 2;
351  }
352  targetParticle.SetDefinitionAndUpdateE( aProton );
353  targetHasChanged = true;
354  }
355  else if( nneg == np+1 )
356  {
357  targetParticle.SetDefinitionAndUpdateE( aProton );
358  targetHasChanged = true;
359  }
360  }
361 
362  SetUpPions(np, nneg, nz, vec, vecLen);
363  for (G4int i = 0; i < vecLen && nvefix > 0; ++i) {
364  if (vec[i]->GetDefinition() == aPiMinus) {
365  if( nvefix >= 1 )vec[i]->SetDefinitionAndUpdateE(aKaonMinus);
366  --nvefix;
367  }
368  }
369 
370  return;
371 }
372 
373  /* end of file */
374