Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-18 22:50:38

0001 //
0002 // ********************************************************************
0003 // * DISCLAIMER                                                       *
0004 // *                                                                  *
0005 // * The following disclaimer summarizes all the specific disclaimers *
0006 // * of contributors to this software. The specific disclaimers,which *
0007 // * govern, are listed with their locations in:                      *
0008 // *   http://cern.ch/geant4/license                                  *
0009 // *                                                                  *
0010 // * Neither the authors of this software system, nor their employing *
0011 // * institutes,nor the agencies providing financial support for this *
0012 // * work  make  any representation or  warranty, express or implied, *
0013 // * regarding  this  software system or assume any liability for its *
0014 // * use.                                                             *
0015 // *                                                                  *
0016 // * This  code  implementation is the  intellectual property  of the *
0017 // * GEANT4 collaboration.                                            *
0018 // * By copying,  distributing  or modifying the Program (or any work *
0019 // * based  on  the Program)  you indicate  your  acceptance of  this *
0020 // * statement, and all its terms.                                    *
0021 // ********************************************************************
0022 //
0023 //
0024 //
0025 // Hadronic Process: Reaction Dynamics
0026 // original by H.P. Wellisch
0027 // modified by J.L. Chuma, TRIUMF, 19-Nov-1996
0028 // Last modified: 27-Mar-1997
0029 // modified by H.P. Wellisch, 24-Apr-97
0030 // H.P. Wellisch, 25.Apr-97: Side of current and target particle taken into account
0031 // H.P. Wellisch, 29.Apr-97: Bug fix in NuclearReaction. (pseudo1 was without energy)
0032 // J.L. Chuma, 30-Apr-97:  Changed return value for GenerateXandPt.  It seems possible
0033 //                         that GenerateXandPt could eliminate some secondaries, but
0034 //                         still finish its calculations, thus we would not want to
0035 //                         then use TwoCluster to again calculate the momenta if vecLen
0036 //                         was less than 6.
0037 // J.L. Chuma, 10-Jun-97:  Modified NuclearReaction.  Was not creating ReactionProduct's
0038 //                         with the new operator, thus they would be meaningless when
0039 //                         going out of scope.
0040 // J.L. Chuma, 20-Jun-97:  Modified GenerateXandPt and TwoCluster to fix units problems
0041 // J.L. Chuma, 23-Jun-97:  Modified ProduceStrangeParticlePairs to fix units problems
0042 // J.L. Chuma, 26-Jun-97:  Modified ProduceStrangeParticlePairs to fix array indices
0043 //                         which were sometimes going out of bounds
0044 // J.L. Chuma, 04-Jul-97:  Many minor modifications to GenerateXandPt and TwoCluster
0045 // J.L. Chuma, 06-Aug-97:  Added original incident particle, before Fermi motion and
0046 //                         evaporation effects are included, needed for self absorption
0047 //                         and corrections for single particle spectra (shower particles)
0048 // logging stopped 1997
0049 // J. Allison, 17-Jun-99:  Replaced a min function to get correct behaviour on DEC.
0050 
0051 #include "SimG4Core/CustomPhysics/interface/FullModelReactionDynamics.h"
0052 #include "G4AntiProton.hh"
0053 #include "G4AntiNeutron.hh"
0054 #include "Randomize.hh"
0055 #include <iostream>
0056 #include "G4ParticleTable.hh"
0057 
0058 using namespace CLHEP;
0059 
0060 G4bool FullModelReactionDynamics::GenerateXandPt(
0061     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
0062     G4int &vecLen,
0063     G4ReactionProduct &modifiedOriginal,      // Fermi motion & evap. effects included
0064     const G4HadProjectile *originalIncident,  // the original incident particle
0065     G4ReactionProduct &currentParticle,
0066     G4ReactionProduct &targetParticle,
0067     const G4Nucleus &targetNucleus,
0068     G4bool &incidentHasChanged,
0069     G4bool &targetHasChanged,
0070     G4bool leadFlag,
0071     G4ReactionProduct &leadingStrangeParticle) {
0072   //
0073   // derived from original FORTRAN code GENXPT by H. Fesefeldt (11-Oct-1987)
0074   //
0075   //  Generation of X- and PT- values for incident, target, and all secondary particles
0076   //  A simple single variable description E D3S/DP3= F(Q) with
0077   //  Q^2 = (M*X)^2 + PT^2 is used. Final state kinematic is produced
0078   //  by an FF-type iterative cascade method
0079   //
0080   //  internal units are GeV
0081   //
0082   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0083 
0084   // Protection in case no secondary has been created; cascades down to two-body.
0085   if (vecLen == 0)
0086     return false;
0087 
0088   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
0089   // G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
0090   // G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
0091   G4ParticleDefinition *aProton = G4Proton::Proton();
0092   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
0093   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
0094   G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
0095   G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
0096   G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
0097   G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
0098   G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
0099 
0100   G4int i, l;
0101   //    G4double forVeryForward = 0.;
0102   G4bool veryForward = false;
0103 
0104   const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
0105   const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
0106   const G4double mOriginal = modifiedOriginal.GetMass() / GeV;
0107   const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
0108   G4double targetMass = targetParticle.GetDefinition()->GetPDGMass() / GeV;
0109   G4double centerofmassEnergy =
0110       std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal);  // GeV
0111   G4double currentMass = currentParticle.GetMass() / GeV;
0112   targetMass = targetParticle.GetMass() / GeV;
0113   //
0114   //  randomize the order of the secondary particles
0115   //  note that the current and target particles are not affected
0116   //
0117   for (i = 0; i < vecLen; ++i) {
0118     G4int itemp = G4int(G4UniformRand() * vecLen);
0119     G4ReactionProduct pTemp = *vec[itemp];
0120     *vec[itemp] = *vec[i];
0121     *vec[i] = pTemp;
0122     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0123   }
0124 
0125   if (currentMass == 0.0 && targetMass == 0.0)  // annihilation
0126   {
0127     // no kinetic energy in target .....
0128     G4double ek = currentParticle.GetKineticEnergy();
0129     G4ThreeVector m = currentParticle.GetMomentum();
0130     currentParticle = *vec[0];
0131     targetParticle = *vec[1];
0132     for (i = 0; i < (vecLen - 2); ++i)
0133       *vec[i] = *vec[i + 2];
0134     G4ReactionProduct *temp = vec[vecLen - 1];
0135     delete temp;
0136     temp = vec[vecLen - 2];
0137     delete temp;
0138     vecLen -= 2;
0139     currentMass = currentParticle.GetMass() / GeV;
0140     targetMass = targetParticle.GetMass() / GeV;
0141     incidentHasChanged = true;
0142     targetHasChanged = true;
0143     currentParticle.SetKineticEnergy(ek);
0144     currentParticle.SetMomentum(m);
0145     veryForward = true;
0146   }
0147   const G4double atomicWeight = targetNucleus.GetN_asInt();
0148   const G4double atomicNumber = targetNucleus.GetZ_asInt();
0149   const G4double protonMass = aProton->GetPDGMass() / MeV;
0150   if ((originalIncident->GetDefinition() == aKaonMinus || originalIncident->GetDefinition() == aKaonZeroL ||
0151        originalIncident->GetDefinition() == aKaonZeroS || originalIncident->GetDefinition() == aKaonPlus) &&
0152       G4UniformRand() >= 0.7) {
0153     G4ReactionProduct temp = currentParticle;
0154     currentParticle = targetParticle;
0155     targetParticle = temp;
0156     incidentHasChanged = true;
0157     targetHasChanged = true;
0158     currentMass = currentParticle.GetMass() / GeV;
0159     targetMass = targetParticle.GetMass() / GeV;
0160   }
0161   const G4double afc = std::min(0.75,
0162                                 0.312 + 0.200 * std::log(std::log(centerofmassEnergy * centerofmassEnergy)) +
0163                                     std::pow(centerofmassEnergy * centerofmassEnergy, 1.5) / 6000.0);
0164 
0165   // PROBLEMET ER HER!!!
0166 
0167   G4double freeEnergy = centerofmassEnergy - currentMass - targetMass;
0168 
0169   if (freeEnergy < 0) {
0170     G4cout << "Free energy < 0!" << G4endl;
0171     G4cout << "E_CMS = " << centerofmassEnergy << " GeV" << G4endl;
0172     G4cout << "m_curr = " << currentMass << " GeV" << G4endl;
0173     G4cout << "m_orig = " << mOriginal << " GeV" << G4endl;
0174     G4cout << "m_targ = " << targetMass << " GeV" << G4endl;
0175     G4cout << "E_free = " << freeEnergy << " GeV" << G4endl;
0176   }
0177 
0178   G4double forwardEnergy = freeEnergy / 2.;
0179   G4int forwardCount = 1;  // number of particles in forward hemisphere
0180 
0181   G4double backwardEnergy = freeEnergy / 2.;
0182   G4int backwardCount = 1;  // number of particles in backward hemisphere
0183   if (veryForward) {
0184     if (currentParticle.GetSide() == -1) {
0185       forwardEnergy += currentMass;
0186       forwardCount--;
0187       backwardEnergy -= currentMass;
0188       backwardCount++;
0189     }
0190     if (targetParticle.GetSide() != -1) {
0191       backwardEnergy += targetMass;
0192       backwardCount--;
0193       forwardEnergy -= targetMass;
0194       forwardCount++;
0195     }
0196   }
0197   for (i = 0; i < vecLen; ++i) {
0198     if (vec[i]->GetSide() == -1) {
0199       ++backwardCount;
0200       backwardEnergy -= vec[i]->GetMass() / GeV;
0201     } else {
0202       ++forwardCount;
0203       forwardEnergy -= vec[i]->GetMass() / GeV;
0204     }
0205   }
0206   //
0207   //  add particles from intranuclear cascade
0208   //  nuclearExcitationCount = number of new secondaries produced by nuclear excitation
0209   //  extraCount = number of nucleons within these new secondaries
0210   //
0211   G4double xtarg;
0212   if (centerofmassEnergy < (2.0 + G4UniformRand()))
0213     xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2.0 * backwardCount + vecLen + 2) / 2.0;
0214   else
0215     xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2.0 * backwardCount);
0216   if (xtarg <= 0.0)
0217     xtarg = 0.01;
0218   G4int nuclearExcitationCount = Poisson(xtarg);
0219   if (atomicWeight < 1.0001)
0220     nuclearExcitationCount = 0;
0221   G4int extraNucleonCount = 0;
0222   G4double extraNucleonMass = 0.0;
0223   if (nuclearExcitationCount > 0) {
0224     const G4double nucsup[] = {1.00, 0.7, 0.5, 0.4, 0.35, 0.3};
0225     const G4double psup[] = {3., 6., 20., 50., 100., 1000.};
0226     G4int momentumBin = 0;
0227     while ((momentumBin < 6) && (modifiedOriginal.GetTotalMomentum() / GeV > psup[momentumBin]))
0228       ++momentumBin;
0229     momentumBin = std::min(5, momentumBin);
0230     //
0231     //  NOTE: in GENXPT, these new particles were given negative codes
0232     //        here I use  NewlyAdded = true  instead
0233     //
0234     for (i = 0; i < nuclearExcitationCount; ++i) {
0235       G4ReactionProduct *pVec = new G4ReactionProduct();
0236       if (G4UniformRand() < nucsup[momentumBin]) {
0237         if (G4UniformRand() > 1.0 - atomicNumber / atomicWeight)
0238           pVec->SetDefinition(aProton);
0239         else
0240           pVec->SetDefinition(aNeutron);
0241         pVec->SetSide(-2);  // -2 means backside nucleon
0242         ++extraNucleonCount;
0243         backwardEnergy += pVec->GetMass() / GeV;
0244         extraNucleonMass += pVec->GetMass() / GeV;
0245       } else {
0246         G4double ran = G4UniformRand();
0247         if (ran < 0.3181)
0248           pVec->SetDefinition(aPiPlus);
0249         else if (ran < 0.6819)
0250           pVec->SetDefinition(aPiZero);
0251         else
0252           pVec->SetDefinition(aPiMinus);
0253         pVec->SetSide(-1);  // backside particle, but not a nucleon
0254       }
0255       pVec->SetNewlyAdded(true);  // true is the same as IPA(i)<0
0256       vec.SetElement(vecLen++, pVec);
0257       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0258       backwardEnergy -= pVec->GetMass() / GeV;
0259       ++backwardCount;
0260     }
0261   }
0262   //
0263   //  assume conservation of kinetic energy in forward & backward hemispheres
0264   //
0265   G4int is, iskip;
0266   while (forwardEnergy <= 0.0)  // must eliminate a particle from the forward side
0267   {
0268     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0269     iskip = G4int(G4UniformRand() * forwardCount) + 1;  // 1 <= iskip <= forwardCount
0270     is = 0;
0271     G4int forwardParticlesLeft = 0;
0272     for (i = (vecLen - 1); i >= 0; --i) {
0273       if (vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled()) {
0274         forwardParticlesLeft = 1;
0275         if (++is == iskip) {
0276           forwardEnergy += vec[i]->GetMass() / GeV;
0277           for (G4int j = i; j < (vecLen - 1); j++)
0278             *vec[j] = *vec[j + 1];  // shift up
0279           --forwardCount;
0280           G4ReactionProduct *temp = vec[vecLen - 1];
0281           delete temp;
0282           if (--vecLen == 0)
0283             return false;  // all the secondaries have been eliminated
0284           break;           // --+
0285         }                  //   |
0286       }                    //   |
0287     }                      // break goes down to here
0288     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0289     if (forwardParticlesLeft == 0) {
0290       forwardEnergy += currentParticle.GetMass() / GeV;
0291       currentParticle.SetDefinitionAndUpdateE(targetParticle.GetDefinition());
0292       targetParticle.SetDefinitionAndUpdateE(vec[0]->GetDefinition());
0293       // above two lines modified 20-oct-97: were just simple equalities
0294       --forwardCount;
0295       for (G4int j = 0; j < (vecLen - 1); ++j)
0296         *vec[j] = *vec[j + 1];
0297       G4ReactionProduct *temp = vec[vecLen - 1];
0298       delete temp;
0299       if (--vecLen == 0)
0300         return false;  // all the secondaries have been eliminated
0301       break;
0302     }
0303   }
0304   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0305   while (backwardEnergy <= 0.0)  // must eliminate a particle from the backward side
0306   {
0307     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0308     iskip = G4int(G4UniformRand() * backwardCount) + 1;  // 1 <= iskip <= backwardCount
0309     is = 0;
0310     G4int backwardParticlesLeft = 0;
0311     for (i = (vecLen - 1); i >= 0; --i) {
0312       if (vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled()) {
0313         backwardParticlesLeft = 1;
0314         if (++is == iskip)  // eliminate the i'th particle
0315         {
0316           if (vec[i]->GetSide() == -2) {
0317             --extraNucleonCount;
0318             extraNucleonMass -= vec[i]->GetMass() / GeV;
0319             backwardEnergy -= vec[i]->GetTotalEnergy() / GeV;
0320           }
0321           backwardEnergy += vec[i]->GetTotalEnergy() / GeV;
0322           for (G4int j = i; j < (vecLen - 1); ++j)
0323             *vec[j] = *vec[j + 1];  // shift up
0324           --backwardCount;
0325           G4ReactionProduct *temp = vec[vecLen - 1];
0326           delete temp;
0327           if (--vecLen == 0)
0328             return false;  // all the secondaries have been eliminated
0329           break;
0330         }
0331       }
0332     }
0333     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0334     if (backwardParticlesLeft == 0) {
0335       backwardEnergy += targetParticle.GetMass() / GeV;
0336       targetParticle = *vec[0];
0337       --backwardCount;
0338       for (G4int j = 0; j < (vecLen - 1); ++j)
0339         *vec[j] = *vec[j + 1];
0340       G4ReactionProduct *temp = vec[vecLen - 1];
0341       delete temp;
0342       if (--vecLen == 0)
0343         return false;  // all the secondaries have been eliminated
0344       break;
0345     }
0346   }
0347   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0348   //
0349   //  define initial state vectors for Lorentz transformations
0350   //  the pseudoParticles have non-standard masses, hence the "pseudo"
0351   //
0352   G4ReactionProduct pseudoParticle[10];
0353   for (i = 0; i < 10; ++i)
0354     pseudoParticle[i].SetZero();
0355 
0356   pseudoParticle[0].SetMass(mOriginal * GeV);
0357   pseudoParticle[0].SetMomentum(0.0, 0.0, pOriginal * GeV);
0358   pseudoParticle[0].SetTotalEnergy(std::sqrt(pOriginal * pOriginal + mOriginal * mOriginal) * GeV);
0359 
0360   pseudoParticle[1].SetMass(protonMass * MeV);  // this could be targetMass
0361   pseudoParticle[1].SetTotalEnergy(protonMass * MeV);
0362 
0363   pseudoParticle[3].SetMass(protonMass * (1 + extraNucleonCount) * MeV);
0364   pseudoParticle[3].SetTotalEnergy(protonMass * (1 + extraNucleonCount) * MeV);
0365 
0366   pseudoParticle[8].SetMomentum(1.0 * GeV, 0.0, 0.0);
0367 
0368   pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
0369   pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
0370 
0371   pseudoParticle[0].Lorentz(pseudoParticle[0], pseudoParticle[2]);
0372   pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[2]);
0373 
0374   G4double dndl[20];
0375   //
0376   //  main loop for 4-momentum generation
0377   //  see Pitha-report (Aachen) for a detailed description of the method
0378   //
0379   G4double aspar, pt, et, x, pp, pp1, rthnve, phinve, rmb, wgt;
0380   G4int innerCounter, outerCounter;
0381   G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
0382 
0383   G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
0384   //
0385   // process the secondary particles in reverse order
0386   // the incident particle is Done after the secondaries
0387   // nucleons, including the target, in the backward hemisphere are also Done later
0388   //
0389   G4double binl[20] = {0.,  0.1,  0.2,  0.3,  0.4,  0.5, 0.6, 0.7,  0.8,  0.9,
0390                        1.0, 1.11, 1.25, 1.43, 1.67, 2.0, 2.5, 3.33, 5.00, 10.00};
0391   G4int backwardNucleonCount = 0;  // number of nucleons in backward hemisphere
0392   G4double totalEnergy, kineticEnergy, vecMass;
0393 
0394   for (i = (vecLen - 1); i >= 0; --i) {
0395     G4double phi = G4UniformRand() * twopi;
0396     if (vec[i]->GetNewlyAdded())  // added from intranuclear cascade
0397     {
0398       if (vec[i]->GetSide() == -2)  //  is a nucleon
0399       {
0400         if (backwardNucleonCount < 18) {
0401           if (vec[i]->GetDefinition() == G4PionMinus::PionMinus() ||
0402               vec[i]->GetDefinition() == G4PionPlus::PionPlus() || vec[i]->GetDefinition() == G4PionZero::PionZero()) {
0403             for (G4int i = 0; i < vecLen; i++)
0404               delete vec[i];
0405             vecLen = 0;
0406             G4ExceptionDescription ed;
0407             ed << "FullModelReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon";
0408             G4Exception("FullModelReactionDynamics::GenerateXandPt", "had064", FatalException, ed);
0409           }
0410           vec[i]->SetSide(-3);
0411           ++backwardNucleonCount;
0412           continue;
0413         }
0414       }
0415     }
0416     //
0417     //  set pt and phi values, they are changed somewhat in the iteration loop
0418     //  set mass parameter for lambda fragmentation model
0419     //
0420     vecMass = vec[i]->GetMass() / GeV;
0421     G4double ran = -std::log(1.0 - G4UniformRand()) / 3.5;
0422     if (vec[i]->GetSide() == -2)  // backward nucleon
0423     {
0424       if (vec[i]->GetDefinition() == aKaonMinus || vec[i]->GetDefinition() == aKaonZeroL ||
0425           vec[i]->GetDefinition() == aKaonZeroS || vec[i]->GetDefinition() == aKaonPlus ||
0426           vec[i]->GetDefinition() == aPiMinus || vec[i]->GetDefinition() == aPiZero ||
0427           vec[i]->GetDefinition() == aPiPlus) {
0428         aspar = 0.75;
0429         pt = std::sqrt(std::pow(ran, 1.7));
0430       } else {         // vec[i] must be a proton, neutron,
0431         aspar = 0.20;  //  lambda, sigma, xsi, or ion
0432         pt = std::sqrt(std::pow(ran, 1.2));
0433       }
0434     } else {  // not a backward nucleon
0435       if (vec[i]->GetDefinition() == aPiMinus || vec[i]->GetDefinition() == aPiZero ||
0436           vec[i]->GetDefinition() == aPiPlus) {
0437         aspar = 0.75;
0438         pt = std::sqrt(std::pow(ran, 1.7));
0439       } else if (vec[i]->GetDefinition() == aKaonMinus || vec[i]->GetDefinition() == aKaonZeroL ||
0440                  vec[i]->GetDefinition() == aKaonZeroS || vec[i]->GetDefinition() == aKaonPlus) {
0441         aspar = 0.70;
0442         pt = std::sqrt(std::pow(ran, 1.7));
0443       } else {         // vec[i] must be a proton, neutron,
0444         aspar = 0.65;  //  lambda, sigma, xsi, or ion
0445         pt = std::sqrt(std::pow(ran, 1.5));
0446       }
0447     }
0448     pt = std::max(0.001, pt);
0449     vec[i]->SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0450     for (G4int j = 0; j < 20; ++j)
0451       binl[j] = j / (19. * pt);
0452     if (vec[i]->GetSide() > 0)
0453       et = pseudoParticle[0].GetTotalEnergy() / GeV;
0454     else
0455       et = pseudoParticle[1].GetTotalEnergy() / GeV;
0456     dndl[0] = 0.0;
0457     //
0458     //   start of outer iteration loop
0459     //
0460     outerCounter = 0;
0461     eliminateThisParticle = true;
0462     resetEnergies = true;
0463     while (++outerCounter < 3) {
0464       for (l = 1; l < 20; ++l) {
0465         x = (binl[l] + binl[l - 1]) / 2.;
0466         pt = std::max(0.001, pt);
0467         if (x > 1.0 / pt)
0468           dndl[l] += dndl[l - 1];  //  changed from just  =  on 02 April 98
0469         else
0470           dndl[l] = et * aspar / std::sqrt(std::pow((1. + aspar * x * aspar * x), 3)) * (binl[l] - binl[l - 1]) /
0471                         std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
0472                     dndl[l - 1];
0473       }
0474       innerCounter = 0;
0475       vec[i]->SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0476       //
0477       //   start of inner iteration loop
0478       //
0479       while (++innerCounter < 7) {
0480         ran = G4UniformRand() * dndl[19];
0481         l = 1;
0482         while ((ran >= dndl[l]) && (l < 20))
0483           l++;
0484         l = std::min(19, l);
0485         x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
0486         if (vec[i]->GetSide() < 0)
0487           x *= -1.;
0488         vec[i]->SetMomentum(x * et * GeV);  // set the z-momentum
0489         totalEnergy = std::sqrt(x * et * x * et + pt * pt + vecMass * vecMass);
0490         vec[i]->SetTotalEnergy(totalEnergy * GeV);
0491         kineticEnergy = vec[i]->GetKineticEnergy() / GeV;
0492         if (vec[i]->GetSide() > 0)  // forward side
0493         {
0494           if ((forwardKinetic + kineticEnergy) < 0.95 * forwardEnergy) {
0495             pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
0496             forwardKinetic += kineticEnergy;
0497             pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0498             pseudoParticle[6].SetMomentum(0.0);  // set the z-momentum
0499             phi = pseudoParticle[6].Angle(pseudoParticle[8]);
0500             if (pseudoParticle[6].GetMomentum().y() / MeV < 0.0)
0501               phi = twopi - phi;
0502             phi += pi + normal() * pi / 12.0;
0503             if (phi > twopi)
0504               phi -= twopi;
0505             if (phi < 0.0)
0506               phi = twopi - phi;
0507             outerCounter = 2;               // leave outer loop
0508             eliminateThisParticle = false;  // don't eliminate this particle
0509             resetEnergies = false;
0510             break;  // leave inner loop
0511           }
0512           if (innerCounter > 5)
0513             break;                        // leave inner loop
0514           if (backwardEnergy >= vecMass)  // switch sides
0515           {
0516             vec[i]->SetSide(-1);
0517             forwardEnergy += vecMass;
0518             backwardEnergy -= vecMass;
0519             ++backwardCount;
0520           }
0521         } else {  // backward side
0522           G4double xxx = 0.95 + 0.05 * extraNucleonCount / 20.0;
0523           if ((backwardKinetic + kineticEnergy) < xxx * backwardEnergy) {
0524             pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
0525             backwardKinetic += kineticEnergy;
0526             pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0527             pseudoParticle[6].SetMomentum(0.0);  // set the z-momentum
0528             phi = pseudoParticle[6].Angle(pseudoParticle[8]);
0529             if (pseudoParticle[6].GetMomentum().y() / MeV < 0.0)
0530               phi = twopi - phi;
0531             phi += pi + normal() * pi / 12.0;
0532             if (phi > twopi)
0533               phi -= twopi;
0534             if (phi < 0.0)
0535               phi = twopi - phi;
0536             outerCounter = 2;               // leave outer loop
0537             eliminateThisParticle = false;  // don't eliminate this particle
0538             resetEnergies = false;
0539             break;  // leave inner loop
0540           }
0541           if (innerCounter > 5)
0542             break;                       // leave inner loop
0543           if (forwardEnergy >= vecMass)  // switch sides
0544           {
0545             vec[i]->SetSide(1);
0546             forwardEnergy -= vecMass;
0547             backwardEnergy += vecMass;
0548             backwardCount--;
0549           }
0550         }
0551         G4ThreeVector momentum = vec[i]->GetMomentum();
0552         vec[i]->SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
0553         pt *= 0.9;
0554         dndl[19] *= 0.9;
0555       }  // closes inner loop
0556       if (resetEnergies) {
0557         //   if we get to here, the inner loop has been Done 6 Times
0558         //   reset the kinetic energies of previously Done particles, if they are lighter
0559         //    than protons and in the forward hemisphere
0560         //   then continue with outer loop
0561         //
0562         forwardKinetic = 0.0;
0563         backwardKinetic = 0.0;
0564         pseudoParticle[4].SetZero();
0565         pseudoParticle[5].SetZero();
0566         for (l = i + 1; l < vecLen; ++l) {
0567           if (vec[l]->GetSide() > 0 || vec[l]->GetDefinition() == aKaonMinus || vec[l]->GetDefinition() == aKaonZeroL ||
0568               vec[l]->GetDefinition() == aKaonZeroS || vec[l]->GetDefinition() == aKaonPlus ||
0569               vec[l]->GetDefinition() == aPiMinus || vec[l]->GetDefinition() == aPiZero ||
0570               vec[l]->GetDefinition() == aPiPlus) {
0571             G4double tempMass = vec[l]->GetMass() / MeV;
0572             totalEnergy = 0.95 * vec[l]->GetTotalEnergy() / MeV + 0.05 * tempMass;
0573             totalEnergy = std::max(tempMass, totalEnergy);
0574             vec[l]->SetTotalEnergy(totalEnergy * MeV);
0575             pp = std::sqrt(std::abs(totalEnergy * totalEnergy - tempMass * tempMass));
0576             pp1 = vec[l]->GetMomentum().mag() / MeV;
0577             if (pp1 < 1.0e-6 * GeV) {
0578               G4double rthnve = pi * G4UniformRand();
0579               G4double phinve = twopi * G4UniformRand();
0580               G4double srth = std::sin(rthnve);
0581               vec[l]->SetMomentum(
0582                   pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0583             } else {
0584               vec[l]->SetMomentum(vec[l]->GetMomentum() * (pp / pp1));
0585             }
0586             G4double px = vec[l]->GetMomentum().x() / MeV;
0587             G4double py = vec[l]->GetMomentum().y() / MeV;
0588             pt = std::max(1.0, std::sqrt(px * px + py * py)) / GeV;
0589             if (vec[l]->GetSide() > 0) {
0590               forwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0591               pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
0592             } else {
0593               backwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0594               pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
0595             }
0596           }
0597         }
0598       }
0599     }  // closes outer loop
0600 
0601     if (eliminateThisParticle && vec[i]->GetMayBeKilled())  // not enough energy, eliminate this particle
0602     {
0603       if (vec[i]->GetSide() > 0) {
0604         --forwardCount;
0605         forwardEnergy += vecMass;
0606       } else {
0607         if (vec[i]->GetSide() == -2) {
0608           --extraNucleonCount;
0609           extraNucleonMass -= vecMass;
0610           backwardEnergy -= vecMass;
0611         }
0612         --backwardCount;
0613         backwardEnergy += vecMass;
0614       }
0615       for (G4int j = i; j < (vecLen - 1); ++j)
0616         *vec[j] = *vec[j + 1];  // shift up
0617       G4ReactionProduct *temp = vec[vecLen - 1];
0618       delete temp;
0619       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0620       if (--vecLen == 0)
0621         return false;  // all the secondaries have been eliminated
0622       pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0623       pseudoParticle[6].SetMomentum(0.0);  // set z-momentum
0624     }
0625   }  // closes main for loop
0626 
0627   //
0628   //  for the incident particle:  it was placed in the forward hemisphere
0629   //   set pt and phi values, they are changed somewhat in the iteration loop
0630   //   set mass parameter for lambda fragmentation model
0631   //
0632   G4double phi = G4UniformRand() * twopi;
0633   G4double ran = -std::log(1.0 - G4UniformRand());
0634   if (currentParticle.GetDefinition() == aPiMinus || currentParticle.GetDefinition() == aPiZero ||
0635       currentParticle.GetDefinition() == aPiPlus) {
0636     aspar = 0.60;
0637     pt = std::sqrt(std::pow(ran / 6.0, 1.7));
0638   } else if (currentParticle.GetDefinition() == aKaonMinus || currentParticle.GetDefinition() == aKaonZeroL ||
0639              currentParticle.GetDefinition() == aKaonZeroS || currentParticle.GetDefinition() == aKaonPlus) {
0640     aspar = 0.50;
0641     pt = std::sqrt(std::pow(ran / 5.0, 1.4));
0642   } else {
0643     aspar = 0.40;
0644     pt = std::sqrt(std::pow(ran / 4.0, 1.2));
0645   }
0646   for (G4int j = 0; j < 20; ++j)
0647     binl[j] = j / (19. * pt);
0648   currentParticle.SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0649   et = pseudoParticle[0].GetTotalEnergy() / GeV;
0650   dndl[0] = 0.0;
0651   vecMass = currentParticle.GetMass() / GeV;
0652   for (l = 1; l < 20; ++l) {
0653     x = (binl[l] + binl[l - 1]) / 2.;
0654     if (x > 1.0 / pt)
0655       dndl[l] += dndl[l - 1];  //  changed from just  =   on 02 April 98
0656     else
0657       dndl[l] = aspar / std::sqrt(std::pow((1. + sqr(aspar * x)), 3)) * (binl[l] - binl[l - 1]) * et /
0658                     std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
0659                 dndl[l - 1];
0660   }
0661   ran = G4UniformRand() * dndl[19];
0662   l = 1;
0663   while ((ran > dndl[l]) && (l < 20))
0664     l++;
0665   l = std::min(19, l);
0666   x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
0667   currentParticle.SetMomentum(x * et * GeV);  // set the z-momentum
0668   if (forwardEnergy < forwardKinetic)
0669     totalEnergy = vecMass + 0.04 * std::fabs(normal());
0670   else
0671     totalEnergy = vecMass + forwardEnergy - forwardKinetic;
0672   currentParticle.SetTotalEnergy(totalEnergy * GeV);
0673   pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0674   pp1 = currentParticle.GetMomentum().mag() / MeV;
0675   if (pp1 < 1.0e-6 * GeV) {
0676     G4double rthnve = pi * G4UniformRand();
0677     G4double phinve = twopi * G4UniformRand();
0678     G4double srth = std::sin(rthnve);
0679     currentParticle.SetMomentum(
0680         pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0681   } else {
0682     currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
0683   }
0684   pseudoParticle[4] = pseudoParticle[4] + currentParticle;
0685   //
0686   // this finishes the current particle
0687   // now for the target particle
0688   //
0689   if (backwardNucleonCount < 18) {
0690     targetParticle.SetSide(-3);
0691     ++backwardNucleonCount;
0692   } else {
0693     //  set pt and phi values, they are changed somewhat in the iteration loop
0694     //  set mass parameter for lambda fragmentation model
0695     //
0696     vecMass = targetParticle.GetMass() / GeV;
0697     ran = -std::log(1.0 - G4UniformRand());
0698     aspar = 0.40;
0699     pt = std::max(0.001, std::sqrt(std::pow(ran / 4.0, 1.2)));
0700     targetParticle.SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0701     for (G4int j = 0; j < 20; ++j)
0702       binl[j] = (j - 1.) / (19. * pt);
0703     et = pseudoParticle[1].GetTotalEnergy() / GeV;
0704     dndl[0] = 0.0;
0705     outerCounter = 0;
0706     resetEnergies = true;
0707     while (++outerCounter < 3)  // start of outer iteration loop
0708     {
0709       for (l = 1; l < 20; ++l) {
0710         x = (binl[l] + binl[l - 1]) / 2.;
0711         if (x > 1.0 / pt)
0712           dndl[l] += dndl[l - 1];  // changed from just  =  on 02 April 98
0713         else
0714           dndl[l] = aspar / std::sqrt(std::pow((1. + aspar * x * aspar * x), 3)) * (binl[l] - binl[l - 1]) * et /
0715                         std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
0716                     dndl[l - 1];
0717       }
0718       innerCounter = 0;
0719       while (++innerCounter < 7)  // start of inner iteration loop
0720       {
0721         l = 1;
0722         ran = G4UniformRand() * dndl[19];
0723         while ((ran >= dndl[l]) && (l < 20))
0724           l++;
0725         l = std::min(19, l);
0726         x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
0727         if (targetParticle.GetSide() < 0)
0728           x *= -1.;
0729         targetParticle.SetMomentum(x * et * GeV);  // set the z-momentum
0730         totalEnergy = std::sqrt(x * et * x * et + pt * pt + vecMass * vecMass);
0731         targetParticle.SetTotalEnergy(totalEnergy * GeV);
0732         if (targetParticle.GetSide() < 0) {
0733           G4double xxx = 0.95 + 0.05 * extraNucleonCount / 20.0;
0734           if ((backwardKinetic + totalEnergy - vecMass) < xxx * backwardEnergy) {
0735             pseudoParticle[5] = pseudoParticle[5] + targetParticle;
0736             backwardKinetic += totalEnergy - vecMass;
0737             pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0738             pseudoParticle[6].SetMomentum(0.0);  // set z-momentum
0739             outerCounter = 2;                    // leave outer loop
0740             resetEnergies = false;
0741             break;  // leave inner loop
0742           }
0743           if (innerCounter > 5)
0744             break;                       // leave inner loop
0745           if (forwardEnergy >= vecMass)  // switch sides
0746           {
0747             targetParticle.SetSide(1);
0748             forwardEnergy -= vecMass;
0749             backwardEnergy += vecMass;
0750             --backwardCount;
0751           }
0752           G4ThreeVector momentum = targetParticle.GetMomentum();
0753           targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
0754           pt *= 0.9;
0755           dndl[19] *= 0.9;
0756         } else  // target has gone to forward side
0757         {
0758           if (forwardEnergy < forwardKinetic)
0759             totalEnergy = vecMass + 0.04 * std::fabs(normal());
0760           else
0761             totalEnergy = vecMass + forwardEnergy - forwardKinetic;
0762           targetParticle.SetTotalEnergy(totalEnergy * GeV);
0763           pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0764           pp1 = targetParticle.GetMomentum().mag() / MeV;
0765           if (pp1 < 1.0e-6 * GeV) {
0766             G4double rthnve = pi * G4UniformRand();
0767             G4double phinve = twopi * G4UniformRand();
0768             G4double srth = std::sin(rthnve);
0769             targetParticle.SetMomentum(
0770                 pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0771           } else
0772             targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
0773 
0774           pseudoParticle[4] = pseudoParticle[4] + targetParticle;
0775           outerCounter = 2;  // leave outer loop
0776           //eliminateThisParticle = false;       // don't eliminate this particle
0777           resetEnergies = false;
0778           break;  // leave inner loop
0779         }
0780       }  // closes inner loop
0781       if (resetEnergies) {
0782         //   if we get to here, the inner loop has been Done 6 Times
0783         //   reset the kinetic energies of previously Done particles, if they are lighter
0784         //    than protons and in the forward hemisphere
0785         //   then continue with outer loop
0786 
0787         forwardKinetic = backwardKinetic = 0.0;
0788         pseudoParticle[4].SetZero();
0789         pseudoParticle[5].SetZero();
0790         for (l = 0; l < vecLen; ++l)  // changed from l=1  on 02 April 98
0791         {
0792           if (vec[l]->GetSide() > 0 || vec[l]->GetDefinition() == aKaonMinus || vec[l]->GetDefinition() == aKaonZeroL ||
0793               vec[l]->GetDefinition() == aKaonZeroS || vec[l]->GetDefinition() == aKaonPlus ||
0794               vec[l]->GetDefinition() == aPiMinus || vec[l]->GetDefinition() == aPiZero ||
0795               vec[l]->GetDefinition() == aPiPlus) {
0796             G4double tempMass = vec[l]->GetMass() / GeV;
0797             totalEnergy = std::max(tempMass, 0.95 * vec[l]->GetTotalEnergy() / GeV + 0.05 * tempMass);
0798             vec[l]->SetTotalEnergy(totalEnergy * GeV);
0799             pp = std::sqrt(std::abs(totalEnergy * totalEnergy - tempMass * tempMass)) * GeV;
0800             pp1 = vec[l]->GetMomentum().mag() / MeV;
0801             if (pp1 < 1.0e-6 * GeV) {
0802               G4double rthnve = pi * G4UniformRand();
0803               G4double phinve = twopi * G4UniformRand();
0804               G4double srth = std::sin(rthnve);
0805               vec[l]->SetMomentum(
0806                   pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0807             } else
0808               vec[l]->SetMomentum(vec[l]->GetMomentum() * (pp / pp1));
0809 
0810             pt = std::max(0.001 * GeV,
0811                           std::sqrt(sqr(vec[l]->GetMomentum().x() / MeV) + sqr(vec[l]->GetMomentum().y() / MeV))) /
0812                  GeV;
0813             if (vec[l]->GetSide() > 0) {
0814               forwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0815               pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
0816             } else {
0817               backwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0818               pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
0819             }
0820           }
0821         }
0822       }
0823     }  // closes outer loop
0824   }
0825   //
0826   //  this finishes the target particle
0827   // backward nucleons produced with a cluster model
0828   //
0829   pseudoParticle[6].Lorentz(pseudoParticle[3], pseudoParticle[2]);
0830   pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
0831   pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
0832   if (backwardNucleonCount == 1)  // target particle is the only backward nucleon
0833   {
0834     G4double ekin = std::min(backwardEnergy - backwardKinetic, centerofmassEnergy / 2.0 - protonMass / GeV);
0835     if (ekin < 0.04)
0836       ekin = 0.04 * std::fabs(normal());
0837     vecMass = targetParticle.GetMass() / GeV;
0838     totalEnergy = ekin + vecMass;
0839     targetParticle.SetTotalEnergy(totalEnergy * GeV);
0840     pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0841     pp1 = pseudoParticle[6].GetMomentum().mag() / MeV;
0842     if (pp1 < 1.0e-6 * GeV) {
0843       rthnve = pi * G4UniformRand();
0844       phinve = twopi * G4UniformRand();
0845       G4double srth = std::sin(rthnve);
0846       targetParticle.SetMomentum(
0847           pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0848     } else {
0849       targetParticle.SetMomentum(pseudoParticle[6].GetMomentum() * (pp / pp1));
0850     }
0851     pseudoParticle[5] = pseudoParticle[5] + targetParticle;
0852   } else  // more than one backward nucleon
0853   {
0854     const G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
0855     const G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
0856     // Replaced the following min function to get correct behaviour on DEC.
0857     G4int tempCount = std::max(1, std::min(5, backwardNucleonCount)) - 1;
0858     //cout << "backwardNucleonCount " << backwardNucleonCount << G4endl;
0859     //cout << "tempCount " << tempCount << G4endl;
0860     G4double rmb0 = 0.0;
0861     if (targetParticle.GetSide() == -3)
0862       rmb0 += targetParticle.GetMass() / GeV;
0863     for (i = 0; i < vecLen; ++i) {
0864       if (vec[i]->GetSide() == -3)
0865         rmb0 += vec[i]->GetMass() / GeV;
0866     }
0867     rmb = rmb0 + std::pow(-std::log(1.0 - G4UniformRand()), cpar[tempCount]) / gpar[tempCount];
0868     totalEnergy = pseudoParticle[6].GetTotalEnergy() / GeV;
0869     vecMass = std::min(rmb, totalEnergy);
0870     pseudoParticle[6].SetMass(vecMass * GeV);
0871     pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0872     pp1 = pseudoParticle[6].GetMomentum().mag() / MeV;
0873     if (pp1 < 1.0e-6 * GeV) {
0874       rthnve = pi * G4UniformRand();
0875       phinve = twopi * G4UniformRand();
0876       G4double srth = std::sin(rthnve);
0877       pseudoParticle[6].SetMomentum(
0878           -pp * srth * std::cos(phinve) * MeV, -pp * srth * std::sin(phinve) * MeV, -pp * std::cos(rthnve) * MeV);
0879     } else
0880       pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp / pp1));
0881 
0882     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;  // tempV contains the backward nucleons
0883     tempV.Initialize(backwardNucleonCount);
0884     G4int tempLen = 0;
0885     if (targetParticle.GetSide() == -3)
0886       tempV.SetElement(tempLen++, &targetParticle);
0887     for (i = 0; i < vecLen; ++i) {
0888       if (vec[i]->GetSide() == -3)
0889         tempV.SetElement(tempLen++, vec[i]);
0890     }
0891     if (tempLen != backwardNucleonCount) {
0892       G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
0893       G4cerr << "tempLen = " << tempLen;
0894       G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
0895       G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
0896       G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
0897       for (i = 0; i < vecLen; ++i)
0898         G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
0899       exit(EXIT_FAILURE);
0900     }
0901     constantCrossSection = true;
0902     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0903     if (tempLen >= 2) {
0904       GenerateNBodyEvent(pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen);
0905       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0906       if (targetParticle.GetSide() == -3) {
0907         targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
0908         // tempV contains the real stuff
0909         pseudoParticle[5] = pseudoParticle[5] + targetParticle;
0910       }
0911       for (i = 0; i < vecLen; ++i) {
0912         if (vec[i]->GetSide() == -3) {
0913           vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
0914           pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
0915         }
0916       }
0917       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0918     }
0919   }
0920   //
0921   //  Lorentz transformation in lab system
0922   //
0923   if (vecLen == 0)
0924     return false;  // all the secondaries have been eliminated
0925   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0926 
0927   G4int numberofFinalStateNucleons = 0;
0928   if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron ||
0929       currentParticle.GetDefinition() == G4SigmaMinus::SigmaMinus() ||
0930       currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
0931       currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
0932       currentParticle.GetDefinition() == G4XiZero::XiZero() ||
0933       currentParticle.GetDefinition() == G4XiMinus::XiMinus() ||
0934       currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
0935       currentParticle.GetDefinition() == G4Lambda::Lambda())
0936     ++numberofFinalStateNucleons;
0937   currentParticle.Lorentz(currentParticle, pseudoParticle[1]);
0938 
0939   if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron ||
0940       targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() ||
0941       targetParticle.GetDefinition() == G4XiMinus::XiMinus() ||
0942       targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
0943       targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
0944       targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
0945       targetParticle.GetDefinition() == G4SigmaMinus::SigmaMinus())
0946     ++numberofFinalStateNucleons;
0947   if (targetParticle.GetDefinition() == G4AntiProton::AntiProton())
0948     --numberofFinalStateNucleons;
0949   if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron())
0950     --numberofFinalStateNucleons;
0951   if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
0952     --numberofFinalStateNucleons;
0953   if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
0954     --numberofFinalStateNucleons;
0955   if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
0956     --numberofFinalStateNucleons;
0957   if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero())
0958     --numberofFinalStateNucleons;
0959   if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus())
0960     --numberofFinalStateNucleons;
0961   if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
0962     --numberofFinalStateNucleons;
0963   if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda())
0964     --numberofFinalStateNucleons;
0965   targetParticle.Lorentz(targetParticle, pseudoParticle[1]);
0966 
0967   for (i = 0; i < vecLen; ++i) {
0968     if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron ||
0969         vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() ||
0970         vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
0971         vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() ||
0972         vec[i]->GetDefinition() == G4SigmaMinus::SigmaMinus())
0973       ++numberofFinalStateNucleons;
0974     if (vec[i]->GetDefinition() == G4AntiProton::AntiProton())
0975       --numberofFinalStateNucleons;
0976     if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron())
0977       --numberofFinalStateNucleons;
0978     if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
0979       --numberofFinalStateNucleons;
0980     if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
0981       --numberofFinalStateNucleons;
0982     if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
0983       --numberofFinalStateNucleons;
0984     if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda())
0985       --numberofFinalStateNucleons;
0986     if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero())
0987       --numberofFinalStateNucleons;
0988     if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus())
0989       --numberofFinalStateNucleons;
0990     if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
0991       --numberofFinalStateNucleons;
0992     vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
0993   }
0994   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
0995   if (veryForward)
0996     numberofFinalStateNucleons++;
0997   numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
0998   //
0999   // leadFlag will be true
1000   //  iff original particle is at least as heavy as K+ and not a proton or neutron AND
1001   //   if
1002   //    incident particle is at least as heavy as K+ and it is not a proton or neutron
1003   //     leadFlag is set to the incident particle
1004   //   or
1005   //    target particle is at least as heavy as K+ and it is not a proton or neutron
1006   //     leadFlag is set to the target particle
1007   //
1008   G4bool leadingStrangeParticleHasChanged = true;
1009   if (leadFlag) {
1010     if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1011       leadingStrangeParticleHasChanged = false;
1012     if (leadingStrangeParticleHasChanged && (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()))
1013       leadingStrangeParticleHasChanged = false;
1014     if (leadingStrangeParticleHasChanged) {
1015       for (i = 0; i < vecLen; i++) {
1016         if (vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition()) {
1017           leadingStrangeParticleHasChanged = false;
1018           break;
1019         }
1020       }
1021     }
1022     if (leadingStrangeParticleHasChanged) {
1023       G4bool leadTest =
1024           (leadingStrangeParticle.GetDefinition() == aKaonMinus ||
1025            leadingStrangeParticle.GetDefinition() == aKaonZeroL ||
1026            leadingStrangeParticle.GetDefinition() == aKaonZeroS ||
1027            leadingStrangeParticle.GetDefinition() == aKaonPlus || leadingStrangeParticle.GetDefinition() == aPiMinus ||
1028            leadingStrangeParticle.GetDefinition() == aPiZero || leadingStrangeParticle.GetDefinition() == aPiPlus);
1029       G4bool targetTest = false;
1030 
1031       // following modified by JLC 22-Oct-97
1032 
1033       if ((leadTest && targetTest) || !(leadTest || targetTest))  // both true or both false
1034       {
1035         targetParticle.SetDefinitionAndUpdateE(leadingStrangeParticle.GetDefinition());
1036         targetHasChanged = true;
1037         // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1038       } else {
1039         currentParticle.SetDefinitionAndUpdateE(leadingStrangeParticle.GetDefinition());
1040         incidentHasChanged = false;
1041         // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1042       }
1043     }
1044   }  // end of if( leadFlag )
1045 
1046   pseudoParticle[3].SetMomentum(0.0, 0.0, pOriginal * GeV);
1047   pseudoParticle[3].SetMass(mOriginal * GeV);
1048   pseudoParticle[3].SetTotalEnergy(std::sqrt(pOriginal * pOriginal + mOriginal * mOriginal) * GeV);
1049 
1050   const G4ParticleDefinition *aOrgDef = modifiedOriginal.GetDefinition();
1051   G4int diff = 0;
1052   if (aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron())
1053     diff = 1;
1054   if (numberofFinalStateNucleons == 1)
1055     diff = 0;
1056   pseudoParticle[4].SetMomentum(0.0, 0.0, 0.0);
1057   pseudoParticle[4].SetMass(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1058   pseudoParticle[4].SetTotalEnergy(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1059 
1060   G4double theoreticalKinetic = pseudoParticle[3].GetTotalEnergy() / MeV + pseudoParticle[4].GetTotalEnergy() / MeV -
1061                                 currentParticle.GetMass() / MeV - targetParticle.GetMass() / MeV;
1062 
1063   G4double simulatedKinetic = currentParticle.GetKineticEnergy() / MeV + targetParticle.GetKineticEnergy() / MeV;
1064 
1065   pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
1066   pseudoParticle[3].Lorentz(pseudoParticle[3], pseudoParticle[5]);
1067   pseudoParticle[4].Lorentz(pseudoParticle[4], pseudoParticle[5]);
1068 
1069   pseudoParticle[7].SetZero();
1070   pseudoParticle[7] = pseudoParticle[7] + currentParticle;
1071   pseudoParticle[7] = pseudoParticle[7] + targetParticle;
1072 
1073   for (i = 0; i < vecLen; ++i) {
1074     pseudoParticle[7] = pseudoParticle[7] + *vec[i];
1075     simulatedKinetic += vec[i]->GetKineticEnergy() / MeV;
1076     theoreticalKinetic -= vec[i]->GetMass() / MeV;
1077   }
1078   if (vecLen <= 16 && vecLen > 0) {
1079     // must create a new set of ReactionProducts here because GenerateNBody will
1080     // modify the momenta for the particles, and we don't want to do this
1081     //
1082     G4ReactionProduct tempR[130];
1083     tempR[0] = currentParticle;
1084     tempR[1] = targetParticle;
1085     for (i = 0; i < vecLen; ++i)
1086       tempR[i + 2] = *vec[i];
1087     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1088     tempV.Initialize(vecLen + 2);
1089     G4int tempLen = 0;
1090     for (i = 0; i < vecLen + 2; ++i)
1091       tempV.SetElement(tempLen++, &tempR[i]);
1092     constantCrossSection = true;
1093 
1094     wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() / MeV + pseudoParticle[4].GetTotalEnergy() / MeV,
1095                              constantCrossSection,
1096                              tempV,
1097                              tempLen);
1098     if (wgt > -.5) {
1099       theoreticalKinetic = 0.0;
1100       for (i = 0; i < tempLen; ++i) {
1101         pseudoParticle[6].Lorentz(*tempV[i], pseudoParticle[4]);
1102         theoreticalKinetic += pseudoParticle[6].GetKineticEnergy() / MeV;
1103       }
1104     }
1105     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1106   }
1107   //
1108   //  Make sure, that the kinetic energies are correct
1109   //
1110   if (simulatedKinetic != 0.0) {
1111     wgt = (theoreticalKinetic) / simulatedKinetic;
1112     theoreticalKinetic = currentParticle.GetKineticEnergy() / MeV * wgt;
1113     simulatedKinetic = theoreticalKinetic;
1114     currentParticle.SetKineticEnergy(theoreticalKinetic * MeV);
1115     pp = currentParticle.GetTotalMomentum() / MeV;
1116     pp1 = currentParticle.GetMomentum().mag() / MeV;
1117     if (pp1 < 1.0e-6 * GeV) {
1118       rthnve = pi * G4UniformRand();
1119       phinve = twopi * G4UniformRand();
1120       currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1121                                   pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1122                                   pp * std::cos(rthnve) * MeV);
1123     } else {
1124       currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1125     }
1126     theoreticalKinetic = targetParticle.GetKineticEnergy() / MeV * wgt;
1127     targetParticle.SetKineticEnergy(theoreticalKinetic * MeV);
1128     simulatedKinetic += theoreticalKinetic;
1129     pp = targetParticle.GetTotalMomentum() / MeV;
1130     pp1 = targetParticle.GetMomentum().mag() / MeV;
1131     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1132     if (pp1 < 1.0e-6 * GeV) {
1133       rthnve = pi * G4UniformRand();
1134       phinve = twopi * G4UniformRand();
1135       targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1136                                  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1137                                  pp * std::cos(rthnve) * MeV);
1138     } else {
1139       targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1140     }
1141     for (i = 0; i < vecLen; ++i) {
1142       theoreticalKinetic = vec[i]->GetKineticEnergy() / MeV * wgt;
1143       simulatedKinetic += theoreticalKinetic;
1144       vec[i]->SetKineticEnergy(theoreticalKinetic * MeV);
1145       pp = vec[i]->GetTotalMomentum() / MeV;
1146       pp1 = vec[i]->GetMomentum().mag() / MeV;
1147       if (pp1 < 1.0e-6 * GeV) {
1148         rthnve = pi * G4UniformRand();
1149         phinve = twopi * G4UniformRand();
1150         vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1151                             pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1152                             pp * std::cos(rthnve) * MeV);
1153       } else
1154         vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
1155     }
1156   }
1157   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1158   Rotate(numberofFinalStateNucleons,
1159          pseudoParticle[3].GetMomentum(),
1160          modifiedOriginal,
1161          originalIncident,
1162          targetNucleus,
1163          currentParticle,
1164          targetParticle,
1165          vec,
1166          vecLen);
1167   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1168   //
1169   // add black track particles
1170   // the total number of particles produced is restricted to 198
1171   // this may have influence on very high energies
1172   //
1173   if (atomicWeight >= 1.5) {
1174     // npnb is number of proton/neutron black track particles
1175     // ndta is the number of deuterons, tritons, and alphas produced
1176     // epnb is the kinetic energy available for proton/neutron black track particles
1177     // edta is the kinetic energy available for deuteron/triton/alpha particles
1178     //
1179     G4double epnb, edta;
1180     G4int npnb = 0;
1181     G4int ndta = 0;
1182 
1183     epnb = targetNucleus.GetPNBlackTrackEnergy();   // was enp1 in fortran code
1184     edta = targetNucleus.GetDTABlackTrackEnergy();  // was enp3 in fortran code
1185     const G4double pnCutOff = 0.001;
1186     const G4double dtaCutOff = 0.001;
1187     const G4double kineticMinimum = 1.e-6;
1188     const G4double kineticFactor = -0.010;
1189     G4double sprob = 0.0;  // sprob = probability of self-absorption in heavy molecules
1190     const G4double ekIncident = originalIncident->GetKineticEnergy() / GeV;
1191     if (ekIncident >= 5.0)
1192       sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0));
1193     if (epnb >= pnCutOff) {
1194       npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta));
1195       if (numberofFinalStateNucleons + npnb > atomicWeight)
1196         npnb = G4int(atomicWeight + 0.00001 - numberofFinalStateNucleons);
1197       npnb = std::min(npnb, 127 - vecLen);
1198     }
1199     if (edta >= dtaCutOff) {
1200       ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta));
1201       ndta = std::min(ndta, 127 - vecLen);
1202     }
1203     G4double spall = numberofFinalStateNucleons;
1204     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1205 
1206     AddBlackTrackParticles(epnb,
1207                            npnb,
1208                            edta,
1209                            ndta,
1210                            sprob,
1211                            kineticMinimum,
1212                            kineticFactor,
1213                            modifiedOriginal,
1214                            spall,
1215                            targetNucleus,
1216                            vec,
1217                            vecLen);
1218 
1219     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1220   }
1221   //
1222   //  calculate time delay for nuclear reactions
1223   //
1224   if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
1225     currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
1226   else
1227     currentParticle.SetTOF(1.0);
1228   return true;
1229   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1230 }
1231 
1232 void FullModelReactionDynamics::SuppressChargedPions(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
1233                                                      G4int &vecLen,
1234                                                      const G4ReactionProduct &modifiedOriginal,
1235                                                      G4ReactionProduct &currentParticle,
1236                                                      G4ReactionProduct &targetParticle,
1237                                                      const G4Nucleus &targetNucleus,
1238                                                      G4bool &incidentHasChanged,
1239                                                      G4bool &targetHasChanged) {
1240   // this code was originally in the fortran code TWOCLU
1241   //
1242   // suppress charged pions, for various reasons
1243   //
1244   const G4double atomicWeight = targetNucleus.GetN_asInt();
1245   const G4double atomicNumber = targetNucleus.GetZ_asInt();
1246   const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / GeV;
1247 
1248   // G4ParticleDefinition *aGamma = G4Gamma::Gamma();
1249   G4ParticleDefinition *aProton = G4Proton::Proton();
1250   G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
1251   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1252   G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
1253   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1254   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1255 
1256   const G4bool antiTest =
1257       modifiedOriginal.GetDefinition() != anAntiProton && modifiedOriginal.GetDefinition() != anAntiNeutron;
1258   if (antiTest &&
1259       (
1260           // currentParticle.GetDefinition() == aGamma ||
1261           currentParticle.GetDefinition() == aPiPlus || currentParticle.GetDefinition() == aPiMinus) &&
1262       (G4UniformRand() <= (10.0 - pOriginal) / 6.0) && (G4UniformRand() <= atomicWeight / 300.0)) {
1263     if (G4UniformRand() > atomicNumber / atomicWeight)
1264       currentParticle.SetDefinitionAndUpdateE(aNeutron);
1265     else
1266       currentParticle.SetDefinitionAndUpdateE(aProton);
1267     incidentHasChanged = true;
1268   }
1269   for (G4int i = 0; i < vecLen; ++i) {
1270     if (antiTest &&
1271         (
1272             // vec[i]->GetDefinition() == aGamma ||
1273             vec[i]->GetDefinition() == aPiPlus || vec[i]->GetDefinition() == aPiMinus) &&
1274         (G4UniformRand() <= (10.0 - pOriginal) / 6.0) && (G4UniformRand() <= atomicWeight / 300.0)) {
1275       if (G4UniformRand() > atomicNumber / atomicWeight)
1276         vec[i]->SetDefinitionAndUpdateE(aNeutron);
1277       else
1278         vec[i]->SetDefinitionAndUpdateE(aProton);
1279       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1280     }
1281   }
1282   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1283 }
1284 
1285 G4bool FullModelReactionDynamics::TwoCluster(
1286     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
1287     G4int &vecLen,
1288     G4ReactionProduct &modifiedOriginal,      // Fermi motion & evap. effects included
1289     const G4HadProjectile *originalIncident,  // the original incident particle
1290     G4ReactionProduct &currentParticle,
1291     G4ReactionProduct &targetParticle,
1292     const G4Nucleus &targetNucleus,
1293     G4bool &incidentHasChanged,
1294     G4bool &targetHasChanged,
1295     G4bool leadFlag,
1296     G4ReactionProduct &leadingStrangeParticle) {
1297   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1298   // derived from original FORTRAN code TWOCLU by H. Fesefeldt (11-Oct-1987)
1299   //
1300   //  Generation of X- and PT- values for incident, target, and all secondary particles
1301   //
1302   //  A simple two cluster model is used.
1303   //  This should be sufficient for low energy interactions.
1304   //
1305 
1306   // + debugging
1307   // raise(SIGSEGV);
1308   // - debugging
1309 
1310   G4int i;
1311   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1312   G4ParticleDefinition *aProton = G4Proton::Proton();
1313   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1314   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1315   G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
1316   G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
1317   const G4double protonMass = aProton->GetPDGMass() / MeV;
1318   const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
1319   const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
1320   const G4double mOriginal = modifiedOriginal.GetMass() / GeV;
1321   const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
1322   G4double targetMass = targetParticle.GetDefinition()->GetPDGMass() / GeV;
1323   G4double centerofmassEnergy =
1324       std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal);  // GeV
1325   G4double currentMass = currentParticle.GetMass() / GeV;
1326   targetMass = targetParticle.GetMass() / GeV;
1327 
1328   if (currentMass == 0.0 && targetMass == 0.0) {
1329     G4double ek = currentParticle.GetKineticEnergy();
1330     G4ThreeVector m = currentParticle.GetMomentum();
1331     currentParticle = *vec[0];
1332     targetParticle = *vec[1];
1333     for (i = 0; i < (vecLen - 2); ++i)
1334       *vec[i] = *vec[i + 2];
1335     if (vecLen < 2) {
1336       for (G4int i = 0; i < vecLen; i++)
1337         delete vec[i];
1338       vecLen = 0;
1339       G4ExceptionDescription ed;
1340       ed << "Negative number of particles";
1341       G4Exception("FullModelReactionDynamics::TwoCluster", "had064", FatalException, ed);
1342     }
1343     delete vec[vecLen - 1];
1344     delete vec[vecLen - 2];
1345     vecLen -= 2;
1346     incidentHasChanged = true;
1347     targetHasChanged = true;
1348     currentParticle.SetKineticEnergy(ek);
1349     currentParticle.SetMomentum(m);
1350   }
1351   const G4double atomicWeight = targetNucleus.GetN_asInt();
1352   const G4double atomicNumber = targetNucleus.GetZ_asInt();
1353   //
1354   // particles have been distributed in forward and backward hemispheres
1355   // in center of mass system of the hadron nucleon interaction
1356   //
1357   // incident is always in forward hemisphere
1358   G4int forwardCount = 1;  // number of particles in forward hemisphere
1359   currentParticle.SetSide(1);
1360   G4double forwardMass = currentParticle.GetMass() / GeV;
1361   G4double cMass = forwardMass;
1362 
1363   // target is always in backward hemisphere
1364   G4int backwardCount = 1;         // number of particles in backward hemisphere
1365   G4int backwardNucleonCount = 1;  // number of nucleons in backward hemisphere
1366   targetParticle.SetSide(-1);
1367   G4double backwardMass = targetParticle.GetMass() / GeV;
1368   G4double bMass = backwardMass;
1369 
1370   for (i = 0; i < vecLen; ++i) {
1371     if (vec[i]->GetSide() < 0)
1372       vec[i]->SetSide(-1);  // added by JLC, 2Jul97
1373     // to take care of the case where vec has been preprocessed by GenerateXandPt
1374     // and some of them have been set to -2 or -3
1375     if (vec[i]->GetSide() == -1) {
1376       ++backwardCount;
1377       backwardMass += vec[i]->GetMass() / GeV;
1378     } else {
1379       ++forwardCount;
1380       forwardMass += vec[i]->GetMass() / GeV;
1381     }
1382   }
1383   //
1384   // nucleons and some pions from intranuclear cascade
1385   //
1386   G4double term1 = std::log(centerofmassEnergy * centerofmassEnergy);
1387   if (term1 < 0)
1388     term1 = 0.0001;  // making sure xtarg<0;
1389   const G4double afc = 0.312 + 0.2 * std::log(term1);
1390   G4double xtarg;
1391   if (centerofmassEnergy < 2.0 + G4UniformRand())  // added +2 below, JLC 4Jul97
1392     xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2 * backwardCount + vecLen + 2) / 2.0;
1393   else
1394     xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2 * backwardCount);
1395   if (xtarg <= 0.0)
1396     xtarg = 0.01;
1397   G4int nuclearExcitationCount = Poisson(xtarg);
1398   if (atomicWeight < 1.0001)
1399     nuclearExcitationCount = 0;
1400   G4int extraNucleonCount = 0;
1401   G4double extraMass = 0.0;
1402   G4double extraNucleonMass = 0.0;
1403   if (nuclearExcitationCount > 0) {
1404     G4int momentumBin = std::min(4, G4int(pOriginal / 3.0));
1405     const G4double nucsup[] = {1.0, 0.8, 0.6, 0.5, 0.4};
1406     //
1407     //  NOTE: in TWOCLU, these new particles were given negative codes
1408     //        here we use  NewlyAdded = true  instead
1409     //
1410     for (i = 0; i < nuclearExcitationCount; ++i) {
1411       G4ReactionProduct *pVec = new G4ReactionProduct();
1412       if (G4UniformRand() < nucsup[momentumBin])  // add proton or neutron
1413       {
1414         if (G4UniformRand() > 1.0 - atomicNumber / atomicWeight)
1415           // HPW: looks like a gheisha bug
1416           pVec->SetDefinition(aProton);
1417         else
1418           pVec->SetDefinition(aNeutron);
1419         ++backwardNucleonCount;
1420         ++extraNucleonCount;
1421         extraNucleonMass += pVec->GetMass() / GeV;
1422       } else {  // add a pion
1423         G4double ran = G4UniformRand();
1424         if (ran < 0.3181)
1425           pVec->SetDefinition(aPiPlus);
1426         else if (ran < 0.6819)
1427           pVec->SetDefinition(aPiZero);
1428         else
1429           pVec->SetDefinition(aPiMinus);
1430       }
1431       pVec->SetSide(-2);  // backside particle
1432       extraMass += pVec->GetMass() / GeV;
1433       pVec->SetNewlyAdded(true);
1434       vec.SetElement(vecLen++, pVec);
1435       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1436     }
1437   }
1438   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1439   G4double forwardEnergy = (centerofmassEnergy - cMass - bMass) / 2.0 + cMass - forwardMass;
1440   G4double backwardEnergy = (centerofmassEnergy - cMass - bMass) / 2.0 + bMass - backwardMass;
1441   G4double eAvailable = centerofmassEnergy - (forwardMass + backwardMass);
1442   G4bool secondaryDeleted;
1443   G4double pMass;
1444   while (eAvailable <= 0.0)  // must eliminate a particle
1445   {
1446     secondaryDeleted = false;
1447     for (i = (vecLen - 1); i >= 0; --i) {
1448       if (vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled()) {
1449         pMass = vec[i]->GetMass() / GeV;
1450         for (G4int j = i; j < (vecLen - 1); ++j)
1451           *vec[j] = *vec[j + 1];  // shift up
1452         --forwardCount;
1453         forwardEnergy += pMass;
1454         forwardMass -= pMass;
1455         secondaryDeleted = true;
1456         break;
1457       } else if (vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled()) {
1458         pMass = vec[i]->GetMass() / GeV;
1459         for (G4int j = i; j < (vecLen - 1); ++j)
1460           *vec[j] = *vec[j + 1];  // shift up
1461         --backwardCount;
1462         backwardEnergy += pMass;
1463         backwardMass -= pMass;
1464         secondaryDeleted = true;
1465         break;
1466       }
1467     }  // breaks go down to here
1468     if (secondaryDeleted) {
1469       G4ReactionProduct *temp = vec[vecLen - 1];
1470       delete temp;
1471       --vecLen;
1472       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1473     } else {
1474       if (vecLen == 0) {
1475         return false;  // all secondaries have been eliminated
1476       }
1477       if (targetParticle.GetSide() == -1) {
1478         pMass = targetParticle.GetMass() / GeV;
1479         targetParticle = *vec[0];
1480         for (G4int j = 0; j < (vecLen - 1); ++j)
1481           *vec[j] = *vec[j + 1];  // shift up
1482         --backwardCount;
1483         backwardEnergy += pMass;
1484         backwardMass -= pMass;
1485         secondaryDeleted = true;
1486       } else if (targetParticle.GetSide() == 1) {
1487         pMass = targetParticle.GetMass() / GeV;
1488         targetParticle = *vec[0];
1489         for (G4int j = 0; j < (vecLen - 1); ++j)
1490           *vec[j] = *vec[j + 1];  // shift up
1491         --forwardCount;
1492         forwardEnergy += pMass;
1493         forwardMass -= pMass;
1494         secondaryDeleted = true;
1495       }
1496       if (secondaryDeleted) {
1497         G4ReactionProduct *temp = vec[vecLen - 1];
1498         delete temp;
1499         --vecLen;
1500         // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1501       } else {
1502         if (currentParticle.GetSide() == -1) {
1503           pMass = currentParticle.GetMass() / GeV;
1504           currentParticle = *vec[0];
1505           for (G4int j = 0; j < (vecLen - 1); ++j)
1506             *vec[j] = *vec[j + 1];  // shift up
1507           --backwardCount;
1508           backwardEnergy += pMass;
1509           backwardMass -= pMass;
1510           secondaryDeleted = true;
1511         } else if (currentParticle.GetSide() == 1) {
1512           pMass = currentParticle.GetMass() / GeV;
1513           currentParticle = *vec[0];
1514           for (G4int j = 0; j < (vecLen - 1); ++j)
1515             *vec[j] = *vec[j + 1];  // shift up
1516           --forwardCount;           //This line can cause infinite loop
1517           forwardEnergy += pMass;
1518           forwardMass -= pMass;
1519           secondaryDeleted = true;
1520         }
1521         if (secondaryDeleted) {
1522           G4ReactionProduct *temp = vec[vecLen - 1];
1523           delete temp;
1524           --vecLen;
1525           // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1526         } else
1527           break;
1528       }
1529     }
1530     eAvailable = centerofmassEnergy - (forwardMass + backwardMass);
1531   }
1532   //
1533   // This is the start of the TwoCluster function
1534   //  Choose masses for the 3 clusters:
1535   //   forward cluster
1536   //   backward meson cluster
1537   //   backward nucleon cluster
1538   //
1539   G4double rmc = 0.0, rmd = 0.0;  // rme = 0.0;
1540   const G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
1541   const G4double gpar[] = {2.6, 2.6, 1.8, 1.30, 1.20};
1542 
1543   if (forwardCount == 0)
1544     return false;
1545 
1546   if (forwardCount == 1)
1547     rmc = forwardMass;
1548   else {
1549     G4int ntc = std::max(1, std::min(5, forwardCount)) - 1;  // check if offset by 1 @@
1550     rmc = forwardMass + std::pow(-std::log(1.0 - G4UniformRand()), cpar[ntc]) / gpar[ntc];
1551   }
1552   if (backwardCount == 1)
1553     rmd = backwardMass;
1554   else {
1555     G4int ntc = std::max(1, std::min(5, backwardCount)) - 1;  // check, if offfset by 1 @@
1556     rmd = backwardMass + std::pow(-std::log(1.0 - G4UniformRand()), cpar[ntc]) / gpar[ntc];
1557   }
1558   while (rmc + rmd > centerofmassEnergy) {
1559     if ((rmc <= forwardMass) && (rmd <= backwardMass)) {
1560       G4double temp = 0.999 * centerofmassEnergy / (rmc + rmd);
1561       rmc *= temp;
1562       rmd *= temp;
1563     } else {
1564       rmc = 0.1 * forwardMass + 0.9 * rmc;
1565       rmd = 0.1 * backwardMass + 0.9 * rmd;
1566     }
1567   }
1568   //
1569   //  Set beam, target of first interaction in centre of mass system
1570   //
1571   G4ReactionProduct pseudoParticle[8];
1572   for (i = 0; i < 8; ++i)
1573     pseudoParticle[i].SetZero();
1574 
1575   pseudoParticle[1].SetMass(mOriginal * GeV);
1576   pseudoParticle[1].SetTotalEnergy(etOriginal * GeV);
1577   pseudoParticle[1].SetMomentum(0.0, 0.0, pOriginal * GeV);
1578 
1579   pseudoParticle[2].SetMass(protonMass * MeV);
1580   pseudoParticle[2].SetTotalEnergy(protonMass * MeV);
1581   pseudoParticle[2].SetMomentum(0.0, 0.0, 0.0);
1582   //
1583   //  transform into centre of mass system
1584   //
1585   pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
1586   pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[0]);
1587   pseudoParticle[2].Lorentz(pseudoParticle[2], pseudoParticle[0]);
1588 
1589   const G4double pfMin = 0.0001;
1590   G4double pf = (centerofmassEnergy * centerofmassEnergy + rmd * rmd - rmc * rmc);
1591   pf *= pf;
1592   pf -= 4 * centerofmassEnergy * centerofmassEnergy * rmd * rmd;
1593   pf = std::sqrt(std::max(pf, pfMin)) / (2.0 * centerofmassEnergy);
1594   //
1595   //  set final state masses and energies in centre of mass system
1596   //
1597   pseudoParticle[3].SetMass(rmc * GeV);
1598   pseudoParticle[3].SetTotalEnergy(std::sqrt(pf * pf + rmc * rmc) * GeV);
1599 
1600   pseudoParticle[4].SetMass(rmd * GeV);
1601   pseudoParticle[4].SetTotalEnergy(std::sqrt(pf * pf + rmd * rmd) * GeV);
1602   //
1603   // set |T| and |TMIN|
1604   //
1605   const G4double bMin = 0.01;
1606   const G4double b1 = 4.0;
1607   const G4double b2 = 1.6;
1608   G4double t = std::log(1.0 - G4UniformRand()) / std::max(bMin, b1 + b2 * std::log(pOriginal));
1609   G4double t1 = pseudoParticle[1].GetTotalEnergy() / GeV - pseudoParticle[3].GetTotalEnergy() / GeV;
1610   G4double pin = pseudoParticle[1].GetMomentum().mag() / GeV;
1611   G4double tacmin = t1 * t1 - (pin - pf) * (pin - pf);
1612   //
1613   // calculate (std::sin(teta/2.)^2 and std::cos(teta), set azimuth angle phi
1614   //
1615   const G4double smallValue = 1.0e-10;
1616   G4double dumnve = 4.0 * pin * pf;
1617   if (dumnve == 0.0)
1618     dumnve = smallValue;
1619   G4double ctet = std::max(-1.0, std::min(1.0, 1.0 + 2.0 * (t - tacmin) / dumnve));
1620   dumnve = std::max(0.0, 1.0 - ctet * ctet);
1621   G4double stet = std::sqrt(dumnve);
1622   G4double phi = G4UniformRand() * twopi;
1623   //
1624   // calculate final state momenta in centre of mass system
1625   //
1626   pseudoParticle[3].SetMomentum(pf * stet * std::sin(phi) * GeV, pf * stet * std::cos(phi) * GeV, pf * ctet * GeV);
1627   pseudoParticle[4].SetMomentum(pseudoParticle[3].GetMomentum() * (-1.0));
1628   //
1629   // simulate backward nucleon cluster in lab. system and transform in cms
1630   //
1631   G4double pp, pp1, rthnve, phinve;
1632   if (nuclearExcitationCount > 0) {
1633     const G4double ga = 1.2;
1634     G4double ekit1 = 0.04;
1635     G4double ekit2 = 0.6;
1636     if (ekOriginal <= 5.0) {
1637       ekit1 *= ekOriginal * ekOriginal / 25.0;
1638       ekit2 *= ekOriginal * ekOriginal / 25.0;
1639     }
1640     const G4double a = (1.0 - ga) / (std::pow(ekit2, (1.0 - ga)) - std::pow(ekit1, (1.0 - ga)));
1641     for (i = 0; i < vecLen; ++i) {
1642       if (vec[i]->GetSide() == -2) {
1643         G4double kineticE =
1644             std::pow((G4UniformRand() * (1.0 - ga) / a + std::pow(ekit1, (1.0 - ga))), (1.0 / (1.0 - ga)));
1645         vec[i]->SetKineticEnergy(kineticE * GeV);
1646         G4double vMass = vec[i]->GetMass() / MeV;
1647         G4double totalE = kineticE + vMass;
1648         pp = std::sqrt(std::abs(totalE * totalE - vMass * vMass));
1649         G4double cost = std::min(1.0, std::max(-1.0, std::log(2.23 * G4UniformRand() + 0.383) / 0.96));
1650         G4double sint = std::sqrt(std::max(0.0, (1.0 - cost * cost)));
1651         phi = twopi * G4UniformRand();
1652         vec[i]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
1653         vec[i]->Lorentz(*vec[i], pseudoParticle[0]);
1654       }
1655     }
1656     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1657   }
1658   //
1659   // fragmentation of forward cluster and backward meson cluster
1660   //
1661   currentParticle.SetMomentum(pseudoParticle[3].GetMomentum());
1662   currentParticle.SetTotalEnergy(pseudoParticle[3].GetTotalEnergy());
1663 
1664   targetParticle.SetMomentum(pseudoParticle[4].GetMomentum());
1665   targetParticle.SetTotalEnergy(pseudoParticle[4].GetTotalEnergy());
1666 
1667   pseudoParticle[5].SetMomentum(pseudoParticle[3].GetMomentum() * (-1.0));
1668   pseudoParticle[5].SetMass(pseudoParticle[3].GetMass());
1669   pseudoParticle[5].SetTotalEnergy(pseudoParticle[3].GetTotalEnergy());
1670 
1671   pseudoParticle[6].SetMomentum(pseudoParticle[4].GetMomentum() * (-1.0));
1672   pseudoParticle[6].SetMass(pseudoParticle[4].GetMass());
1673   pseudoParticle[6].SetTotalEnergy(pseudoParticle[4].GetTotalEnergy());
1674 
1675   G4double wgt;
1676   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1677   if (forwardCount > 1)  // tempV will contain the forward particles
1678   {
1679     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1680     tempV.Initialize(forwardCount);
1681     G4bool constantCrossSection = true;
1682     G4int tempLen = 0;
1683     if (currentParticle.GetSide() == 1)
1684       tempV.SetElement(tempLen++, &currentParticle);
1685     if (targetParticle.GetSide() == 1)
1686       tempV.SetElement(tempLen++, &targetParticle);
1687     for (i = 0; i < vecLen; ++i) {
1688       if (vec[i]->GetSide() == 1) {
1689         if (tempLen < 18)
1690           tempV.SetElement(tempLen++, vec[i]);
1691         else {
1692           vec[i]->SetSide(-1);
1693           continue;
1694         }
1695       }
1696     }
1697     if (tempLen >= 2) {
1698       GenerateNBodyEvent(pseudoParticle[3].GetMass() / MeV, constantCrossSection, tempV, tempLen);
1699       if (currentParticle.GetSide() == 1)
1700         currentParticle.Lorentz(currentParticle, pseudoParticle[5]);
1701       if (targetParticle.GetSide() == 1)
1702         targetParticle.Lorentz(targetParticle, pseudoParticle[5]);
1703       for (i = 0; i < vecLen; ++i) {
1704         if (vec[i]->GetSide() == 1)
1705           vec[i]->Lorentz(*vec[i], pseudoParticle[5]);
1706       }
1707     }
1708   }
1709   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1710   if (backwardCount > 1)  //  tempV will contain the backward particles,
1711   {                       //  but not those created from the intranuclear cascade
1712     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1713     tempV.Initialize(backwardCount);
1714     G4bool constantCrossSection = true;
1715     G4int tempLen = 0;
1716     if (currentParticle.GetSide() == -1)
1717       tempV.SetElement(tempLen++, &currentParticle);
1718     if (targetParticle.GetSide() == -1)
1719       tempV.SetElement(tempLen++, &targetParticle);
1720     for (i = 0; i < vecLen; ++i) {
1721       if (vec[i]->GetSide() == -1) {
1722         if (tempLen < 18)
1723           tempV.SetElement(tempLen++, vec[i]);
1724         else {
1725           vec[i]->SetSide(-2);
1726           vec[i]->SetKineticEnergy(0.0);
1727           vec[i]->SetMomentum(0.0, 0.0, 0.0);
1728           continue;
1729         }
1730       }
1731     }
1732     if (tempLen >= 2) {
1733       GenerateNBodyEvent(pseudoParticle[4].GetMass() / MeV, constantCrossSection, tempV, tempLen);
1734       if (currentParticle.GetSide() == -1)
1735         currentParticle.Lorentz(currentParticle, pseudoParticle[6]);
1736       if (targetParticle.GetSide() == -1)
1737         targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
1738       for (i = 0; i < vecLen; ++i) {
1739         if (vec[i]->GetSide() == -1)
1740           vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
1741       }
1742     }
1743   }
1744   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1745   //
1746   // Lorentz transformation in lab system
1747   //
1748   G4int numberofFinalStateNucleons = 0;
1749   if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron ||
1750       currentParticle.GetDefinition() == aSigmaMinus || currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
1751       currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
1752       currentParticle.GetDefinition() == G4XiZero::XiZero() ||
1753       currentParticle.GetDefinition() == G4XiMinus::XiMinus() ||
1754       currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1755       currentParticle.GetDefinition() == G4Lambda::Lambda())
1756     ++numberofFinalStateNucleons;
1757   currentParticle.Lorentz(currentParticle, pseudoParticle[2]);
1758   if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron ||
1759       targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() ||
1760       targetParticle.GetDefinition() == G4XiMinus::XiMinus() ||
1761       targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1762       targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
1763       targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() || targetParticle.GetDefinition() == aSigmaMinus)
1764     ++numberofFinalStateNucleons;
1765   if (targetParticle.GetDefinition() == G4AntiProton::AntiProton())
1766     --numberofFinalStateNucleons;
1767   if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron())
1768     --numberofFinalStateNucleons;
1769   if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
1770     --numberofFinalStateNucleons;
1771   if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
1772     --numberofFinalStateNucleons;
1773   if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
1774     --numberofFinalStateNucleons;
1775   if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero())
1776     --numberofFinalStateNucleons;
1777   if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus())
1778     --numberofFinalStateNucleons;
1779   if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
1780     --numberofFinalStateNucleons;
1781   if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda())
1782     --numberofFinalStateNucleons;
1783   targetParticle.Lorentz(targetParticle, pseudoParticle[2]);
1784   for (i = 0; i < vecLen; ++i) {
1785     if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron ||
1786         vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() ||
1787         vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1788         vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() ||
1789         vec[i]->GetDefinition() == aSigmaMinus)
1790       ++numberofFinalStateNucleons;
1791     if (vec[i]->GetDefinition() == G4AntiProton::AntiProton())
1792       --numberofFinalStateNucleons;
1793     if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron())
1794       --numberofFinalStateNucleons;
1795     if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
1796       --numberofFinalStateNucleons;
1797     if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
1798       --numberofFinalStateNucleons;
1799     if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
1800       --numberofFinalStateNucleons;
1801     if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda())
1802       --numberofFinalStateNucleons;
1803     if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero())
1804       --numberofFinalStateNucleons;
1805     if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus())
1806       --numberofFinalStateNucleons;
1807     if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
1808       --numberofFinalStateNucleons;
1809     vec[i]->Lorentz(*vec[i], pseudoParticle[2]);
1810   }
1811   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1812   numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
1813   //
1814   // sometimes the leading strange particle is lost, set it back
1815   //
1816   G4bool dum = true;
1817   if (leadFlag) {
1818     // leadFlag will be true
1819     //  iff original particle is at least as heavy as K+ and not a proton or neutron AND
1820     //   if
1821     //    incident particle is at least as heavy as K+ and it is not a proton or neutron
1822     //     leadFlag is set to the incident particle
1823     //   or
1824     //    target particle is at least as heavy as K+ and it is not a proton or neutron
1825     //     leadFlag is set to the target particle
1826     //
1827     if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1828       dum = false;
1829     else if (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1830       dum = false;
1831     else {
1832       for (i = 0; i < vecLen; ++i) {
1833         if (vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition()) {
1834           dum = false;
1835           break;
1836         }
1837       }
1838     }
1839     if (dum) {
1840       G4double leadMass = leadingStrangeParticle.GetMass() / MeV;
1841       G4double ekin;
1842       if (((leadMass < protonMass) && (targetParticle.GetMass() / MeV < protonMass)) ||
1843           ((leadMass >= protonMass) && (targetParticle.GetMass() / MeV >= protonMass))) {
1844         ekin = targetParticle.GetKineticEnergy() / GeV;
1845         pp1 = targetParticle.GetMomentum().mag() / MeV;  // old momentum
1846         targetParticle.SetDefinition(leadingStrangeParticle.GetDefinition());
1847         targetParticle.SetKineticEnergy(ekin * GeV);
1848         pp = targetParticle.GetTotalMomentum() / MeV;  // new momentum
1849         if (pp1 < 1.0e-3) {
1850           rthnve = pi * G4UniformRand();
1851           phinve = twopi * G4UniformRand();
1852           targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1853                                      pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1854                                      pp * std::cos(rthnve) * MeV);
1855         } else
1856           targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1857 
1858         targetHasChanged = true;
1859       } else {
1860         ekin = currentParticle.GetKineticEnergy() / GeV;
1861         pp1 = currentParticle.GetMomentum().mag() / MeV;
1862         currentParticle.SetDefinition(leadingStrangeParticle.GetDefinition());
1863         currentParticle.SetKineticEnergy(ekin * GeV);
1864         pp = currentParticle.GetTotalMomentum() / MeV;
1865         if (pp1 < 1.0e-3) {
1866           rthnve = pi * G4UniformRand();
1867           phinve = twopi * G4UniformRand();
1868           currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1869                                       pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1870                                       pp * std::cos(rthnve) * MeV);
1871         } else
1872           currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1873 
1874         incidentHasChanged = true;
1875       }
1876     }
1877   }  // end of if( leadFlag )
1878   //
1879   //  for various reasons, the energy balance is not sufficient,
1880   //  check that,  energy balance, angle of final system, etc.
1881   //
1882   pseudoParticle[4].SetMass(mOriginal * GeV);
1883   pseudoParticle[4].SetTotalEnergy(etOriginal * GeV);
1884   pseudoParticle[4].SetMomentum(0.0, 0.0, pOriginal * GeV);
1885 
1886   const G4ParticleDefinition *aOrgDef = modifiedOriginal.GetDefinition();
1887   G4int diff = 0;
1888   if (aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron())
1889     diff = 1;
1890   if (numberofFinalStateNucleons == 1)
1891     diff = 0;
1892   pseudoParticle[5].SetMomentum(0.0, 0.0, 0.0);
1893   pseudoParticle[5].SetMass(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1894   pseudoParticle[5].SetTotalEnergy(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1895 
1896   //    G4double ekin0 = pseudoParticle[4].GetKineticEnergy()/GeV;
1897   G4double theoreticalKinetic = pseudoParticle[4].GetTotalEnergy() / GeV + pseudoParticle[5].GetTotalEnergy() / GeV;
1898 
1899   pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
1900   pseudoParticle[4].Lorentz(pseudoParticle[4], pseudoParticle[6]);
1901   pseudoParticle[5].Lorentz(pseudoParticle[5], pseudoParticle[6]);
1902 
1903   if (vecLen < 16) {
1904     G4ReactionProduct tempR[130];
1905     //G4ReactionProduct *tempR = new G4ReactionProduct [vecLen+2];
1906     tempR[0] = currentParticle;
1907     tempR[1] = targetParticle;
1908     for (i = 0; i < vecLen; ++i)
1909       tempR[i + 2] = *vec[i];
1910 
1911     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1912     tempV.Initialize(vecLen + 2);
1913     G4bool constantCrossSection = true;
1914     G4int tempLen = 0;
1915     for (i = 0; i < vecLen + 2; ++i)
1916       tempV.SetElement(tempLen++, &tempR[i]);
1917 
1918     if (tempLen >= 2) {
1919       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1920       GenerateNBodyEvent(pseudoParticle[4].GetTotalEnergy() / MeV + pseudoParticle[5].GetTotalEnergy() / MeV,
1921                          constantCrossSection,
1922                          tempV,
1923                          tempLen);
1924       theoreticalKinetic = 0.0;
1925       for (i = 0; i < vecLen + 2; ++i) {
1926         pseudoParticle[7].SetMomentum(tempV[i]->GetMomentum());
1927         pseudoParticle[7].SetMass(tempV[i]->GetMass());
1928         pseudoParticle[7].SetTotalEnergy(tempV[i]->GetTotalEnergy());
1929         pseudoParticle[7].Lorentz(pseudoParticle[7], pseudoParticle[5]);
1930         theoreticalKinetic += pseudoParticle[7].GetKineticEnergy() / GeV;
1931       }
1932     }
1933     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1934     //delete [] tempR;
1935   } else {
1936     theoreticalKinetic -= (currentParticle.GetMass() / GeV + targetParticle.GetMass() / GeV);
1937     for (i = 0; i < vecLen; ++i)
1938       theoreticalKinetic -= vec[i]->GetMass() / GeV;
1939   }
1940   G4double simulatedKinetic = currentParticle.GetKineticEnergy() / GeV + targetParticle.GetKineticEnergy() / GeV;
1941   for (i = 0; i < vecLen; ++i)
1942     simulatedKinetic += vec[i]->GetKineticEnergy() / GeV;
1943   //
1944   // make sure that kinetic energies are correct
1945   // the backward nucleon cluster is not produced within proper kinematics!!!
1946   //
1947 
1948   if (simulatedKinetic != 0.0) {
1949     wgt = (theoreticalKinetic) / simulatedKinetic;
1950     currentParticle.SetKineticEnergy(wgt * currentParticle.GetKineticEnergy());
1951     pp = currentParticle.GetTotalMomentum() / MeV;
1952     pp1 = currentParticle.GetMomentum().mag() / MeV;
1953     if (pp1 < 0.001 * MeV) {
1954       rthnve = pi * G4UniformRand();
1955       phinve = twopi * G4UniformRand();
1956       currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1957                                   pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1958                                   pp * std::cos(rthnve) * MeV);
1959     } else
1960       currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1961 
1962     targetParticle.SetKineticEnergy(wgt * targetParticle.GetKineticEnergy());
1963     pp = targetParticle.GetTotalMomentum() / MeV;
1964     pp1 = targetParticle.GetMomentum().mag() / MeV;
1965     if (pp1 < 0.001 * MeV) {
1966       rthnve = pi * G4UniformRand();
1967       phinve = twopi * G4UniformRand();
1968       targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1969                                  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1970                                  pp * std::cos(rthnve) * MeV);
1971     } else
1972       targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1973 
1974     for (i = 0; i < vecLen; ++i) {
1975       vec[i]->SetKineticEnergy(wgt * vec[i]->GetKineticEnergy());
1976       pp = vec[i]->GetTotalMomentum() / MeV;
1977       pp1 = vec[i]->GetMomentum().mag() / MeV;
1978       if (pp1 < 0.001) {
1979         rthnve = pi * G4UniformRand();
1980         phinve = twopi * G4UniformRand();
1981         vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1982                             pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1983                             pp * std::cos(rthnve) * MeV);
1984       } else
1985         vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
1986     }
1987   }
1988   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1989   Rotate(numberofFinalStateNucleons,
1990          pseudoParticle[4].GetMomentum(),
1991          modifiedOriginal,
1992          originalIncident,
1993          targetNucleus,
1994          currentParticle,
1995          targetParticle,
1996          vec,
1997          vecLen);
1998   //
1999   //  add black track particles
2000   //  the total number of particles produced is restricted to 198
2001   //  this may have influence on very high energies
2002   //
2003   if (atomicWeight >= 1.5) {
2004     // npnb is number of proton/neutron black track particles
2005     // ndta is the number of deuterons, tritons, and alphas produced
2006     // epnb is the kinetic energy available for proton/neutron black track particles
2007     // edta is the kinetic energy available for deuteron/triton/alpha particles
2008     //
2009     G4double epnb, edta;
2010     G4int npnb = 0;
2011     G4int ndta = 0;
2012 
2013     epnb = targetNucleus.GetPNBlackTrackEnergy();   // was enp1 in fortran code
2014     edta = targetNucleus.GetDTABlackTrackEnergy();  // was enp3 in fortran code
2015     const G4double pnCutOff = 0.001;                // GeV
2016     const G4double dtaCutOff = 0.001;               // GeV
2017     const G4double kineticMinimum = 1.e-6;
2018     const G4double kineticFactor = -0.005;
2019 
2020     G4double sprob = 0.0;  // sprob = probability of self-absorption in heavy molecules
2021     const G4double ekIncident = originalIncident->GetKineticEnergy() / GeV;
2022     if (ekIncident >= 5.0)
2023       sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0));
2024 
2025     if (epnb >= pnCutOff) {
2026       npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta));
2027       if (numberofFinalStateNucleons + npnb > atomicWeight)
2028         npnb = G4int(atomicWeight - numberofFinalStateNucleons);
2029       npnb = std::min(npnb, 127 - vecLen);
2030     }
2031     if (edta >= dtaCutOff) {
2032       ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta));
2033       ndta = std::min(ndta, 127 - vecLen);
2034     }
2035     G4double spall = numberofFinalStateNucleons;
2036     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2037 
2038     AddBlackTrackParticles(epnb,
2039                            npnb,
2040                            edta,
2041                            ndta,
2042                            sprob,
2043                            kineticMinimum,
2044                            kineticFactor,
2045                            modifiedOriginal,
2046                            spall,
2047                            targetNucleus,
2048                            vec,
2049                            vecLen);
2050 
2051     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2052   }
2053   //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
2054   //  MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2055   //
2056   //  calculate time delay for nuclear reactions
2057   //
2058   if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
2059     currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
2060   else
2061     currentParticle.SetTOF(1.0);
2062 
2063   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2064   return true;
2065 }
2066 
2067 void FullModelReactionDynamics::TwoBody(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2068                                         G4int &vecLen,
2069                                         G4ReactionProduct &modifiedOriginal,
2070                                         const G4DynamicParticle * /*originalTarget*/,
2071                                         G4ReactionProduct &currentParticle,
2072                                         G4ReactionProduct &targetParticle,
2073                                         const G4Nucleus &targetNucleus,
2074                                         G4bool & /* targetHasChanged*/) {
2075   //    G4cout<<"TwoBody called"<<G4endl;
2076   //
2077   // derived from original FORTRAN code TWOB by H. Fesefeldt (15-Sep-1987)
2078   //
2079   // Generation of momenta for elastic and quasi-elastic 2 body reactions
2080   //
2081   // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used.
2082   // The b values are parametrizations from experimental data.
2083   // Not available values are taken from those of similar reactions.
2084   //
2085   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2086   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2087   G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2088   G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
2089   G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
2090   G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
2091   G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
2092 
2093   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2094   static const G4double expxu = 82.;     // upper bound for arg. of exp
2095   static const G4double expxl = -expxu;  // lower bound for arg. of exp
2096 
2097   const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
2098   const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
2099   G4double currentMass = currentParticle.GetMass() / GeV;
2100   G4double targetMass = targetParticle.GetMass() / GeV;
2101   const G4double atomicWeight = targetNucleus.GetN_asInt();
2102   //    G4cout<<"Atomic weight is found to be: "<<atomicWeight<<G4endl;
2103   G4double etCurrent = currentParticle.GetTotalEnergy() / GeV;
2104   G4double pCurrent = currentParticle.GetTotalMomentum() / GeV;
2105 
2106   G4double cmEnergy =
2107       std::sqrt(currentMass * currentMass + targetMass * targetMass + 2.0 * targetMass * etCurrent);  // in GeV
2108 
2109   //if( (pOriginal < 0.1) ||
2110   //    (centerofmassEnergy < 0.01) ) // 2-body scattering not possible
2111   // Continue with original particle, but spend the nuclear evaporation energy
2112   //  targetParticle.SetMass( 0.0 );  // flag that the target doesn't exist
2113   //else                           // Two-body scattering is possible
2114 
2115   if ((pCurrent < 0.1) || (cmEnergy < 0.01))  // 2-body scattering not possible
2116   {
2117     targetParticle.SetMass(0.0);  // flag that the target particle doesn't exist
2118   } else {
2119     // moved this if-block to a later stage, i.e. to the assignment of the scattering angle
2120     // @@@@@ double-check.
2121     //      if( targetParticle.GetDefinition() == aKaonMinus ||
2122     //          targetParticle.GetDefinition() == aKaonZeroL ||
2123     //          targetParticle.GetDefinition() == aKaonZeroS ||
2124     //          targetParticle.GetDefinition() == aKaonPlus ||
2125     //          targetParticle.GetDefinition() == aPiMinus ||
2126     //          targetParticle.GetDefinition() == aPiZero ||
2127     //          targetParticle.GetDefinition() == aPiPlus )
2128     //      {
2129     //        if( G4UniformRand() < 0.5 )
2130     //          targetParticle.SetDefinitionAndUpdateE( aNeutron );
2131     //        else
2132     //          targetParticle.SetDefinitionAndUpdateE( aProton );
2133     //        targetHasChanged = true;
2134     //        targetMass = targetParticle.GetMass()/GeV;
2135     //      }
2136     //
2137     // Set masses and momenta for final state particles
2138     //
2139     /*
2140       G4cout<<"Check 0"<<G4endl;
2141       G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()/GeV<<G4endl;
2142       G4cout<<"target mass: "<<targetParticle.GetMass()/GeV<<G4endl;
2143       G4cout<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
2144       G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()/GeV<<G4endl;
2145       G4cout<<"current mass: "<<currentParticle.GetMass()/GeV<<G4endl;
2146       G4cout<<currentParticle.GetDefinition()->GetParticleName()<<G4endl;
2147       */
2148     G4double pf = cmEnergy * cmEnergy + targetMass * targetMass - currentMass * currentMass;
2149     pf = pf * pf - 4 * cmEnergy * cmEnergy * targetMass * targetMass;
2150     //      G4cout << "pf: " << pf<< G4endl;
2151 
2152     if (pf <= 0.)  // 0.001 )
2153     {
2154       for (G4int i = 0; i < vecLen; i++)
2155         delete vec[i];
2156       vecLen = 0;
2157       throw G4HadronicException(__FILE__, __LINE__, "FullModelReactionDynamics::TwoBody: pf is too small ");
2158     }
2159 
2160     pf = std::sqrt(pf) / (2.0 * cmEnergy);
2161     //
2162     // Set beam and target in centre of mass system
2163     //
2164     G4ReactionProduct pseudoParticle[3];
2165     pseudoParticle[0].SetMass(currentMass * GeV);
2166     pseudoParticle[0].SetTotalEnergy(etCurrent * GeV);
2167     pseudoParticle[0].SetMomentum(0.0, 0.0, pCurrent * GeV);
2168 
2169     pseudoParticle[1].SetMomentum(0.0, 0.0, 0.0);
2170     pseudoParticle[1].SetMass(targetMass * GeV);
2171     pseudoParticle[1].SetKineticEnergy(0.0);
2172     //
2173     // Transform into centre of mass system
2174     //
2175     pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
2176     pseudoParticle[0].Lorentz(pseudoParticle[0], pseudoParticle[2]);
2177     pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[2]);
2178     //
2179     // Set final state masses and energies in centre of mass system
2180     //
2181     currentParticle.SetTotalEnergy(std::sqrt(pf * pf + currentMass * currentMass) * GeV);
2182     targetParticle.SetTotalEnergy(std::sqrt(pf * pf + targetMass * targetMass) * GeV);
2183     //
2184     // Set |t| and |tmin|
2185     //
2186     const G4double cb = 0.01;
2187     const G4double b1 = 4.225;
2188     const G4double b2 = 1.795;
2189     //
2190     // Calculate slope b for elastic scattering on proton/neutron
2191     //
2192     G4double b = std::max(cb, b1 + b2 * std::log(pOriginal));
2193     //      G4double b = std::max( cb, b1+b2*std::log(ptemp) );
2194     G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag() / GeV;
2195 
2196     G4double exindt = -1.0;
2197     exindt += std::exp(std::max(-btrang, expxl));
2198     //
2199     // Calculate sqr(std::sin(teta/2.) and std::cos(teta), set azimuth angle phi
2200     //
2201     G4double ctet = 1.0 + 2 * std::log(1.0 + G4UniformRand() * exindt) / btrang;
2202     if (std::fabs(ctet) > 1.0)
2203       ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
2204     G4double stet = std::sqrt((1.0 - ctet) * (1.0 + ctet));
2205     G4double phi = twopi * G4UniformRand();
2206     //
2207     // Calculate final state momenta in centre of mass system
2208     //
2209     if (targetParticle.GetDefinition() == aKaonMinus || targetParticle.GetDefinition() == aKaonZeroL ||
2210         targetParticle.GetDefinition() == aKaonZeroS || targetParticle.GetDefinition() == aKaonPlus ||
2211         targetParticle.GetDefinition() == aPiMinus || targetParticle.GetDefinition() == aPiZero ||
2212         targetParticle.GetDefinition() == aPiPlus) {
2213       currentParticle.SetMomentum(-pf * stet * std::sin(phi) * GeV, -pf * stet * std::cos(phi) * GeV, -pf * ctet * GeV);
2214     } else {
2215       currentParticle.SetMomentum(pf * stet * std::sin(phi) * GeV, pf * stet * std::cos(phi) * GeV, pf * ctet * GeV);
2216     }
2217     targetParticle.SetMomentum(currentParticle.GetMomentum() * (-1.0));
2218     //
2219     // Transform into lab system
2220     //
2221     currentParticle.Lorentz(currentParticle, pseudoParticle[1]);
2222     targetParticle.Lorentz(targetParticle, pseudoParticle[1]);
2223 
2224     /*
2225       G4cout<<"Check 1"<<G4endl;
2226       G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
2227       G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
2228       G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
2229       G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
2230       */
2231 
2232     Defs1(modifiedOriginal, currentParticle, targetParticle, vec, vecLen);
2233 
2234     G4double pp, pp1, ekin;
2235     if (atomicWeight >= 1.5) {
2236       const G4double cfa = 0.025 * ((atomicWeight - 1.) / 120.) * std::exp(-(atomicWeight - 1.) / 120.);
2237       pp1 = currentParticle.GetMomentum().mag() / MeV;
2238       if (pp1 >= 1.0) {
2239         ekin = currentParticle.GetKineticEnergy() / MeV - cfa * (1.0 + 0.5 * normal()) * GeV;
2240         ekin = std::max(0.0001 * GeV, ekin);
2241         currentParticle.SetKineticEnergy(ekin * MeV);
2242         pp = currentParticle.GetTotalMomentum() / MeV;
2243         currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2244       }
2245       pp1 = targetParticle.GetMomentum().mag() / MeV;
2246       if (pp1 >= 1.0) {
2247         ekin = targetParticle.GetKineticEnergy() / MeV - cfa * (1.0 + normal() / 2.) * GeV;
2248         ekin = std::max(0.0001 * GeV, ekin);
2249         targetParticle.SetKineticEnergy(ekin * MeV);
2250         pp = targetParticle.GetTotalMomentum() / MeV;
2251         targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
2252       }
2253     }
2254   }
2255   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2256   if (atomicWeight >= 1.5) {
2257     // Add black track particles
2258     //  the procedure is somewhat different than in TwoCluster and GenerateXandPt.
2259     //  The reason is that we have to also simulate the nuclear reactions
2260     //  at low energies like a(h,p)b, a(h,p p)b, a(h,n)b etc.
2261     //
2262     // npnb is number of proton/neutron black track particles
2263     // ndta is the number of deuterons, tritons, and alphas produced
2264     // epnb is the kinetic energy available for proton/neutron black track particles
2265     // edta is the kinetic energy available for deuteron/triton/alpha particles
2266     //
2267     G4double epnb, edta;
2268     G4int npnb = 0, ndta = 0;
2269 
2270     epnb = targetNucleus.GetPNBlackTrackEnergy();   // was enp1 in fortran code
2271     edta = targetNucleus.GetDTABlackTrackEnergy();  // was enp3 in fortran code
2272     const G4double pnCutOff = 0.0001;               // GeV
2273     const G4double dtaCutOff = 0.0001;              // GeV
2274     const G4double kineticMinimum = 0.0001;
2275     const G4double kineticFactor = -0.010;
2276     G4double sprob = 0.0;  // sprob = probability of self-absorption in heavy molecules
2277     if (epnb >= pnCutOff) {
2278       npnb = Poisson(epnb / 0.02);
2279       /*
2280     G4cout<<"A couple of Poisson numbers:"<<G4endl;
2281     for (int n=0;n!=10;n++) G4cout<<Poisson(epnb/0.02)<<", ";
2282     G4cout<<G4endl;
2283     */
2284       if (npnb > atomicWeight)
2285         npnb = G4int(atomicWeight);
2286       if ((epnb > pnCutOff) && (npnb <= 0))
2287         npnb = 1;
2288       npnb = std::min(npnb, 127 - vecLen);
2289     }
2290     if (edta >= dtaCutOff) {
2291       ndta = G4int(2.0 * std::log(atomicWeight));
2292       ndta = std::min(ndta, 127 - vecLen);
2293     }
2294     G4double spall = 0.0;
2295     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2296 
2297     /*
2298       G4cout<<"Check 2"<<G4endl;
2299       G4cout<<"target E_kin: "<<targetParticle.GetKineticEnergy()<<G4endl;
2300       G4cout<<"target mass: "<<targetParticle.GetMass()<<G4endl;
2301       G4cout<<"current E_kin: "<<currentParticle.GetKineticEnergy()<<G4endl;
2302       G4cout<<"current mass: "<<currentParticle.GetMass()<<G4endl;
2303 
2304       G4cout<<"------------------------------------------------------------------------"<<G4endl;
2305       G4cout<<"Atomic weight: "<<atomicWeight<<G4endl;
2306       G4cout<<"number of proton/neutron black track particles: "<<npnb<<G4endl;
2307       G4cout<<"number of deuterons, tritons, and alphas produced: "<<ndta <<G4endl;
2308       G4cout<<"kinetic energy available for proton/neutron black track particles: "<<epnb/GeV<<" GeV" <<G4endl;
2309       G4cout<<"kinetic energy available for deuteron/triton/alpha particles: "<<edta/GeV <<" GeV"<<G4endl;
2310       G4cout<<"------------------------------------------------------------------------"<<G4endl;
2311       */
2312 
2313     AddBlackTrackParticles(epnb,
2314                            npnb,
2315                            edta,
2316                            ndta,
2317                            sprob,
2318                            kineticMinimum,
2319                            kineticFactor,
2320                            modifiedOriginal,
2321                            spall,
2322                            targetNucleus,
2323                            vec,
2324                            vecLen);
2325 
2326     // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2327   }
2328   //
2329   //  calculate time delay for nuclear reactions
2330   //
2331   if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
2332     currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
2333   else
2334     currentParticle.SetTOF(1.0);
2335   return;
2336 }
2337 
2338 G4double FullModelReactionDynamics::GenerateNBodyEvent(const G4double totalEnergy,  // MeV
2339                                                        const G4bool constantCrossSection,
2340                                                        G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2341                                                        G4int &vecLen) {
2342   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2343   // derived from original FORTRAN code PHASP by H. Fesefeldt (02-Dec-1986)
2344   // Returns the weight of the event
2345   //
2346   G4int i;
2347   const G4double expxu = 82.;     // upper bound for arg. of exp
2348   const G4double expxl = -expxu;  // lower bound for arg. of exp
2349   if (vecLen < 2) {
2350     G4cerr << "*** Error in FullModelReactionDynamics::GenerateNBodyEvent" << G4endl;
2351     G4cerr << "    number of particles < 2" << G4endl;
2352     G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
2353     return -1.0;
2354   }
2355   G4double mass[18];    // mass of each particle
2356   G4double energy[18];  // total energy of each particle
2357   G4double pcm[3][18];  // pcm is an array with 3 rows and vecLen columns
2358   G4double totalMass = 0.0;
2359   G4double extraMass = 0;
2360   G4double sm[18];
2361 
2362   for (i = 0; i < vecLen; ++i) {
2363     mass[i] = vec[i]->GetMass() / GeV;
2364     if (vec[i]->GetSide() == -2)
2365       extraMass += vec[i]->GetMass() / GeV;
2366     vec[i]->SetMomentum(0.0, 0.0, 0.0);
2367     pcm[0][i] = 0.0;      // x-momentum of i-th particle
2368     pcm[1][i] = 0.0;      // y-momentum of i-th particle
2369     pcm[2][i] = 0.0;      // z-momentum of i-th particle
2370     energy[i] = mass[i];  // total energy of i-th particle
2371     totalMass += mass[i];
2372     sm[i] = totalMass;
2373   }
2374   G4double totalE = totalEnergy / GeV;
2375   if (totalMass > totalE) {
2376     return -1.0;
2377   }
2378   G4double kineticEnergy = totalE - totalMass;
2379   G4double emm[18];
2380   //G4double *emm = new G4double [vecLen];
2381   emm[0] = mass[0];
2382   emm[vecLen - 1] = totalE;
2383   if (vecLen > 2)  // the random numbers are sorted
2384   {
2385     G4double ran[18];
2386     for (i = 0; i < vecLen; ++i)
2387       ran[i] = G4UniformRand();
2388     for (i = 0; i < vecLen - 2; ++i) {
2389       for (G4int j = vecLen - 2; j > i; --j) {
2390         if (ran[i] > ran[j]) {
2391           G4double temp = ran[i];
2392           ran[i] = ran[j];
2393           ran[j] = temp;
2394         }
2395       }
2396     }
2397     for (i = 1; i < vecLen - 1; ++i)
2398       emm[i] = ran[i - 1] * kineticEnergy + sm[i];
2399   }
2400   //   Weight is the sum of logarithms of terms instead of the product of terms
2401   G4bool lzero = true;
2402   G4double wtmax = 0.0;
2403   if (constantCrossSection)  // this is KGENEV=1 in PHASP
2404   {
2405     G4double emmax = kineticEnergy + mass[0];
2406     G4double emmin = 0.0;
2407     for (i = 1; i < vecLen; ++i) {
2408       emmin += mass[i - 1];
2409       emmax += mass[i];
2410       G4double wtfc = 0.0;
2411       if (emmax * emmax > 0.0) {
2412         G4double arg = emmax * emmax +
2413                        (emmin * emmin - mass[i] * mass[i]) * (emmin * emmin - mass[i] * mass[i]) / (emmax * emmax) -
2414                        2.0 * (emmin * emmin + mass[i] * mass[i]);
2415         if (arg > 0.0)
2416           wtfc = 0.5 * std::sqrt(arg);
2417       }
2418       if (wtfc == 0.0) {
2419         lzero = false;
2420         break;
2421       }
2422       wtmax += std::log(wtfc);
2423     }
2424     if (lzero)
2425       wtmax = -wtmax;
2426     else
2427       wtmax = expxu;
2428   } else {
2429     //   ffq(n) = pi*(2*pi)^(n-2)/(n-2)!
2430     const G4double ffq[18] = {0.,
2431                               3.141592,
2432                               19.73921,
2433                               62.01255,
2434                               129.8788,
2435                               204.0131,
2436                               256.3704,
2437                               268.4705,
2438                               240.9780,
2439                               189.2637,
2440                               132.1308,
2441                               83.0202,
2442                               47.4210,
2443                               24.8295,
2444                               12.0006,
2445                               5.3858,
2446                               2.2560,
2447                               0.8859};
2448     wtmax = std::log(std::pow(kineticEnergy, vecLen - 2) * ffq[vecLen - 1] / totalE);
2449   }
2450   lzero = true;
2451   G4double pd[50];
2452   //G4double *pd = new G4double [vecLen-1];
2453   for (i = 0; i < vecLen - 1; ++i) {
2454     pd[i] = 0.0;
2455     if (emm[i + 1] * emm[i + 1] > 0.0) {
2456       G4double arg = emm[i + 1] * emm[i + 1] +
2457                      (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) * (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) /
2458                          (emm[i + 1] * emm[i + 1]) -
2459                      2.0 * (emm[i] * emm[i] + mass[i + 1] * mass[i + 1]);
2460       if (arg > 0.0)
2461         pd[i] = 0.5 * std::sqrt(arg);
2462     }
2463     if (pd[i] <= 0.0)  //  changed from  ==  on 02 April 98
2464       lzero = false;
2465     else
2466       wtmax += std::log(pd[i]);
2467   }
2468   G4double weight = 0.0;  // weight is returned by GenerateNBodyEvent
2469   if (lzero)
2470     weight = std::exp(std::max(std::min(wtmax, expxu), expxl));
2471 
2472   G4double bang, cb, sb, s0, s1, s2, c, s, esys, a, b, gama, beta;
2473   pcm[0][0] = 0.0;
2474   pcm[1][0] = pd[0];
2475   pcm[2][0] = 0.0;
2476   for (i = 1; i < vecLen; ++i) {
2477     pcm[0][i] = 0.0;
2478     pcm[1][i] = -pd[i - 1];
2479     pcm[2][i] = 0.0;
2480     bang = twopi * G4UniformRand();
2481     cb = std::cos(bang);
2482     sb = std::sin(bang);
2483     c = 2.0 * G4UniformRand() - 1.0;
2484     s = std::sqrt(std::fabs(1.0 - c * c));
2485     if (i < vecLen - 1) {
2486       esys = std::sqrt(pd[i] * pd[i] + emm[i] * emm[i]);
2487       beta = pd[i] / esys;
2488       gama = esys / emm[i];
2489       for (G4int j = 0; j <= i; ++j) {
2490         s0 = pcm[0][j];
2491         s1 = pcm[1][j];
2492         s2 = pcm[2][j];
2493         energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2494         a = s0 * c - s1 * s;  //  rotation
2495         pcm[1][j] = s0 * s + s1 * c;
2496         b = pcm[2][j];
2497         pcm[0][j] = a * cb - b * sb;
2498         pcm[2][j] = a * sb + b * cb;
2499         pcm[1][j] = gama * (pcm[1][j] + beta * energy[j]);
2500       }
2501     } else {
2502       for (G4int j = 0; j <= i; ++j) {
2503         s0 = pcm[0][j];
2504         s1 = pcm[1][j];
2505         s2 = pcm[2][j];
2506         energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2507         a = s0 * c - s1 * s;  //  rotation
2508         pcm[1][j] = s0 * s + s1 * c;
2509         b = pcm[2][j];
2510         pcm[0][j] = a * cb - b * sb;
2511         pcm[2][j] = a * sb + b * cb;
2512       }
2513     }
2514   }
2515   for (i = 0; i < vecLen; ++i) {
2516     vec[i]->SetMomentum(pcm[0][i] * GeV, pcm[1][i] * GeV, pcm[2][i] * GeV);
2517     vec[i]->SetTotalEnergy(energy[i] * GeV);
2518   }
2519   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2520   return weight;
2521 }
2522 
2523 G4double FullModelReactionDynamics::normal() {
2524   G4double ran = -6.0;
2525   for (G4int i = 0; i < 12; ++i)
2526     ran += G4UniformRand();
2527   return ran;
2528 }
2529 
2530 G4int FullModelReactionDynamics::Poisson(G4double x)  // generation of poisson distribution
2531 {
2532   G4int iran;
2533   G4double ran;
2534 
2535   if (x > 9.9)  // use normal distribution with sigma^2 = <x>
2536     iran = static_cast<G4int>(std::max(0.0, x + normal() * std::sqrt(x)));
2537   else {
2538     G4int mm = G4int(5.0 * x);
2539     if (mm <= 0)  // for very small x try iran=1,2,3
2540     {
2541       G4double p1 = x * std::exp(-x);
2542       G4double p2 = x * p1 / 2.0;
2543       G4double p3 = x * p2 / 3.0;
2544       ran = G4UniformRand();
2545       if (ran < p3)
2546         iran = 3;
2547       else if (ran < p2)  // this is original Geisha, it should be ran < p2+p3
2548         iran = 2;
2549       else if (ran < p1)  // should be ran < p1+p2+p3
2550         iran = 1;
2551       else
2552         iran = 0;
2553     } else {
2554       iran = 0;
2555       G4double r = std::exp(-x);
2556       ran = G4UniformRand();
2557       if (ran > r) {
2558         G4double rrr;
2559         G4double rr = r;
2560         for (G4int i = 1; i <= mm; ++i) {
2561           iran++;
2562           if (i > 5)  // Stirling's formula for large numbers
2563             rrr = std::exp(i * std::log(x) - (i + 0.5) * std::log((G4double)i) + i - 0.9189385);
2564           else
2565             rrr = std::pow(x, i) / Factorial(i);
2566           rr += r * rrr;
2567           if (ran <= rr)
2568             break;
2569         }
2570       }
2571     }
2572   }
2573   return iran;
2574 }
2575 
2576 G4int FullModelReactionDynamics::Factorial(G4int n) {  // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
2577   G4int m = std::min(n, 10);
2578   G4int result = 1;
2579   if (m <= 1)
2580     return result;
2581   for (G4int i = 2; i <= m; ++i)
2582     result *= i;
2583   return result;
2584 }
2585 
2586 void FullModelReactionDynamics::Defs1(const G4ReactionProduct &modifiedOriginal,
2587                                       G4ReactionProduct &currentParticle,
2588                                       G4ReactionProduct &targetParticle,
2589                                       G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2590                                       G4int &vecLen) {
2591   const G4double pjx = modifiedOriginal.GetMomentum().x() / MeV;
2592   const G4double pjy = modifiedOriginal.GetMomentum().y() / MeV;
2593   const G4double pjz = modifiedOriginal.GetMomentum().z() / MeV;
2594   const G4double p = modifiedOriginal.GetMomentum().mag() / MeV;
2595   if (pjx * pjx + pjy * pjy > 0.0) {
2596     G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
2597     cost = pjz / p;
2598     sint = 0.5 * (std::sqrt(std::abs((1.0 - cost) * (1.0 + cost))) + std::sqrt(pjx * pjx + pjy * pjy) / p);
2599     if (pjy < 0.0)
2600       ph = 3 * halfpi;
2601     else
2602       ph = halfpi;
2603     if (std::abs(pjx) > 0.001 * MeV)
2604       ph = std::atan2(pjy, pjx);
2605     cosp = std::cos(ph);
2606     sinp = std::sin(ph);
2607     pix = currentParticle.GetMomentum().x() / MeV;
2608     piy = currentParticle.GetMomentum().y() / MeV;
2609     piz = currentParticle.GetMomentum().z() / MeV;
2610     currentParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2611                                 cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2612                                 -sint * pix * MeV + cost * piz * MeV);
2613     pix = targetParticle.GetMomentum().x() / MeV;
2614     piy = targetParticle.GetMomentum().y() / MeV;
2615     piz = targetParticle.GetMomentum().z() / MeV;
2616     targetParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2617                                cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2618                                -sint * pix * MeV + cost * piz * MeV);
2619     for (G4int i = 0; i < vecLen; ++i) {
2620       pix = vec[i]->GetMomentum().x() / MeV;
2621       piy = vec[i]->GetMomentum().y() / MeV;
2622       piz = vec[i]->GetMomentum().z() / MeV;
2623       vec[i]->SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2624                           cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2625                           -sint * pix * MeV + cost * piz * MeV);
2626     }
2627   } else {
2628     if (pjz < 0.0) {
2629       currentParticle.SetMomentum(-currentParticle.GetMomentum().z());
2630       targetParticle.SetMomentum(-targetParticle.GetMomentum().z());
2631       for (G4int i = 0; i < vecLen; ++i)
2632         vec[i]->SetMomentum(-vec[i]->GetMomentum().z());
2633     }
2634   }
2635 }
2636 
2637 void FullModelReactionDynamics::Rotate(
2638     const G4double numberofFinalStateNucleons,
2639     const G4ThreeVector &temp,
2640     const G4ReactionProduct &modifiedOriginal,  // Fermi motion & evap. effect included
2641     const G4HadProjectile *originalIncident,    // original incident particle
2642     const G4Nucleus &targetNucleus,
2643     G4ReactionProduct &currentParticle,
2644     G4ReactionProduct &targetParticle,
2645     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2646     G4int &vecLen) {
2647   // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
2648   //
2649   //   Rotate in direction of z-axis, this does disturb in some way our
2650   //    inclusive distributions, but it is necessary for momentum conservation
2651   //
2652   const G4double atomicWeight = targetNucleus.GetN_asInt();
2653   const G4double logWeight = std::log(atomicWeight);
2654 
2655   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2656   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2657   G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2658 
2659   G4int i;
2660   G4ThreeVector pseudoParticle[4];
2661   for (i = 0; i < 4; ++i)
2662     pseudoParticle[i].set(0, 0, 0);
2663   pseudoParticle[0] = currentParticle.GetMomentum() + targetParticle.GetMomentum();
2664   for (i = 0; i < vecLen; ++i)
2665     pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
2666   //
2667   //  Some smearing in transverse direction from Fermi motion
2668   //
2669   G4float pp, pp1;
2670   G4double alekw, p, rthnve, phinve;
2671   G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
2672 
2673   r1 = twopi * G4UniformRand();
2674   r2 = G4UniformRand();
2675   a1 = std::sqrt(-2.0 * std::log(r2));
2676   ran1 = a1 * std::sin(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2677   ran2 = a1 * std::cos(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2678   G4ThreeVector fermi(ran1, ran2, 0);
2679 
2680   pseudoParticle[0] = pseudoParticle[0] + fermi;  // all particles + fermi
2681   pseudoParticle[2] = temp;                       // original in cms system
2682   pseudoParticle[3] = pseudoParticle[0];
2683 
2684   pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
2685   G4double rotation = 2. * pi * G4UniformRand();
2686   pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
2687   pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
2688   for (G4int ii = 1; ii <= 3; ii++) {
2689     p = pseudoParticle[ii].mag();
2690     if (p == 0.0)
2691       pseudoParticle[ii] = G4ThreeVector(0.0, 0.0, 0.0);
2692     else
2693       pseudoParticle[ii] = pseudoParticle[ii] * (1. / p);
2694   }
2695 
2696   pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
2697   pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
2698   pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
2699   currentParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2700 
2701   pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
2702   pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
2703   pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
2704   targetParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2705 
2706   for (i = 0; i < vecLen; ++i) {
2707     pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
2708     pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
2709     pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
2710     vec[i]->SetMomentum(pxTemp, pyTemp, pzTemp);
2711   }
2712   //
2713   //  Rotate in direction of primary particle, subtract binding energies
2714   //   and make some further corrections if required
2715   //
2716   Defs1(modifiedOriginal, currentParticle, targetParticle, vec, vecLen);
2717   G4double ekin;
2718   G4double dekin = 0.0;
2719   G4double ek1 = 0.0;
2720   G4int npions = 0;
2721   if (atomicWeight >= 1.5)  // self-absorption in heavy molecules
2722   {
2723     // corrections for single particle spectra (shower particles)
2724     //
2725     const G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00};
2726     const G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65};
2727     alekw = std::log(originalIncident->GetKineticEnergy() / GeV);
2728     exh = 1.0;
2729     if (alekw > alem[0])  //   get energy bin
2730     {
2731       exh = val0[6];
2732       for (G4int j = 1; j < 7; ++j) {
2733         if (alekw < alem[j])  // use linear interpolation/extrapolation
2734         {
2735           G4double rcnve = (val0[j] - val0[j - 1]) / (alem[j] - alem[j - 1]);
2736           exh = rcnve * alekw + val0[j - 1] - rcnve * alem[j - 1];
2737           break;
2738         }
2739       }
2740       exh = 1.0 - exh;
2741     }
2742     const G4double cfa = 0.025 * ((atomicWeight - 1.) / 120.) * std::exp(-(atomicWeight - 1.) / 120.);
2743     ekin = currentParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2744     ekin = std::max(1.0e-6, ekin);
2745     xxh = 1.0;
2746     dekin += ekin * (1.0 - xxh);
2747     ekin *= xxh;
2748     currentParticle.SetKineticEnergy(ekin * GeV);
2749     pp = currentParticle.GetTotalMomentum() / MeV;
2750     pp1 = currentParticle.GetMomentum().mag() / MeV;
2751     if (pp1 < 0.001 * MeV) {
2752       rthnve = pi * G4UniformRand();
2753       phinve = twopi * G4UniformRand();
2754       currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2755                                   pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2756                                   pp * std::cos(rthnve) * MeV);
2757     } else
2758       currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2759     ekin = targetParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2760     ekin = std::max(1.0e-6, ekin);
2761     xxh = 1.0;
2762     if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2763         (targetParticle.GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2764       xxh = exh;
2765     dekin += ekin * (1.0 - xxh);
2766     ekin *= xxh;
2767     if ((targetParticle.GetDefinition() == aPiPlus) || (targetParticle.GetDefinition() == aPiZero) ||
2768         (targetParticle.GetDefinition() == aPiMinus)) {
2769       ++npions;
2770       ek1 += ekin;
2771     }
2772     targetParticle.SetKineticEnergy(ekin * GeV);
2773     pp = targetParticle.GetTotalMomentum() / MeV;
2774     pp1 = targetParticle.GetMomentum().mag() / MeV;
2775     if (pp1 < 0.001 * MeV) {
2776       rthnve = pi * G4UniformRand();
2777       phinve = twopi * G4UniformRand();
2778       targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2779                                  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2780                                  pp * std::cos(rthnve) * MeV);
2781     } else
2782       targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
2783     for (i = 0; i < vecLen; ++i) {
2784       ekin = vec[i]->GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2785       ekin = std::max(1.0e-6, ekin);
2786       xxh = 1.0;
2787       if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2788           (vec[i]->GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2789         xxh = exh;
2790       dekin += ekin * (1.0 - xxh);
2791       ekin *= xxh;
2792       if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2793           (vec[i]->GetDefinition() == aPiMinus)) {
2794         ++npions;
2795         ek1 += ekin;
2796       }
2797       vec[i]->SetKineticEnergy(ekin * GeV);
2798       pp = vec[i]->GetTotalMomentum() / MeV;
2799       pp1 = vec[i]->GetMomentum().mag() / MeV;
2800       if (pp1 < 0.001 * MeV) {
2801         rthnve = pi * G4UniformRand();
2802         phinve = twopi * G4UniformRand();
2803         vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2804                             pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2805                             pp * std::cos(rthnve) * MeV);
2806       } else
2807         vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2808     }
2809   }
2810   if ((ek1 != 0.0) && (npions > 0)) {
2811     dekin = 1.0 + dekin / ek1;
2812     //
2813     //  first do the incident particle
2814     //
2815     if ((currentParticle.GetDefinition() == aPiPlus) || (currentParticle.GetDefinition() == aPiZero) ||
2816         (currentParticle.GetDefinition() == aPiMinus)) {
2817       currentParticle.SetKineticEnergy(std::max(0.001 * MeV, dekin * currentParticle.GetKineticEnergy()));
2818       pp = currentParticle.GetTotalMomentum() / MeV;
2819       pp1 = currentParticle.GetMomentum().mag() / MeV;
2820       if (pp1 < 0.001) {
2821         rthnve = pi * G4UniformRand();
2822         phinve = twopi * G4UniformRand();
2823         currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2824                                     pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2825                                     pp * std::cos(rthnve) * MeV);
2826       } else
2827         currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2828     }
2829     for (i = 0; i < vecLen; ++i) {
2830       if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2831           (vec[i]->GetDefinition() == aPiMinus)) {
2832         vec[i]->SetKineticEnergy(std::max(0.001 * MeV, dekin * vec[i]->GetKineticEnergy()));
2833         pp = vec[i]->GetTotalMomentum() / MeV;
2834         pp1 = vec[i]->GetMomentum().mag() / MeV;
2835         if (pp1 < 0.001) {
2836           rthnve = pi * G4UniformRand();
2837           phinve = twopi * G4UniformRand();
2838           vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2839                               pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2840                               pp * std::cos(rthnve) * MeV);
2841         } else
2842           vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2843       }
2844     }
2845   }
2846 }
2847 
2848 void FullModelReactionDynamics::AddBlackTrackParticles(const G4double epnb,  // GeV
2849                                                        const G4int npnb,
2850                                                        const G4double edta,  // GeV
2851                                                        const G4int ndta,
2852                                                        const G4double sprob,
2853                                                        const G4double kineticMinimum,  // GeV
2854                                                        const G4double kineticFactor,   // GeV
2855                                                        const G4ReactionProduct &modifiedOriginal,
2856                                                        G4double spall,
2857                                                        const G4Nucleus &targetNucleus,
2858                                                        G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2859                                                        G4int &vecLen) {
2860   // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
2861   //
2862   // npnb is number of proton/neutron black track particles
2863   // ndta is the number of deuterons, tritons, and alphas produced
2864   // epnb is the kinetic energy available for proton/neutron black track particles
2865   // edta is the kinetic energy available for deuteron/triton/alpha particles
2866   //
2867 
2868   G4ParticleDefinition *aProton = G4Proton::Proton();
2869   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
2870   G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
2871   G4ParticleDefinition *aTriton = G4Triton::Triton();
2872   G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
2873 
2874   const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / MeV;
2875   G4double atomicWeight = targetNucleus.GetN_asInt();
2876   G4double atomicNumber = targetNucleus.GetZ_asInt();
2877 
2878   const G4double ika1 = 3.6;
2879   const G4double ika2 = 35.56;
2880   const G4double ika3 = 6.45;
2881   const G4double sp1 = 1.066;
2882 
2883   G4int i;
2884   G4double pp;
2885   // G4double totalQ = 0;
2886   G4double kinCreated = 0;
2887   G4double cfa = 0.025 * ((atomicWeight - 1.0) / 120.0) * std::exp(-(atomicWeight - 1.0) / 120.0);
2888   if (npnb > 0)  // first add protons and neutrons
2889   {
2890     G4double backwardKinetic = 0.0;
2891     G4int local_npnb = npnb;
2892     for (i = 0; i < npnb; ++i)
2893       if (G4UniformRand() < sprob)
2894         local_npnb--;
2895     G4double ekin = epnb / local_npnb;
2896 
2897     for (i = 0; i < local_npnb; ++i) {
2898       G4ReactionProduct *p1 = new G4ReactionProduct();
2899       if (backwardKinetic > epnb) {
2900         delete p1;
2901         break;
2902       }
2903       G4double ran = G4UniformRand();
2904       G4double kinetic = -ekin * std::log(ran) - cfa * (1.0 + 0.5 * normal());
2905       if (kinetic < 0.0)
2906         kinetic = -0.010 * std::log(ran);
2907       backwardKinetic += kinetic;
2908       if (backwardKinetic > epnb)
2909         kinetic = std::max(kineticMinimum, epnb - (backwardKinetic - kinetic));
2910       if (G4UniformRand() > (1.0 - atomicNumber / atomicWeight))
2911         p1->SetDefinition(aProton);
2912       else
2913         p1->SetDefinition(aNeutron);
2914       vec.SetElement(vecLen, p1);
2915       ++spall;
2916       G4double cost = G4UniformRand() * 2.0 - 1.0;
2917       G4double sint = std::sqrt(std::fabs(1.0 - cost * cost));
2918       G4double phi = twopi * G4UniformRand();
2919       vec[vecLen]->SetNewlyAdded(true);
2920       vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2921       kinCreated += kinetic;
2922       pp = vec[vecLen]->GetTotalMomentum() / MeV;
2923       vec[vecLen]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2924       vecLen++;
2925       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2926     }
2927     if ((atomicWeight >= 10.0) && (ekOriginal <= 2.0 * GeV)) {
2928       G4double ekw = ekOriginal / GeV;
2929       G4int ika, kk = 0;
2930       if (ekw > 1.0)
2931         ekw *= ekw;
2932       ekw = std::max(0.1, ekw);
2933       ika = G4int(ika1 * std::exp((atomicNumber * atomicNumber / atomicWeight - ika2) / ika3) / ekw);
2934       if (ika > 0) {
2935         for (i = (vecLen - 1); i >= 0; --i) {
2936           if ((vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded()) {
2937             vec[i]->SetDefinitionAndUpdateE(aNeutron);  // modified 22-Oct-97
2938             if (++kk > ika)
2939               break;
2940           }
2941         }
2942       }
2943     }
2944   }
2945   if (ndta > 0)  //  now, try to add deuterons, tritons and alphas
2946   {
2947     G4double backwardKinetic = 0.0;
2948     G4int local_ndta = ndta;
2949     for (i = 0; i < ndta; ++i)
2950       if (G4UniformRand() < sprob)
2951         local_ndta--;
2952     G4double ekin = edta / local_ndta;
2953 
2954     for (i = 0; i < local_ndta; ++i) {
2955       G4ReactionProduct *p2 = new G4ReactionProduct();
2956       if (backwardKinetic > edta) {
2957         delete p2;
2958         break;
2959       }
2960       G4double ran = G4UniformRand();
2961       G4double kinetic = -ekin * std::log(ran) - cfa * (1. + 0.5 * normal());
2962       if (kinetic < 0.0)
2963         kinetic = kineticFactor * std::log(ran);
2964       backwardKinetic += kinetic;
2965       if (backwardKinetic > edta)
2966         kinetic = edta - (backwardKinetic - kinetic);
2967       if (kinetic < 0.0)
2968         kinetic = kineticMinimum;
2969       G4double cost = 2.0 * G4UniformRand() - 1.0;
2970       G4double sint = std::sqrt(std::max(0.0, (1.0 - cost * cost)));
2971       G4double phi = twopi * G4UniformRand();
2972       ran = G4UniformRand();
2973       if (ran <= 0.60)
2974         p2->SetDefinition(aDeuteron);
2975       else if (ran <= 0.90)
2976         p2->SetDefinition(aTriton);
2977       else
2978         p2->SetDefinition(anAlpha);
2979       spall += p2->GetMass() / GeV * sp1;
2980       if (spall > atomicWeight) {
2981         delete p2;
2982         break;
2983       }
2984       vec.SetElement(vecLen, p2);
2985       vec[vecLen]->SetNewlyAdded(true);
2986       vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2987       kinCreated += kinetic;
2988       pp = vec[vecLen]->GetTotalMomentum() / MeV;
2989       vec[vecLen++]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2990       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2991     }
2992   }
2993   // G4double delta = epnb+edta - kinCreated;
2994 }
2995 
2996 void FullModelReactionDynamics::MomentumCheck(const G4ReactionProduct &modifiedOriginal,
2997                                               G4ReactionProduct &currentParticle,
2998                                               G4ReactionProduct &targetParticle,
2999                                               G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
3000                                               G4int &vecLen) {
3001   const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / MeV;
3002   G4double testMomentum = currentParticle.GetMomentum().mag() / MeV;
3003   G4double pMass;
3004   if (testMomentum >= pOriginal) {
3005     pMass = currentParticle.GetMass() / MeV;
3006     currentParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3007     currentParticle.SetMomentum(currentParticle.GetMomentum() * (pOriginal / testMomentum));
3008   }
3009   testMomentum = targetParticle.GetMomentum().mag() / MeV;
3010   if (testMomentum >= pOriginal) {
3011     pMass = targetParticle.GetMass() / MeV;
3012     targetParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3013     targetParticle.SetMomentum(targetParticle.GetMomentum() * (pOriginal / testMomentum));
3014   }
3015   for (G4int i = 0; i < vecLen; ++i) {
3016     testMomentum = vec[i]->GetMomentum().mag() / MeV;
3017     if (testMomentum >= pOriginal) {
3018       pMass = vec[i]->GetMass() / MeV;
3019       vec[i]->SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3020       vec[i]->SetMomentum(vec[i]->GetMomentum() * (pOriginal / testMomentum));
3021     }
3022   }
3023 }
3024 
3025 void FullModelReactionDynamics::ProduceStrangeParticlePairs(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
3026                                                             G4int &vecLen,
3027                                                             const G4ReactionProduct &modifiedOriginal,
3028                                                             const G4DynamicParticle *originalTarget,
3029                                                             G4ReactionProduct &currentParticle,
3030                                                             G4ReactionProduct &targetParticle,
3031                                                             G4bool &incidentHasChanged,
3032                                                             G4bool &targetHasChanged) {
3033   // derived from original FORTRAN code STPAIR by H. Fesefeldt (16-Dec-1987)
3034   //
3035   // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
3036   //                            K+ Y0, K0 Y+,  K0 Y-
3037   // For antibaryon induced reactions half of the cross sections KB YB
3038   // pairs are produced.  Charge is not conserved, no experimental data available
3039   // for exclusive reactions, therefore some average behaviour assumed.
3040   // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
3041   //
3042   if (vecLen == 0)
3043     return;
3044   //
3045   // the following protects against annihilation processes
3046   //
3047   if (currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0)
3048     return;
3049 
3050   const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
3051   const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass() / GeV;
3052   G4double targetMass = originalTarget->GetDefinition()->GetPDGMass() / GeV;
3053   G4double centerofmassEnergy =
3054       std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal);  // GeV
3055   G4double currentMass = currentParticle.GetMass() / GeV;
3056   G4double availableEnergy = centerofmassEnergy - (targetMass + currentMass);
3057   if (availableEnergy <= 1.0)
3058     return;
3059 
3060   G4ParticleDefinition *aProton = G4Proton::Proton();
3061   G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
3062   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3063   G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
3064   G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
3065   G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
3066   G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
3067   G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
3068   G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
3069   G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
3070   G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
3071   G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
3072   G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
3073   G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
3074   G4ParticleDefinition *aLambda = G4Lambda::Lambda();
3075   G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
3076 
3077   const G4double protonMass = aProton->GetPDGMass() / GeV;
3078   const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass() / GeV;
3079   //
3080   // determine the center of mass energy bin
3081   //
3082   const G4double avrs[] = {3., 4., 5., 6., 7., 8., 9., 10., 20., 30., 40., 50.};
3083 
3084   G4int ibin, i3, i4;
3085   G4double avk, avy, avn, ran;
3086   G4int i = 1;
3087   while ((i < 12) && (centerofmassEnergy > avrs[i]))
3088     ++i;
3089   if (i == 12)
3090     ibin = 11;
3091   else
3092     ibin = i;
3093   //
3094   // the fortran code chooses a random replacement of produced kaons
3095   //  but does not take into account charge conservation
3096   //
3097   if (vecLen == 1)  // we know that vecLen > 0
3098   {
3099     i3 = 0;
3100     i4 = 1;  // note that we will be adding a new secondary particle in this case only
3101   } else     // otherwise  0 <= i3,i4 < vecLen
3102   {
3103     G4double ran = G4UniformRand();
3104     while (ran == 1.0)
3105       ran = G4UniformRand();
3106     i4 = i3 = G4int(vecLen * ran);
3107     while (i3 == i4) {
3108       ran = G4UniformRand();
3109       while (ran == 1.0)
3110         ran = G4UniformRand();
3111       i4 = G4int(vecLen * ran);
3112     }
3113   }
3114   //
3115   // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
3116   //
3117   const G4double avkkb[] = {0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075, 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573};
3118   const G4double avky[] = {0.005, 0.03, 0.064, 0.095, 0.115, 0.13, 0.145, 0.155, 0.20, 0.205, 0.210, 0.212};
3119   const G4double avnnb[] = {0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02, 0.04, 0.05, 0.12, 0.15, 0.18, 0.20};
3120 
3121   avk = (std::log(avkkb[ibin]) - std::log(avkkb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3122             (avrs[ibin] - avrs[ibin - 1]) +
3123         std::log(avkkb[ibin - 1]);
3124   avk = std::exp(avk);
3125 
3126   avy = (std::log(avky[ibin]) - std::log(avky[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3127             (avrs[ibin] - avrs[ibin - 1]) +
3128         std::log(avky[ibin - 1]);
3129   avy = std::exp(avy);
3130 
3131   avn = (std::log(avnnb[ibin]) - std::log(avnnb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3132             (avrs[ibin] - avrs[ibin - 1]) +
3133         std::log(avnnb[ibin - 1]);
3134   avn = std::exp(avn);
3135 
3136   if (avk + avy + avn <= 0.0)
3137     return;
3138 
3139   if (currentMass < protonMass)
3140     avy /= 2.0;
3141   if (targetMass < protonMass)
3142     avy = 0.0;
3143   avy += avk + avn;
3144   avk += avn;
3145   ran = G4UniformRand();
3146   if (ran < avn) {
3147     if (availableEnergy < 2.0)
3148       return;
3149     if (vecLen == 1)  // add a new secondary
3150     {
3151       G4ReactionProduct *p1 = new G4ReactionProduct;
3152       if (G4UniformRand() < 0.5) {
3153         vec[0]->SetDefinition(aNeutron);
3154         p1->SetDefinition(anAntiNeutron);
3155         (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3156         vec[0]->SetMayBeKilled(false);
3157         p1->SetMayBeKilled(false);
3158       } else {
3159         vec[0]->SetDefinition(aProton);
3160         p1->SetDefinition(anAntiProton);
3161         (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3162         vec[0]->SetMayBeKilled(false);
3163         p1->SetMayBeKilled(false);
3164       }
3165       vec.SetElement(vecLen++, p1);
3166       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3167     } else {  // replace two secondaries
3168       if (G4UniformRand() < 0.5) {
3169         vec[i3]->SetDefinition(aNeutron);
3170         vec[i4]->SetDefinition(anAntiNeutron);
3171         vec[i3]->SetMayBeKilled(false);
3172         vec[i4]->SetMayBeKilled(false);
3173       } else {
3174         vec[i3]->SetDefinition(aProton);
3175         vec[i4]->SetDefinition(anAntiProton);
3176         vec[i3]->SetMayBeKilled(false);
3177         vec[i4]->SetMayBeKilled(false);
3178       }
3179     }
3180   } else if (ran < avk) {
3181     if (availableEnergy < 1.0)
3182       return;
3183 
3184     const G4double kkb[] = {0.2500, 0.3750, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8750, 1.000};
3185     const G4int ipakkb1[] = {10, 10, 10, 11, 11, 12, 12, 11, 12};
3186     const G4int ipakkb2[] = {13, 11, 12, 11, 12, 11, 12, 13, 13};
3187     ran = G4UniformRand();
3188     i = 0;
3189     while ((i < 9) && (ran >= kkb[i]))
3190       ++i;
3191     if (i == 9)
3192       return;
3193     //
3194     // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
3195     // charge       +  -   +  0   +  0   0  0   0  0   0  0   0  0   0  -   0  -
3196     //
3197     switch (ipakkb1[i]) {
3198       case 10:
3199         vec[i3]->SetDefinition(aKaonPlus);
3200         vec[i3]->SetMayBeKilled(false);
3201         break;
3202       case 11:
3203         vec[i3]->SetDefinition(aKaonZS);
3204         vec[i3]->SetMayBeKilled(false);
3205         break;
3206       case 12:
3207         vec[i3]->SetDefinition(aKaonZL);
3208         vec[i3]->SetMayBeKilled(false);
3209         break;
3210     }
3211     if (vecLen == 1)  // add a secondary
3212     {
3213       G4ReactionProduct *p1 = new G4ReactionProduct;
3214       switch (ipakkb2[i]) {
3215         case 11:
3216           p1->SetDefinition(aKaonZS);
3217           p1->SetMayBeKilled(false);
3218           break;
3219         case 12:
3220           p1->SetDefinition(aKaonZL);
3221           p1->SetMayBeKilled(false);
3222           break;
3223         case 13:
3224           p1->SetDefinition(aKaonMinus);
3225           p1->SetMayBeKilled(false);
3226           break;
3227       }
3228       (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3229       vec.SetElement(vecLen++, p1);
3230       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3231     } else  // replace
3232     {
3233       switch (ipakkb2[i]) {
3234         case 11:
3235           vec[i4]->SetDefinition(aKaonZS);
3236           vec[i4]->SetMayBeKilled(false);
3237           break;
3238         case 12:
3239           vec[i4]->SetDefinition(aKaonZL);
3240           vec[i4]->SetMayBeKilled(false);
3241           break;
3242         case 13:
3243           vec[i4]->SetDefinition(aKaonMinus);
3244           vec[i4]->SetMayBeKilled(false);
3245           break;
3246       }
3247     }
3248   } else if (ran < avy) {
3249     if (availableEnergy < 1.6)
3250       return;
3251 
3252     const G4double ky[] = {0.200, 0.300, 0.400, 0.550, 0.625, 0.700, 0.800, 0.850, 0.900, 0.950, 0.975, 1.000};
3253     const G4int ipaky1[] = {18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22};
3254     const G4int ipaky2[] = {10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12};
3255     const G4int ipakyb1[] = {19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25};
3256     const G4int ipakyb2[] = {13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11};
3257     ran = G4UniformRand();
3258     i = 0;
3259     while ((i < 12) && (ran > ky[i]))
3260       ++i;
3261     if (i == 12)
3262       return;
3263     if ((currentMass < protonMass) || (G4UniformRand() < 0.5)) {
3264       // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
3265       //             0  +   0  0   0  0   +  +   +  0   +  0
3266       //
3267       //             21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
3268       //             0  +   0  0   0  0   -  +   -  0   -  0
3269       switch (ipaky1[i]) {
3270         case 18:
3271           targetParticle.SetDefinition(aLambda);
3272           break;
3273         case 20:
3274           targetParticle.SetDefinition(aSigmaPlus);
3275           break;
3276         case 21:
3277           targetParticle.SetDefinition(aSigmaZero);
3278           break;
3279         case 22:
3280           targetParticle.SetDefinition(aSigmaMinus);
3281           break;
3282       }
3283       targetHasChanged = true;
3284       switch (ipaky2[i]) {
3285         case 10:
3286           vec[i3]->SetDefinition(aKaonPlus);
3287           vec[i3]->SetMayBeKilled(false);
3288           break;
3289         case 11:
3290           vec[i3]->SetDefinition(aKaonZS);
3291           vec[i3]->SetMayBeKilled(false);
3292           break;
3293         case 12:
3294           vec[i3]->SetDefinition(aKaonZL);
3295           vec[i3]->SetMayBeKilled(false);
3296           break;
3297       }
3298     } else  // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
3299     {
3300       if ((currentParticle.GetDefinition() == anAntiProton) || (currentParticle.GetDefinition() == anAntiNeutron) ||
3301           (currentParticle.GetDefinition() == anAntiLambda) || (currentMass > sigmaMinusMass)) {
3302         switch (ipakyb1[i]) {
3303           case 19:
3304             currentParticle.SetDefinitionAndUpdateE(anAntiLambda);
3305             break;
3306           case 23:
3307             currentParticle.SetDefinitionAndUpdateE(anAntiSigmaPlus);
3308             break;
3309           case 24:
3310             currentParticle.SetDefinitionAndUpdateE(anAntiSigmaZero);
3311             break;
3312           case 25:
3313             currentParticle.SetDefinitionAndUpdateE(anAntiSigmaMinus);
3314             break;
3315         }
3316         incidentHasChanged = true;
3317         switch (ipakyb2[i]) {
3318           case 11:
3319             vec[i3]->SetDefinition(aKaonZS);
3320             vec[i3]->SetMayBeKilled(false);
3321             break;
3322           case 12:
3323             vec[i3]->SetDefinition(aKaonZL);
3324             vec[i3]->SetMayBeKilled(false);
3325             break;
3326           case 13:
3327             vec[i3]->SetDefinition(aKaonMinus);
3328             vec[i3]->SetMayBeKilled(false);
3329             break;
3330         }
3331       } else {
3332         switch (ipaky1[i]) {
3333           case 18:
3334             currentParticle.SetDefinitionAndUpdateE(aLambda);
3335             break;
3336           case 20:
3337             currentParticle.SetDefinitionAndUpdateE(aSigmaPlus);
3338             break;
3339           case 21:
3340             currentParticle.SetDefinitionAndUpdateE(aSigmaZero);
3341             break;
3342           case 22:
3343             currentParticle.SetDefinitionAndUpdateE(aSigmaMinus);
3344             break;
3345         }
3346         incidentHasChanged = true;
3347         switch (ipaky2[i]) {
3348           case 10:
3349             vec[i3]->SetDefinition(aKaonPlus);
3350             vec[i3]->SetMayBeKilled(false);
3351             break;
3352           case 11:
3353             vec[i3]->SetDefinition(aKaonZS);
3354             vec[i3]->SetMayBeKilled(false);
3355             break;
3356           case 12:
3357             vec[i3]->SetDefinition(aKaonZL);
3358             vec[i3]->SetMayBeKilled(false);
3359             break;
3360         }
3361       }
3362     }
3363   } else
3364     return;
3365   //
3366   //  check the available energy
3367   //   if there is not enough energy for kkb/ky pair production
3368   //   then reduce the number of secondary particles
3369   //  NOTE:
3370   //        the number of secondaries may have been changed
3371   //        the incident and/or target particles may have changed
3372   //        charge conservation is ignored (as well as strangness conservation)
3373   //
3374   currentMass = currentParticle.GetMass() / GeV;
3375   targetMass = targetParticle.GetMass() / GeV;
3376 
3377   G4double energyCheck = centerofmassEnergy - (currentMass + targetMass);
3378   for (i = 0; i < vecLen; ++i) {
3379     energyCheck -= vec[i]->GetMass() / GeV;
3380     if (energyCheck < 0.0)  // chop off the secondary List
3381     {
3382       vecLen = std::max(0, --i);  // looks like a memory leak @@@@@@@@@@@@
3383       G4int j;
3384       for (j = i; j < vecLen; j++)
3385         delete vec[j];
3386       break;
3387     }
3388   }
3389   return;
3390 }
3391 
3392 void FullModelReactionDynamics::NuclearReaction(G4FastVector<G4ReactionProduct, 4> &vec,
3393                                                 G4int &vecLen,
3394                                                 const G4HadProjectile *originalIncident,
3395                                                 const G4Nucleus &targetNucleus,
3396                                                 const G4double theAtomicMass,
3397                                                 const G4double *mass) {
3398   // derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
3399   //
3400   // Nuclear reaction kinematics at low energies
3401   //
3402   G4ParticleDefinition *aGamma = G4Gamma::Gamma();
3403   G4ParticleDefinition *aProton = G4Proton::Proton();
3404   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3405   G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
3406   G4ParticleDefinition *aTriton = G4Triton::Triton();
3407   G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
3408 
3409   const G4double aProtonMass = aProton->GetPDGMass() / MeV;
3410   const G4double aNeutronMass = aNeutron->GetPDGMass() / MeV;
3411   const G4double aDeuteronMass = aDeuteron->GetPDGMass() / MeV;
3412   const G4double aTritonMass = aTriton->GetPDGMass() / MeV;
3413   const G4double anAlphaMass = anAlpha->GetPDGMass() / MeV;
3414 
3415   G4ReactionProduct currentParticle;
3416   currentParticle = *originalIncident;
3417   //
3418   // Set beam particle, take kinetic energy of current particle as the
3419   // fundamental quantity.  Due to the difficult kinematic, all masses have to
3420   // be assigned the best measured values
3421   //
3422   G4double p = currentParticle.GetTotalMomentum();
3423   G4double pp = currentParticle.GetMomentum().mag();
3424   if (pp <= 0.001 * MeV) {
3425     G4double phinve = twopi * G4UniformRand();
3426     G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3427     currentParticle.SetMomentum(
3428         p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3429   } else
3430     currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3431   //
3432   // calculate Q-value of reactions
3433   //
3434   G4double currentKinetic = currentParticle.GetKineticEnergy() / MeV;
3435   G4double currentMass = currentParticle.GetDefinition()->GetPDGMass() / MeV;
3436   G4double qv = currentKinetic + theAtomicMass + currentMass;
3437 
3438   G4double qval[9];
3439   qval[0] = qv - mass[0];
3440   qval[1] = qv - mass[1] - aNeutronMass;
3441   qval[2] = qv - mass[2] - aProtonMass;
3442   qval[3] = qv - mass[3] - aDeuteronMass;
3443   qval[4] = qv - mass[4] - aTritonMass;
3444   qval[5] = qv - mass[5] - anAlphaMass;
3445   qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
3446   qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
3447   qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
3448 
3449   if (currentParticle.GetDefinition() == aNeutron) {
3450     const G4double A = targetNucleus.GetN_asInt();  // atomic weight
3451     if (G4UniformRand() > ((A - 1.0) / 230.0) * ((A - 1.0) / 230.0))
3452       qval[0] = 0.0;
3453     if (G4UniformRand() >= currentKinetic / 7.9254 * A)
3454       qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
3455   } else
3456     qval[0] = 0.0;
3457 
3458   G4int i;
3459   qv = 0.0;
3460   for (i = 0; i < 9; ++i) {
3461     if (mass[i] < 500.0 * MeV)
3462       qval[i] = 0.0;
3463     if (qval[i] < 0.0)
3464       qval[i] = 0.0;
3465     qv += qval[i];
3466   }
3467   G4double qv1 = 0.0;
3468   G4double ran = G4UniformRand();
3469   G4int index;
3470   for (index = 0; index < 9; ++index) {
3471     if (qval[index] > 0.0) {
3472       qv1 += qval[index] / qv;
3473       if (ran <= qv1)
3474         break;
3475     }
3476   }
3477   if (index == 9)  // loop continued to the end
3478   {
3479     throw G4HadronicException(
3480         __FILE__,
3481         __LINE__,
3482         "FullModelReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
3483   }
3484   G4double ke = currentParticle.GetKineticEnergy() / GeV;
3485   G4int nt = 2;
3486   if ((index >= 6) || (G4UniformRand() < std::min(0.5, ke * 10.0)))
3487     nt = 3;
3488 
3489   G4ReactionProduct **v = new G4ReactionProduct *[3];
3490   v[0] = new G4ReactionProduct;
3491   v[1] = new G4ReactionProduct;
3492   v[2] = new G4ReactionProduct;
3493 
3494   v[0]->SetMass(mass[index] * MeV);
3495   switch (index) {
3496     case 0:
3497       v[1]->SetDefinition(aGamma);
3498       v[2]->SetDefinition(aGamma);
3499       break;
3500     case 1:
3501       v[1]->SetDefinition(aNeutron);
3502       v[2]->SetDefinition(aGamma);
3503       break;
3504     case 2:
3505       v[1]->SetDefinition(aProton);
3506       v[2]->SetDefinition(aGamma);
3507       break;
3508     case 3:
3509       v[1]->SetDefinition(aDeuteron);
3510       v[2]->SetDefinition(aGamma);
3511       break;
3512     case 4:
3513       v[1]->SetDefinition(aTriton);
3514       v[2]->SetDefinition(aGamma);
3515       break;
3516     case 5:
3517       v[1]->SetDefinition(anAlpha);
3518       v[2]->SetDefinition(aGamma);
3519       break;
3520     case 6:
3521       v[1]->SetDefinition(aNeutron);
3522       v[2]->SetDefinition(aNeutron);
3523       break;
3524     case 7:
3525       v[1]->SetDefinition(aNeutron);
3526       v[2]->SetDefinition(aProton);
3527       break;
3528     case 8:
3529       v[1]->SetDefinition(aProton);
3530       v[2]->SetDefinition(aProton);
3531       break;
3532   }
3533   //
3534   // calculate centre of mass energy
3535   //
3536   G4ReactionProduct pseudo1;
3537   pseudo1.SetMass(theAtomicMass * MeV);
3538   pseudo1.SetTotalEnergy(theAtomicMass * MeV);
3539   G4ReactionProduct pseudo2 = currentParticle + pseudo1;
3540   pseudo2.SetMomentum(pseudo2.GetMomentum() * (-1.0));
3541   //
3542   // use phase space routine in centre of mass system
3543   //
3544   G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
3545   tempV.Initialize(nt);
3546   G4int tempLen = 0;
3547   tempV.SetElement(tempLen++, v[0]);
3548   tempV.SetElement(tempLen++, v[1]);
3549   if (nt == 3)
3550     tempV.SetElement(tempLen++, v[2]);
3551   G4bool constantCrossSection = true;
3552   GenerateNBodyEvent(pseudo2.GetMass() / MeV, constantCrossSection, tempV, tempLen);
3553   v[0]->Lorentz(*v[0], pseudo2);
3554   v[1]->Lorentz(*v[1], pseudo2);
3555   if (nt == 3)
3556     v[2]->Lorentz(*v[2], pseudo2);
3557 
3558   G4bool particleIsDefined = false;
3559   if (v[0]->GetMass() / MeV - aProtonMass < 0.1) {
3560     v[0]->SetDefinition(aProton);
3561     particleIsDefined = true;
3562   } else if (v[0]->GetMass() / MeV - aNeutronMass < 0.1) {
3563     v[0]->SetDefinition(aNeutron);
3564     particleIsDefined = true;
3565   } else if (v[0]->GetMass() / MeV - aDeuteronMass < 0.1) {
3566     v[0]->SetDefinition(aDeuteron);
3567     particleIsDefined = true;
3568   } else if (v[0]->GetMass() / MeV - aTritonMass < 0.1) {
3569     v[0]->SetDefinition(aTriton);
3570     particleIsDefined = true;
3571   } else if (v[0]->GetMass() / MeV - anAlphaMass < 0.1) {
3572     v[0]->SetDefinition(anAlpha);
3573     particleIsDefined = true;
3574   }
3575   currentParticle.SetKineticEnergy(std::max(0.001, currentParticle.GetKineticEnergy() / MeV));
3576   p = currentParticle.GetTotalMomentum();
3577   pp = currentParticle.GetMomentum().mag();
3578   if (pp <= 0.001 * MeV) {
3579     G4double phinve = twopi * G4UniformRand();
3580     G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3581     currentParticle.SetMomentum(
3582         p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3583   } else
3584     currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3585 
3586   if (particleIsDefined) {
3587     v[0]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[0]->GetKineticEnergy() / MeV));
3588     p = v[0]->GetTotalMomentum();
3589     pp = v[0]->GetMomentum().mag();
3590     if (pp <= 0.001 * MeV) {
3591       G4double phinve = twopi * G4UniformRand();
3592       G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3593       v[0]->SetMomentum(
3594           p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3595     } else
3596       v[0]->SetMomentum(v[0]->GetMomentum() * (p / pp));
3597   }
3598   if ((v[1]->GetDefinition() == aDeuteron) || (v[1]->GetDefinition() == aTriton) || (v[1]->GetDefinition() == anAlpha))
3599     v[1]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[1]->GetKineticEnergy() / MeV));
3600   else
3601     v[1]->SetKineticEnergy(std::max(0.001, v[1]->GetKineticEnergy() / MeV));
3602 
3603   p = v[1]->GetTotalMomentum();
3604   pp = v[1]->GetMomentum().mag();
3605   if (pp <= 0.001 * MeV) {
3606     G4double phinve = twopi * G4UniformRand();
3607     G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3608     v[1]->SetMomentum(
3609         p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3610   } else
3611     v[1]->SetMomentum(v[1]->GetMomentum() * (p / pp));
3612 
3613   if (nt == 3) {
3614     if ((v[2]->GetDefinition() == aDeuteron) || (v[2]->GetDefinition() == aTriton) ||
3615         (v[2]->GetDefinition() == anAlpha))
3616       v[2]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[2]->GetKineticEnergy() / MeV));
3617     else
3618       v[2]->SetKineticEnergy(std::max(0.001, v[2]->GetKineticEnergy() / MeV));
3619 
3620     p = v[2]->GetTotalMomentum();
3621     pp = v[2]->GetMomentum().mag();
3622     if (pp <= 0.001 * MeV) {
3623       G4double phinve = twopi * G4UniformRand();
3624       G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3625       v[2]->SetMomentum(
3626           p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3627     } else
3628       v[2]->SetMomentum(v[2]->GetMomentum() * (p / pp));
3629   }
3630   G4int del;
3631   for (del = 0; del < vecLen; del++)
3632     delete vec[del];
3633   vecLen = 0;
3634   if (particleIsDefined) {
3635     vec.SetElement(vecLen++, v[0]);
3636   } else {
3637     delete v[0];
3638   }
3639   vec.SetElement(vecLen++, v[1]);
3640   if (nt == 3) {
3641     vec.SetElement(vecLen++, v[2]);
3642   } else {
3643     delete v[2];
3644   }
3645   delete[] v;
3646   return;
3647 }
3648 
3649 /* end of file */