ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiOmegaMinusInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiOmegaMinusInelastic.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, CASAOM, simply called the
29 // routine for the OmegaMinus particle. Hence, the Cascade function
30 // below is just a copy of the Cascade from the OmegaMinus 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 << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
59  G4cout << "target material = " << targetMaterial->GetName() << ", ";
60  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
61  << G4endl;
62  }
63  //
64  // Fermi motion and evaporation
65  // As of Geant3, the Fermi energy calculation had not been Done
66  //
67  G4double ek = originalIncident->GetKineticEnergy()/MeV;
68  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
69  G4ReactionProduct modifiedOriginal;
70  modifiedOriginal = *originalIncident;
71 
72  G4double tkin = targetNucleus.Cinema( ek );
73  ek += tkin;
74  modifiedOriginal.SetKineticEnergy( ek*MeV );
75  G4double et = ek + amas;
76  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
77  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
78  if( pp > 0.0 )
79  {
80  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
81  modifiedOriginal.SetMomentum( momentum * (p/pp) );
82  }
83  //
84  // calculate black track energies
85  //
86  tkin = targetNucleus.EvaporationEffects( ek );
87  ek -= tkin;
88  modifiedOriginal.SetKineticEnergy( ek*MeV );
89  et = ek + amas;
90  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
91  pp = modifiedOriginal.GetMomentum().mag()/MeV;
92  if( pp > 0.0 )
93  {
94  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
95  modifiedOriginal.SetMomentum( momentum * (p/pp) );
96  }
97  G4ReactionProduct currentParticle = modifiedOriginal;
98  G4ReactionProduct targetParticle;
99  targetParticle = *originalTarget;
100  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
101  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
102  G4bool incidentHasChanged = false;
103  G4bool targetHasChanged = false;
104  G4bool quasiElastic = false;
105  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
106  G4int vecLen = 0;
107  vec.Initialize( 0 );
108 
109  const G4double cutOff = 0.1;
110  const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
111 
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 
131 void
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 CASOM
142  // AntiOmegaMinus undergoes interaction with nucleon within a nucleus. Check if it is
143  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
144  // occurs and input particle is degraded in energy. No other particles are produced.
145  // If reaction is possible, find the correct number of pions/protons/neutrons
146  // produced using an interpolation to multiplicity data. Replace some pions or
147  // protons/neutrons by kaons or strange baryons according to the average
148  // multiplicity per Inelastic reaction.
149 
150  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
151  const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
152  const G4double targetMass = targetParticle.GetMass()/MeV;
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()/MeV) {
158  // not energetically possible to produce pion(s)
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  G4double n, anpn;
233  GetNormalizationConstant( availableEnergy, n, anpn );
235  G4double dum, excs = 0.0;
236  G4int nvefix = 0;
237  if( targetParticle.GetDefinition() == aProton )
238  {
239  counter = -1;
240  for( np=0; np<numSec/3 && ran>=excs; ++np )
241  {
242  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
243  {
244  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
245  {
246  if( ++counter < numMul )
247  {
248  nt = np+nneg+nz;
249  if( nt>0 && nt<=numSec )
250  {
251  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
252  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
253  if( std::fabs(dum) < 1.0 )
254  {
255  if( test >= 1.0e-10 )excs += dum*test;
256  }
257  else
258  excs += dum*test;
259  }
260  }
261  }
262  }
263  }
264  if( ran >= excs ) // 3 previous loops continued to the end
265  {
266  quasiElastic = true;
267  return;
268  }
269  np--; nneg--; nz--;
270  //
271  // number of secondary mesons determined by kno distribution
272  // check for total charge of final state mesons to determine
273  // the kind of baryons to be produced, taking into account
274  // charge and strangeness conservation
275  //
276  if( np < nneg )
277  {
278  if( np+1 == nneg )
279  {
280  currentParticle.SetDefinitionAndUpdateE( aXiZero );
281  incidentHasChanged = true;
282  nvefix = 1;
283  }
284  else // charge mismatch
285  {
286  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
287  incidentHasChanged = true;
288  nvefix = 2;
289  }
290  }
291  else if( np > nneg )
292  {
293  targetParticle.SetDefinitionAndUpdateE( aNeutron );
294  targetHasChanged = true;
295  }
296  }
297  else // target must be a neutron
298  {
299  counter = -1;
300  for( np=0; np<numSec/3 && ran>=excs; ++np )
301  {
302  for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
303  {
304  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
305  {
306  if( ++counter < numMul )
307  {
308  nt = np+nneg+nz;
309  if( nt>0 && nt<=numSec )
310  {
311  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
312  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
313  if( std::fabs(dum) < 1.0 )
314  {
315  if( test >= 1.0e-10 )excs += dum*test;
316  }
317  else
318  excs += dum*test;
319  }
320  }
321  }
322  }
323  }
324  if( ran >= excs ) // 3 previous loops continued to the end
325  {
326  quasiElastic = true;
327  return;
328  }
329  np--; nneg--; nz--;
330  if( np+1 < nneg )
331  {
332  if( np+2 == nneg )
333  {
334  currentParticle.SetDefinitionAndUpdateE( aXiZero );
335  incidentHasChanged = true;
336  nvefix = 1;
337  }
338  else // charge mismatch
339  {
340  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
341  incidentHasChanged = true;
342  nvefix = 2;
343  }
344  targetParticle.SetDefinitionAndUpdateE( aProton );
345  targetHasChanged = true;
346  }
347  else if( np+1 == nneg )
348  {
349  targetParticle.SetDefinitionAndUpdateE( aProton );
350  targetHasChanged = true;
351  }
352  }
353 
354  SetUpPions(np, nneg, nz, vec, vecLen);
355  for (G4int i = 0; i < vecLen && nvefix > 0; ++i) {
356  if (vec[i]->GetDefinition() == G4PionMinus::PionMinus() ) {
357  // correct the strangeness by replacing a pi- by a kaon-
358  if( nvefix >= 1 )vec[i]->SetDefinitionAndUpdateE( aKaonMinus );
359  --nvefix;
360  }
361  }
362 
363  return;
364 }
365 
366  /* end of file */
367