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