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