Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEOmegaMinusInelastic.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 
29 // G4 Process: Gheisha High Energy Collision model.
30 // This includes the high energy cascading model, the two-body-resonance model
31 // and the low energy two-body model. Not included are the low energy stuff
32 // like nuclear reactions, nuclear fission without any cascading and all
33 // processes for particles at rest.
34 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
35 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
36 
37 #include "G4HEOmegaMinusInelastic.hh"
38 #include "globals.hh"
39 #include "G4ios.hh"
40 #include "G4PhysicalConstants.hh"
41 
42 void G4HEOmegaMinusInelastic::ModelDescription(std::ostream& outFile) const
43 {
44  outFile << "G4HEOmegaMinusInelastic is one of the High Energy\n"
45  << "Parameterized (HEP) models used to implement inelastic\n"
46  << "Omega- scattering from nuclei. It is a re-engineered\n"
47  << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
48  << "initial collision products into backward- and forward-going\n"
49  << "clusters which are then decayed into final state hadrons.\n"
50  << "The model does not conserve energy on an event-by-event\n"
51  << "basis. It may be applied to Omega- with initial energies\n"
52  << "above 20 GeV.\n";
53 }
54 
55 
56 G4HadFinalState*
57 G4HEOmegaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
58  G4Nucleus& targetNucleus)
59 {
60  G4HEVector* pv = new G4HEVector[MAXPART];
61  const G4HadProjectile* aParticle = &aTrack;
62  const G4double A = targetNucleus.GetA_asInt();
63  const G4double Z = targetNucleus.GetZ_asInt();
64  G4HEVector incidentParticle(aParticle);
65 
66  G4double atomicNumber = Z;
67  G4double atomicWeight = A;
68 
69  G4int incidentCode = incidentParticle.getCode();
70  G4double incidentMass = incidentParticle.getMass();
71  G4double incidentTotalEnergy = incidentParticle.getEnergy();
72 
73  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
74  // DHW 19 may 2011: variable set but not used
75 
76  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
77 
78  if (incidentKineticEnergy < 1.)
79  G4cout << "GHEOmegaMinusInelastic: incident energy < 1 GeV" << G4endl;
80 
81  if (verboseLevel > 1) {
82  G4cout << "G4HEOmegaMinusInelastic::ApplyYourself" << G4endl;
83  G4cout << "incident particle " << incidentParticle.getName()
84  << "mass " << incidentMass
85  << "kinetic energy " << incidentKineticEnergy
86  << G4endl;
87  G4cout << "target material with (A,Z) = ("
88  << atomicWeight << "," << atomicNumber << ")" << G4endl;
89  }
90 
91  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
92  atomicWeight, atomicNumber);
93  if (verboseLevel > 1)
94  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
95 
96  incidentKineticEnergy -= inelasticity;
97 
98  G4double excitationEnergyGNP = 0.;
99  G4double excitationEnergyDTA = 0.;
100 
101  G4double excitation = NuclearExcitation(incidentKineticEnergy,
102  atomicWeight, atomicNumber,
103  excitationEnergyGNP,
104  excitationEnergyDTA);
105  if (verboseLevel > 1)
106  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
107  << excitationEnergyDTA << G4endl;
108 
109  incidentKineticEnergy -= excitation;
110  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
111  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
112  // *(incidentTotalEnergy+incidentMass));
113  // DHW 19 May 2011: variable set but not used
114 
115  G4HEVector targetParticle;
116  if (G4UniformRand() < atomicNumber/atomicWeight) {
117  targetParticle.setDefinition("Proton");
118  } else {
119  targetParticle.setDefinition("Neutron");
120  }
121 
122  G4double targetMass = targetParticle.getMass();
123  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
124  + targetMass*targetMass
125  + 2.0*targetMass*incidentTotalEnergy);
126  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
127 
128  G4bool inElastic = true;
129  vecLength = 0;
130 
131  if (verboseLevel > 1)
132  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
133  << incidentCode << G4endl;
134 
135  G4bool successful = false;
136 
137  FirstIntInCasOmegaMinus(inElastic, availableEnergy, pv, vecLength,
138  incidentParticle, targetParticle, atomicWeight);
139 
140  if (verboseLevel > 1)
141  G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
142 
143  if ((vecLength > 0) && (availableEnergy > 1.))
144  StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
145  pv, vecLength,
146  incidentParticle, targetParticle);
147 
148  HighEnergyCascading(successful, pv, vecLength,
149  excitationEnergyGNP, excitationEnergyDTA,
150  incidentParticle, targetParticle,
151  atomicWeight, atomicNumber);
152  if (!successful)
153  HighEnergyClusterProduction(successful, pv, vecLength,
154  excitationEnergyGNP, excitationEnergyDTA,
155  incidentParticle, targetParticle,
156  atomicWeight, atomicNumber);
157  if (!successful)
158  MediumEnergyCascading(successful, pv, vecLength,
159  excitationEnergyGNP, excitationEnergyDTA,
160  incidentParticle, targetParticle,
161  atomicWeight, atomicNumber);
162 
163  if (!successful)
164  MediumEnergyClusterProduction(successful, pv, vecLength,
165  excitationEnergyGNP, excitationEnergyDTA,
166  incidentParticle, targetParticle,
167  atomicWeight, atomicNumber);
168  if (!successful)
169  QuasiElasticScattering(successful, pv, vecLength,
170  excitationEnergyGNP, excitationEnergyDTA,
171  incidentParticle, targetParticle,
172  atomicWeight, atomicNumber);
173  if (!successful)
174  ElasticScattering(successful, pv, vecLength,
175  incidentParticle,
176  atomicWeight, atomicNumber);
177 
178  if (!successful)
179  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
180  << G4endl;
181 
182  FillParticleChange(pv, vecLength);
183  delete [] pv;
184  theParticleChange.SetStatusChange(stopAndKill);
185  return &theParticleChange;
186 }
187 
188 
189 void
190 G4HEOmegaMinusInelastic::FirstIntInCasOmegaMinus(G4bool& inElastic,
191  const G4double availableEnergy,
192  G4HEVector pv[],
193  G4int& vecLen,
194  const G4HEVector& incidentParticle,
195  const G4HEVector& targetParticle,
196  const G4double atomicWeight)
197 
198 // Xi0 undergoes interaction with nucleon within a nucleus. Check if it is
199 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
200 // occurs and input particle is degraded in energy. No other particles are produced.
201 // If reaction is possible, find the correct number of pions/protons/neutrons
202 // produced using an interpolation to multiplicity data. Replace some pions or
203 // protons/neutrons by kaons or strange baryons according to the average
204 // multiplicity per inelastic reaction.
205 {
206  static const G4double expxu = 82.; // upper bound for arg. of exp
207  static const G4double expxl = -expxu; // lower bound for arg. of exp
208 
209  static const G4double protb = 0.7;
210  static const G4double neutb = 0.7;
211  static const G4double c = 1.25;
212 
213  static const G4int numMul = 1200;
214  static const G4int numSec = 60;
215 
216  G4int protonCode = Proton.getCode();
217 
218  G4int targetCode = targetParticle.getCode();
219  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
220 
221  static G4bool first = true;
222  static G4double protmul[numMul], protnorm[numSec]; // proton constants
223  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
224 
225  // misc. local variables
226  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
227 
228  G4int i, counter, nt, npos, nneg, nzero;
229 
230  if (first) {
231  // compute normalization constants, this will only be done once
232  first = false;
233  for( i=0; i<numMul; i++ )protmul[i] = 0.0;
234  for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
235  counter = -1;
236  for( npos=0; npos<(numSec/3); npos++ )
237  {
238  for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ )
239  {
240  for( nzero=0; nzero<numSec/3; nzero++ )
241  {
242  if( ++counter < numMul )
243  {
244  nt = npos+nneg+nzero;
245  if( (nt>0) && (nt<=numSec) )
246  {
247  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
248  protnorm[nt-1] += protmul[counter];
249  }
250  }
251  }
252  }
253  }
254  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
255  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
256  counter = -1;
257  for( npos=0; npos<numSec/3; npos++ )
258  {
259  for( nneg=npos; nneg<=(npos+1); nneg++ )
260  {
261  for( nzero=0; nzero<numSec/3; nzero++ )
262  {
263  if( ++counter < numMul )
264  {
265  nt = npos+nneg+nzero;
266  if( (nt>0) && (nt<=numSec) )
267  {
268  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
269  neutnorm[nt-1] += neutmul[counter];
270  }
271  }
272  }
273  }
274  }
275  for( i=0; i<numSec; i++ )
276  {
277  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
278  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
279  }
280  } // end of initialization
281 
282 
283  // initialize the first two places
284  // the same as beam and target
285  pv[0] = incidentParticle;
286  pv[1] = targetParticle;
287  vecLen = 2;
288 
289  if( !inElastic )
290  { // quasi-elastic scattering, no pions produced
291  G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
292  G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
293  if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
294  {
295  G4double ran = G4UniformRand();
296  if( targetCode == protonCode)
297  {
298  if (ran < 0.2)
299  {
300  pv[0] = XiZero;
301  pv[1] = SigmaZero;
302  }
303  else if (ran < 0.4)
304  {
305  pv[0] = SigmaZero;
306  pv[1] = XiZero;
307  }
308  else if (ran < 0.6)
309  {
310  pv[0] = XiZero;
311  pv[1] = Lambda;
312  }
313  else if (ran < 0.8)
314  {
315  pv[0] = Lambda;
316  pv[1] = XiZero;
317  }
318  else
319  {
320  pv[0] = Proton;
321  pv[1] = OmegaMinus;
322  }
323  }
324  else
325  {
326  if (ran < 0.2)
327  {
328  pv[0] = Neutron;
329  pv[1] = OmegaMinus;
330  }
331  else if (ran < 0.4)
332  {
333  pv[0] = XiZero;
334  pv[1] = SigmaMinus;
335  }
336  else if (ran < 0.6)
337  {
338  pv[0] = SigmaMinus;
339  pv[1] = XiZero;
340  }
341  else if (ran < 0.8)
342  {
343  pv[0] = XiMinus;
344  pv[1] = Lambda;
345  }
346  else
347  {
348  pv[0] = Lambda;
349  pv[1] = XiMinus;
350  }
351  }
352  }
353  return;
354  }
355  else if (availableEnergy <= PionPlus.getMass())
356  return;
357 
358  // inelastic scattering
359  npos = 0; nneg = 0; nzero = 0;
360 
361  // number of total particles vs. centre of mass Energy - 2*proton mass
362  G4double aleab = std::log(availableEnergy);
363  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
364  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
365 
366  // normalization constant for kno-distribution.
367  // calculate first the sum of all constants, check for numerical problems.
368  G4double test, dum, anpn = 0.0;
369 
370  for (nt=1; nt<=numSec; nt++) {
371  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
372  dum = pi*nt/(2.0*n*n);
373  if (std::fabs(dum) < 1.0) {
374  if( test >= 1.0e-10 )anpn += dum*test;
375  } else {
376  anpn += dum*test;
377  }
378  }
379 
380  G4double ran = G4UniformRand();
381  G4double excs = 0.0;
382  if (targetCode == protonCode) {
383  counter = -1;
384  for (npos = 0; npos < numSec/3; npos++) {
385  for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ )
386  {
387  for( nzero=0; nzero<numSec/3; nzero++ )
388  {
389  if( ++counter < numMul )
390  {
391  nt = npos+nneg+nzero;
392  if ( (nt>0) && (nt<=numSec) ) {
393  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
394  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
395  if (std::fabs(dum) < 1.0) {
396  if( test >= 1.0e-10 )excs += dum*test;
397  } else {
398  excs += dum*test;
399  }
400  if (ran < excs) goto outOfLoop; //----------------------->
401  }
402  }
403  }
404  }
405  }
406 
407  // 3 previous loops continued to the end
408  inElastic = false; // quasi-elastic scattering
409  return;
410  }
411  else
412  { // target must be a neutron
413  counter = -1;
414  for( npos=0; npos<numSec/3; npos++ )
415  {
416  for( nneg=npos; nneg<=(npos+1); nneg++ )
417  {
418  for( nzero=0; nzero<numSec/3; nzero++ )
419  {
420  if( ++counter < numMul )
421  {
422  nt = npos+nneg+nzero;
423  if ( (nt>=1) && (nt<=numSec) ) {
424  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
425  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
426  if (std::fabs(dum) < 1.0) {
427  if( test >= 1.0e-10 )excs += dum*test;
428  } else {
429  excs += dum*test;
430  }
431  if (ran < excs) goto outOfLoop; // -------------------------->
432  }
433  }
434  }
435  }
436  }
437  // 3 previous loops continued to the end
438  inElastic = false; // quasi-elastic scattering.
439  return;
440  }
441 
442  outOfLoop: // <------------------------------------------------------------------------
443 
444  // in the following we do not consider strangeness transfer at high
445  // multiplicity events. YK combinations are added in
446  // StrangeParticlePairProduction
447  ran = G4UniformRand();
448  if (targetCode == protonCode) {
449  if( npos == nneg)
450  {
451  }
452  else
453  {
454  pv[1] = Neutron;
455  }
456  } else {
457  if (npos == nneg)
458  {
459  }
460  else
461  {
462  pv[1] = Proton;
463  }
464  }
465 
466  nt = npos + nneg + nzero;
467  while (nt > 0) {
468  G4double rnd = G4UniformRand();
469  if (rnd < (G4double)npos/nt) {
470  if (npos > 0) {
471  pv[vecLen++] = PionPlus;
472  npos--;
473  }
474  } else if (rnd < (G4double)(npos+nneg)/nt) {
475  if (nneg > 0) {
476  pv[vecLen++] = PionMinus;
477  nneg--;
478  }
479  } else {
480  if (nzero > 0) {
481  pv[vecLen++] = PionZero;
482  nzero--;
483  }
484  }
485  nt = npos + nneg + nzero;
486  }
487 
488  if (verboseLevel > 1) {
489  G4cout << "Particles produced: " ;
490  G4cout << pv[0].getName() << " " ;
491  G4cout << pv[1].getName() << " " ;
492  for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
493  G4cout << G4endl;
494  }
495  return;
496 }
497 
Generated on Sat May 25 2013 14:33:56 for Geant4 by   doxygen 1.8.4