Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4HadronElastic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/coherent_elastic/src/G4HadronElastic.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4HadronElastic.cc (Version 8.1.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4HadronElastic.cc,v 1.30 2006/10/31 18:36:04 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-08-01-patch-02 $
 26 //                                                 28 //
 27 // Geant4 Header : G4HadronElastic             << 
 28 //                                                 29 //
 29 // Author : V.Ivanchenko 29 June 2009 (redesig <<  30 // Physics model class G4HadronElastic (derived from G4LElastic)
 30 //                                             <<  31 //
                                                   >>  32 //
                                                   >>  33 // G4 Model: Low-energy Elastic scattering with 4-momentum balance
                                                   >>  34 // F.W. Jones, TRIUMF, 04-JUN-96
                                                   >>  35 // Uses  G4ElasticHadrNucleusHE and G4VQCrossSection
                                                   >>  36 //
                                                   >>  37 //
                                                   >>  38 // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange.
                                                   >>  39 // 09-Set-05 V.Ivanchenko HARP version of the model: fix scattering
                                                   >>  40 //           on hydrogen, use relativistic Lorentz transformation
                                                   >>  41 // 24-Nov-05 V.Ivanchenko sample cost in center of mass reference system
                                                   >>  42 // 03-Dec-05 V.Ivanchenko add protection to initial momentum 20 MeV/c in
                                                   >>  43 //           center of mass system (before it was in lab system)
                                                   >>  44 //           below model is not valid
                                                   >>  45 // 14-Dec-05 V.Ivanchenko change protection to cos(theta) < -1 and
                                                   >>  46 //           rename the class
                                                   >>  47 // 13-Apr-06 V.Ivanchenko move to coherent_elastic subdirectory; remove
                                                   >>  48 //           charge exchange; remove limitation on incident momentum;
                                                   >>  49 //           add s-wave regim below some momentum        
                                                   >>  50 // 24-Apr-06 V.Ivanchenko add neutron scattering on hydrogen from CHIPS
                                                   >>  51 // 07-Jun-06 V.Ivanchenko fix problem of rotation
                                                   >>  52 // 25-Jul-06 V.Ivanchenko add 19 MeV low energy, below which S-wave is sampled
                                                   >>  53 // 02-Aug-06 V.Ivanchenko introduce energy cut on the aria of S-wave for pions
                                                   >>  54 // 24-Aug-06 V.Ivanchenko switch on G4ElasticHadrNucleusHE
                                                   >>  55 // 31-Aug-06 V.Ivanchenko do not sample sacttering for particles with kinetic 
                                                   >>  56 //                        energy below 10 keV
                                                   >>  57 //
 31                                                    58 
 32 #include "G4HadronElastic.hh"                      59 #include "G4HadronElastic.hh"
 33 #include "G4SystemOfUnits.hh"                  << 
 34 #include "G4ParticleTable.hh"                      60 #include "G4ParticleTable.hh"
 35 #include "G4ParticleDefinition.hh"                 61 #include "G4ParticleDefinition.hh"
 36 #include "G4IonTable.hh"                           62 #include "G4IonTable.hh"
                                                   >>  63 #include "G4QElasticCrossSection.hh"
                                                   >>  64 #include "G4VQCrossSection.hh"
                                                   >>  65 #include "G4ElasticHadrNucleusHE.hh"
 37 #include "Randomize.hh"                            66 #include "Randomize.hh"
 38 #include "G4Proton.hh"                             67 #include "G4Proton.hh"
 39 #include "G4Neutron.hh"                            68 #include "G4Neutron.hh"
 40 #include "G4Deuteron.hh"                           69 #include "G4Deuteron.hh"
 41 #include "G4Alpha.hh"                              70 #include "G4Alpha.hh"
 42 #include "G4Pow.hh"                            <<  71 #include "G4PionPlus.hh"
 43 #include "G4Exp.hh"                            <<  72 #include "G4PionMinus.hh"
 44 #include "G4Log.hh"                            << 
 45 #include "G4HadronicParameters.hh"             << 
 46 #include "G4PhysicsModelCatalog.hh"            << 
 47                                                << 
 48                                                    73 
 49 G4HadronElastic::G4HadronElastic(const G4Strin <<  74 G4HadronElastic::G4HadronElastic(G4double plow, G4double elim, G4double ehigh) 
 50   : G4HadronicInteraction(name), secID(-1)     <<  75   : G4HadronicInteraction()
 51 {                                                  76 {
 52   SetMinEnergy( 0.0*GeV );                         77   SetMinEnergy( 0.0*GeV );
 53   SetMaxEnergy( G4HadronicParameters::Instance <<  78   SetMaxEnergy( 100.*TeV );
 54   lowestEnergyLimit= 1.e-6*eV;                 <<  79   verboseLevel= 0;
 55   pLocalTmax  = 0.0;                           <<  80   ekinIon     = elim;
 56   nwarn = 0;                                   <<  81   ekinhigh    = ehigh;
                                                   >>  82   ekinlow     = 19.*MeV;
                                                   >>  83   ekinpi      = MeV;
                                                   >>  84   plablow     = plow;
                                                   >>  85 
                                                   >>  86   qCManager   = G4QElasticCrossSection::GetPointer();
                                                   >>  87   hElastic    = new G4ElasticHadrNucleusHE();
 57                                                    88 
 58   theProton   = G4Proton::Proton();                89   theProton   = G4Proton::Proton();
 59   theNeutron  = G4Neutron::Neutron();              90   theNeutron  = G4Neutron::Neutron();
 60   theDeuteron = G4Deuteron::Deuteron();            91   theDeuteron = G4Deuteron::Deuteron();
 61   theAlpha    = G4Alpha::Alpha();                  92   theAlpha    = G4Alpha::Alpha();
 62                                                <<  93   thePionPlus = G4PionPlus::PionPlus();
 63   secID = G4PhysicsModelCatalog::GetModelID( " <<  94   thePionMinus= G4PionMinus::PionMinus();
 64 }                                                  95 }
 65                                                    96 
 66 G4HadronElastic::~G4HadronElastic()                97 G4HadronElastic::~G4HadronElastic()
 67 {}                                             <<  98 {
                                                   >>  99   delete hElastic;
                                                   >> 100 }
 68                                                   101 
                                                   >> 102 G4VQCrossSection* G4HadronElastic::GetCS()
                                                   >> 103 {
                                                   >> 104   return qCManager;
                                                   >> 105 }
 69                                                   106 
 70 void G4HadronElastic::ModelDescription(std::os << 107 G4ElasticHadrNucleusHE* G4HadronElastic::GetHElastic()
 71 {                                                 108 {
 72   outFile << "G4HadronElastic is the base clas << 109   return hElastic;
 73           << "elastic scattering models except << 
 74           << "By default it uses the Gheisha t << 
 75     << "transfer parameterization.  The model  << 
 76     << "as opposed to the original Gheisha mod << 
 77     << "This model may be used for all long-li << 
 78     << "incident energies but fit the data onl << 
 79 }                                                 110 }
 80                                                   111 
 81 G4HadFinalState* G4HadronElastic::ApplyYoursel    112 G4HadFinalState* G4HadronElastic::ApplyYourself(
 82      const G4HadProjectile& aTrack, G4Nucleus&    113      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
 83 {                                                 114 {
 84   theParticleChange.Clear();                      115   theParticleChange.Clear();
 85                                                   116 
 86   const G4HadProjectile* aParticle = &aTrack;     117   const G4HadProjectile* aParticle = &aTrack;
 87   G4double ekin = aParticle->GetKineticEnergy(    118   G4double ekin = aParticle->GetKineticEnergy();
 88                                                << 119   if(ekin <= 10.0*keV) {
 89   // no scattering below the limit             << 
 90   if(ekin <= lowestEnergyLimit) {              << 
 91     theParticleChange.SetEnergyChange(ekin);      120     theParticleChange.SetEnergyChange(ekin);
 92     theParticleChange.SetMomentumChange(0.,0., << 121     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
 93     return &theParticleChange;                    122     return &theParticleChange;
 94   }                                               123   }
 95                                                   124 
 96   G4int A = targetNucleus.GetA_asInt();        << 125   G4double aTarget = targetNucleus.GetN();
 97   G4int Z = targetNucleus.GetZ_asInt();        << 126   G4double zTarget = targetNucleus.GetZ();
                                                   >> 127 
                                                   >> 128   G4double plab = aParticle->GetTotalMomentum();
                                                   >> 129   if (verboseLevel >1) 
                                                   >> 130     G4cout << "G4HadronElastic::DoIt: Incident particle plab=" 
                                                   >> 131      << plab/GeV << " GeV/c " 
                                                   >> 132      << " ekin(MeV) = " << ekin/MeV << "  " 
                                                   >> 133      << aParticle->GetDefinition()->GetParticleName() << G4endl;
 98                                                   134 
 99   // Scattered particle referred to axis of in    135   // Scattered particle referred to axis of incident particle
100   const G4ParticleDefinition* theParticle = aP    136   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
101   G4double m1 = theParticle->GetPDGMass();        137   G4double m1 = theParticle->GetPDGMass();
102   G4double plab = std::sqrt(ekin*(ekin + 2.0*m << 
103                                                   138 
104   if (verboseLevel>1) {                        << 139   G4int Z = static_cast<G4int>(zTarget);
105     G4cout << "G4HadronElastic: "              << 140   G4int A = static_cast<G4int>(aTarget);
106      << aParticle->GetDefinition()->GetParticl << 141   G4int N = A - Z;
107      << " Plab(GeV/c)= " << plab/GeV           << 142   G4int projPDG = theParticle->GetPDGEncoding();
108      << " Ekin(MeV) = " << ekin/MeV            << 143   if (verboseLevel>1) 
109      << " scattered off Z= " << Z              << 144     G4cout << "G4HadronElastic for " << theParticle->GetParticleName()
110      << " A= " << A                            << 145      << " PDGcode= " << projPDG << " on nucleus Z= " << Z 
                                                   >> 146      << " A= " << A << " N= " << N 
111      << G4endl;                                   147      << G4endl;
112   }                                            << 
113                                                   148 
114   G4double mass2 = G4NucleiProperties::GetNucl << 149   G4ParticleDefinition * theDef = 0;
115   G4double e1 = m1 + ekin;                     << 150 
116   G4LorentzVector lv(0.0,0.0,plab,e1+mass2);   << 151   if(Z == 1 && A == 1)       theDef = G4Proton::Proton();
                                                   >> 152   else if (Z == 1 && A == 2) theDef = G4Deuteron::Deuteron();
                                                   >> 153   else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
                                                   >> 154   else if (Z == 2 && A == 3) theDef = G4He3::He3();
                                                   >> 155   else if (Z == 2 && A == 4) theDef = G4Alpha::Alpha();
                                                   >> 156   else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
                                                   >> 157  
                                                   >> 158   G4double m2 = theDef->GetPDGMass();
                                                   >> 159   G4LorentzVector lv1 = aParticle->Get4Momentum();
                                                   >> 160   G4LorentzVector lv(0.0,0.0,0.0,m2);   
                                                   >> 161   lv += lv1;
                                                   >> 162 
117   G4ThreeVector bst = lv.boostVector();           163   G4ThreeVector bst = lv.boostVector();
118   G4double momentumCMS = plab*mass2/std::sqrt( << 164   lv1.boost(-bst);
119                                                   165 
120   pLocalTmax = 4.0*momentumCMS*momentumCMS;    << 166   G4ThreeVector p1 = lv1.vect();
                                                   >> 167   G4double ptot = p1.mag();
                                                   >> 168   G4double tmax = 4.0*ptot*ptot;
                                                   >> 169   G4double t = 0.0;
                                                   >> 170 
                                                   >> 171   // Choose generator
                                                   >> 172   G4ElasticGenerator gtype = fLElastic;
                                                   >> 173 
                                                   >> 174   // CHIPS Elastic
                                                   >> 175   if ((theParticle == theProton || theParticle == theNeutron) && 
                                                   >> 176       Z <= 2 && ekin > ekinlow) {
                                                   >> 177     gtype = fQElastic;
                                                   >> 178   } else {
121                                                   179 
122   // Sampling in CM system                     << 180     // High energy elastic 
123   G4double t = SampleInvariantT(theParticle, p << 181     if(ekin >= ekinhigh && A != 2 && A != 3 && A < 238) gtype = fHElastic;
124                                                   182 
125   if(t < 0.0 || t > pLocalTmax) {              << 183     // Simplification - s-wave
126     // For the very rare cases where cos(theta << 184     else if((theParticle == thePionPlus || 
127     // print some debugging information via a  << 185        theParticle == thePionMinus) && 
128     // using the default algorithm             << 186       (ekin <= ekinpi) ) 
129 #ifdef G4VERBOSE                               << 187       gtype = fSWave;
130     if(nwarn < 2) {                            << 188     else if(plab <= plablow) gtype = fSWave;
131       G4ExceptionDescription ed;               << 189   }
132       ed << GetModelName() << " wrong sampling << 190 
133    << " for " << aParticle->GetDefinition()->G << 191   // Sample t
134    << " ekin=" << ekin << " MeV"               << 192   if(gtype == fQElastic) {
135    << " off (Z,A)=(" << Z << "," << A << ") -  << 193     if (verboseLevel >1) 
136       G4Exception( "G4HadronElastic::ApplyYour << 194       G4cout << "G4HadronElastic: Z= " << Z << " N= " 
137       ++nwarn;                                 << 195        << N << " pdg= " <<  projPDG
                                                   >> 196        << " mom(GeV)= " << plab/GeV << "  " << qCManager << G4endl; 
                                                   >> 197     if(Z == 1 && N == 2) N = 1;
                                                   >> 198     else if(Z == 2 && N == 1) N = 2;
                                                   >> 199     G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
                                                   >> 200     if(cs > 0.0) t = qCManager->GetExchangeT(Z,N,projPDG);
                                                   >> 201     else gtype = fLElastic;
                                                   >> 202   }
                                                   >> 203 
                                                   >> 204   if(gtype == fLElastic) {
                                                   >> 205     t = GeV*GeV*SampleT(ptot,m1,m2,aTarget);
                                                   >> 206     if(t > tmax) gtype = fSWave;
                                                   >> 207   }
                                                   >> 208 
                                                   >> 209   if(gtype == fHElastic) {
                                                   >> 210     t = hElastic->SampleT(theParticle,plab,Z,A);
                                                   >> 211     if(t > tmax) gtype = fSWave;
                                                   >> 212   }
                                                   >> 213 
                                                   >> 214   // NaN finder
                                                   >> 215   if(!(t < 0.0 || t >= 0.0)) {
                                                   >> 216     if (verboseLevel > 0) {
                                                   >> 217       G4cout << "G4HadronElastic:WARNING: Z= " << Z << " N= " 
                                                   >> 218        << N << " pdg= " <<  projPDG
                                                   >> 219        << " mom(GeV)= " << plab/GeV 
                                                   >> 220        << " the model type " << gtype;
                                                   >> 221       if(gtype ==  fQElastic) G4cout << " CHIPS ";
                                                   >> 222       else if(gtype ==  fLElastic) G4cout << " LElastic ";
                                                   >> 223       else if(gtype ==  fHElastic) G4cout << " HElastic ";
                                                   >> 224       G4cout << " S-wave will be sampled" 
                                                   >> 225        << G4endl; 
138     }                                             226     }
139 #endif                                         << 227     gtype = fSWave;
140     t = G4HadronElastic::SampleInvariantT(theP << 
141   }                                               228   }
142                                                   229 
143   G4double phi  = G4UniformRand()*CLHEP::twopi << 230   if(gtype == fSWave) t = G4UniformRand()*tmax;
144   G4double cost = 1. - 2.0*t/pLocalTmax;       << 
145                                                   231 
146   if (cost > 1.0) { cost = 1.0; }              << 232   if(verboseLevel>1)
147   else if(cost < -1.0) { cost = -1.0; }        << 233     G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax 
                                                   >> 234      << " ptot= " << ptot << G4endl;
148                                                   235 
                                                   >> 236   // Sampling in CM system
                                                   >> 237   G4double phi  = G4UniformRand()*twopi;
                                                   >> 238   G4double cost = 1. - 2.0*t/tmax;
                                                   >> 239   if(std::abs(cost) > 1.0) cost = -1.0 + 2.0*G4UniformRand();
149   G4double sint = std::sqrt((1.0-cost)*(1.0+co    240   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
                                                   >> 241   
                                                   >> 242   if (verboseLevel>1) 
                                                   >> 243     G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
150                                                   244 
151   if (verboseLevel>1) {                        << 245   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
152     G4cout << " t= " << t << " tmax(GeV^2)= "  << 246   v1 *= ptot;
153      << " Pcms(GeV)= " << momentumCMS/GeV << " << 247   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
154      << " sin(t)=" << sint << G4endl;          << 
155   }                                            << 
156   G4LorentzVector nlv1(momentumCMS*sint*std::c << 
157            momentumCMS*sint*std::sin(phi),     << 
158                        momentumCMS*cost,       << 
159            std::sqrt(momentumCMS*momentumCMS + << 
160                                                   248 
161   nlv1.boost(bst);                                249   nlv1.boost(bst); 
162                                                   250 
163   G4double eFinal = nlv1.e() - m1;                251   G4double eFinal = nlv1.e() - m1;
164   if (verboseLevel > 1) {                      << 252   if (verboseLevel > 1) 
165     G4cout <<"G4HadronElastic: m= " << m1 << " << 253     G4cout << "Scattered: "
166      << " 4-M Final: " << nlv1                 << 254      << nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal 
                                                   >> 255      << " Proj: 4-mom " << lv1 
                                                   >> 256      <<G4endl;
                                                   >> 257   if(eFinal < 0.0) {
                                                   >> 258     G4cout << "G4HadronElastic WARNING ekin= " << eFinal
                                                   >> 259      << " after scattering of " 
                                                   >> 260      << aParticle->GetDefinition()->GetParticleName()
                                                   >> 261      << " p(GeV/c)= " << plab
                                                   >> 262      << " on " << theDef->GetParticleName()
167      << G4endl;                                   263      << G4endl;
                                                   >> 264     eFinal = 0.0;
                                                   >> 265     nlv1.setE(m1);
168   }                                               266   }
169                                                   267 
170   if(eFinal <= 0.0) {                          << 268   theParticleChange.SetMomentumChange(nlv1.vect().unit());
171     theParticleChange.SetMomentumChange(0.0,0. << 269   theParticleChange.SetEnergyChange(eFinal);
172     theParticleChange.SetEnergyChange(0.0);    << 270   
173   } else {                                     << 271   G4LorentzVector nlv0 = lv - nlv1;
174     theParticleChange.SetMomentumChange(nlv1.v << 272   G4double erec =  nlv0.e() - m2;
175     theParticleChange.SetEnergyChange(eFinal); << 273   if (verboseLevel > 1) 
176   }                                            << 274     G4cout << "Recoil: "
177   lv -= nlv1;                                  << 275      << nlv0<<" m= " << m2 << " ekin(MeV)= " << erec 
178   G4double erec =  std::max(lv.e() - mass2, 0. << 276      <<G4endl;
179   if (verboseLevel > 1) {                      << 277 
180     G4cout << "Recoil: " <<" m= " << mass2 <<  << 278   if(erec > ekinIon) {
181      << " 4-mom: " << lv                       << 279     G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
182      << G4endl;                                << 280     theParticleChange.AddSecondary(aSec);
183   }                                            << 
184                                                << 
185   // the recoil is created if kinetic energy a << 
186   if(erec > GetRecoilEnergyThreshold()) {      << 
187     G4ParticleDefinition * theDef = nullptr;   << 
188     if(Z == 1 && A == 1)       { theDef = theP << 
189     else if (Z == 1 && A == 2) { theDef = theD << 
190     else if (Z == 1 && A == 3) { theDef = G4Tr << 
191     else if (Z == 2 && A == 3) { theDef = G4He << 
192     else if (Z == 2 && A == 4) { theDef = theA << 
193     else {                                     << 
194       theDef =                                 << 
195   G4ParticleTable::GetParticleTable()->GetIonT << 
196     }                                          << 
197     G4DynamicParticle * aSec = new G4DynamicPa << 
198     theParticleChange.AddSecondary(aSec, secID << 
199   } else {                                        281   } else {
                                                   >> 282     if(erec < 0.0) erec = 0.0;
200     theParticleChange.SetLocalEnergyDeposit(er    283     theParticleChange.SetLocalEnergyDeposit(erec);
201   }                                               284   }
202                                                   285 
203   return &theParticleChange;                      286   return &theParticleChange;
204 }                                                 287 }
205                                                   288 
206 // sample momentum transfer in the CMS system  << 
207 G4double                                          289 G4double 
208 G4HadronElastic::SampleInvariantT(const G4Part << 290 G4HadronElastic::SampleT(G4double, G4double, G4double, G4double atno2)
209           G4double mom, G4int, G4int A)        << 
210 {                                                 291 {
211   const G4double plabLowLimit = 400.0*CLHEP::M << 292   // G4cout << "Entering elastic scattering 2"<<G4endl;
212   const G4double GeV2 = GeV*GeV;               << 293   // Compute the direction of elastic scattering.
213   const G4double z07in13 = std::pow(0.7, 0.333 << 294   // It is planned to replace this code with a method based on
214   const G4double numLimit = 18.;               << 295   // parameterized functions and a Monte Carlo method to invert the CDF.
215                                                << 296 
216   G4int pdg = std::abs(part->GetPDGEncoding()) << 297   G4double ran = G4UniformRand();
217   G4double tmax = pLocalTmax/GeV2;             << 298   G4double aa, bb, cc, dd, rr;
218                                                << 299   if (atno2 <= 62.) {
219   G4double aa, bb, cc, dd;                     << 300     aa = std::pow(atno2, 1.63);
220   G4Pow* g4pow = G4Pow::GetInstance();         << 301     bb = 14.5*std::pow(atno2, 0.66);
221   if (A <= 62) {                               << 302     cc = 1.4*std::pow(atno2, 0.33);
222     if (pdg == 211){ //Pions                   << 303     dd = 10.;
223       if(mom >= plabLowLimit){     //High ener << 304   } else {
224   bb = 14.5*g4pow->Z23(A);/*14.5*/             << 305     aa = std::pow(atno2, 1.33);
225   dd = 10.;                                    << 306     bb = 60.*std::pow(atno2, 0.33);
226   cc = 0.075*g4pow->Z13(A)/dd;//1.4            << 307     cc = 0.4*std::pow(atno2, 0.40);
227   //aa = g4pow->powZ(A, 1.93)/bb;//1.63        << 308     dd = 10.;
228   aa = (A*A)/bb;//1.63                         << 309   }
229       } else {                       //Low ene << 310   aa = aa/bb;
230   bb = 29.*z07in13*z07in13*g4pow->Z23(A);      << 311   cc = cc/dd;
231   dd = 15.;                                    << 312   rr = (aa + cc)*ran;
232   cc = 0.04*g4pow->Z13(A)/dd;//1.4             << 313   if (verboseLevel > 1) {
233   aa = g4pow->powZ(A, 1.63)/bb;//1.63          << 314     G4cout << "DoIt: aa,bb,cc,dd,rr" << G4endl;
234       }                                        << 315     G4cout << aa << " " << bb << " " << cc << " " << dd << " " << rr << G4endl;
235     } else { //Other particles                 << 316   }
236       bb = 14.5*g4pow->Z23(A);                 << 317   G4double t1 = -std::log(ran)/bb;
237       dd = 20.;                                << 318   G4double t2 = -std::log(ran)/dd;
238       aa = (A*A)/bb;//1.63                     << 319   if (verboseLevel > 1) {
239       cc = 1.4*g4pow->Z13(A)/dd;               << 320     G4cout << "t1,Fctcos " << t1 << " " << Fctcos(t1, aa, bb, cc, dd, rr) << G4endl;
240     }                                          << 321     G4cout << "t2,Fctcos " << t2 << " " << Fctcos(t2, aa, bb, cc, dd, rr) << G4endl;
241       //===========================            << 322   }
242   } else { //(A>62)                            << 323   G4double eps = 0.001;
243     if (pdg == 211) {                          << 324   G4int ind1 = 10;
244       if(mom >= plabLowLimit){ //high          << 325   G4double t = 0.0;
245   bb = 60.*z07in13*g4pow->Z13(A);//60          << 326   G4int ier1;
246   dd = 30.;                                    << 327   ier1 = Rtmi(&t, t1, t2, eps, ind1,
247   aa = 0.5*(A*A)/bb;//1.33                     << 328         aa, bb, cc, dd, rr);
248   cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4     -- << 329   if (verboseLevel > 1) {
249       } else { //low                           << 330     G4cout << "From Rtmi, ier1=" << ier1 << G4endl;
250   bb = 120.*z07in13*g4pow->Z13(A);//60         << 331     G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << G4endl;
251   dd = 30.;                                    << 
252   aa = 2.*g4pow->powZ(A,1.33)/bb;              << 
253   cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4     -- << 
254       }                                        << 
255     } else {                                   << 
256       bb = 60.*g4pow->Z13(A);                  << 
257       dd = 25.;                                << 
258       aa = g4pow->powZ(A,1.33)/bb;//1.33       << 
259       cc = 0.2*g4pow->powZ(A,0.4)/dd;//1:0.4   << 
260     }                                          << 
261   }                                               332   }
262   G4double q1 = 1.0 - G4Exp(-std::min(bb*tmax, << 333   if (ier1 != 0) t = 0.25*(3.*t1 + t2);
263   G4double q2 = 1.0 - G4Exp(-std::min(dd*tmax, << 334   if (verboseLevel > 1) {
264   G4double s1 = q1*aa;                         << 335       G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << 
265   G4double s2 = q2*cc;                         << 336               G4endl;
266   if((s1 + s2)*G4UniformRand() < s2) {         << 
267     q1 = q2;                                   << 
268     bb = dd;                                   << 
269   }                                               337   }
270   return -GeV2*G4Log(1.0 - G4UniformRand()*q1) << 338   return t;
271 }                                                 339 }
272                                                   340 
273 ////////////////////////////////////////////// << 341 // The following is a "translation" of a root-finding routine
274 //                                             << 342 // from GEANT3.21/GHEISHA.  Some of the labelled block structure has
275 // Cofs for s-,c-,b-particles ds/dt slopes     << 343 // been retained for clarity.  This routine will not be needed after
276                                                << 344 // the planned revisions to DoIt().
277 G4double G4HadronElastic::GetSlopeCof(const G4 << 345 
                                                   >> 346 G4int
                                                   >> 347 G4HadronElastic::Rtmi(G4double* x, G4double xli, G4double xri, G4double eps, 
                                                   >> 348           G4int iend, 
                                                   >> 349           G4double aa, G4double bb, G4double cc, G4double dd, 
                                                   >> 350           G4double rr)
278 {                                                 351 {
279   // The input parameter "pdg" should be the a << 352    G4int ier = 0;
280   // (i.e. the same value for a particle and i << 353    G4double xl = xli;
                                                   >> 354    G4double xr = xri;
                                                   >> 355    *x = xl;
                                                   >> 356    G4double tol = *x;
                                                   >> 357    G4double f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 358    if (f == 0.) return ier;
                                                   >> 359    G4double fl, fr;
                                                   >> 360    fl = f;
                                                   >> 361    *x = xr;
                                                   >> 362    tol = *x;
                                                   >> 363    f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 364    if (f == 0.) return ier;
                                                   >> 365    fr = f;
                                                   >> 366 
                                                   >> 367 // Error return in case of wrong input data
                                                   >> 368    if (fl*fr >= 0.) {
                                                   >> 369       ier = 2;
                                                   >> 370       return ier;
                                                   >> 371    }
                                                   >> 372 
                                                   >> 373 // Basic assumption fl*fr less than 0 is satisfied.
                                                   >> 374 // Generate tolerance for function values.
                                                   >> 375    G4int i = 0;
                                                   >> 376    G4double tolf = 100.*eps;
                                                   >> 377 
                                                   >> 378 // Start iteration loop
                                                   >> 379 label4:
                                                   >> 380    i++;
                                                   >> 381 
                                                   >> 382 // Start bisection loop
                                                   >> 383    for (G4int k = 1; k <= iend; k++) {
                                                   >> 384       *x = 0.5*(xl + xr);
                                                   >> 385       tol = *x;
                                                   >> 386       f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 387       if (f == 0.) return 0;
                                                   >> 388       if (f*fr < 0.) {      // Interchange xl and xr in order to get the
                                                   >> 389          tol = xl;          // same Sign in f and fr
                                                   >> 390          xl = xr;
                                                   >> 391          xr = tol;
                                                   >> 392          tol = fl;
                                                   >> 393          fl = fr;
                                                   >> 394          fr = tol;
                                                   >> 395       }
                                                   >> 396       tol = f - fl;
                                                   >> 397       G4double a = f*tol;
                                                   >> 398       a = a + a;
                                                   >> 399       if (a < fr*(fr - fl) && i <= iend) goto label17;
                                                   >> 400       xr = *x;
                                                   >> 401       fr = f;
                                                   >> 402 
                                                   >> 403 // Test on satisfactory accuracy in bisection loop
                                                   >> 404       tol = eps;
                                                   >> 405       a = std::abs(xr);
                                                   >> 406       if (a > 1.) tol = tol*a;
                                                   >> 407       if (std::abs(xr - xl) <= tol && std::abs(fr - fl) <= tolf) goto label14;
                                                   >> 408    }
                                                   >> 409 // End of bisection loop
                                                   >> 410 
                                                   >> 411 // No convergence after iend iteration steps followed by iend
                                                   >> 412 // successive steps of bisection or steadily increasing function
                                                   >> 413 // values at right bounds.  Error return.
                                                   >> 414    ier = 1;
                                                   >> 415 
                                                   >> 416 label14:
                                                   >> 417    if (std::abs(fr) > std::abs(fl)) {
                                                   >> 418       *x = xl;
                                                   >> 419       f = fl;
                                                   >> 420    }
                                                   >> 421    return ier;
                                                   >> 422 
                                                   >> 423 // Computation of iterated x-value by inverse parabolic interp
                                                   >> 424 label17:
                                                   >> 425    G4double a = fr - f;
                                                   >> 426    G4double dx = (*x - xl)*fl*(1. + f*(a - tol)/(a*(fr - fl)))/tol;
                                                   >> 427    G4double xm = *x;
                                                   >> 428    G4double fm = f;
                                                   >> 429    *x = xl - dx;
                                                   >> 430    tol = *x;
                                                   >> 431    f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 432    if (f == 0.) return ier;
                                                   >> 433 
                                                   >> 434 // Test on satisfactory accuracy in iteration loop
                                                   >> 435    tol = eps;
                                                   >> 436    a = std::abs(*x);
                                                   >> 437    if (a > 1) tol = tol*a;
                                                   >> 438    if (std::abs(dx) <= tol && std::abs(f) <= tolf) return ier;
                                                   >> 439 
                                                   >> 440 // Preparation of next bisection loop
                                                   >> 441    if (f*fl < 0.) {
                                                   >> 442       xr = *x;
                                                   >> 443       fr = f;
                                                   >> 444    }
                                                   >> 445    else {
                                                   >> 446       xl = *x;
                                                   >> 447       fl = f;
                                                   >> 448       xr = xm;
                                                   >> 449       fr = fm;
                                                   >> 450    }
                                                   >> 451    goto label4;
                                                   >> 452 }
281                                                   453 
282   G4double coeff = 1.0;                        << 
283                                                   454 
284   // heavy barions                             << 455 // Test function for root-finder
285                                                   456 
286   static const G4double  lBarCof1S  = 0.88;    << 457 G4double
287   static const G4double  lBarCof2S  = 0.76;    << 458 G4HadronElastic::Fctcos(G4double t, 
288   static const G4double  lBarCof3S  = 0.64;    << 459       G4double aa, G4double bb, G4double cc, G4double dd, 
289   static const G4double  lBarCof1C  = 0.784378 << 460       G4double rr)
290   static const G4double  lBarCofSC  = 0.664378 << 461 {
291   static const G4double  lBarCof2SC = 0.544378 << 462    const G4double expxl = -82.;
292   static const G4double  lBarCof1B  = 0.740659 << 463    const G4double expxu = 82.;
293   static const G4double  lBarCofSB  = 0.620659 << 464 
294   static const G4double  lBarCof2SB = 0.500659 << 465    G4double test1 = -bb*t;
295                                                << 466    if (test1 > expxu) test1 = expxu;
296   if( pdg == 3122 || pdg == 3222 ||  pdg == 31 << 467    if (test1 < expxl) test1 = expxl;
297   {                                            << 468 
298     coeff = lBarCof1S; // Lambda, Sigma+, Sigm << 469    G4double test2 = -dd*t;
299                                                << 470    if (test2 > expxu) test2 = expxu;
300   } else if( pdg == 3322 || pdg == 3312   )    << 471    if (test2 < expxl) test2 = expxl;
301   {                                            << 472 
302     coeff = lBarCof2S; // Xi-, Xi0             << 473    return aa*std::exp(test1) + cc*std::exp(test2) - rr;
303   }                                            << 
304   else if( pdg == 3324)                        << 
305   {                                            << 
306     coeff = lBarCof3S; // Omega                << 
307   }                                            << 
308   else if( pdg == 4122 ||  pdg == 4212 ||   pd << 
309   {                                            << 
310     coeff = lBarCof1C; // LambdaC+, SigmaC+, S << 
311   }                                            << 
312   else if( pdg == 4332 )                       << 
313   {                                            << 
314     coeff = lBarCof2SC; // OmegaC              << 
315   }                                            << 
316   else if( pdg == 4232 || pdg == 4132 )        << 
317   {                                            << 
318     coeff = lBarCofSC; // XiC+, XiC0           << 
319   }                                            << 
320   else if( pdg == 5122 || pdg == 5222 || pdg = << 
321   {                                            << 
322     coeff = lBarCof1B; // LambdaB, SigmaB+, Si << 
323   }                                            << 
324   else if( pdg == 5332 )                       << 
325   {                                            << 
326     coeff = lBarCof2SB; // OmegaB-             << 
327   }                                            << 
328   else if( pdg == 5132 || pdg == 5232 ) // XiB << 
329   {                                            << 
330     coeff = lBarCofSB;                         << 
331   }                                            << 
332   // heavy mesons Kaons?                       << 
333   static const G4double lMesCof1S = 0.82; // K << 
334   static const G4double llMesCof1C = 0.676568; << 
335   static const G4double llMesCof1B = 0.610989; << 
336   static const G4double llMesCof2C = 0.353135; << 
337   static const G4double llMesCof2B = 0.221978; << 
338   static const G4double llMesCofSC = 0.496568; << 
339   static const G4double llMesCofSB = 0.430989; << 
340   static const G4double llMesCofCB = 0.287557; << 
341   static const G4double llMesCofEtaP = 0.88;   << 
342   static const G4double llMesCofEta = 0.76;    << 
343                                                << 
344   if( pdg == 321 || pdg == 311 || pdg == 310 ) << 
345   {                                            << 
346     coeff = lMesCof1S; //K+-0                  << 
347   }                                            << 
348   else if( pdg == 511 ||  pdg == 521  )        << 
349   {                                            << 
350     coeff = llMesCof1B; // BMeson0, BMeson+    << 
351   }                                            << 
352   else if(pdg == 421 ||  pdg == 411 )          << 
353   {                                            << 
354     coeff = llMesCof1C; // DMeson+, DMeson0    << 
355   }                                            << 
356   else if( pdg == 531  )                       << 
357   {                                            << 
358     coeff = llMesCofSB; // BSMeson0            << 
359   }                                            << 
360   else if( pdg == 541 )                        << 
361   {                                            << 
362     coeff = llMesCofCB; // BCMeson+-           << 
363   }                                            << 
364   else if(pdg == 431 )                         << 
365   {                                            << 
366     coeff = llMesCofSC; // DSMeson+-           << 
367   }                                            << 
368   else if(pdg == 441 || pdg == 443 )           << 
369   {                                            << 
370     coeff = llMesCof2C; // Etac, JPsi          << 
371   }                                            << 
372   else if(pdg == 553 )                         << 
373   {                                            << 
374     coeff = llMesCof2B; // Upsilon             << 
375   }                                            << 
376   else if(pdg == 221 )                         << 
377   {                                            << 
378     coeff = llMesCofEta; // Eta                << 
379   }                                            << 
380   else if(pdg == 331 )                         << 
381   {                                            << 
382     coeff = llMesCofEtaP; // Eta'              << 
383   }                                            << 
384   return coeff;                                << 
385 }                                                 474 }
386                                                   475 
387                                                   476 
                                                   >> 477 void
                                                   >> 478 G4HadronElastic::Defs1(G4double p, G4double px, G4double py, G4double pz, 
                                                   >> 479            G4double pxinc, G4double pyinc, G4double pzinc, 
                                                   >> 480            G4double* pxnew, G4double* pynew, G4double* pznew)
                                                   >> 481 {
                                                   >> 482 // Transform scattered particle to reflect direction of incident particle
                                                   >> 483    G4double pt2 = pxinc*pxinc + pyinc*pyinc;
                                                   >> 484    if (pt2 > 0.) {
                                                   >> 485       G4double cost = pzinc/p;
                                                   >> 486       G4double sint1 = std::sqrt(std::abs((1. - cost )*(1.+cost)));
                                                   >> 487       G4double sint2 = std::sqrt(pt2)/p;
                                                   >> 488       G4double sint = 0.5*(sint1 + sint2);
                                                   >> 489       G4double ph = pi*0.5;
                                                   >> 490       if (pyinc < 0.) ph = pi*1.5;
                                                   >> 491       if (std::abs(pxinc) > 1.e-6) ph = std::atan2(pyinc, pxinc);
                                                   >> 492       G4double cosp = std::cos(ph);
                                                   >> 493       G4double sinp = std::sin(ph);
                                                   >> 494       if (verboseLevel > 1) {
                                                   >> 495          G4cout << "cost sint " << cost << " " << sint << G4endl;
                                                   >> 496          G4cout << "cosp sinp " << cosp << " " << sinp << G4endl;
                                                   >> 497       }
                                                   >> 498       *pxnew = cost*cosp*px - sinp*py + sint*cosp*pz;
                                                   >> 499       *pynew = cost*sinp*px + cosp*py + sint*sinp*pz;
                                                   >> 500       *pznew =     -sint*px                 +cost*pz;
                                                   >> 501    }
                                                   >> 502    else {
                                                   >> 503        G4double cost=pzinc/p;
                                                   >> 504        *pxnew = cost*px;
                                                   >> 505        *pynew = py;
                                                   >> 506        *pznew = cost*pz;
                                                   >> 507    }
                                                   >> 508 }
388                                                   509