Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ElementaryParticleCollider.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 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100316 D. Wright (restored by M. Kelsey) -- Replace original (incorrect)
30 // pp, pn, nn 2-body to 2-body scattering angular distributions
31 // with a new parametrization of elastic scattering data using
32 // the sum of two exponentials.
33 // 20100319 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
34 // eliminate some unnecessary std::pow()
35 // 20100407 M. Kelsey -- Replace std::vector<>::resize(0) with ::clear()
36 // Eliminate return-by-value std::vector<> by creating buffers
37 // buffers for particles, momenta, and particle types.
38 // The following functions now return void and are non-const:
39 // ::generateSCMfinalState()
40 // ::generateMomModules()
41 // ::generateStrangeChannelPartTypes()
42 // ::generateSCMpionAbsorption()
43 // 20100408 M. Kelsey -- Follow changes to G4*Sampler to pass particle_kinds
44 // as input buffer.
45 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
46 // 20100428 M. Kelsey -- Use G4InuclParticleNames enum
47 // 20100429 M. Kelsey -- Change "photon()" to "isPhoton()"
48 // 20100507 M. Kelsey -- Rationalize multiplicity returns to be actual value
49 // 20100511 M. Kelsey -- Replace G4PionSampler and G4NucleonSampler with new
50 // pi-N and N-N classes, reorganize if-cascades
51 // 20100511 M. Kelsey -- Eliminate three residual random-angles blocks.
52 // 20100511 M. Kelsey -- Bug fix: pi-N two-body final states not correctly
53 // tested for charge-exchange case.
54 // 20100512 M. Kelsey -- Rationalize multiplicity returns to be actual value
55 // 20100512 M. Kelsey -- Add some additional debugging messages for 2-to-2
56 // 20100512 M. Kelsey -- Replace "if (is==)" cascades with switch blocks.
57 // Use G4CascadeInterpolator for angular distributions.
58 // 20100517 M. Kelsey -- Inherit from common base class, make arrays static
59 // 20100519 M. Kelsey -- Use G4InteractionCase to compute "is" values.
60 // 20100625 M. Kelsey -- Two bugs in n-body momentum, last particle recoil
61 // 20100713 M. Kelsey -- Bump collide start message up to verbose > 1
62 // 20100714 M. Kelsey -- Move conservation checking to base class
63 // 20100714 M. Kelsey -- Add sanity check for two-body final state, to ensure
64 // that final state total mass is below etot_scm; also compute
65 // kinematics without "rescaling" (which led to non-conservation)
66 // 20100726 M. Kelsey -- Move remaining std::vector<> buffers to .hh file
67 // 20100804 M. Kelsey -- Add printing of final-state tables, protected by
68 // G4CASCADE_DEBUG_SAMPLER preprocessor flag
69 // 20101019 M. Kelsey -- CoVerity report: check dynamic_cast<> for null
70 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
71 // 20110718 V. Uzhinskiy -- Drop IL,IM reset for multiplicity-three collisions
72 // 20110718 M. Kelsey -- Use enum names in switch blocks (c.f. G4NucleiModel)
73 // 20110720 M. Kelsey -- Follow interface change for cross-section tables,
74 // eliminating switch blocks.
75 // 20110806 M. Kelsey -- Pre-allocate buffers to avoid memory churn
76 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
77 // 20110926 M. Kelsey -- Protect sampleCMcosFor2to2 from unphysical arguments
78 // 20111003 M. Kelsey -- Prepare for gamma-N interactions by checking for
79 // final-state tables instead of particle "isPhoton()"
80 // 20111007 M. Kelsey -- Add gamma-N final-state tables to printFinalState
81 // 20111107 M. Kelsey -- In sampleCMmomentumFor2to2(), hide message about
82 // unrecognized gamma-N initial state behind verbosity.
83 // 20111216 M. Kelsey -- Add diagnostics to generateMomModulesFor2toMany(),
84 // defer allocation of buffer in generateSCMfinalState() so that
85 // multiplicity failures return zero output, and can be trapped.
86 // 20120308 M. Kelsey -- Allow photons to interact with dibaryons (see
87 // changes in G4NucleiModel).
88 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
89 
90 #include "G4ElementaryParticleCollider.hh"
91 #include "G4CascadeChannel.hh"
92 #include "G4CascadeChannelTables.hh"
93 #include "G4CascadeInterpolator.hh"
94 #include "G4CollisionOutput.hh"
95 #include "G4InuclParticleNames.hh"
96 #include "G4InuclSpecialFunctions.hh"
97 #include "G4LorentzConvertor.hh"
98 #include "G4ParticleLargerEkin.hh"
99 #include "Randomize.hh"
100 #include <algorithm>
101 #include <cfloat>
102 #include <vector>
103 
104 using namespace G4InuclParticleNames;
105 using namespace G4InuclSpecialFunctions;
106 
107 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
108 
109 
110 G4ElementaryParticleCollider::G4ElementaryParticleCollider()
111  : G4CascadeColliderBase("G4ElementaryParticleCollider") {}
112 
113 
114 void
115 G4ElementaryParticleCollider::collide(G4InuclParticle* bullet,
116  G4InuclParticle* target,
117  G4CollisionOutput& output)
118 {
119  if (verboseLevel > 1)
120  G4cout << " >>> G4ElementaryParticleCollider::collide" << G4endl;
121 
122  if (!useEPCollider(bullet,target)) { // Sanity check
123  G4cerr << " ElementaryParticleCollider -> can collide only particle with particle "
124  << G4endl;
125  return;
126  }
127 
128 #ifdef G4CASCADE_DEBUG_SAMPLER
129  static G4bool doPrintTables = true; // Once and only once per job
130  if (doPrintTables) {
131  printFinalStateTables(); // For diagnostic reporting
132  doPrintTables = false;
133  }
134 #endif
135 
136  interCase.set(bullet, target); // To identify kind of collision
137 
138  if (verboseLevel > 1) G4cout << *bullet << G4endl << *target << G4endl;
139 
140  G4InuclElementaryParticle* particle1 =
141  dynamic_cast<G4InuclElementaryParticle*>(bullet);
142  G4InuclElementaryParticle* particle2 =
143  dynamic_cast<G4InuclElementaryParticle*>(target);
144 
145  if (!particle1 || !particle2) { // Redundant with useEPCollider()
146  G4cerr << " ElementaryParticleCollider -> can only collide hadrons"
147  << G4endl;
148  return;
149  }
150 
151  // Check for available interaction, or pion+dibaryon special case
152  if (!G4CascadeChannelTables::GetTable(interCase.hadrons()) &&
153  !particle1->quasi_deutron() && !particle2->quasi_deutron()) {
154  G4cerr << " ElementaryParticleCollider -> cannot collide "
155  << particle1->getDefinition()->GetParticleName() << " with "
156  << particle2->getDefinition()->GetParticleName() << G4endl;
157  return;
158  }
159  // Generate nucleon or pion collision with nucleon
160  // or pion with quasi-deuteron
161 
162  if (particle1->nucleon() || particle2->nucleon()) { // ok
163  G4LorentzConvertor convertToSCM;
164  if(particle2->nucleon()) {
165  convertToSCM.setBullet(particle1);
166  convertToSCM.setTarget(particle2);
167  } else {
168  convertToSCM.setBullet(particle2);
169  convertToSCM.setTarget(particle1);
170  };
171 
172  convertToSCM.setVerbose(verboseLevel);
173 
174  convertToSCM.toTheCenterOfMass();
175  G4double ekin = convertToSCM.getKinEnergyInTheTRS();
176  G4double etot_scm = convertToSCM.getTotalSCMEnergy();
177  G4double pscm = convertToSCM.getSCMMomentum();
178 
179  generateSCMfinalState(ekin, etot_scm, pscm, particle1, particle2,
180  &convertToSCM);
181 
182  if (particles.empty()) { // No final state possible, pass bullet through
183  if (verboseLevel) {
184  G4cerr << " ElementaryParticleCollider -> failed to collide "
185  << particle1->getMomModule() << " GeV/c "
186  << particle1->getDefinition()->GetParticleName() << " with "
187  << particle2->getDefinition()->GetParticleName() << G4endl;
188  }
189  } else { // convert back to Lab
190  G4LorentzVector mom; // Buffer to avoid memory churn
191  particleIterator ipart;
192  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
193  mom = convertToSCM.backToTheLab(ipart->getMomentum());
194  ipart->setMomentum(mom);
195  };
196 
197  // Check conservation in multibody final state
198  if (verboseLevel && !validateOutput(bullet, target, particles)) {
199  G4cout << " incoming particles: \n" << *particle1 << G4endl
200  << *particle2 << G4endl
201  << " outgoing particles: " << G4endl;
202  for(ipart = particles.begin(); ipart != particles.end(); ipart++)
203  G4cout << *ipart << G4endl;
204 
205  G4cout << " <<< Non-conservation in G4ElementaryParticleCollider"
206  << G4endl;
207  }
208 
209  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
210  output.addOutgoingParticles(particles);
211  }
212  } else { // neither particle is nucleon: pion on quasideuteron
213  if (particle1->quasi_deutron() || particle2->quasi_deutron()) {
214  if (particle1->pion() || particle2->pion() ||
215  particle1->isPhoton() || particle2->isPhoton()) {
216  G4LorentzConvertor convertToSCM;
217  if(particle2->quasi_deutron()) { // Quasideuteron is target
218  convertToSCM.setBullet(particle1);
219  convertToSCM.setTarget(particle2);
220  } else {
221  convertToSCM.setBullet(particle2);
222  convertToSCM.setTarget(particle1);
223  };
224  convertToSCM.toTheCenterOfMass();
225  G4double etot_scm = convertToSCM.getTotalSCMEnergy();
226 
227  generateSCMpionAbsorption(etot_scm, particle1, particle2);
228 
229  if (particles.empty()) { // Failed to generate final state
230  if (verboseLevel) {
231  G4cerr << " ElementaryParticleCollider -> failed to collide "
232  << particle1->getMomModule() << " GeV/c "
233  << particle1->getDefinition()->GetParticleName() << " with "
234  << particle2->getDefinition()->GetParticleName() << G4endl;
235  }
236  } else { // convert back to Lab
237  G4LorentzVector mom; // Buffer to avoid memory churn
238  particleIterator ipart;
239  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
240  mom = convertToSCM.backToTheLab(ipart->getMomentum());
241  ipart->setMomentum(mom);
242  };
243 
244  validateOutput(bullet, target, particles); // Check conservation
245 
246  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
247  output.addOutgoingParticles(particles);
248  };
249  } else {
250  G4cerr << " ElementaryParticleCollider -> can only collide pions with dibaryons "
251  << G4endl;
252  };
253  } else {
254  G4cerr << " ElementaryParticleCollider -> can only collide something with nucleon or dibaryon "
255  << G4endl;
256  };
257  };
258 }
259 
260 
261 G4int
262 G4ElementaryParticleCollider::generateMultiplicity(G4int is,
263  G4double ekin) const
264 {
265  G4int mul = 0;
266 
267  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);
268 
269  if (xsecTable) mul = xsecTable->getMultiplicity(ekin);
270  else {
271  G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
272  << is << " - multiplicity not generated " << G4endl;
273  }
274 
275  if(verboseLevel > 3){
276  G4cout << " G4ElementaryParticleCollider::generateMultiplicity: "
277  << " multiplicity = " << mul << G4endl;
278  }
279 
280  return mul;
281 }
282 
283 
284 void
285 G4ElementaryParticleCollider::generateSCMfinalState(G4double ekin,
286  G4double etot_scm,
287  G4double pscm,
288  G4InuclElementaryParticle* particle1,
289  G4InuclElementaryParticle* particle2,
290  G4LorentzConvertor* toSCM) {
291  if (verboseLevel > 3) {
292  G4cout << " >>> G4ElementaryParticleCollider::generateSCMfinalState"
293  << G4endl;
294  }
295 
296  const G4double ang_cut = 0.9999;
297  const G4double difr_const = 0.3678794;
298  const G4int itry_max = 10;
299  G4InuclElementaryParticle dummy;
300 
301  G4int type1 = particle1->type();
302  G4int type2 = particle2->type();
303 
304  G4int is = type1 * type2;
305 
306  if(verboseLevel > 3){
307  G4cout << " is " << is << G4endl;
308  }
309 
310  G4int multiplicity = 0;
311  G4bool generate = true;
312 
313  while (generate) {
314  particles.clear(); // Initialize buffers for this event
315  particle_kinds.clear();
316 
317  // Generate list of final-state particles
318  multiplicity = generateMultiplicity(is, ekin);
319 
320  generateOutgoingPartTypes(is, multiplicity, ekin);
321  if (particle_kinds.empty()) {
322  if (verboseLevel > 3) {
323  G4cout << " generateOutgoingPartTypes failed mult " << multiplicity
324  << G4endl;
325  }
326  continue;
327  }
328 
329  if (multiplicity == 2) {
330  // Identify charge or strangeness exchange (non-elastic scatter)
331  G4int finaltype = particle_kinds[0]*particle_kinds[1];
332  G4int kw = (finaltype != is) ? 2 : 1;
333 
334  G4double pmod = pscm; // Elastic scattering preserves CM momentum
335 
336  if (kw == 2) { // Non-elastic needs new CM momentum value
337  G4double mone = dummy.getParticleMass(particle_kinds[0]);
338  G4double mtwo = dummy.getParticleMass(particle_kinds[1]);
339 
340  if (etot_scm < mone+mtwo) { // Can't produce final state
341  if (verboseLevel > 2) {
342  G4cerr << " bad final state " << particle_kinds[0]
343  << " , " << particle_kinds[1] << " etot_scm " << etot_scm
344  << " < mone+mtwo " << mone+mtwo << " , but ekin " << ekin
345  << G4endl;
346  }
347  continue;
348  }
349 
350  G4double ecm_sq = etot_scm*etot_scm;
351  G4double msumsq = mone+mtwo; msumsq *= msumsq;
352  G4double mdifsq = mone-mtwo; mdifsq *= mdifsq;
353 
354  G4double a = (ecm_sq - msumsq) * (ecm_sq - mdifsq);
355 
356  pmod = std::sqrt(a)/(2.*etot_scm);
357  }
358 
359  G4LorentzVector mom = sampleCMmomentumFor2to2(is, kw, ekin, pmod);
360 
361  if (verboseLevel > 3) {
362  G4cout << " Particle kinds = " << particle_kinds[0] << " , "
363  << particle_kinds[1] << G4endl
364  << " pscm " << pscm << " pmod " << pmod << G4endl
365  << " before rotation px " << mom.x() << " py " << mom.y()
366  << " pz " << mom.z() << G4endl;
367  }
368 
369  mom = toSCM->rotate(mom);
370 
371  if (verboseLevel > 3){
372  G4cout << " after rotation px " << mom.x() << " py " << mom.y() <<
373  " pz " << mom.z() << G4endl;
374  }
375  G4LorentzVector mom1(-mom.vect(), mom.e());
376 
377  particles.resize(multiplicity); // Preallocate buffer
378  particles[0].fill(mom, particle_kinds[0], G4InuclParticle::EPCollider);
379  particles[1].fill(mom1, particle_kinds[1], G4InuclParticle::EPCollider);
380  generate = false;
381  } else { // 2 -> many
382  G4int itry = 0;
383  G4bool bad = true;
384  G4int knd_last = particle_kinds[multiplicity - 1];
385  G4double mass_last = dummy.getParticleMass(knd_last);
386 
387  if (verboseLevel > 3){
388  G4cout << " knd_last " << knd_last << " mass " << mass_last << G4endl;
389  }
390 
391  while (bad && itry < itry_max) {
392  itry++;
393 
394  if (verboseLevel > 3){
395  G4cout << " itry in generateSCMfinalState " << itry << G4endl;
396  }
397 
398  generateMomModules(multiplicity, is, ekin, etot_scm);
399  if (G4int(modules.size()) != multiplicity) {
400  if (verboseLevel > 3) {
401  G4cerr << " generateMomModule failed at mult " << multiplicity
402  << " ekin " << ekin << " etot_scm " << etot_scm << G4endl;
403  }
404  continue;
405  }
406 
407  if (multiplicity == 3) {
408  G4LorentzVector mom3 =
409  particleSCMmomentumFor2to3(is, knd_last, ekin, modules[2]);
410 
411  mom3 = toSCM->rotate(mom3);
412 
413  // generate the momentum of first
414  G4double ct = -0.5 * (modules[2] * modules[2] +
415  modules[0] * modules[0] -
416  modules[1] * modules[1]) /
417  modules[2] / modules[0];
418 
419  if(std::fabs(ct) < ang_cut) {
420 
421  if(verboseLevel > 2){
422  G4cout << " ok for mult " << multiplicity << G4endl;
423  }
424 
425  G4LorentzVector mom1 = generateWithFixedTheta(ct, modules[0]);
426  mom1 = toSCM->rotate(mom3, mom1);
427 
428  // Second particle recoils off 1 & 3
429  G4LorentzVector mom2(etot_scm);
430  mom2 -= mom1+mom3;
431 
432  bad = false;
433  generate = false;
434 
435  particles.resize(multiplicity); // Preallocate buffer
436 
437  particles[0].fill(mom1, particle_kinds[0], G4InuclParticle::EPCollider);
438  particles[1].fill(mom2, particle_kinds[1], G4InuclParticle::EPCollider);
439  particles[2].fill(mom3, particle_kinds[2], G4InuclParticle::EPCollider);
440  };
441  } else { // multiplicity > 3
442  // generate first mult - 2 momentums
443  G4LorentzVector tot_mom;
444  scm_momentums.clear();
445 
446  for (G4int i = 0; i < multiplicity - 2; i++) {
447  G4double p0 = particle_kinds[i] < 3 ? 0.36 : 0.25;
448  G4double alf = 1.0 / p0 / (p0 - (modules[i] + p0) *
449  std::exp(-modules[i] / p0));
450  G4double st = 2.0;
451  G4int itry1 = 0;
452 
453  while (std::fabs(st) > ang_cut && itry1 < itry_max) {
454  itry1++;
455  G4double s1 = modules[i] * inuclRndm();
456  G4double s2 = alf * difr_const * p0 * inuclRndm();
457 
458  if(verboseLevel > 3){
459  G4cout << " s1 * alf * std::exp(-s1 / p0) "
460  << s1 * alf * std::exp(-s1 / p0)
461  << " s2 " << s2 << G4endl;
462  }
463 
464  if(s1 * alf * std::exp(-s1 / p0) > s2) st = s1 / modules[i];
465  };
466 
467  if(verboseLevel > 3){
468  G4cout << " itry1 " << itry1 << " i " << i << " st " << st
469  << G4endl;
470  }
471 
472  if(itry1 == itry_max) {
473  if(verboseLevel > 2){
474  G4cout << " high energy angles generation: itry1 " << itry1
475  << G4endl;
476  }
477 
478  st = 0.5 * inuclRndm();
479  };
480 
481  G4double ct = std::sqrt(1.0 - st * st);
482  if (inuclRndm() > 0.5) ct = -ct;
483 
484  G4LorentzVector mom = generateWithFixedTheta(ct,modules[i]);
485 
486  tot_mom += mom;
487 
488  scm_momentums.push_back(mom);
489  };
490 
491  // handle last two
492  G4double tot_mod = tot_mom.rho();
493  G4double ct = -0.5 * (tot_mod * tot_mod +
494  modules[multiplicity - 2] * modules[multiplicity - 2] -
495  modules[multiplicity - 1] * modules[multiplicity - 1]) / tot_mod /
496  modules[multiplicity - 2];
497 
498  if (verboseLevel > 2){
499  G4cout << " ct last " << ct << G4endl;
500  }
501 
502  if (std::fabs(ct) < ang_cut) {
503  G4int i(0);
504  for (i = 0; i < multiplicity - 2; i++)
505  scm_momentums[i] = toSCM->rotate(scm_momentums[i]);
506 
507  tot_mom = toSCM->rotate(tot_mom);
508 
509  G4LorentzVector mom =
510  generateWithFixedTheta(ct, modules[multiplicity - 2]);
511 
512  mom = toSCM->rotate(tot_mom, mom);
513  scm_momentums.push_back(mom);
514 
515  // Last particle recoils off everything else
516  G4LorentzVector mom1(etot_scm);
517  mom1 -= mom+tot_mom;
518 
519  scm_momentums.push_back(mom1);
520  bad = false;
521  generate = false;
522 
523  if (verboseLevel > 2){
524  G4cout << " ok for mult " << multiplicity << G4endl;
525  }
526 
527  particles.resize(multiplicity); // Preallocate buffer
528  for (i = 0; i < multiplicity; i++) {
529  particles[i].fill(scm_momentums[i], particle_kinds[i],
530  G4InuclParticle::EPCollider);
531  }
532  } // |ct| < ang_cut
533  } // multiplicity > 3
534  } // while (bad && itry < itry_max)
535 
536  if (itry == itry_max) {
537  if (verboseLevel > 2) {
538  G4cout << " cannot generate the distr. for mult " << multiplicity
539  << G4endl;
540  }
541  break;
542  }
543  } // multiplicity > 2
544  } // while (generate)
545 
546  if (verboseLevel > 3) {
547  G4cout << " <<< G4ElementaryParticleCollider::generateSCMfinalState"
548  << G4endl;
549  }
550 
551  return; // Particles buffer filled
552 }
553 
554 void
555 G4ElementaryParticleCollider::generateMomModules(G4int mult,
556  G4int is,
557  G4double ekin,
558  G4double etot_cm) {
559  if (verboseLevel > 3) {
560  G4cout << " >>> G4ElementaryParticleCollider::generateMomModules"
561  << G4endl;
562  }
563 
564  if (verboseLevel > 2){
565  G4cout << " mult " << mult << " is " << is << " ekin " << ekin
566  << " etot_cm " << etot_cm << G4endl;
567  }
568 
569  const G4int itry_max = 10;
570  const G4double small = 1.e-10;
571  G4InuclElementaryParticle dummy;
572  G4int itry = 0;
573 
574  modules.clear(); // Initialize buffer for this attempt
575  modules.resize(mult,0.);
576 
577  masses2.clear();
578  masses2.resize(mult,0.); // Allows direct [i] setting
579 
580  for (G4int i = 0; i < mult; i++) {
581  G4double mass = dummy.getParticleMass(particle_kinds[i]);
582  masses2[i] = mass * mass;
583  };
584 
585  G4double mass_last = std::sqrt(masses2[mult - 1]);
586 
587  if (verboseLevel > 3){
588  G4cout << " knd_last " << particle_kinds[mult - 1] << " mass_last "
589  << mass_last << G4endl;
590  }
591 
592  while (itry < itry_max) {
593  itry++;
594  if(verboseLevel > 3){
595  G4cout << " itry in generateMomModules " << itry << G4endl;
596  }
597 
598  G4int ilast = -1;
599  G4double eleft = etot_cm;
600 
601  for (G4int i = 0; i < mult - 1; i++) {
602  G4double pmod =
603  getMomModuleFor2toMany(is, mult, particle_kinds[i], ekin);
604 
605  if (pmod < small) break;
606  eleft -= std::sqrt(pmod * pmod + masses2[i]);
607 
608  if (verboseLevel > 3){
609  G4cout << " kp " << particle_kinds[i] << " pmod " << pmod
610  << " mass2 " << masses2[i] << " eleft " << eleft
611  << "\n x1 " << eleft - mass_last << G4endl;
612  }
613 
614  if (eleft <= mass_last) break;
615  ilast++;
616  modules[i] = pmod;
617  };
618 
619  if (ilast == mult - 2) {
620  G4double plast = eleft * eleft - masses2[mult - 1];
621  if (verboseLevel > 2){
622  G4cout << " plast ** 2 " << plast << G4endl;
623  }
624 
625  if (plast > small) {
626  plast = std::sqrt(plast);
627  modules[mult - 1] = plast;
628 
629  if (mult == 3) {
630  if (satisfyTriangle(modules)) return;
631  } else return;
632  }
633  }
634  } // while (itry < itry_max)
635 
636  if (verboseLevel > 2)
637  G4cerr << " Unable to generate momenta for multiplicity " << mult << G4endl;
638 
639  modules.clear(); // Something went wrong, throw away partial
640  return;
641 }
642 
643 
644 G4bool
645 G4ElementaryParticleCollider::satisfyTriangle(
646  const std::vector<G4double>& pmod) const
647 {
648  if (verboseLevel > 3) {
649  G4cout << " >>> G4ElementaryParticleCollider::satisfyTriangle"
650  << G4endl;
651  }
652 
653  G4bool good = ( (pmod.size() != 3) ||
654  !(std::fabs(pmod[1] - pmod[2]) > pmod[0] ||
655  pmod[0] > pmod[1] + pmod[2] ||
656  std::fabs(pmod[0] - pmod[2]) > pmod[1] ||
657  pmod[1] > pmod[0] + pmod[2] ||
658  std::fabs(pmod[0] - pmod[1]) > pmod[2] ||
659  pmod[2] > pmod[1] + pmod[0]));
660 
661  return good;
662 }
663 
664 
665 void
666 G4ElementaryParticleCollider::generateOutgoingPartTypes(G4int is, G4int mult,
667  G4double ekin)
668 {
669  particle_kinds.clear(); // Initialize buffer for generation
670 
671  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);
672 
673  if (xsecTable)
674  xsecTable->getOutgoingParticleTypes(particle_kinds, mult, ekin);
675  else {
676  G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
677  << is << " - outgoing kinds not generated " << G4endl;
678  }
679 
680  return;
681 }
682 
683 
684 G4double
685 G4ElementaryParticleCollider::getMomModuleFor2toMany(G4int is, G4int /*mult*/,
686  G4int knd,
687  G4double ekin) const
688 {
689  if (verboseLevel > 2) {
690  G4cout << " >>> G4ElementaryParticleCollider::getMomModuleFor2toMany "
691  << " is " << is << " knd " << knd << " ekin " << ekin << G4endl;
692  }
693 
694  G4double S = inuclRndm();
695  G4double PS = 0.0;
696  G4double PR = 0.0;
697  G4double PQ = 0.0;
698  G4int KM = 2;
699  G4int IL = 4;
700  G4int JK = 4;
701  G4int JM = 2;
702  G4int IM = 3;
703 
704  if (is == 1 || is == 2 || is == 4) KM = 1;
705  if (knd == 1 || knd == 2) JK = 0;
706 
707  if (verboseLevel > 3) {
708  G4cout << " S " << S << " KM " << KM << " IL " << IL << " JK " << JK
709  << " JM " << JM << " IM " << IM << G4endl;
710  }
711 
712  for(G4int i = 0; i < 4; i++) {
713  G4double V = 0.0;
714  for(G4int k = 0; k < 4; k++) {
715  if (verboseLevel > 3) {
716  G4cout << " k " << k << " : rmn[k+JK][i+IL][KM-1] "
717  << rmn[k+JK][i+IL][KM-1] << " ekin^k " << std::pow(ekin, k)
718  << G4endl;
719  }
720 
721  V += rmn[k + JK][i + IL][KM - 1] * std::pow(ekin, k);
722  }
723 
724  if (verboseLevel > 3) {
725  G4cout << " i " << i << " : V " << V << " S^i " << std::pow(S, i)
726  << G4endl;
727  }
728 
729  PR += V * std::pow(S, i);
730  PQ += V;
731  }
732 
733  if (verboseLevel > 3) G4cout << " PR " << PR << " PQ " << PQ << G4endl;
734 
735  if (knd == 1 || knd == 2) JM = 1;
736  if (verboseLevel > 3) G4cout << " JM " << JM << G4endl;
737 
738  for(G4int im = 0; im < 3; im++) {
739  if (verboseLevel >3) {
740  G4cout << " im " << im << " : rmn[8+IM+im][7+JM][KM-1] "
741  << rmn[8+IM+im][7+JM][KM-1] << " ekin^im " << std::pow(ekin, im)
742  << G4endl;
743  }
744  PS += rmn[8 + IM + im][7 + JM][KM - 1] * std::pow(ekin, im);
745  }
746 
747  G4double PRA = PS * std::sqrt(S) * (PR + (1 - PQ) * (S*S*S*S));
748 
749  if (verboseLevel > 3)
750  G4cout << " PS " << PS << " PRA = PS*sqrt(S)*(PR+(1-PQ)*S^4) " << PRA
751  << G4endl;
752 
753  return std::fabs(PRA);
754 }
755 
756 
757 G4LorentzVector
758 G4ElementaryParticleCollider::particleSCMmomentumFor2to3(
759  G4int is,
760  G4int knd,
761  G4double ekin,
762  G4double pmod) const
763 {
764  if (verboseLevel > 3) {
765  G4cout << " >>> G4ElementaryParticleCollider::particleSCMmomentumFor2to3"
766  << G4endl;
767  }
768 
769  const G4int itry_max = 100;
770  G4double ct = 2.0;
771  G4int K = 3;
772  G4int J = 1;
773 
774  if(is == 1 || is == 2 || is == 4) K = 1;
775 
776  if(knd == 1 || knd == 2) J = 0;
777 
778  G4int itry = 0;
779 
780  while(std::fabs(ct) > 1.0 && itry < itry_max) {
781  itry++;
782  G4double S = inuclRndm();
783  G4double U = 0.0;
784  G4double W = 0.0;
785 
786  for(G4int l = 0; l < 4; l++) {
787  G4double V = 0.0;
788 
789  for(G4int im = 0; im < 4; im++) {
790  V += abn[im][l][K+J-1] * std::pow(ekin, im);
791  };
792 
793  U += V;
794  W += V * std::pow(S, l);
795  };
796  ct = 2.0 * std::sqrt(S) * (W + (1.0 - U) * (S*S*S*S)) - 1.0;
797  };
798 
799  if(itry == itry_max) {
800 
801  if(verboseLevel > 2){
802  G4cout << " particleSCMmomentumFor2to3 -> itry = itry_max " << itry << G4endl;
803  }
804 
805  ct = 2.0 * inuclRndm() - 1.0;
806  };
807 
808  return generateWithFixedTheta(ct, pmod);
809 }
810 
811 
812 G4LorentzVector
813 G4ElementaryParticleCollider::sampleCMmomentumFor2to2(G4int is, G4int kw,
814  G4double ekin,
815  G4double pscm) const
816 {
817  if (verboseLevel > 3)
818  G4cout << " >>> G4ElementaryParticleCollider::sampleCMmomentumFor2to2"
819  << " is " << is << " kw " << kw << " ekin " << ekin << " pscm "
820  << pscm << G4endl;
821 
822  G4double pA=0.0, pC=0.0, pCos=0.0, pFrac=0.0; // Angular parameters
823 
824  // Arrays below are parameters for two-exponential sampling of angular
825  // distributions of two-body scattering in the specified channels
826 
827  if (is == 1 || is == 2 || is == 4 ||
828  is == 21 || is == 23 || is == 25 || is == 27 || is ==29 || is == 31 ||
829  is == 42 || is == 46 || is == 50 || is == 54 || is ==58 || is == 62) {
830  // nucleon-nucleon or hyperon-nucleon
831  if (verboseLevel > 3) G4cout << " nucleon/hyperon elastic" << G4endl;
832 
833  static const G4double nnke[9] = { 0.0, 0.44, 0.59, 0.80, 1.00, 2.24, 4.40, 6.15, 10.00};
834  static const G4double nnA[9] = { 0.0, 0.34, 2.51, 4.59, 4.2, 5.61, 6.38, 7.93, 8.7};
835  static const G4double nnC[9] = { 0.0, 0.0, 1.21, 1.54, 1.88, 1.24, 1.91, 4.04, 8.7};
836  static const G4double nnCos[9] = {-1.0, -1.0, 0.4226, 0.4226, 0.4384, 0.7193, 0.8788, 0.9164, 0.95};
837  static const G4double nnFrac[9] = {1.0, 1.0, 0.4898, 0.7243, 0.7990, 0.8892, 0.8493, 0.9583, 1.0};
838 
839  static G4CascadeInterpolator<9> interp(nnke); // Only need one!
840  pA = interp.interpolate(ekin, nnA);
841  pC = interp.interpolate(ekin, nnC);
842  pCos = interp.interpolate(ekin, nnCos);
843  pFrac = interp.interpolate(ekin, nnFrac);
844 
845  } else if (kw == 2 && (is == 9 || is == 18)) {
846  // gamma p -> pi+ n, gamma p -> pi0 p, gamma p -> K Y (and isospin variants)
847  // for now and due to lack of better data, use the gamma p -> pi+ n angular
848  // distribution for all of these channels
849  if (verboseLevel > 3)
850  G4cout << " gamma-nucleon inelastic with 2-body final state" << G4endl;
851 
852  static const G4double gnke[10] = {0.0, 0.11, 0.22, 0.26, 0.30, 0.34, 0.42, 0.59, 0.79, 10.0};
853  static const G4double gnA[10] = {0.0, 0.0, 5.16, 5.55, 5.33, 7.40, 13.55, 13.44, 13.31, 7.3};
854  static const G4double gnC[10] = {0.0, -10.33, -5.44, -5.92, -4.27, -0.66, 1.37, 1.07, 0.52, 7.3};
855  static const G4double gnCos[10] = {1.0, 1.0, 0.906, 0.940, 0.940, 0.906, 0.906, 0.91, 0.91, 0.94};
856  static const G4double gnFrac[10] = {0.0, 0.0, 0.028, 0.012, 0.014, 0.044, 0.087, 0.122, 0.16, 1.0};
857 
858  static G4CascadeInterpolator<10> interp(gnke);
859  pA = interp.interpolate(ekin, gnA);
860  pC = interp.interpolate(ekin, gnC);
861  pCos = interp.interpolate(ekin, gnCos);
862  pFrac = interp.interpolate(ekin, gnFrac);
863 
864  } else if (kw == 2) {
865  // pi- p -> pi0 n, pi0 p -> pi+ n, pi- p -> K Y, pi0 p -> K Y (and isospin variants)
866  // includes charge and strangeness exchange
867  if (verboseLevel > 3)
868  G4cout << " pion-nucleon inelastic with 2-body final state " << G4endl;
869 
870  static const G4double qxke[10] = {0.0, 0.062, 0.12, 0.217, 0.533, 0.873, 1.34, 2.86, 5.86, 10.0};
871  static const G4double qxA[10] = {0.0, 0.0, 2.48, 7.93, 10.0, 9.78, 5.08, 8.13, 8.13, 8.13};
872  static const G4double qxC[10] = {0.0, -39.58, -12.55, -4.38, 1.81, -1.99, -0.33, 1.2, 1.43, 8.13};
873  static const G4double qxCos[10] = {1.0, 1.0, 0.604, -0.033, 0.25, 0.55, 0.65, 0.80, 0.916, 0.916};
874  static const G4double qxFrac[10] = {0.0, 0.0, 0.1156, 0.5832, 0.8125, 0.3357, 0.3269, 0.7765, 0.8633, 1.0};
875 
876  static G4CascadeInterpolator<10> interp(qxke); // Only need one!
877  pA = interp.interpolate(ekin, qxA);
878  pC = interp.interpolate(ekin, qxC);
879  pCos = interp.interpolate(ekin, qxCos);
880  pFrac = interp.interpolate(ekin, qxFrac);
881 
882  } else if (is == 3 || is == 7 || is == 9 || is == 11 || is == 17 ||
883  is == 10 || is == 14 || is == 18 || is == 26 || is == 30) {
884  // pi+p, pi0p, gammap, k+p, k0bp, pi-n, pi0n, gamman, k-n, or k0n
885  if (verboseLevel > 3) G4cout << " meson-nucleon elastic (1)" << G4endl;
886 
887  static const G4double hn1ke[10] = {0.0, 0.062, 0.12, 0.217, 0.533, 0.873, 1.34, 2.86, 5.86, 10.0};
888  static const G4double hn1A[10] = {0.0, 0.0, 27.58, 19.83, 6.46, 4.59, 6.47, 6.68, 6.43, 6.7};
889  static const G4double hn1C[10] = {0.0, -26.4, -30.55, -19.42, -5.05, -5.24, -1.00, 2.14, 2.9, 6.7};
890  static const G4double hn1Cos[10] = {1.0, 1.0, 0.174, -0.174, -0.7, -0.295, 0.5, 0.732, 0.837, 0.89};
891  static const G4double hn1Frac[10] = {0.0, 0.0, 0.2980, 0.7196, 0.9812, 0.8363, 0.5602, 0.9601, 0.9901, 1.0};
892 
893  static G4CascadeInterpolator<10> interp(hn1ke); // Only need one!
894  pA = interp.interpolate(ekin, hn1A);
895  pC = interp.interpolate(ekin, hn1C);
896  pCos = interp.interpolate(ekin, hn1Cos);
897  pFrac = interp.interpolate(ekin, hn1Frac);
898 
899  } else if (is == 5 || is == 6 || is == 13 || is == 34 || is == 22 ||
900  is == 15) {
901  // pi-p, pi+n, k-p, k0bn, k+n, or k0p
902  if (verboseLevel > 3) G4cout << " meson-nucleon elastic (2)" << G4endl;
903 
904  static const G4double hn2ke[10] = {0.0, 0.062, 0.12, 0.217, 0.533, 0.873, 1.34, 2.86, 5.86, 10.0};
905  static const G4double hn2A[10] = {0.0, 27.08, 19.32, 9.92, 7.74, 9.86, 5.51, 7.25, 7.23, 7.3};
906  static const G4double hn2C[10] = {0.0, 0.0, -19.49, -15.78, -9.78, -2.74, -1.16, 2.31, 2.96, 7.3};
907  static const G4double hn2Cos[10] = {-1.0, -1.0, -0.235, -0.259, -0.276, 0.336, 0.250, 0.732, 0.875, 0.9};
908  static const G4double hn2Frac[10] = {1.0, 1.0, 0.6918, 0.6419, 0.7821, 0.6542, 0.8382, 0.9722, 0.9784, 1.0};
909 
910  static G4CascadeInterpolator<10> interp(hn2ke); // Only need one!
911  pA = interp.interpolate(ekin, hn2A);
912  pC = interp.interpolate(ekin, hn2C);
913  pCos = interp.interpolate(ekin, hn2Cos);
914  pFrac = interp.interpolate(ekin, hn2Frac);
915 
916  } else {
917  if (verboseLevel)
918  G4cerr << " G4ElementaryParticleCollider::sampleCMmomentumFor2to2:"
919  << " interaction is=" << is << " not recognized " << G4endl;
920  }
921 
922  // Bound parameters by their physical ranges
923  pCos = std::max(-1.,std::min(pCos,1.));
924  pFrac = std::max(0.,std::min(pFrac,1.));
925 
926  // Use parameters determined above to get polar angle
927  G4double ct = sampleCMcosFor2to2(pscm, pFrac, pA, pC, pCos);
928 
929  return generateWithFixedTheta(ct, pscm);
930 }
931 
932 
933 G4double
934 G4ElementaryParticleCollider::sampleCMcosFor2to2(G4double pscm, G4double pFrac,
935  G4double pA, G4double pC,
936  G4double pCos) const
937 {
938  if (verboseLevel>3) {
939  G4cout << " sampleCMcosFor2to2: pscm " << pscm << " pFrac " << pFrac
940  << " pA " << pA << " pC " << pC << " pCos " << pCos << G4endl;
941  }
942 
943  G4bool smallAngle = (G4UniformRand() < pFrac); // 0 < theta < theta0
944 
945  G4double term1 = 2.0 * pscm*pscm * (smallAngle ? pA : pC);
946 
947  if (std::abs(term1) < 1e-7) return 1.0; // No actual scattering here!
948  if (term1 > DBL_MAX_EXP) return 1.0;
949 
950  G4double term2 = std::exp(-2.0*term1);
951  G4double randScale = (std::exp(-term1*(1.0 - pCos)) - term2)/(1.0 - term2);
952 
953  G4double randVal;
954  if (smallAngle) randVal = (1.0 - randScale)*G4UniformRand() + randScale;
955  else randVal = randScale*G4UniformRand();
956 
957  G4double costheta = 1.0 + std::log(randVal*(1.0 - term2) + term2)/term1;
958 
959  if (verboseLevel>3) {
960  G4cout << " term1 " << term1 << " term2 " << term2 << " randVal "
961  << randVal << " => costheta " << costheta << G4endl;
962  }
963 
964  return costheta;
965 }
966 
967 
968 void
969 G4ElementaryParticleCollider::generateSCMpionAbsorption(G4double etot_scm,
970  G4InuclElementaryParticle* particle1,
971  G4InuclElementaryParticle* particle2) {
972  if (verboseLevel > 3)
973  G4cout << " >>> G4ElementaryParticleCollider::generateSCMpionAbsorption"
974  << G4endl;
975 
976  // generate nucleons momenta for pion or photon absorption
977  // the nucleon distribution assumed to be isotropic in SCM
978 
979  particles.clear(); // Initialize buffers for this event
980  particles.resize(2);
981 
982  particle_kinds.clear();
983 
984  G4int type1 = particle1->type();
985  G4int type2 = particle2->type();
986 
987  // generate kinds
988 
989  if (type1 == pionPlus) {
990  if (type2 == diproton) { // pi+ + PP -> ?
991  G4cerr << " pion absorption: pi+ + PP -> ? " << G4endl;
992  return;
993  } else if (type2 == unboundPN) { // pi+ + PN -> PP
994  particle_kinds.push_back(proton);
995  particle_kinds.push_back(proton);
996  } else if (type2 == dineutron) { // pi+ + NN -> PN
997  particle_kinds.push_back(proton);
998  particle_kinds.push_back(neutron);
999  }
1000  } else if (type1 == pionMinus) {
1001  if (type2 == diproton) { // pi- + PP -> PN
1002  particle_kinds.push_back(proton);
1003  particle_kinds.push_back(neutron);
1004  } else if (type2 == unboundPN) { // pi- + PN -> NN
1005  particle_kinds.push_back(neutron);
1006  particle_kinds.push_back(neutron);
1007  } else if (type2 == dineutron) { // pi- + NN -> ?
1008  G4cerr << " pion absorption: pi- + NN -> ? " << G4endl;
1009  return;
1010  }
1011  } else if (type1 == pionZero || type1 == photon) {
1012  if (type2 == diproton) { // pi0/gamma + PP -> PP
1013  particle_kinds.push_back(proton);
1014  particle_kinds.push_back(proton);
1015  } else if (type2 == unboundPN) { // pi0/gamma + PN -> PN
1016  particle_kinds.push_back(proton);
1017  particle_kinds.push_back(neutron);
1018  } else if (type2 == dineutron) { // pi0/gamma + NN -> ?
1019  particle_kinds.push_back(neutron);
1020  particle_kinds.push_back(neutron);
1021  }
1022  }
1023 
1024  G4InuclElementaryParticle dummy;
1025 
1026  G4double mone = dummy.getParticleMass(particle_kinds[0]);
1027  G4double m1sq = mone*mone;
1028 
1029  G4double mtwo = dummy.getParticleMass(particle_kinds[1]);
1030  G4double m2sq = mtwo*mtwo;
1031 
1032  G4double a = 0.5 * (etot_scm * etot_scm - m1sq - m2sq);
1033 
1034  G4double pmod = std::sqrt((a * a - m1sq * m2sq) / (m1sq + m2sq + 2.0 * a));
1035  G4LorentzVector mom1 = generateWithRandomAngles(pmod, mone);
1036  G4LorentzVector mom2;
1037  mom2.setVectM(-mom1.vect(), mtwo);
1038 
1039  particles[0].fill(mom1, particle_kinds[0], G4InuclParticle::EPCollider);
1040  particles[1].fill(mom2, particle_kinds[1], G4InuclParticle::EPCollider);
1041 
1042  return;
1043 }
1044 
1045 
1046 // Dump lookup tables for N-body final states
1047 
1048 void G4ElementaryParticleCollider::
1049 printFinalStateTables(std::ostream& os) const {
1050  G4CascadeChannelTables::PrintTable(pro*pro, os);
1051  G4CascadeChannelTables::PrintTable(neu*pro, os);
1052  G4CascadeChannelTables::PrintTable(neu*neu, os);
1053  G4CascadeChannelTables::PrintTable(kmi*neu, os);
1054  G4CascadeChannelTables::PrintTable(kmi*pro, os);
1055  G4CascadeChannelTables::PrintTable(kpl*neu, os);
1056  G4CascadeChannelTables::PrintTable(kpl*pro, os);
1057  G4CascadeChannelTables::PrintTable(k0b*neu, os);
1058  G4CascadeChannelTables::PrintTable(k0b*pro, os);
1059  G4CascadeChannelTables::PrintTable(k0*neu, os);
1060  G4CascadeChannelTables::PrintTable(k0*pro, os);
1061  G4CascadeChannelTables::PrintTable(lam*neu, os);
1062  G4CascadeChannelTables::PrintTable(lam*pro, os);
1063  G4CascadeChannelTables::PrintTable(pim*neu, os);
1064  G4CascadeChannelTables::PrintTable(pim*pro, os);
1065  G4CascadeChannelTables::PrintTable(pip*neu, os);
1066  G4CascadeChannelTables::PrintTable(pip*pro, os);
1067  G4CascadeChannelTables::PrintTable(pi0*neu, os);
1068  G4CascadeChannelTables::PrintTable(pi0*pro, os);
1069  G4CascadeChannelTables::PrintTable(sm*neu, os);
1070  G4CascadeChannelTables::PrintTable(sm*pro, os);
1071  G4CascadeChannelTables::PrintTable(sp*neu, os);
1072  G4CascadeChannelTables::PrintTable(sp*pro, os);
1073  G4CascadeChannelTables::PrintTable(s0*neu, os);
1074  G4CascadeChannelTables::PrintTable(s0*pro, os);
1075  G4CascadeChannelTables::PrintTable(xim*neu, os);
1076  G4CascadeChannelTables::PrintTable(xim*pro, os);
1077  G4CascadeChannelTables::PrintTable(xi0*neu, os);
1078  G4CascadeChannelTables::PrintTable(xi0*pro, os);
1079 
1080  os << " * * * PRELIMINARY -- GAMMA-NUCLEON TABLES * * *" << G4endl;
1081  G4CascadeChannelTables::PrintTable(gam*neu, os);
1082  G4CascadeChannelTables::PrintTable(gam*pro, os);
1083 }
1084 
1085 
1086 // Parameter array for momentum calculation of many body final states
1087 const G4double G4ElementaryParticleCollider::rmn[14][10][2] = {
1088  {{0.5028, 0.6305}, {3.1442, -3.7333}, {-7.8172, 13.464}, {8.1667, -18.594},
1089  {1.6208, 1.9439}, {-4.3139, -4.6268}, {12.291, 9.7879}, {-15.288, -9.6074},
1090  { 0.0, 0.0}, { 0.0, 0.0}},
1091 
1092  {{0.9348, 2.1801}, {-10.59, 1.5163}, { 29.227, -16.38}, {-34.55, 27.944},
1093  {-0.2009, -0.3464}, {1.3641, 1.1093}, {-3.403, -1.9313}, { 3.8559, 1.7064},
1094  { 0.0, 0.0}, { 0.0, 0.0}},
1095 
1096  {{-0.0967, -1.2886}, {4.7335, -2.457}, {-14.298, 15.129}, {17.685, -23.295},
1097  { 0.0126, 0.0271}, {-0.0835, -0.1164}, { 0.186, 0.2697}, {-0.2004, -0.3185},
1098  { 0.0, 0.0}, { 0.0, 0.0}},
1099 
1100  {{-0.025, 0.2091}, {-0.6248, 0.5228}, { 2.0282, -2.8687}, {-2.5895, 4.2688},
1101  {-0.0002, -0.0007}, {0.0014, 0.0051}, {-0.0024, -0.015}, { 0.0022, 0.0196},
1102  { 0.0, 0.0}, { 0.0, 0.0}},
1103 
1104  {{1.1965, 0.9336}, {-0.8289,-1.8181}, { 1.0426, 5.5157}, { -1.909,-8.5216},
1105  { 1.2419, 1.8693}, {-4.3633, -5.5678}, { 13.743, 14.795}, {-18.592, -16.903},
1106  { 0.0, 0.0}, { 0.0, 0.0}},
1107 
1108  {{0.287, 1.7811}, {-4.9065,-8.2927}, { 16.264, 20.607}, {-19.904,-20.827},
1109  {-0.244, -0.4996}, {1.3158, 1.7874}, {-3.5691, -4.133}, { 4.3867, 3.8393},
1110  { 0.0, 0.0}, { 0.0, 0.0}},
1111 
1112  {{-0.2449, -1.5264}, { 2.9191, 6.8433}, {-9.5776, -16.067}, { 11.938, 16.845},
1113  {0.0157, 0.0462}, {-0.0826, -0.1854}, { 0.2143, 0.4531}, {-0.2585, -0.4627},
1114  { 0.0, 0.0}, { 0.0, 0.0}},
1115 
1116  {{0.0373, 0.2713}, {-0.422, -1.1944}, { 1.3883, 2.7495}, {-1.7476,-2.9045},
1117  {-0.0003, -0.0013}, {0.0014, 0.0058}, {-0.0034,-0.0146}, { 0.0039, 0.0156},
1118  { 0.0, 0.0}, { 0.0, 0.0}},
1119 
1120  {{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1121  { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1122  { 0.1451, 0.0929},{ 0.1538, 0.1303}},
1123 
1124  {{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1125  { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1126  { 0.4652, 0.5389},{ 0.2744, 0.4071}},
1127 
1128  {{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1129  { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1130  { -0.033, -0.0545},{-0.0146, -0.0288}},
1131 
1132  {{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1133  { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1134  { 0.6296, 0.1491},{ 0.8381, 0.1802}},
1135 
1136  {{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1137  { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1138  { 0.1787, 0.385},{ 0.0086, 0.3302}},
1139 
1140  {{ 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1141  { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0}, { 0.0, 0.0},
1142  {-0.0026, -0.0128},{ 0.0033, -0.0094}}
1143 };
1144 
1145 const G4double G4ElementaryParticleCollider::abn[4][4][4] = {
1146  {{0.0856, 0.0716, 0.1729, 0.0376}, {5.0390, 3.0960, 7.1080, 1.4331},
1147  {-13.782, -11.125, -17.961, -3.1350}, {14.661, 18.130, 16.403, 6.4864}},
1148  {{0.0543, 0.0926, -0.1450, 0.2383}, {-9.2324, -3.2186, -13.032, 1.8253},
1149  {36.397, 20.273, 41.781, 1.7648}, {-42.962, -33.245, -40.799, -16.735}},
1150  {{-0.0511, -0.0515, 0.0454, -0.1541}, {4.6003, 0.8989, 8.3515, -1.5201},
1151  {-20.534, -7.5084, -30.260, -1.5692}, {27.731, 13.188, 32.882, 17.185}},
1152  {{0.0075, 0.0058, -0.0048, 0.0250}, {-0.6253, -0.0017, -1.4095, 0.3059},
1153  {2.9159, 0.7022, 5.3505, 0.3252}, {-4.1101, -1.4854, -6.0946, -3.5277}}
1154 };
Generated on Sat May 25 2013 14:33:44 for Geant4 by   doxygen 1.8.4