G4ElementaryParticleCollider.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// $Id$
//
// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
// 20100316  D. Wright (restored by M. Kelsey) -- Replace original (incorrect)
//              pp, pn, nn 2-body to 2-body scattering angular distributions
//              with a new parametrization of elastic scattering data using
//              the sum of two exponentials.
// 20100319  M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
//              eliminate some unnecessary std::pow()
// 20100407  M. Kelsey -- Replace std::vector<>::resize(0) with ::clear()
//              Eliminate return-by-value std::vector<> by creating buffers
//              buffers for particles, momenta, and particle types.
//              The following functions now return void and are non-const:
//                ::generateSCMfinalState()
//                ::generateMomModules()
//                ::generateStrangeChannelPartTypes()
//                ::generateSCMpionAbsorption()
// 20100408  M. Kelsey -- Follow changes to G4*Sampler to pass particle_kinds
//              as input buffer.
// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
// 20100428  M. Kelsey -- Use G4InuclParticleNames enum
// 20100429  M. Kelsey -- Change "photon()" to "isPhoton()"
// 20100507  M. Kelsey -- Rationalize multiplicity returns to be actual value
// 20100511  M. Kelsey -- Replace G4PionSampler and G4NucleonSampler with new
//              pi-N and N-N classes, reorganize if-cascades 
// 20100511  M. Kelsey -- Eliminate three residual random-angles blocks.
// 20100511  M. Kelsey -- Bug fix: pi-N two-body final states not correctly
//              tested for charge-exchange case.
// 20100512  M. Kelsey -- Rationalize multiplicity returns to be actual value
// 20100512  M. Kelsey -- Add some additional debugging messages for 2-to-2
// 20100512  M. Kelsey -- Replace "if (is==)" cascades with switch blocks.
//              Use G4CascadeInterpolator for angular distributions.
// 20100517  M. Kelsey -- Inherit from common base class, make arrays static
// 20100519  M. Kelsey -- Use G4InteractionCase to compute "is" values.
// 20100625  M. Kelsey -- Two bugs in n-body momentum, last particle recoil
// 20100713  M. Kelsey -- Bump collide start message up to verbose > 1
// 20100714  M. Kelsey -- Move conservation checking to base class
// 20100714  M. Kelsey -- Add sanity check for two-body final state, to ensure
//              that final state total mass is below etot_scm; also compute
//              kinematics without "rescaling" (which led to non-conservation)
// 20100726  M. Kelsey -- Move remaining std::vector<> buffers to .hh file
// 20100804  M. Kelsey -- Add printing of final-state tables, protected by
//              G4CASCADE_DEBUG_SAMPLER preprocessor flag
// 20101019  M. Kelsey -- CoVerity report: check dynamic_cast<> for null
// 20110214  M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
// 20110718  V. Uzhinskiy -- Drop IL,IM reset for multiplicity-three collisions
// 20110718  M. Kelsey -- Use enum names in switch blocks (c.f. G4NucleiModel)
// 20110720  M. Kelsey -- Follow interface change for cross-section tables,
//              eliminating switch blocks.
// 20110806  M. Kelsey -- Pre-allocate buffers to avoid memory churn
// 20110922  M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
// 20110926  M. Kelsey -- Protect sampleCMcosFor2to2 from unphysical arguments
// 20111003  M. Kelsey -- Prepare for gamma-N interactions by checking for
//              final-state tables instead of particle "isPhoton()"
// 20111007  M. Kelsey -- Add gamma-N final-state tables to printFinalState
// 20111107  M. Kelsey -- In sampleCMmomentumFor2to2(), hide message about
//              unrecognized gamma-N initial state behind verbosity.
// 20111216  M. Kelsey -- Add diagnostics to generateMomModulesFor2toMany(),
//              defer allocation of buffer in generateSCMfinalState() so that
//              multiplicity failures return zero output, and can be trapped.
// 20120308  M. Kelsey -- Allow photons to interact with dibaryons (see
//              changes in G4NucleiModel).
// 20120608  M. Kelsey -- Fix variable-name "shadowing" compiler warnings.

#include "G4ElementaryParticleCollider.hh"
#include "G4CascadeChannel.hh"
#include "G4CascadeChannelTables.hh"
#include "G4CascadeInterpolator.hh"
#include "G4CollisionOutput.hh"
#include "G4InuclParticleNames.hh"
#include "G4InuclSpecialFunctions.hh"
#include "G4LorentzConvertor.hh"
#include "G4ParticleLargerEkin.hh"
#include "Randomize.hh"
#include <algorithm>
#include <cfloat>
#include <vector>

using namespace G4InuclParticleNames;
using namespace G4InuclSpecialFunctions;

typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;


G4ElementaryParticleCollider::G4ElementaryParticleCollider()
  : G4CascadeColliderBase("G4ElementaryParticleCollider") {}


void
G4ElementaryParticleCollider::collide(G4InuclParticle* bullet,
                                      G4InuclParticle* target,
                                      G4CollisionOutput& output) 
{
  if (verboseLevel > 1)
    G4cout << " >>> G4ElementaryParticleCollider::collide" << G4endl;

  if (!useEPCollider(bullet,target)) {          // Sanity check
    G4cerr << " ElementaryParticleCollider -> can collide only particle with particle " 
           << G4endl;
    return;
  }

#ifdef G4CASCADE_DEBUG_SAMPLER
  static G4bool doPrintTables = true;   // Once and only once per job
  if (doPrintTables) {
    printFinalStateTables();            // For diagnostic reporting
    doPrintTables = false;
  }
#endif

  interCase.set(bullet, target);        // To identify kind of collision

  if (verboseLevel > 1) G4cout << *bullet << G4endl << *target << G4endl;

  G4InuclElementaryParticle* particle1 =
    dynamic_cast<G4InuclElementaryParticle*>(bullet);
  G4InuclElementaryParticle* particle2 =        
    dynamic_cast<G4InuclElementaryParticle*>(target);

  if (!particle1 || !particle2) {       // Redundant with useEPCollider()
    G4cerr << " ElementaryParticleCollider -> can only collide hadrons"
           << G4endl;
    return;
  }

  // Check for available interaction, or pion+dibaryon special case
  if (!G4CascadeChannelTables::GetTable(interCase.hadrons()) &&
      !particle1->quasi_deutron() && !particle2->quasi_deutron()) {
    G4cerr << " ElementaryParticleCollider -> cannot collide " 
           << particle1->getDefinition()->GetParticleName() << " with "
           << particle2->getDefinition()->GetParticleName() << G4endl;
    return;
  }
  // Generate nucleon or pion collision with nucleon
  // or pion with quasi-deuteron

  if (particle1->nucleon() || particle2->nucleon()) { // ok
    G4LorentzConvertor convertToSCM;
    if(particle2->nucleon()) {
      convertToSCM.setBullet(particle1);
      convertToSCM.setTarget(particle2);
    } else {
      convertToSCM.setBullet(particle2);
      convertToSCM.setTarget(particle1);
    };

    convertToSCM.setVerbose(verboseLevel);

    convertToSCM.toTheCenterOfMass();
    G4double ekin = convertToSCM.getKinEnergyInTheTRS();
    G4double etot_scm = convertToSCM.getTotalSCMEnergy();
    G4double pscm = convertToSCM.getSCMMomentum();

    generateSCMfinalState(ekin, etot_scm, pscm, particle1, particle2,
                          &convertToSCM);

    if (particles.empty()) {    // No final state possible, pass bullet through
      if (verboseLevel) {
        G4cerr << " ElementaryParticleCollider -> failed to collide " 
               << particle1->getMomModule() << " GeV/c " 
               << particle1->getDefinition()->GetParticleName() << " with "
               << particle2->getDefinition()->GetParticleName() << G4endl;
      }
    } else {                     // convert back to Lab
      G4LorentzVector mom;              // Buffer to avoid memory churn
      particleIterator ipart;
      for(ipart = particles.begin(); ipart != particles.end(); ipart++) {       
        mom = convertToSCM.backToTheLab(ipart->getMomentum());
        ipart->setMomentum(mom); 
      };

      // Check conservation in multibody final state
      if (verboseLevel && !validateOutput(bullet, target, particles)) {
        G4cout << " incoming particles: \n" << *particle1 << G4endl
               << *particle2 << G4endl
               << " outgoing particles: " << G4endl;
        for(ipart = particles.begin(); ipart != particles.end(); ipart++)
          G4cout << *ipart << G4endl;

        G4cout << " <<< Non-conservation in G4ElementaryParticleCollider"
               << G4endl;
      }

      std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
      output.addOutgoingParticles(particles);
    }
  } else {      // neither particle is nucleon: pion on quasideuteron
    if (particle1->quasi_deutron() || particle2->quasi_deutron()) {
      if (particle1->pion() || particle2->pion() ||
          particle1->isPhoton() || particle2->isPhoton()) {
        G4LorentzConvertor convertToSCM;
        if(particle2->quasi_deutron()) {        // Quasideuteron is target
          convertToSCM.setBullet(particle1);
          convertToSCM.setTarget(particle2);
        } else {
          convertToSCM.setBullet(particle2);
          convertToSCM.setTarget(particle1);
        }; 
        convertToSCM.toTheCenterOfMass(); 
        G4double etot_scm = convertToSCM.getTotalSCMEnergy();
        
        generateSCMpionAbsorption(etot_scm, particle1, particle2);

        if (particles.empty()) {        // Failed to generate final state
          if (verboseLevel) {
            G4cerr << " ElementaryParticleCollider -> failed to collide " 
                   << particle1->getMomModule() << " GeV/c " 
                   << particle1->getDefinition()->GetParticleName() << " with "
                   << particle2->getDefinition()->GetParticleName() << G4endl;
          }
        } else {                        // convert back to Lab
          G4LorentzVector mom;  // Buffer to avoid memory churn
          particleIterator ipart;
          for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
            mom = convertToSCM.backToTheLab(ipart->getMomentum());
            ipart->setMomentum(mom); 
          };

          validateOutput(bullet, target, particles);    // Check conservation

          std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
          output.addOutgoingParticles(particles);
        };
      } else {
        G4cerr << " ElementaryParticleCollider -> can only collide pions with dibaryons " 
               << G4endl;
      };
    } else {
      G4cerr << " ElementaryParticleCollider -> can only collide something with nucleon or dibaryon " 
             << G4endl;
    };
  };  
}


G4int 
G4ElementaryParticleCollider::generateMultiplicity(G4int is, 
                                                   G4double ekin) const 
{
  G4int mul = 0;

  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);

  if (xsecTable) mul = xsecTable->getMultiplicity(ekin);
  else {
    G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
           << is << " - multiplicity not generated " << G4endl;
  }

  if(verboseLevel > 3){
    G4cout << " G4ElementaryParticleCollider::generateMultiplicity: "  
           << " multiplicity = " << mul << G4endl; 
  }

  return mul;
}

 
void
G4ElementaryParticleCollider::generateSCMfinalState(G4double ekin, 
                                     G4double etot_scm, 
                                     G4double pscm,
                                     G4InuclElementaryParticle* particle1,
                                     G4InuclElementaryParticle* particle2, 
                                     G4LorentzConvertor* toSCM) {
  if (verboseLevel > 3) {
    G4cout << " >>> G4ElementaryParticleCollider::generateSCMfinalState" 
           << G4endl;
  }

  const G4double ang_cut = 0.9999;
  const G4double difr_const = 0.3678794;   
  const G4int itry_max = 10;
  G4InuclElementaryParticle dummy;

  G4int type1 = particle1->type();
  G4int type2 = particle2->type();

  G4int is = type1 * type2;

  if(verboseLevel > 3){
    G4cout << " is " << is << G4endl;
  }

  G4int multiplicity = 0;
  G4bool generate = true;

  while (generate) {
    particles.clear();          // Initialize buffers for this event
    particle_kinds.clear();

    // Generate list of final-state particles
    multiplicity = generateMultiplicity(is, ekin);

    generateOutgoingPartTypes(is, multiplicity, ekin);
    if (particle_kinds.empty()) {
      if (verboseLevel > 3) {
        G4cout << " generateOutgoingPartTypes failed mult " << multiplicity
               << G4endl;
      }
      continue;
    }

    if (multiplicity == 2) {
      // Identify charge or strangeness exchange (non-elastic scatter)
      G4int finaltype = particle_kinds[0]*particle_kinds[1];
      G4int kw = (finaltype != is) ? 2 : 1;

      G4double pmod = pscm;     // Elastic scattering preserves CM momentum

      if (kw == 2) {            // Non-elastic needs new CM momentum value
        G4double mone = dummy.getParticleMass(particle_kinds[0]);
        G4double mtwo = dummy.getParticleMass(particle_kinds[1]);

        if (etot_scm < mone+mtwo) {             // Can't produce final state
          if (verboseLevel > 2) {
            G4cerr << " bad final state " << particle_kinds[0]
                   << " , " << particle_kinds[1] << " etot_scm " << etot_scm
                   << " < mone+mtwo " << mone+mtwo << " , but ekin " << ekin
                   << G4endl;
          }
          continue;
        }

        G4double ecm_sq = etot_scm*etot_scm;
        G4double msumsq = mone+mtwo; msumsq *= msumsq;
        G4double mdifsq = mone-mtwo; mdifsq *= mdifsq;

        G4double a = (ecm_sq - msumsq) * (ecm_sq - mdifsq);

        pmod = std::sqrt(a)/(2.*etot_scm);
      }

      G4LorentzVector mom = sampleCMmomentumFor2to2(is, kw, ekin, pmod);

      if (verboseLevel > 3) {
        G4cout << " Particle kinds = " << particle_kinds[0] << " , "
               << particle_kinds[1] << G4endl
               << " pscm " << pscm << " pmod " << pmod << G4endl
               << " before rotation px " << mom.x() << " py " << mom.y()
               << " pz " << mom.z() << G4endl;
      }

      mom = toSCM->rotate(mom); 

      if (verboseLevel > 3){
        G4cout << " after rotation px " << mom.x() << " py " << mom.y() <<
          " pz " << mom.z() << G4endl;
      }
      G4LorentzVector mom1(-mom.vect(), mom.e());

      particles.resize(multiplicity);           // Preallocate buffer
      particles[0].fill(mom, particle_kinds[0], G4InuclParticle::EPCollider);
      particles[1].fill(mom1, particle_kinds[1], G4InuclParticle::EPCollider);
      generate = false;
    } else {                     // 2 -> many
      G4int itry = 0;
      G4bool bad = true;
      G4int knd_last = particle_kinds[multiplicity - 1];
      G4double mass_last = dummy.getParticleMass(knd_last);

      if (verboseLevel > 3){
        G4cout << " knd_last " << knd_last << " mass " << mass_last << G4endl;
      }

      while (bad && itry < itry_max) {
        itry++;

        if (verboseLevel > 3){
          G4cout << " itry in generateSCMfinalState " << itry << G4endl;
        }

        generateMomModules(multiplicity, is, ekin, etot_scm);
        if (G4int(modules.size()) != multiplicity) {
          if (verboseLevel > 3) {
            G4cerr << " generateMomModule failed at mult " << multiplicity
                   << " ekin " << ekin << " etot_scm " << etot_scm << G4endl;
          }
          continue;
        }

        if (multiplicity == 3) { 
          G4LorentzVector mom3 = 
            particleSCMmomentumFor2to3(is, knd_last, ekin, modules[2]);
          
          mom3 = toSCM->rotate(mom3);
          
          // generate the momentum of first
          G4double ct = -0.5 * (modules[2] * modules[2] + 
                                modules[0] * modules[0] - 
                                modules[1] * modules[1]) /
            modules[2] / modules[0];   
          
          if(std::fabs(ct) < ang_cut) {
            
            if(verboseLevel > 2){
              G4cout << " ok for mult " << multiplicity << G4endl;
            }
            
            G4LorentzVector mom1 = generateWithFixedTheta(ct, modules[0]);
            mom1 = toSCM->rotate(mom3, mom1);

            // Second particle recoils off 1 & 3
            G4LorentzVector mom2(etot_scm);
            mom2 -= mom1+mom3;
            
            bad = false;
            generate = false;

            particles.resize(multiplicity);             // Preallocate buffer
   
            particles[0].fill(mom1, particle_kinds[0], G4InuclParticle::EPCollider);
            particles[1].fill(mom2, particle_kinds[1], G4InuclParticle::EPCollider);
            particles[2].fill(mom3, particle_kinds[2], G4InuclParticle::EPCollider);
          };
        } else { // multiplicity > 3
          // generate first mult - 2 momentums
          G4LorentzVector tot_mom;
          scm_momentums.clear();
          
          for (G4int i = 0; i < multiplicity - 2; i++) {
            G4double p0 = particle_kinds[i] < 3 ? 0.36 : 0.25;
            G4double alf = 1.0 / p0 / (p0 - (modules[i] + p0) *
                                       std::exp(-modules[i] / p0));
            G4double st = 2.0;
            G4int itry1 = 0;
            
            while (std::fabs(st) > ang_cut && itry1 < itry_max) {
              itry1++;
              G4double s1 = modules[i] * inuclRndm();
              G4double s2 = alf * difr_const * p0 * inuclRndm();
              
              if(verboseLevel > 3){
                G4cout << " s1 * alf * std::exp(-s1 / p0) "
                       << s1 * alf * std::exp(-s1 / p0) 
                       << " s2 " << s2 << G4endl;
              }
              
              if(s1 * alf * std::exp(-s1 / p0) > s2) st = s1 / modules[i];
            };
            
            if(verboseLevel > 3){
              G4cout << " itry1 " << itry1 << " i " << i << " st " << st
                     << G4endl;
            }
            
            if(itry1 == itry_max) {
              if(verboseLevel > 2){
                G4cout << " high energy angles generation: itry1 " << itry1
                       << G4endl;
              }
              
              st = 0.5 * inuclRndm();
            };

            G4double ct = std::sqrt(1.0 - st * st);
            if (inuclRndm() > 0.5) ct = -ct;
            
            G4LorentzVector mom = generateWithFixedTheta(ct,modules[i]);

            tot_mom += mom;
            
            scm_momentums.push_back(mom);
          }; 
          
          // handle last two
          G4double tot_mod = tot_mom.rho(); 
          G4double ct = -0.5 * (tot_mod * tot_mod + 
                                modules[multiplicity - 2] * modules[multiplicity - 2] -
                                modules[multiplicity - 1] * modules[multiplicity - 1]) / tot_mod /
            modules[multiplicity - 2];  
          
          if (verboseLevel > 2){
            G4cout << " ct last " << ct << G4endl;
          }            
          
          if (std::fabs(ct) < ang_cut) {
            G4int i(0);
            for (i = 0; i < multiplicity - 2; i++) 
              scm_momentums[i] = toSCM->rotate(scm_momentums[i]);
            
            tot_mom = toSCM->rotate(tot_mom);  
            
            G4LorentzVector mom = 
              generateWithFixedTheta(ct, modules[multiplicity - 2]);
            
            mom = toSCM->rotate(tot_mom, mom);
            scm_momentums.push_back(mom);

            // Last particle recoils off everything else
            G4LorentzVector mom1(etot_scm);
            mom1 -= mom+tot_mom;
            
            scm_momentums.push_back(mom1);  
            bad = false;
            generate = false;
            
            if (verboseLevel > 2){
              G4cout << " ok for mult " << multiplicity << G4endl;
            }

            particles.resize(multiplicity);             // Preallocate buffer
            for (i = 0; i < multiplicity; i++) {
              particles[i].fill(scm_momentums[i], particle_kinds[i],
                                G4InuclParticle::EPCollider);
            }
          }     // |ct| < ang_cut
        }       // multiplicity > 3
      }         // while (bad && itry < itry_max)

      if (itry == itry_max) {
        if (verboseLevel > 2) {
          G4cout << " cannot generate the distr. for mult " << multiplicity
                 << G4endl;
        }
        break;
      }
    }   // multiplicity > 2
  }     // while (generate) 

  if (verboseLevel > 3) {
    G4cout << " <<< G4ElementaryParticleCollider::generateSCMfinalState"
           << G4endl;
  }

  return;       // Particles buffer filled
}

void 
G4ElementaryParticleCollider::generateMomModules(G4int mult, 
                                                 G4int is, 
                                                 G4double ekin, 
                                                 G4double etot_cm) {
  if (verboseLevel > 3) {
    G4cout << " >>> G4ElementaryParticleCollider::generateMomModules" 
           << G4endl;
  }

  if (verboseLevel > 2){
    G4cout << " mult " << mult << " is " << is << " ekin " << ekin
           << " etot_cm " << etot_cm << G4endl;
  }

  const G4int itry_max = 10;
  const G4double small = 1.e-10;
  G4InuclElementaryParticle dummy;
  G4int itry = 0;

  modules.clear();                      // Initialize buffer for this attempt
  modules.resize(mult,0.);

  masses2.clear();
  masses2.resize(mult,0.);              // Allows direct [i] setting

  for (G4int i = 0; i < mult; i++) {
    G4double mass = dummy.getParticleMass(particle_kinds[i]);
    masses2[i] = mass * mass;
  };

  G4double mass_last = std::sqrt(masses2[mult - 1]);

  if (verboseLevel > 3){
    G4cout << " knd_last " << particle_kinds[mult - 1] << " mass_last " 
           << mass_last << G4endl;
  }

  while (itry < itry_max) {
    itry++;
    if(verboseLevel > 3){
      G4cout << " itry in generateMomModules " << itry << G4endl;
    }

    G4int ilast = -1;
    G4double eleft = etot_cm;

    for (G4int i = 0; i < mult - 1; i++) {
      G4double pmod = 
        getMomModuleFor2toMany(is, mult, particle_kinds[i], ekin);

      if (pmod < small) break;
      eleft -= std::sqrt(pmod * pmod + masses2[i]);

      if (verboseLevel > 3){
        G4cout << " kp " << particle_kinds[i] << " pmod " << pmod
               << " mass2 " << masses2[i] << " eleft " << eleft
               << "\n x1 " << eleft - mass_last << G4endl;
      }

      if (eleft <= mass_last) break;
      ilast++;
      modules[i] = pmod;
    };

    if (ilast == mult - 2) {
      G4double plast = eleft * eleft - masses2[mult - 1];
      if (verboseLevel > 2){
        G4cout << " plast ** 2 " << plast << G4endl;
      }

      if (plast > small) {
        plast = std::sqrt(plast);
        modules[mult - 1] = plast;      

        if (mult == 3) { 
          if (satisfyTriangle(modules)) return;
        } else return;
      }
    }
  }     // while (itry < itry_max)

  if (verboseLevel > 2)
    G4cerr << " Unable to generate momenta for multiplicity " << mult << G4endl;

  modules.clear();              // Something went wrong, throw away partial
  return;    
}


G4bool 
G4ElementaryParticleCollider::satisfyTriangle(
                        const std::vector<G4double>& pmod) const 
{
  if (verboseLevel > 3) {
    G4cout << " >>> G4ElementaryParticleCollider::satisfyTriangle" 
           << G4endl;
  }

  G4bool good = ( (pmod.size() != 3) ||
                  !(std::fabs(pmod[1] - pmod[2]) > pmod[0] || 
                    pmod[0] > pmod[1] + pmod[2] ||
                    std::fabs(pmod[0] - pmod[2]) > pmod[1] ||
                    pmod[1] > pmod[0] + pmod[2] ||
                    std::fabs(pmod[0] - pmod[1]) > pmod[2] ||
                    pmod[2] > pmod[1] + pmod[0]));

  return good;
}


void 
G4ElementaryParticleCollider::generateOutgoingPartTypes(G4int is, G4int mult,
                                                        G4double ekin)
{
  particle_kinds.clear();       // Initialize buffer for generation

  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(is);

  if (xsecTable)
    xsecTable->getOutgoingParticleTypes(particle_kinds, mult, ekin);
  else {
    G4cerr << " G4ElementaryParticleCollider: Unknown interaction channel "
           << is << " - outgoing kinds not generated " << G4endl;
  }

  return;
}


G4double 
G4ElementaryParticleCollider::getMomModuleFor2toMany(G4int is, G4int /*mult*/, 
                                                     G4int knd, 
                                                     G4double ekin) const 
{
  if (verboseLevel > 2) {
    G4cout << " >>> G4ElementaryParticleCollider::getMomModuleFor2toMany "
           << " is " << is << " knd " << knd << " ekin " << ekin << G4endl;
  }

  G4double S = inuclRndm();
  G4double PS = 0.0;
  G4double PR = 0.0;
  G4double PQ = 0.0;
  G4int KM = 2;
  G4int IL = 4;
  G4int JK = 4;
  G4int JM = 2;
  G4int IM = 3;

  if (is == 1 || is == 2 || is == 4) KM = 1;
  if (knd == 1 || knd == 2) JK = 0;

  if (verboseLevel > 3) {
    G4cout << " S " << S << " KM " << KM << " IL " << IL << " JK " << JK
           << " JM " << JM << " IM " << IM << G4endl;
  }

  for(G4int i = 0; i < 4; i++) {
    G4double V = 0.0;
    for(G4int k = 0; k < 4; k++) {
      if (verboseLevel > 3) {
        G4cout << " k " << k << " : rmn[k+JK][i+IL][KM-1] "
               << rmn[k+JK][i+IL][KM-1] << " ekin^k " << std::pow(ekin, k)
               << G4endl;
      }

      V += rmn[k + JK][i + IL][KM - 1] * std::pow(ekin, k);
    }

    if (verboseLevel > 3) {
      G4cout << " i " << i << " : V " << V << " S^i " << std::pow(S, i)
             << G4endl;
    }

    PR += V * std::pow(S, i);
    PQ += V;
  }

  if (verboseLevel > 3) G4cout << " PR " << PR << " PQ " << PQ << G4endl;

  if (knd == 1 || knd == 2) JM = 1;
  if (verboseLevel > 3) G4cout << " JM " << JM << G4endl;

  for(G4int im = 0; im < 3; im++) {
    if (verboseLevel >3) {
      G4cout << " im " << im << " : rmn[8+IM+im][7+JM][KM-1] "
             << rmn[8+IM+im][7+JM][KM-1] << " ekin^im " << std::pow(ekin, im)
             << G4endl;
    }
    PS += rmn[8 + IM + im][7 + JM][KM - 1] * std::pow(ekin, im);
  }
  
  G4double PRA = PS * std::sqrt(S) * (PR + (1 - PQ) * (S*S*S*S));

  if (verboseLevel > 3) 
    G4cout << " PS " << PS << " PRA = PS*sqrt(S)*(PR+(1-PQ)*S^4) " << PRA
           << G4endl;

  return std::fabs(PRA);
}


G4LorentzVector 
G4ElementaryParticleCollider::particleSCMmomentumFor2to3(
                           G4int is, 
                           G4int knd, 
                           G4double ekin, 
                           G4double pmod) const 
{
  if (verboseLevel > 3) {
    G4cout << " >>> G4ElementaryParticleCollider::particleSCMmomentumFor2to3" 
           << G4endl;
  }

  const G4int itry_max = 100;
  G4double ct = 2.0;
  G4int K = 3;
  G4int J = 1;

  if(is == 1 || is == 2 || is == 4) K = 1;

  if(knd == 1 || knd == 2) J = 0;

  G4int itry = 0;

  while(std::fabs(ct) > 1.0 && itry < itry_max) {
    itry++;
    G4double S = inuclRndm();
    G4double U = 0.0;
    G4double W = 0.0;

    for(G4int l = 0; l < 4; l++) {
      G4double V = 0.0;

      for(G4int im = 0; im < 4; im++) {
        V += abn[im][l][K+J-1] * std::pow(ekin, im);
      };

      U += V;
      W += V * std::pow(S, l);
    };  
    ct = 2.0 * std::sqrt(S) * (W + (1.0 - U) * (S*S*S*S)) - 1.0;
  };

  if(itry == itry_max) {

    if(verboseLevel > 2){
      G4cout << " particleSCMmomentumFor2to3 -> itry = itry_max " << itry << G4endl;
    }

    ct = 2.0 * inuclRndm() - 1.0;
  };

  return generateWithFixedTheta(ct, pmod);
}


G4LorentzVector 
G4ElementaryParticleCollider::sampleCMmomentumFor2to2(G4int is, G4int kw, 
                                                      G4double ekin, 
                                                      G4double pscm) const 
{
  if (verboseLevel > 3)
    G4cout << " >>> G4ElementaryParticleCollider::sampleCMmomentumFor2to2" 
           << " is " << is << " kw " << kw << " ekin " << ekin << " pscm "
           << pscm << G4endl;

  G4double pA=0.0, pC=0.0, pCos=0.0, pFrac=0.0;         // Angular parameters

  // Arrays below are parameters for two-exponential sampling of angular
  // distributions of two-body scattering in the specified channels

  if (is == 1 || is == 2 || is == 4 ||
      is == 21 || is == 23 || is == 25 || is == 27 || is ==29 || is == 31 ||
      is == 42 || is == 46 || is == 50 || is == 54 || is ==58 || is == 62) {
    // nucleon-nucleon or hyperon-nucleon
    if (verboseLevel > 3) G4cout << " nucleon/hyperon elastic" << G4endl;

    static const G4double nnke[9] =  { 0.0,   0.44, 0.59,   0.80,   1.00,   2.24,   4.40,   6.15,  10.00};
    static const G4double nnA[9] =   { 0.0,   0.34, 2.51,   4.59,   4.2,    5.61,   6.38,   7.93,   8.7};
    static const G4double nnC[9] =   { 0.0,   0.0,  1.21,   1.54,   1.88,   1.24,   1.91,   4.04,   8.7};
    static const G4double nnCos[9] = {-1.0,  -1.0, 0.4226, 0.4226, 0.4384, 0.7193, 0.8788, 0.9164,  0.95};
    static const G4double nnFrac[9] = {1.0,   1.0, 0.4898, 0.7243, 0.7990, 0.8892, 0.8493, 0.9583,  1.0};

    static G4CascadeInterpolator<9> interp(nnke);       // Only need one!
    pA = interp.interpolate(ekin, nnA);
    pC = interp.interpolate(ekin, nnC);
    pCos = interp.interpolate(ekin, nnCos);
    pFrac = interp.interpolate(ekin, nnFrac);

  } else if (kw == 2 && (is == 9 || is == 18)) {
    // gamma p -> pi+ n, gamma p -> pi0 p, gamma p -> K Y (and isospin variants)
    // for now and due to lack of better data, use the gamma p -> pi+ n angular
    // distribution for all of these channels
    if (verboseLevel > 3)
      G4cout << " gamma-nucleon inelastic with 2-body final state" << G4endl;

    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};
    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};
    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};
    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};
    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};

    static G4CascadeInterpolator<10> interp(gnke);
    pA = interp.interpolate(ekin, gnA);
    pC = interp.interpolate(ekin, gnC);
    pCos = interp.interpolate(ekin, gnCos);
    pFrac = interp.interpolate(ekin, gnFrac);

  } else if (kw == 2) {
    // pi- p -> pi0 n, pi0 p -> pi+ n, pi- p -> K Y, pi0 p -> K Y (and isospin variants)
    // includes charge and strangeness exchange  
    if (verboseLevel > 3)
      G4cout << " pion-nucleon inelastic with 2-body final state " << G4endl;

    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};
    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};
    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};
    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};
    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};

    static G4CascadeInterpolator<10> interp(qxke);      // Only need one!
    pA = interp.interpolate(ekin, qxA);
    pC = interp.interpolate(ekin, qxC);
    pCos = interp.interpolate(ekin, qxCos);
    pFrac = interp.interpolate(ekin, qxFrac);

  } else if (is == 3 || is == 7 || is == 9 || is == 11 || is == 17 ||
             is == 10 || is == 14 || is == 18 || is == 26 || is == 30) {
    // pi+p, pi0p, gammap, k+p, k0bp, pi-n, pi0n, gamman, k-n, or k0n
    if (verboseLevel > 3) G4cout << " meson-nucleon elastic (1)" << G4endl;

    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};
    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};
    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};
    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};
    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};

    static G4CascadeInterpolator<10> interp(hn1ke);     // Only need one!
    pA = interp.interpolate(ekin, hn1A);
    pC = interp.interpolate(ekin, hn1C);
    pCos = interp.interpolate(ekin, hn1Cos);
    pFrac = interp.interpolate(ekin, hn1Frac);

  } else if (is == 5 || is == 6 || is == 13 || is == 34 || is == 22 ||
             is == 15) {
    // pi-p, pi+n, k-p, k0bn, k+n, or k0p
    if (verboseLevel > 3) G4cout << " meson-nucleon elastic (2)" << G4endl;

    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};
    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};
    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};
    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};
    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};

    static G4CascadeInterpolator<10> interp(hn2ke);     // Only need one!
    pA = interp.interpolate(ekin, hn2A);
    pC = interp.interpolate(ekin, hn2C);
    pCos = interp.interpolate(ekin, hn2Cos);
    pFrac = interp.interpolate(ekin, hn2Frac);

  } else {
    if (verboseLevel)
      G4cerr << " G4ElementaryParticleCollider::sampleCMmomentumFor2to2:"
             << " interaction is=" << is << " not recognized " << G4endl;
  } 

  // Bound parameters by their physical ranges
  pCos = std::max(-1.,std::min(pCos,1.));
  pFrac = std::max(0.,std::min(pFrac,1.));

  // Use parameters determined above to get polar angle
  G4double ct = sampleCMcosFor2to2(pscm, pFrac, pA, pC, pCos);

  return generateWithFixedTheta(ct, pscm);
}


G4double
G4ElementaryParticleCollider::sampleCMcosFor2to2(G4double pscm, G4double pFrac,
                                                 G4double pA, G4double pC,
                                                 G4double pCos) const 
{
  if (verboseLevel>3) {
    G4cout << " sampleCMcosFor2to2: pscm " << pscm << " pFrac " << pFrac
           << " pA " << pA << " pC " << pC << " pCos " << pCos << G4endl;
  }

  G4bool smallAngle = (G4UniformRand() < pFrac);        // 0 < theta < theta0

  G4double term1 = 2.0 * pscm*pscm * (smallAngle ? pA : pC);

  if (std::abs(term1) < 1e-7) return 1.0;       // No actual scattering here!
  if (term1 > DBL_MAX_EXP) return 1.0;

  G4double term2 = std::exp(-2.0*term1);
  G4double randScale = (std::exp(-term1*(1.0 - pCos)) - term2)/(1.0 - term2);

  G4double randVal;
  if (smallAngle) randVal = (1.0 - randScale)*G4UniformRand() + randScale;
  else randVal = randScale*G4UniformRand();

  G4double costheta = 1.0 + std::log(randVal*(1.0 - term2) + term2)/term1;

  if (verboseLevel>3) {
    G4cout << " term1 " << term1 << " term2 " << term2 << " randVal "
           << randVal << " => costheta " << costheta << G4endl;
  }

  return costheta;
}


void
G4ElementaryParticleCollider::generateSCMpionAbsorption(G4double etot_scm,
                                     G4InuclElementaryParticle* particle1,
                                     G4InuclElementaryParticle* particle2) {
  if (verboseLevel > 3)
    G4cout << " >>> G4ElementaryParticleCollider::generateSCMpionAbsorption" 
           << G4endl;

  // generate nucleons momenta for pion or photon absorption
  // the nucleon distribution assumed to be isotropic in SCM

  particles.clear();            // Initialize buffers for this event
  particles.resize(2);

  particle_kinds.clear();

  G4int type1 = particle1->type();
  G4int type2 = particle2->type();

  // generate kinds

  if (type1 == pionPlus) {
    if (type2 == diproton) {            // pi+ + PP -> ? 
      G4cerr << " pion absorption: pi+ + PP -> ? " << G4endl;
      return;
    } else if (type2 == unboundPN) {    // pi+ + PN -> PP
      particle_kinds.push_back(proton);
      particle_kinds.push_back(proton);
    } else if (type2 == dineutron) {    // pi+ + NN -> PN
      particle_kinds.push_back(proton);
      particle_kinds.push_back(neutron);
    }
  } else if (type1 == pionMinus) { 
    if (type2 == diproton) {            // pi- + PP -> PN 
      particle_kinds.push_back(proton);
      particle_kinds.push_back(neutron);
    } else if (type2 == unboundPN) {     // pi- + PN -> NN
      particle_kinds.push_back(neutron);
      particle_kinds.push_back(neutron);
    } else if (type2 == dineutron) {    // pi- + NN -> ?
      G4cerr << " pion absorption: pi- + NN -> ? " << G4endl;
      return;
    }
  } else if (type1 == pionZero || type1 == photon) {
    if (type2 == diproton) {            // pi0/gamma + PP -> PP 
      particle_kinds.push_back(proton);
      particle_kinds.push_back(proton);
    } else if (type2 == unboundPN) {    // pi0/gamma + PN -> PN
      particle_kinds.push_back(proton);
      particle_kinds.push_back(neutron);
    } else if (type2 == dineutron) {    // pi0/gamma + NN -> ?
      particle_kinds.push_back(neutron);
      particle_kinds.push_back(neutron);
    }
  }
    
  G4InuclElementaryParticle dummy;

  G4double mone = dummy.getParticleMass(particle_kinds[0]);
  G4double m1sq = mone*mone;

  G4double mtwo = dummy.getParticleMass(particle_kinds[1]);
  G4double m2sq = mtwo*mtwo;

  G4double a = 0.5 * (etot_scm * etot_scm - m1sq - m2sq);

  G4double pmod = std::sqrt((a * a - m1sq * m2sq) / (m1sq + m2sq + 2.0 * a));
  G4LorentzVector mom1 = generateWithRandomAngles(pmod, mone);
  G4LorentzVector mom2;
  mom2.setVectM(-mom1.vect(), mtwo);

  particles[0].fill(mom1, particle_kinds[0], G4InuclParticle::EPCollider);
  particles[1].fill(mom2, particle_kinds[1], G4InuclParticle::EPCollider);

  return;
}


// Dump lookup tables for N-body final states

void G4ElementaryParticleCollider::
printFinalStateTables(std::ostream& os) const {
  G4CascadeChannelTables::PrintTable(pro*pro, os);
  G4CascadeChannelTables::PrintTable(neu*pro, os);
  G4CascadeChannelTables::PrintTable(neu*neu, os);
  G4CascadeChannelTables::PrintTable(kmi*neu, os);
  G4CascadeChannelTables::PrintTable(kmi*pro, os);
  G4CascadeChannelTables::PrintTable(kpl*neu, os);
  G4CascadeChannelTables::PrintTable(kpl*pro, os);
  G4CascadeChannelTables::PrintTable(k0b*neu, os);
  G4CascadeChannelTables::PrintTable(k0b*pro, os);
  G4CascadeChannelTables::PrintTable(k0*neu, os);
  G4CascadeChannelTables::PrintTable(k0*pro, os);
  G4CascadeChannelTables::PrintTable(lam*neu, os);
  G4CascadeChannelTables::PrintTable(lam*pro, os);
  G4CascadeChannelTables::PrintTable(pim*neu, os);
  G4CascadeChannelTables::PrintTable(pim*pro, os);
  G4CascadeChannelTables::PrintTable(pip*neu, os);
  G4CascadeChannelTables::PrintTable(pip*pro, os);
  G4CascadeChannelTables::PrintTable(pi0*neu, os);
  G4CascadeChannelTables::PrintTable(pi0*pro, os);
  G4CascadeChannelTables::PrintTable(sm*neu, os);
  G4CascadeChannelTables::PrintTable(sm*pro, os);
  G4CascadeChannelTables::PrintTable(sp*neu, os);
  G4CascadeChannelTables::PrintTable(sp*pro, os);
  G4CascadeChannelTables::PrintTable(s0*neu, os);
  G4CascadeChannelTables::PrintTable(s0*pro, os);
  G4CascadeChannelTables::PrintTable(xim*neu, os);
  G4CascadeChannelTables::PrintTable(xim*pro, os);
  G4CascadeChannelTables::PrintTable(xi0*neu, os);
  G4CascadeChannelTables::PrintTable(xi0*pro, os);

  os << " * * * PRELIMINARY -- GAMMA-NUCLEON TABLES * * *" << G4endl;
  G4CascadeChannelTables::PrintTable(gam*neu, os);
  G4CascadeChannelTables::PrintTable(gam*pro, os);
}


// Parameter array for momentum calculation of many body final states
const G4double G4ElementaryParticleCollider::rmn[14][10][2] = {
  {{0.5028,   0.6305}, {3.1442, -3.7333}, {-7.8172,  13.464}, {8.1667, -18.594}, 
   {1.6208,   1.9439}, {-4.3139, -4.6268}, {12.291,  9.7879}, {-15.288, -9.6074}, 
   {   0.0,     0.0}, {   0.0,      0.0}},     

  {{0.9348,   2.1801}, {-10.59,  1.5163}, { 29.227,  -16.38}, {-34.55,  27.944}, 
   {-0.2009, -0.3464}, {1.3641,   1.1093}, {-3.403, -1.9313}, { 3.8559,  1.7064}, 
   {   0.0,     0.0}, {    0.0,     0.0}},    
  
  {{-0.0967, -1.2886}, {4.7335,  -2.457}, {-14.298,  15.129}, {17.685, -23.295}, 
   { 0.0126,  0.0271}, {-0.0835, -0.1164}, { 0.186,  0.2697}, {-0.2004, -0.3185}, 
   {   0.0,     0.0}, {    0.0,     0.0}},    
  
  {{-0.025,   0.2091}, {-0.6248, 0.5228}, { 2.0282, -2.8687}, {-2.5895, 4.2688}, 
   {-0.0002, -0.0007}, {0.0014,   0.0051}, {-0.0024, -0.015}, { 0.0022,  0.0196}, 
   {    0.0,    0.0}, {    0.0,     0.0}},     
  
  {{1.1965,   0.9336}, {-0.8289,-1.8181}, { 1.0426,  5.5157}, { -1.909,-8.5216}, 
   { 1.2419,  1.8693}, {-4.3633, -5.5678}, { 13.743, 14.795}, {-18.592, -16.903}, 
   {    0.0,    0.0}, {    0.0,     0.0}},     
  
  {{0.287,    1.7811}, {-4.9065,-8.2927}, { 16.264,  20.607}, {-19.904,-20.827}, 
   {-0.244,  -0.4996}, {1.3158,   1.7874}, {-3.5691, -4.133}, { 4.3867,  3.8393}, 
   {    0.0,    0.0}, {   0.0,      0.0}}, 
  
  {{-0.2449, -1.5264}, { 2.9191, 6.8433}, {-9.5776, -16.067}, { 11.938, 16.845}, 
   {0.0157,   0.0462}, {-0.0826, -0.1854}, { 0.2143, 0.4531}, {-0.2585, -0.4627}, 
   {    0.0,    0.0}, {   0.0,      0.0}},
  
  {{0.0373,   0.2713}, {-0.422, -1.1944}, { 1.3883,  2.7495}, {-1.7476,-2.9045}, 
   {-0.0003, -0.0013}, {0.0014,   0.0058}, {-0.0034,-0.0146}, { 0.0039,  0.0156}, 
   {    0.0,    0.0}, {    0.0,     0.0}},     
  
  {{   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, {    0.0,     0.0},
   {    0.0,     0.0}, {   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, 
   { 0.1451,  0.0929},{ 0.1538,  0.1303}},  
  
  {{   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, {    0.0,     0.0},
   {    0.0,     0.0}, {   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, 
   { 0.4652,  0.5389},{ 0.2744,  0.4071}},  
  
  {{   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, {    0.0,     0.0},
   {    0.0,     0.0}, {   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0},
   { -0.033, -0.0545},{-0.0146, -0.0288}},  
  
  {{   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, {    0.0,     0.0},
   {    0.0,     0.0}, {   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0},
   { 0.6296,  0.1491},{ 0.8381,  0.1802}},  
  
  {{   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, {    0.0,     0.0},
   {    0.0,     0.0}, {   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0},
   { 0.1787,   0.385},{ 0.0086,  0.3302}},  
  
  {{   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0}, {    0.0,     0.0},
   {    0.0,     0.0}, {   0.0,      0.0}, {    0.0,    0.0}, {    0.0,     0.0},
   {-0.0026, -0.0128},{ 0.0033, -0.0094}}
};

const G4double G4ElementaryParticleCollider::abn[4][4][4] = {
  {{0.0856,  0.0716,  0.1729,  0.0376},  {5.0390,  3.0960,  7.1080,  1.4331},
   {-13.782, -11.125, -17.961, -3.1350},  {14.661,  18.130,  16.403,  6.4864}},
  {{0.0543,  0.0926, -0.1450,  0.2383}, {-9.2324, -3.2186, -13.032,  1.8253},
   {36.397,  20.273,  41.781,  1.7648}, {-42.962, -33.245, -40.799, -16.735}},
  {{-0.0511, -0.0515,  0.0454, -0.1541}, {4.6003,  0.8989,  8.3515, -1.5201},
   {-20.534, -7.5084, -30.260, -1.5692},  {27.731,  13.188,  32.882,  17.185}},
  {{0.0075,  0.0058, -0.0048,  0.0250}, {-0.6253, -0.0017, -1.4095,  0.3059},
   {2.9159,  0.7022,  5.3505,  0.3252}, {-4.1101, -1.4854, -6.0946, -3.5277}} 
};

---