Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:04:28

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