Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-23 16:00:31 

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] = {0.0};
2285   //G4double *pd = new G4double [vecLen-1];
2286   for (i = 0; i < vecLen - 1; ++i) {
2287     if (emm[i + 1] * emm[i + 1] > 0.0) {
2288       G4double arg = emm[i + 1] * emm[i + 1] +
2289                      (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) * (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) /
2290                          (emm[i + 1] * emm[i + 1]) -
2291                      2.0 * (emm[i] * emm[i] + mass[i + 1] * mass[i + 1]);
2292       if (arg > 0.0)
2293         pd[i] = 0.5 * std::sqrt(arg);
2294     }
2295     if (pd[i] <= 0.0)  //  changed from  ==  on 02 April 98
2296       lzero = false;
2297     else
2298       wtmax += std::log(pd[i]);
2299   }
2300   G4double weight = 0.0;  // weight is returned by GenerateNBodyEvent
2301   if (lzero)
2302     weight = std::exp(std::max(std::min(wtmax, expxu), expxl));
2303 
2304   G4double bang, cb, sb, s0, s1, s2, c, s, esys, a, b, gama, beta;
2305   pcm[0][0] = 0.0;
2306   pcm[1][0] = pd[0];
2307   pcm[2][0] = 0.0;
2308   for (i = 1; i < vecLen; ++i) {
2309     pcm[0][i] = 0.0;
2310     pcm[1][i] = -pd[i - 1];
2311     pcm[2][i] = 0.0;
2312     bang = twopi * G4UniformRand();
2313     cb = std::cos(bang);
2314     sb = std::sin(bang);
2315     c = 2.0 * G4UniformRand() - 1.0;
2316     s = std::sqrt(std::fabs(1.0 - c * c));
2317     if (i < vecLen - 1) {
2318       esys = std::sqrt(pd[i] * pd[i] + emm[i] * emm[i]);
2319       beta = pd[i] / esys;
2320       gama = esys / emm[i];
2321       for (G4int j = 0; j <= i; ++j) {
2322         s0 = pcm[0][j];
2323         s1 = pcm[1][j];
2324         s2 = pcm[2][j];
2325         energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2326         a = s0 * c - s1 * s;  //  rotation
2327         pcm[1][j] = s0 * s + s1 * c;
2328         b = pcm[2][j];
2329         pcm[0][j] = a * cb - b * sb;
2330         pcm[2][j] = a * sb + b * cb;
2331         pcm[1][j] = gama * (pcm[1][j] + beta * energy[j]);
2332       }
2333     } else {
2334       for (G4int j = 0; j <= i; ++j) {
2335         s0 = pcm[0][j];
2336         s1 = pcm[1][j];
2337         s2 = pcm[2][j];
2338         energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2339         a = s0 * c - s1 * s;  //  rotation
2340         pcm[1][j] = s0 * s + s1 * c;
2341         b = pcm[2][j];
2342         pcm[0][j] = a * cb - b * sb;
2343         pcm[2][j] = a * sb + b * cb;
2344       }
2345     }
2346   }
2347   for (i = 0; i < vecLen; ++i) {
2348     vec[i]->SetMomentum(pcm[0][i] * GeV, pcm[1][i] * GeV, pcm[2][i] * GeV);
2349     vec[i]->SetTotalEnergy(energy[i] * GeV);
2350   }
2351   // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2352   return weight;
2353 }
2354 
2355 G4double FullModelReactionDynamics::normal() {
2356   G4double ran = -6.0;
2357   for (G4int i = 0; i < 12; ++i)
2358     ran += G4UniformRand();
2359   return ran;
2360 }
2361 
2362 G4int FullModelReactionDynamics::Factorial(G4int n) {  // calculates factorial( n ) = n*(n-1)*(n-2)*...*1
2363   G4int m = std::min(n, 10);
2364   G4int result = 1;
2365   if (m <= 1)
2366     return result;
2367   for (G4int i = 2; i <= m; ++i)
2368     result *= i;
2369   return result;
2370 }
2371 
2372 void FullModelReactionDynamics::Defs1(const G4ReactionProduct &modifiedOriginal,
2373                                       G4ReactionProduct &currentParticle,
2374                                       G4ReactionProduct &targetParticle,
2375                                       G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2376                                       G4int &vecLen) {
2377   const G4double pjx = modifiedOriginal.GetMomentum().x() / MeV;
2378   const G4double pjy = modifiedOriginal.GetMomentum().y() / MeV;
2379   const G4double pjz = modifiedOriginal.GetMomentum().z() / MeV;
2380   const G4double p = modifiedOriginal.GetMomentum().mag() / MeV;
2381   if (pjx * pjx + pjy * pjy > 0.0) {
2382     G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
2383     cost = pjz / p;
2384     sint = 0.5 * (std::sqrt(std::abs((1.0 - cost) * (1.0 + cost))) + std::sqrt(pjx * pjx + pjy * pjy) / p);
2385     if (pjy < 0.0)
2386       ph = 3 * halfpi;
2387     else
2388       ph = halfpi;
2389     if (std::abs(pjx) > 0.001 * MeV)
2390       ph = std::atan2(pjy, pjx);
2391     cosp = std::cos(ph);
2392     sinp = std::sin(ph);
2393     pix = currentParticle.GetMomentum().x() / MeV;
2394     piy = currentParticle.GetMomentum().y() / MeV;
2395     piz = currentParticle.GetMomentum().z() / MeV;
2396     currentParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2397                                 cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2398                                 -sint * pix * MeV + cost * piz * MeV);
2399     pix = targetParticle.GetMomentum().x() / MeV;
2400     piy = targetParticle.GetMomentum().y() / MeV;
2401     piz = targetParticle.GetMomentum().z() / MeV;
2402     targetParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2403                                cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2404                                -sint * pix * MeV + cost * piz * MeV);
2405     for (G4int i = 0; i < vecLen; ++i) {
2406       pix = vec[i]->GetMomentum().x() / MeV;
2407       piy = vec[i]->GetMomentum().y() / MeV;
2408       piz = vec[i]->GetMomentum().z() / MeV;
2409       vec[i]->SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2410                           cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2411                           -sint * pix * MeV + cost * piz * MeV);
2412     }
2413   } else {
2414     if (pjz < 0.0) {
2415       currentParticle.SetMomentum(-currentParticle.GetMomentum().z());
2416       targetParticle.SetMomentum(-targetParticle.GetMomentum().z());
2417       for (G4int i = 0; i < vecLen; ++i)
2418         vec[i]->SetMomentum(-vec[i]->GetMomentum().z());
2419     }
2420   }
2421 }
2422 
2423 void FullModelReactionDynamics::Rotate(
2424     const G4double numberofFinalStateNucleons,
2425     const G4ThreeVector &temp,
2426     const G4ReactionProduct &modifiedOriginal,  // Fermi motion & evap. effect included
2427     const G4HadProjectile *originalIncident,    // original incident particle
2428     const G4Nucleus &targetNucleus,
2429     G4ReactionProduct &currentParticle,
2430     G4ReactionProduct &targetParticle,
2431     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2432     G4int &vecLen) {
2433   // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
2434   //
2435   //   Rotate in direction of z-axis, this does disturb in some way our
2436   //    inclusive distributions, but it is necessary for momentum conservation
2437   //
2438   const G4double atomicWeight = targetNucleus.GetN_asInt();
2439   const G4double logWeight = std::log(atomicWeight);
2440 
2441   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2442   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2443   G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2444 
2445   G4int i;
2446   G4ThreeVector pseudoParticle[4];
2447   for (i = 0; i < 4; ++i)
2448     pseudoParticle[i].set(0, 0, 0);
2449   pseudoParticle[0] = currentParticle.GetMomentum() + targetParticle.GetMomentum();
2450   for (i = 0; i < vecLen; ++i)
2451     pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
2452   //
2453   //  Some smearing in transverse direction from Fermi motion
2454   //
2455   G4float pp, pp1;
2456   G4double alekw, p, rthnve, phinve;
2457   G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
2458 
2459   r1 = twopi * G4UniformRand();
2460   r2 = G4UniformRand();
2461   a1 = std::sqrt(-2.0 * std::log(r2));
2462   ran1 = a1 * std::sin(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2463   ran2 = a1 * std::cos(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2464   G4ThreeVector fermi(ran1, ran2, 0);
2465 
2466   pseudoParticle[0] = pseudoParticle[0] + fermi;  // all particles + fermi
2467   pseudoParticle[2] = temp;                       // original in cms system
2468   pseudoParticle[3] = pseudoParticle[0];
2469 
2470   pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
2471   G4double rotation = 2. * pi * G4UniformRand();
2472   pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
2473   pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
2474   for (G4int ii = 1; ii <= 3; ii++) {
2475     p = pseudoParticle[ii].mag();
2476     if (p == 0.0)
2477       pseudoParticle[ii] = G4ThreeVector(0.0, 0.0, 0.0);
2478     else
2479       pseudoParticle[ii] = pseudoParticle[ii] * (1. / p);
2480   }
2481 
2482   pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
2483   pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
2484   pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
2485   currentParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2486 
2487   pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
2488   pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
2489   pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
2490   targetParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2491 
2492   for (i = 0; i < vecLen; ++i) {
2493     pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
2494     pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
2495     pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
2496     vec[i]->SetMomentum(pxTemp, pyTemp, pzTemp);
2497   }
2498   //
2499   //  Rotate in direction of primary particle, subtract binding energies
2500   //   and make some further corrections if required
2501   //
2502   Defs1(modifiedOriginal, currentParticle, targetParticle, vec, vecLen);
2503   G4double ekin;
2504   G4double dekin = 0.0;
2505   G4double ek1 = 0.0;
2506   G4int npions = 0;
2507   if (atomicWeight >= 1.5)  // self-absorption in heavy molecules
2508   {
2509     // corrections for single particle spectra (shower particles)
2510     //
2511     const G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00};
2512     const G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65};
2513     alekw = std::log(originalIncident->GetKineticEnergy() / GeV);
2514     exh = 1.0;
2515     if (alekw > alem[0])  //   get energy bin
2516     {
2517       exh = val0[6];
2518       for (G4int j = 1; j < 7; ++j) {
2519         if (alekw < alem[j])  // use linear interpolation/extrapolation
2520         {
2521           G4double rcnve = (val0[j] - val0[j - 1]) / (alem[j] - alem[j - 1]);
2522           exh = rcnve * alekw + val0[j - 1] - rcnve * alem[j - 1];
2523           break;
2524         }
2525       }
2526       exh = 1.0 - exh;
2527     }
2528     const G4double cfa = 0.025 * ((atomicWeight - 1.) / 120.) * std::exp(-(atomicWeight - 1.) / 120.);
2529     ekin = currentParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2530     ekin = std::max(1.0e-6, ekin);
2531     xxh = 1.0;
2532     dekin += ekin * (1.0 - xxh);
2533     ekin *= xxh;
2534     currentParticle.SetKineticEnergy(ekin * GeV);
2535     pp = currentParticle.GetTotalMomentum() / MeV;
2536     pp1 = currentParticle.GetMomentum().mag() / MeV;
2537     if (pp1 < 0.001 * MeV) {
2538       rthnve = pi * G4UniformRand();
2539       phinve = twopi * G4UniformRand();
2540       currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2541                                   pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2542                                   pp * std::cos(rthnve) * MeV);
2543     } else
2544       currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2545     ekin = targetParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2546     ekin = std::max(1.0e-6, ekin);
2547     xxh = 1.0;
2548     if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2549         (targetParticle.GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2550       xxh = exh;
2551     dekin += ekin * (1.0 - xxh);
2552     ekin *= xxh;
2553     if ((targetParticle.GetDefinition() == aPiPlus) || (targetParticle.GetDefinition() == aPiZero) ||
2554         (targetParticle.GetDefinition() == aPiMinus)) {
2555       ++npions;
2556       ek1 += ekin;
2557     }
2558     targetParticle.SetKineticEnergy(ekin * GeV);
2559     pp = targetParticle.GetTotalMomentum() / MeV;
2560     pp1 = targetParticle.GetMomentum().mag() / MeV;
2561     if (pp1 < 0.001 * MeV) {
2562       rthnve = pi * G4UniformRand();
2563       phinve = twopi * G4UniformRand();
2564       targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2565                                  pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2566                                  pp * std::cos(rthnve) * MeV);
2567     } else
2568       targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
2569     for (i = 0; i < vecLen; ++i) {
2570       ekin = vec[i]->GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2571       ekin = std::max(1.0e-6, ekin);
2572       xxh = 1.0;
2573       if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2574           (vec[i]->GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2575         xxh = exh;
2576       dekin += ekin * (1.0 - xxh);
2577       ekin *= xxh;
2578       if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2579           (vec[i]->GetDefinition() == aPiMinus)) {
2580         ++npions;
2581         ek1 += ekin;
2582       }
2583       vec[i]->SetKineticEnergy(ekin * GeV);
2584       pp = vec[i]->GetTotalMomentum() / MeV;
2585       pp1 = vec[i]->GetMomentum().mag() / MeV;
2586       if (pp1 < 0.001 * MeV) {
2587         rthnve = pi * G4UniformRand();
2588         phinve = twopi * G4UniformRand();
2589         vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2590                             pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2591                             pp * std::cos(rthnve) * MeV);
2592       } else
2593         vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2594     }
2595   }
2596   if ((ek1 != 0.0) && (npions > 0)) {
2597     dekin = 1.0 + dekin / ek1;
2598     //
2599     //  first do the incident particle
2600     //
2601     if ((currentParticle.GetDefinition() == aPiPlus) || (currentParticle.GetDefinition() == aPiZero) ||
2602         (currentParticle.GetDefinition() == aPiMinus)) {
2603       currentParticle.SetKineticEnergy(std::max(0.001 * MeV, dekin * currentParticle.GetKineticEnergy()));
2604       pp = currentParticle.GetTotalMomentum() / MeV;
2605       pp1 = currentParticle.GetMomentum().mag() / MeV;
2606       if (pp1 < 0.001) {
2607         rthnve = pi * G4UniformRand();
2608         phinve = twopi * G4UniformRand();
2609         currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2610                                     pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2611                                     pp * std::cos(rthnve) * MeV);
2612       } else
2613         currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2614     }
2615     for (i = 0; i < vecLen; ++i) {
2616       if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2617           (vec[i]->GetDefinition() == aPiMinus)) {
2618         vec[i]->SetKineticEnergy(std::max(0.001 * MeV, dekin * vec[i]->GetKineticEnergy()));
2619         pp = vec[i]->GetTotalMomentum() / MeV;
2620         pp1 = vec[i]->GetMomentum().mag() / MeV;
2621         if (pp1 < 0.001) {
2622           rthnve = pi * G4UniformRand();
2623           phinve = twopi * G4UniformRand();
2624           vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2625                               pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2626                               pp * std::cos(rthnve) * MeV);
2627         } else
2628           vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2629       }
2630     }
2631   }
2632 }
2633 
2634 void FullModelReactionDynamics::AddBlackTrackParticles(const G4double epnb,  // GeV
2635                                                        const G4int npnb,
2636                                                        const G4double edta,  // GeV
2637                                                        const G4int ndta,
2638                                                        const G4double sprob,
2639                                                        const G4double kineticMinimum,  // GeV
2640                                                        const G4double kineticFactor,   // GeV
2641                                                        const G4ReactionProduct &modifiedOriginal,
2642                                                        G4double spall,
2643                                                        const G4Nucleus &targetNucleus,
2644                                                        G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2645                                                        G4int &vecLen) {
2646   // derived from original FORTRAN code in GENXPT and TWOCLU by H. Fesefeldt
2647   //
2648   // npnb is number of proton/neutron black track particles
2649   // ndta is the number of deuterons, tritons, and alphas produced
2650   // epnb is the kinetic energy available for proton/neutron black track particles
2651   // edta is the kinetic energy available for deuteron/triton/alpha particles
2652   //
2653 
2654   G4ParticleDefinition *aProton = G4Proton::Proton();
2655   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
2656   G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
2657   G4ParticleDefinition *aTriton = G4Triton::Triton();
2658   G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
2659 
2660   const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / MeV;
2661   G4double atomicWeight = targetNucleus.GetN_asInt();
2662   G4double atomicNumber = targetNucleus.GetZ_asInt();
2663 
2664   const G4double ika1 = 3.6;
2665   const G4double ika2 = 35.56;
2666   const G4double ika3 = 6.45;
2667   const G4double sp1 = 1.066;
2668 
2669   G4int i;
2670   G4double pp;
2671   // G4double totalQ = 0;
2672   //G4double kinCreated = 0;
2673   G4double cfa = 0.025 * ((atomicWeight - 1.0) / 120.0) * std::exp(-(atomicWeight - 1.0) / 120.0);
2674   if (npnb > 0)  // first add protons and neutrons
2675   {
2676     G4double backwardKinetic = 0.0;
2677     G4int local_npnb = npnb;
2678     for (i = 0; i < npnb; ++i)
2679       if (G4UniformRand() < sprob)
2680         local_npnb--;
2681     G4double ekin = epnb / local_npnb;
2682 
2683     for (i = 0; i < local_npnb; ++i) {
2684       G4ReactionProduct *p1 = new G4ReactionProduct();
2685       if (backwardKinetic > epnb) {
2686         delete p1;
2687         break;
2688       }
2689       G4double ran = G4UniformRand();
2690       G4double kinetic = -ekin * std::log(ran) - cfa * (1.0 + 0.5 * normal());
2691       if (kinetic < 0.0)
2692         kinetic = -0.010 * std::log(ran);
2693       backwardKinetic += kinetic;
2694       if (backwardKinetic > epnb)
2695         kinetic = std::max(kineticMinimum, epnb - (backwardKinetic - kinetic));
2696       if (G4UniformRand() > (1.0 - atomicNumber / atomicWeight))
2697         p1->SetDefinition(aProton);
2698       else
2699         p1->SetDefinition(aNeutron);
2700       vec.SetElement(vecLen, p1);
2701       ++spall;
2702       G4double cost = G4UniformRand() * 2.0 - 1.0;
2703       G4double sint = std::sqrt(std::fabs(1.0 - cost * cost));
2704       G4double phi = twopi * G4UniformRand();
2705       vec[vecLen]->SetNewlyAdded(true);
2706       vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2707       //kinCreated += kinetic;
2708       pp = vec[vecLen]->GetTotalMomentum() / MeV;
2709       vec[vecLen]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2710       vecLen++;
2711       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2712     }
2713     if ((atomicWeight >= 10.0) && (ekOriginal <= 2.0 * GeV)) {
2714       G4double ekw = ekOriginal / GeV;
2715       G4int ika, kk = 0;
2716       if (ekw > 1.0)
2717         ekw *= ekw;
2718       ekw = std::max(0.1, ekw);
2719       ika = G4int(ika1 * std::exp((atomicNumber * atomicNumber / atomicWeight - ika2) / ika3) / ekw);
2720       if (ika > 0) {
2721         for (i = (vecLen - 1); i >= 0; --i) {
2722           if ((vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded()) {
2723             vec[i]->SetDefinitionAndUpdateE(aNeutron);  // modified 22-Oct-97
2724             if (++kk > ika)
2725               break;
2726           }
2727         }
2728       }
2729     }
2730   }
2731   if (ndta > 0)  //  now, try to add deuterons, tritons and alphas
2732   {
2733     G4double backwardKinetic = 0.0;
2734     G4int local_ndta = ndta;
2735     for (i = 0; i < ndta; ++i)
2736       if (G4UniformRand() < sprob)
2737         local_ndta--;
2738     G4double ekin = edta / local_ndta;
2739 
2740     for (i = 0; i < local_ndta; ++i) {
2741       G4ReactionProduct *p2 = new G4ReactionProduct();
2742       if (backwardKinetic > edta) {
2743         delete p2;
2744         break;
2745       }
2746       G4double ran = G4UniformRand();
2747       G4double kinetic = -ekin * std::log(ran) - cfa * (1. + 0.5 * normal());
2748       if (kinetic < 0.0)
2749         kinetic = kineticFactor * std::log(ran);
2750       backwardKinetic += kinetic;
2751       if (backwardKinetic > edta)
2752         kinetic = edta - (backwardKinetic - kinetic);
2753       if (kinetic < 0.0)
2754         kinetic = kineticMinimum;
2755       G4double cost = 2.0 * G4UniformRand() - 1.0;
2756       G4double sint = std::sqrt(std::max(0.0, (1.0 - cost * cost)));
2757       G4double phi = twopi * G4UniformRand();
2758       ran = G4UniformRand();
2759       if (ran <= 0.60)
2760         p2->SetDefinition(aDeuteron);
2761       else if (ran <= 0.90)
2762         p2->SetDefinition(aTriton);
2763       else
2764         p2->SetDefinition(anAlpha);
2765       spall += p2->GetMass() / GeV * sp1;
2766       if (spall > atomicWeight) {
2767         delete p2;
2768         break;
2769       }
2770       vec.SetElement(vecLen, p2);
2771       vec[vecLen]->SetNewlyAdded(true);
2772       vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2773       //kinCreated += kinetic;
2774       pp = vec[vecLen]->GetTotalMomentum() / MeV;
2775       vec[vecLen++]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2776       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2777     }
2778   }
2779   // G4double delta = epnb+edta - kinCreated;
2780 }
2781 
2782 void FullModelReactionDynamics::MomentumCheck(const G4ReactionProduct &modifiedOriginal,
2783                                               G4ReactionProduct &currentParticle,
2784                                               G4ReactionProduct &targetParticle,
2785                                               G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2786                                               G4int &vecLen) {
2787   const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / MeV;
2788   G4double testMomentum = currentParticle.GetMomentum().mag() / MeV;
2789   G4double pMass;
2790   if (testMomentum >= pOriginal) {
2791     pMass = currentParticle.GetMass() / MeV;
2792     currentParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
2793     currentParticle.SetMomentum(currentParticle.GetMomentum() * (pOriginal / testMomentum));
2794   }
2795   testMomentum = targetParticle.GetMomentum().mag() / MeV;
2796   if (testMomentum >= pOriginal) {
2797     pMass = targetParticle.GetMass() / MeV;
2798     targetParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
2799     targetParticle.SetMomentum(targetParticle.GetMomentum() * (pOriginal / testMomentum));
2800   }
2801   for (G4int i = 0; i < vecLen; ++i) {
2802     testMomentum = vec[i]->GetMomentum().mag() / MeV;
2803     if (testMomentum >= pOriginal) {
2804       pMass = vec[i]->GetMass() / MeV;
2805       vec[i]->SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
2806       vec[i]->SetMomentum(vec[i]->GetMomentum() * (pOriginal / testMomentum));
2807     }
2808   }
2809 }
2810 
2811 void FullModelReactionDynamics::ProduceStrangeParticlePairs(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2812                                                             G4int &vecLen,
2813                                                             const G4ReactionProduct &modifiedOriginal,
2814                                                             const G4DynamicParticle *originalTarget,
2815                                                             G4ReactionProduct &currentParticle,
2816                                                             G4ReactionProduct &targetParticle,
2817                                                             G4bool &incidentHasChanged,
2818                                                             G4bool &targetHasChanged) {
2819   // derived from original FORTRAN code STPAIR by H. Fesefeldt (16-Dec-1987)
2820   //
2821   // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
2822   //                            K+ Y0, K0 Y+,  K0 Y-
2823   // For antibaryon induced reactions half of the cross sections KB YB
2824   // pairs are produced.  Charge is not conserved, no experimental data available
2825   // for exclusive reactions, therefore some average behaviour assumed.
2826   // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
2827   //
2828   if (vecLen == 0)
2829     return;
2830   //
2831   // the following protects against annihilation processes
2832   //
2833   if (currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0)
2834     return;
2835 
2836   const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
2837   const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass() / GeV;
2838   G4double targetMass = originalTarget->GetDefinition()->GetPDGMass() / GeV;
2839   G4double centerofmassEnergy =
2840       std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal);  // GeV
2841   G4double currentMass = currentParticle.GetMass() / GeV;
2842   G4double availableEnergy = centerofmassEnergy - (targetMass + currentMass);
2843   if (availableEnergy <= 1.0)
2844     return;
2845 
2846   G4ParticleDefinition *aProton = G4Proton::Proton();
2847   G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
2848   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
2849   G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
2850   G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
2851   G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
2852   G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
2853   G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
2854   G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
2855   G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
2856   G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
2857   G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
2858   G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
2859   G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
2860   G4ParticleDefinition *aLambda = G4Lambda::Lambda();
2861   G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
2862 
2863   const G4double protonMass = aProton->GetPDGMass() / GeV;
2864   const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass() / GeV;
2865   //
2866   // determine the center of mass energy bin
2867   //
2868   const G4double avrs[] = {3., 4., 5., 6., 7., 8., 9., 10., 20., 30., 40., 50.};
2869 
2870   G4int ibin, i3, i4;
2871   G4double avk, avy, avn, ran;
2872   G4int i = 1;
2873   while ((i < 12) && (centerofmassEnergy > avrs[i]))
2874     ++i;
2875   if (i == 12)
2876     ibin = 11;
2877   else
2878     ibin = i;
2879   //
2880   // the fortran code chooses a random replacement of produced kaons
2881   //  but does not take into account charge conservation
2882   //
2883   if (vecLen == 1)  // we know that vecLen > 0
2884   {
2885     i3 = 0;
2886     i4 = 1;  // note that we will be adding a new secondary particle in this case only
2887   } else     // otherwise  0 <= i3,i4 < vecLen
2888   {
2889     G4double ran = G4UniformRand();
2890     while (ran == 1.0)
2891       ran = G4UniformRand();
2892     i4 = i3 = G4int(vecLen * ran);
2893     while (i3 == i4) {
2894       ran = G4UniformRand();
2895       while (ran == 1.0)
2896         ran = G4UniformRand();
2897       i4 = G4int(vecLen * ran);
2898     }
2899   }
2900   //
2901   // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
2902   //
2903   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};
2904   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};
2905   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};
2906 
2907   avk = (std::log(avkkb[ibin]) - std::log(avkkb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
2908             (avrs[ibin] - avrs[ibin - 1]) +
2909         std::log(avkkb[ibin - 1]);
2910   avk = std::exp(avk);
2911 
2912   avy = (std::log(avky[ibin]) - std::log(avky[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
2913             (avrs[ibin] - avrs[ibin - 1]) +
2914         std::log(avky[ibin - 1]);
2915   avy = std::exp(avy);
2916 
2917   avn = (std::log(avnnb[ibin]) - std::log(avnnb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
2918             (avrs[ibin] - avrs[ibin - 1]) +
2919         std::log(avnnb[ibin - 1]);
2920   avn = std::exp(avn);
2921 
2922   if (avk + avy + avn <= 0.0)
2923     return;
2924 
2925   if (currentMass < protonMass)
2926     avy /= 2.0;
2927   if (targetMass < protonMass)
2928     avy = 0.0;
2929   avy += avk + avn;
2930   avk += avn;
2931   ran = G4UniformRand();
2932   if (ran < avn) {
2933     if (availableEnergy < 2.0)
2934       return;
2935     if (vecLen == 1)  // add a new secondary
2936     {
2937       G4ReactionProduct *p1 = new G4ReactionProduct;
2938       if (G4UniformRand() < 0.5) {
2939         vec[0]->SetDefinition(aNeutron);
2940         p1->SetDefinition(anAntiNeutron);
2941         (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
2942         vec[0]->SetMayBeKilled(false);
2943         p1->SetMayBeKilled(false);
2944       } else {
2945         vec[0]->SetDefinition(aProton);
2946         p1->SetDefinition(anAntiProton);
2947         (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
2948         vec[0]->SetMayBeKilled(false);
2949         p1->SetMayBeKilled(false);
2950       }
2951       vec.SetElement(vecLen++, p1);
2952       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
2953     } else {  // replace two secondaries
2954       if (G4UniformRand() < 0.5) {
2955         vec[i3]->SetDefinition(aNeutron);
2956         vec[i4]->SetDefinition(anAntiNeutron);
2957         vec[i3]->SetMayBeKilled(false);
2958         vec[i4]->SetMayBeKilled(false);
2959       } else {
2960         vec[i3]->SetDefinition(aProton);
2961         vec[i4]->SetDefinition(anAntiProton);
2962         vec[i3]->SetMayBeKilled(false);
2963         vec[i4]->SetMayBeKilled(false);
2964       }
2965     }
2966   } else if (ran < avk) {
2967     if (availableEnergy < 1.0)
2968       return;
2969 
2970     const G4double kkb[] = {0.2500, 0.3750, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8750, 1.000};
2971     const G4int ipakkb1[] = {10, 10, 10, 11, 11, 12, 12, 11, 12};
2972     const G4int ipakkb2[] = {13, 11, 12, 11, 12, 11, 12, 13, 13};
2973     ran = G4UniformRand();
2974     i = 0;
2975     while ((i < 9) && (ran >= kkb[i]))
2976       ++i;
2977     if (i == 9)
2978       return;
2979     //
2980     // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
2981     // charge       +  -   +  0   +  0   0  0   0  0   0  0   0  0   0  -   0  -
2982     //
2983     switch (ipakkb1[i]) {
2984       case 10:
2985         vec[i3]->SetDefinition(aKaonPlus);
2986         vec[i3]->SetMayBeKilled(false);
2987         break;
2988       case 11:
2989         vec[i3]->SetDefinition(aKaonZS);
2990         vec[i3]->SetMayBeKilled(false);
2991         break;
2992       case 12:
2993         vec[i3]->SetDefinition(aKaonZL);
2994         vec[i3]->SetMayBeKilled(false);
2995         break;
2996     }
2997     if (vecLen == 1)  // add a secondary
2998     {
2999       G4ReactionProduct *p1 = new G4ReactionProduct;
3000       switch (ipakkb2[i]) {
3001         case 11:
3002           p1->SetDefinition(aKaonZS);
3003           p1->SetMayBeKilled(false);
3004           break;
3005         case 12:
3006           p1->SetDefinition(aKaonZL);
3007           p1->SetMayBeKilled(false);
3008           break;
3009         case 13:
3010           p1->SetDefinition(aKaonMinus);
3011           p1->SetMayBeKilled(false);
3012           break;
3013       }
3014       (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3015       vec.SetElement(vecLen++, p1);
3016       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
3017     } else  // replace
3018     {
3019       switch (ipakkb2[i]) {
3020         case 11:
3021           vec[i4]->SetDefinition(aKaonZS);
3022           vec[i4]->SetMayBeKilled(false);
3023           break;
3024         case 12:
3025           vec[i4]->SetDefinition(aKaonZL);
3026           vec[i4]->SetMayBeKilled(false);
3027           break;
3028         case 13:
3029           vec[i4]->SetDefinition(aKaonMinus);
3030           vec[i4]->SetMayBeKilled(false);
3031           break;
3032       }
3033     }
3034   } else if (ran < avy) {
3035     if (availableEnergy < 1.6)
3036       return;
3037 
3038     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};
3039     const G4int ipaky1[] = {18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22};
3040     const G4int ipaky2[] = {10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12};
3041     const G4int ipakyb1[] = {19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25};
3042     const G4int ipakyb2[] = {13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11};
3043     ran = G4UniformRand();
3044     i = 0;
3045     while ((i < 12) && (ran > ky[i]))
3046       ++i;
3047     if (i == 12)
3048       return;
3049     if ((currentMass < protonMass) || (G4UniformRand() < 0.5)) {
3050       switch (ipaky1[i]) {
3051         case 18:
3052           targetParticle.SetDefinition(aLambda);
3053           break;
3054         case 20:
3055           targetParticle.SetDefinition(aSigmaPlus);
3056           break;
3057         case 21:
3058           targetParticle.SetDefinition(aSigmaZero);
3059           break;
3060         case 22:
3061           targetParticle.SetDefinition(aSigmaMinus);
3062           break;
3063       }
3064       targetHasChanged = true;
3065       switch (ipaky2[i]) {
3066         case 10:
3067           vec[i3]->SetDefinition(aKaonPlus);
3068           vec[i3]->SetMayBeKilled(false);
3069           break;
3070         case 11:
3071           vec[i3]->SetDefinition(aKaonZS);
3072           vec[i3]->SetMayBeKilled(false);
3073           break;
3074         case 12:
3075           vec[i3]->SetDefinition(aKaonZL);
3076           vec[i3]->SetMayBeKilled(false);
3077           break;
3078       }
3079     } else  // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
3080     {
3081       if ((currentParticle.GetDefinition() == anAntiProton) || (currentParticle.GetDefinition() == anAntiNeutron) ||
3082           (currentParticle.GetDefinition() == anAntiLambda) || (currentMass > sigmaMinusMass)) {
3083         switch (ipakyb1[i]) {
3084           case 19:
3085             currentParticle.SetDefinitionAndUpdateE(anAntiLambda);
3086             break;
3087           case 23:
3088             currentParticle.SetDefinitionAndUpdateE(anAntiSigmaPlus);
3089             break;
3090           case 24:
3091             currentParticle.SetDefinitionAndUpdateE(anAntiSigmaZero);
3092             break;
3093           case 25:
3094             currentParticle.SetDefinitionAndUpdateE(anAntiSigmaMinus);
3095             break;
3096         }
3097         incidentHasChanged = true;
3098         switch (ipakyb2[i]) {
3099           case 11:
3100             vec[i3]->SetDefinition(aKaonZS);
3101             vec[i3]->SetMayBeKilled(false);
3102             break;
3103           case 12:
3104             vec[i3]->SetDefinition(aKaonZL);
3105             vec[i3]->SetMayBeKilled(false);
3106             break;
3107           case 13:
3108             vec[i3]->SetDefinition(aKaonMinus);
3109             vec[i3]->SetMayBeKilled(false);
3110             break;
3111         }
3112       } else {
3113         switch (ipaky1[i]) {
3114           case 18:
3115             currentParticle.SetDefinitionAndUpdateE(aLambda);
3116             break;
3117           case 20:
3118             currentParticle.SetDefinitionAndUpdateE(aSigmaPlus);
3119             break;
3120           case 21:
3121             currentParticle.SetDefinitionAndUpdateE(aSigmaZero);
3122             break;
3123           case 22:
3124             currentParticle.SetDefinitionAndUpdateE(aSigmaMinus);
3125             break;
3126         }
3127         incidentHasChanged = true;
3128         switch (ipaky2[i]) {
3129           case 10:
3130             vec[i3]->SetDefinition(aKaonPlus);
3131             vec[i3]->SetMayBeKilled(false);
3132             break;
3133           case 11:
3134             vec[i3]->SetDefinition(aKaonZS);
3135             vec[i3]->SetMayBeKilled(false);
3136             break;
3137           case 12:
3138             vec[i3]->SetDefinition(aKaonZL);
3139             vec[i3]->SetMayBeKilled(false);
3140             break;
3141         }
3142       }
3143     }
3144   } else
3145     return;
3146   //
3147   //  check the available energy
3148   //   if there is not enough energy for kkb/ky pair production
3149   //   then reduce the number of secondary particles
3150   //  NOTE:
3151   //        the number of secondaries may have been changed
3152   //        the incident and/or target particles may have changed
3153   //        charge conservation is ignored (as well as strangness conservation)
3154   //
3155   currentMass = currentParticle.GetMass() / GeV;
3156   targetMass = targetParticle.GetMass() / GeV;
3157 
3158   G4double energyCheck = centerofmassEnergy - (currentMass + targetMass);
3159   for (i = 0; i < vecLen; ++i) {
3160     energyCheck -= vec[i]->GetMass() / GeV;
3161     if (energyCheck < 0.0)  // chop off the secondary List
3162     {
3163       vecLen = std::max(0, --i);  // looks like a memory leak @@@@@@@@@@@@
3164       G4int j;
3165       for (j = i; j < vecLen; j++)
3166         delete vec[j];
3167       break;
3168     }
3169   }
3170   return;
3171 }
3172 
3173 void FullModelReactionDynamics::NuclearReaction(G4FastVector<G4ReactionProduct, 4> &vec,
3174                                                 G4int &vecLen,
3175                                                 const G4HadProjectile *originalIncident,
3176                                                 const G4Nucleus &targetNucleus,
3177                                                 const G4double theAtomicMass,
3178                                                 const G4double *mass) {
3179   // derived from original FORTRAN code NUCREC by H. Fesefeldt (12-Feb-1987)
3180   //
3181   // Nuclear reaction kinematics at low energies
3182   //
3183   G4ParticleDefinition *aGamma = G4Gamma::Gamma();
3184   G4ParticleDefinition *aProton = G4Proton::Proton();
3185   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3186   G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
3187   G4ParticleDefinition *aTriton = G4Triton::Triton();
3188   G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
3189 
3190   const G4double aProtonMass = aProton->GetPDGMass() / MeV;
3191   const G4double aNeutronMass = aNeutron->GetPDGMass() / MeV;
3192   const G4double aDeuteronMass = aDeuteron->GetPDGMass() / MeV;
3193   const G4double aTritonMass = aTriton->GetPDGMass() / MeV;
3194   const G4double anAlphaMass = anAlpha->GetPDGMass() / MeV;
3195 
3196   G4ReactionProduct currentParticle;
3197   currentParticle = *originalIncident;
3198   //
3199   // Set beam particle, take kinetic energy of current particle as the
3200   // fundamental quantity.  Due to the difficult kinematic, all masses have to
3201   // be assigned the best measured values
3202   //
3203   G4double p = currentParticle.GetTotalMomentum();
3204   G4double pp = currentParticle.GetMomentum().mag();
3205   if (pp <= 0.001 * MeV) {
3206     G4double phinve = twopi * G4UniformRand();
3207     G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3208     currentParticle.SetMomentum(
3209         p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3210   } else
3211     currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3212   //
3213   // calculate Q-value of reactions
3214   //
3215   G4double currentKinetic = currentParticle.GetKineticEnergy() / MeV;
3216   G4double currentMass = currentParticle.GetDefinition()->GetPDGMass() / MeV;
3217   G4double qv = currentKinetic + theAtomicMass + currentMass;
3218 
3219   G4double qval[9];
3220   qval[0] = qv - mass[0];
3221   qval[1] = qv - mass[1] - aNeutronMass;
3222   qval[2] = qv - mass[2] - aProtonMass;
3223   qval[3] = qv - mass[3] - aDeuteronMass;
3224   qval[4] = qv - mass[4] - aTritonMass;
3225   qval[5] = qv - mass[5] - anAlphaMass;
3226   qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
3227   qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
3228   qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
3229 
3230   if (currentParticle.GetDefinition() == aNeutron) {
3231     const G4double A = targetNucleus.GetN_asInt();  // atomic weight
3232     if (G4UniformRand() > ((A - 1.0) / 230.0) * ((A - 1.0) / 230.0))
3233       qval[0] = 0.0;
3234     if (G4UniformRand() >= currentKinetic / 7.9254 * A)
3235       qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
3236   } else
3237     qval[0] = 0.0;
3238 
3239   G4int i;
3240   qv = 0.0;
3241   for (i = 0; i < 9; ++i) {
3242     if (mass[i] < 500.0 * MeV)
3243       qval[i] = 0.0;
3244     if (qval[i] < 0.0)
3245       qval[i] = 0.0;
3246     qv += qval[i];
3247   }
3248   G4double qv1 = 0.0;
3249   G4double ran = G4UniformRand();
3250   G4int index;
3251   for (index = 0; index < 9; ++index) {
3252     if (qval[index] > 0.0) {
3253       qv1 += qval[index] / qv;
3254       if (ran <= qv1)
3255         break;
3256     }
3257   }
3258   if (index == 9)  // loop continued to the end
3259   {
3260     throw G4HadronicException(
3261         __FILE__,
3262         __LINE__,
3263         "FullModelReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
3264   }
3265   G4double ke = currentParticle.GetKineticEnergy() / GeV;
3266   G4int nt = 2;
3267   if ((index >= 6) || (G4UniformRand() < std::min(0.5, ke * 10.0)))
3268     nt = 3;
3269 
3270   G4ReactionProduct **v = new G4ReactionProduct *[3];
3271   v[0] = new G4ReactionProduct;
3272   v[1] = new G4ReactionProduct;
3273   v[2] = new G4ReactionProduct;
3274 
3275   v[0]->SetMass(mass[index] * MeV);
3276   switch (index) {
3277     case 0:
3278       v[1]->SetDefinition(aGamma);
3279       v[2]->SetDefinition(aGamma);
3280       break;
3281     case 1:
3282       v[1]->SetDefinition(aNeutron);
3283       v[2]->SetDefinition(aGamma);
3284       break;
3285     case 2:
3286       v[1]->SetDefinition(aProton);
3287       v[2]->SetDefinition(aGamma);
3288       break;
3289     case 3:
3290       v[1]->SetDefinition(aDeuteron);
3291       v[2]->SetDefinition(aGamma);
3292       break;
3293     case 4:
3294       v[1]->SetDefinition(aTriton);
3295       v[2]->SetDefinition(aGamma);
3296       break;
3297     case 5:
3298       v[1]->SetDefinition(anAlpha);
3299       v[2]->SetDefinition(aGamma);
3300       break;
3301     case 6:
3302       v[1]->SetDefinition(aNeutron);
3303       v[2]->SetDefinition(aNeutron);
3304       break;
3305     case 7:
3306       v[1]->SetDefinition(aNeutron);
3307       v[2]->SetDefinition(aProton);
3308       break;
3309     case 8:
3310       v[1]->SetDefinition(aProton);
3311       v[2]->SetDefinition(aProton);
3312       break;
3313   }
3314   //
3315   // calculate centre of mass energy
3316   //
3317   G4ReactionProduct pseudo1;
3318   pseudo1.SetMass(theAtomicMass * MeV);
3319   pseudo1.SetTotalEnergy(theAtomicMass * MeV);
3320   G4ReactionProduct pseudo2 = currentParticle + pseudo1;
3321   pseudo2.SetMomentum(pseudo2.GetMomentum() * (-1.0));
3322   //
3323   // use phase space routine in centre of mass system
3324   //
3325   G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
3326   tempV.Initialize(nt);
3327   G4int tempLen = 0;
3328   tempV.SetElement(tempLen++, v[0]);
3329   tempV.SetElement(tempLen++, v[1]);
3330   if (nt == 3)
3331     tempV.SetElement(tempLen++, v[2]);
3332   G4bool constantCrossSection = true;
3333   GenerateNBodyEvent(pseudo2.GetMass() / MeV, constantCrossSection, tempV, tempLen);
3334   v[0]->Lorentz(*v[0], pseudo2);
3335   v[1]->Lorentz(*v[1], pseudo2);
3336   if (nt == 3)
3337     v[2]->Lorentz(*v[2], pseudo2);
3338 
3339   G4bool particleIsDefined = false;
3340   if (v[0]->GetMass() / MeV - aProtonMass < 0.1) {
3341     v[0]->SetDefinition(aProton);
3342     particleIsDefined = true;
3343   } else if (v[0]->GetMass() / MeV - aNeutronMass < 0.1) {
3344     v[0]->SetDefinition(aNeutron);
3345     particleIsDefined = true;
3346   } else if (v[0]->GetMass() / MeV - aDeuteronMass < 0.1) {
3347     v[0]->SetDefinition(aDeuteron);
3348     particleIsDefined = true;
3349   } else if (v[0]->GetMass() / MeV - aTritonMass < 0.1) {
3350     v[0]->SetDefinition(aTriton);
3351     particleIsDefined = true;
3352   } else if (v[0]->GetMass() / MeV - anAlphaMass < 0.1) {
3353     v[0]->SetDefinition(anAlpha);
3354     particleIsDefined = true;
3355   }
3356   currentParticle.SetKineticEnergy(std::max(0.001, currentParticle.GetKineticEnergy() / MeV));
3357   p = currentParticle.GetTotalMomentum();
3358   pp = currentParticle.GetMomentum().mag();
3359   if (pp <= 0.001 * MeV) {
3360     G4double phinve = twopi * G4UniformRand();
3361     G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3362     currentParticle.SetMomentum(
3363         p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3364   } else
3365     currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3366 
3367   if (particleIsDefined) {
3368     v[0]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[0]->GetKineticEnergy() / MeV));
3369     p = v[0]->GetTotalMomentum();
3370     pp = v[0]->GetMomentum().mag();
3371     if (pp <= 0.001 * MeV) {
3372       G4double phinve = twopi * G4UniformRand();
3373       G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3374       v[0]->SetMomentum(
3375           p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3376     } else
3377       v[0]->SetMomentum(v[0]->GetMomentum() * (p / pp));
3378   }
3379   if ((v[1]->GetDefinition() == aDeuteron) || (v[1]->GetDefinition() == aTriton) || (v[1]->GetDefinition() == anAlpha))
3380     v[1]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[1]->GetKineticEnergy() / MeV));
3381   else
3382     v[1]->SetKineticEnergy(std::max(0.001, v[1]->GetKineticEnergy() / MeV));
3383 
3384   p = v[1]->GetTotalMomentum();
3385   pp = v[1]->GetMomentum().mag();
3386   if (pp <= 0.001 * MeV) {
3387     G4double phinve = twopi * G4UniformRand();
3388     G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3389     v[1]->SetMomentum(
3390         p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3391   } else
3392     v[1]->SetMomentum(v[1]->GetMomentum() * (p / pp));
3393 
3394   if (nt == 3) {
3395     if ((v[2]->GetDefinition() == aDeuteron) || (v[2]->GetDefinition() == aTriton) ||
3396         (v[2]->GetDefinition() == anAlpha))
3397       v[2]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[2]->GetKineticEnergy() / MeV));
3398     else
3399       v[2]->SetKineticEnergy(std::max(0.001, v[2]->GetKineticEnergy() / MeV));
3400 
3401     p = v[2]->GetTotalMomentum();
3402     pp = v[2]->GetMomentum().mag();
3403     if (pp <= 0.001 * MeV) {
3404       G4double phinve = twopi * G4UniformRand();
3405       G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3406       v[2]->SetMomentum(
3407           p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3408     } else
3409       v[2]->SetMomentum(v[2]->GetMomentum() * (p / pp));
3410   }
3411   G4int del;
3412   for (del = 0; del < vecLen; del++)
3413     delete vec[del];
3414   vecLen = 0;
3415   if (particleIsDefined) {
3416     vec.SetElement(vecLen++, v[0]);
3417   } else {
3418     delete v[0];
3419   }
3420   vec.SetElement(vecLen++, v[1]);
3421   if (nt == 3) {
3422     vec.SetElement(vecLen++, v[2]);
3423   } else {
3424     delete v[2];
3425   }
3426   delete[] v;
3427   return;
3428 }
3429 
3430 /* end of file */
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!