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.3.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.48 2007/05/05 18:45:23 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-08-03-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 // 16-Nov-06 V.Ivanchenko Simplify logic of choosing of the model for sampling
                                                   >>  58 // 30-Mar-07 V.Ivanchenko lowEnergyLimitQ=0, lowEnergyLimitHE = 1.0*GeV,
                                                   >>  59 //                        lowestEnergyLimit= 0
                                                   >>  60 // 04-May-07 V.Ivanchenko do not use HE model for hydrogen target to avoid NaN;
                                                   >>  61 //                        use QElastic for p, n incident for any energy for 
                                                   >>  62 //                        p and He targets only  
                                                   >>  63 //
 31                                                    64 
 32 #include "G4HadronElastic.hh"                      65 #include "G4HadronElastic.hh"
 33 #include "G4SystemOfUnits.hh"                  << 
 34 #include "G4ParticleTable.hh"                      66 #include "G4ParticleTable.hh"
 35 #include "G4ParticleDefinition.hh"                 67 #include "G4ParticleDefinition.hh"
 36 #include "G4IonTable.hh"                           68 #include "G4IonTable.hh"
                                                   >>  69 #include "G4QElasticCrossSection.hh"
                                                   >>  70 #include "G4VQCrossSection.hh"
                                                   >>  71 #include "G4ElasticHadrNucleusHE.hh"
 37 #include "Randomize.hh"                            72 #include "Randomize.hh"
 38 #include "G4Proton.hh"                             73 #include "G4Proton.hh"
 39 #include "G4Neutron.hh"                            74 #include "G4Neutron.hh"
 40 #include "G4Deuteron.hh"                           75 #include "G4Deuteron.hh"
 41 #include "G4Alpha.hh"                              76 #include "G4Alpha.hh"
 42 #include "G4Pow.hh"                            <<  77 #include "G4PionPlus.hh"
 43 #include "G4Exp.hh"                            <<  78 #include "G4PionMinus.hh"
 44 #include "G4Log.hh"                            <<  79 #include "G4NistManager.hh"
 45 #include "G4HadronicParameters.hh"             << 
 46 #include "G4PhysicsModelCatalog.hh"            << 
 47                                                << 
 48                                                    80 
 49 G4HadronElastic::G4HadronElastic(const G4Strin <<  81 G4HadronElastic::G4HadronElastic() 
 50   : G4HadronicInteraction(name), secID(-1)     <<  82   : G4HadronicInteraction()
 51 {                                                  83 {
 52   SetMinEnergy( 0.0*GeV );                         84   SetMinEnergy( 0.0*GeV );
 53   SetMaxEnergy( G4HadronicParameters::Instance <<  85   SetMaxEnergy( 100.*TeV );
 54   lowestEnergyLimit= 1.e-6*eV;                 <<  86   verboseLevel= 0;
 55   pLocalTmax  = 0.0;                           <<  87   lowEnergyRecoilLimit = 100.*keV;  
 56   nwarn = 0;                                   <<  88   lowEnergyLimitQ  = 0.0*GeV;  
                                                   >>  89   lowEnergyLimitHE = 1.0*GeV;  
                                                   >>  90   lowestEnergyLimit= 0.0*keV;  
                                                   >>  91   plabLowLimit     = 20.0*MeV;
                                                   >>  92 
                                                   >>  93   qCManager   = G4QElasticCrossSection::GetPointer();
                                                   >>  94   nistManager = G4NistManager::Instance();
                                                   >>  95   hElastic    = new G4ElasticHadrNucleusHE();
 57                                                    96 
 58   theProton   = G4Proton::Proton();                97   theProton   = G4Proton::Proton();
 59   theNeutron  = G4Neutron::Neutron();              98   theNeutron  = G4Neutron::Neutron();
 60   theDeuteron = G4Deuteron::Deuteron();            99   theDeuteron = G4Deuteron::Deuteron();
 61   theAlpha    = G4Alpha::Alpha();                 100   theAlpha    = G4Alpha::Alpha();
 62                                                << 101   thePionPlus = G4PionPlus::PionPlus();
 63   secID = G4PhysicsModelCatalog::GetModelID( " << 102   thePionMinus= G4PionMinus::PionMinus();
 64 }                                                 103 }
 65                                                   104 
 66 G4HadronElastic::~G4HadronElastic()               105 G4HadronElastic::~G4HadronElastic()
 67 {}                                             << 106 {
                                                   >> 107   delete hElastic;
                                                   >> 108 }
 68                                                   109 
                                                   >> 110 G4VQCrossSection* G4HadronElastic::GetCS()
                                                   >> 111 {
                                                   >> 112   return qCManager;
                                                   >> 113 }
 69                                                   114 
 70 void G4HadronElastic::ModelDescription(std::os << 115 G4ElasticHadrNucleusHE* G4HadronElastic::GetHElastic()
 71 {                                                 116 {
 72   outFile << "G4HadronElastic is the base clas << 117   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 }                                                 118 }
 80                                                   119 
 81 G4HadFinalState* G4HadronElastic::ApplyYoursel    120 G4HadFinalState* G4HadronElastic::ApplyYourself(
 82      const G4HadProjectile& aTrack, G4Nucleus&    121      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
 83 {                                                 122 {
 84   theParticleChange.Clear();                      123   theParticleChange.Clear();
 85                                                   124 
 86   const G4HadProjectile* aParticle = &aTrack;     125   const G4HadProjectile* aParticle = &aTrack;
 87   G4double ekin = aParticle->GetKineticEnergy(    126   G4double ekin = aParticle->GetKineticEnergy();
 88                                                << 
 89   // no scattering below the limit             << 
 90   if(ekin <= lowestEnergyLimit) {                 127   if(ekin <= lowestEnergyLimit) {
 91     theParticleChange.SetEnergyChange(ekin);      128     theParticleChange.SetEnergyChange(ekin);
 92     theParticleChange.SetMomentumChange(0.,0., << 129     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
 93     return &theParticleChange;                    130     return &theParticleChange;
 94   }                                               131   }
 95                                                   132 
 96   G4int A = targetNucleus.GetA_asInt();        << 133   G4double aTarget = targetNucleus.GetN();
 97   G4int Z = targetNucleus.GetZ_asInt();        << 134   G4double zTarget = targetNucleus.GetZ();
                                                   >> 135 
                                                   >> 136   G4double plab = aParticle->GetTotalMomentum();
                                                   >> 137   if (verboseLevel >1) 
                                                   >> 138     G4cout << "G4HadronElastic::DoIt: Incident particle plab=" 
                                                   >> 139      << plab/GeV << " GeV/c " 
                                                   >> 140      << " ekin(MeV) = " << ekin/MeV << "  " 
                                                   >> 141      << aParticle->GetDefinition()->GetParticleName() << G4endl;
 98                                                   142 
 99   // Scattered particle referred to axis of in    143   // Scattered particle referred to axis of incident particle
100   const G4ParticleDefinition* theParticle = aP    144   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
101   G4double m1 = theParticle->GetPDGMass();        145   G4double m1 = theParticle->GetPDGMass();
102   G4double plab = std::sqrt(ekin*(ekin + 2.0*m << 
103                                                   146 
104   if (verboseLevel>1) {                        << 147   G4int Z = static_cast<G4int>(zTarget+0.5);
105     G4cout << "G4HadronElastic: "              << 148   G4int A = static_cast<G4int>(aTarget+0.5);
106      << aParticle->GetDefinition()->GetParticl << 149   G4int N = A - Z;
107      << " Plab(GeV/c)= " << plab/GeV           << 150   G4int projPDG = theParticle->GetPDGEncoding();
108      << " Ekin(MeV) = " << ekin/MeV            << 151   if (verboseLevel>1) 
109      << " scattered off Z= " << Z              << 152     G4cout << "G4HadronElastic for " << theParticle->GetParticleName()
110      << " A= " << A                            << 153      << " PDGcode= " << projPDG << " on nucleus Z= " << Z 
                                                   >> 154      << " A= " << A << " N= " << N 
111      << G4endl;                                   155      << G4endl;
112   }                                            << 
113                                                   156 
114   G4double mass2 = G4NucleiProperties::GetNucl << 157   G4ParticleDefinition * theDef = 0;
115   G4double e1 = m1 + ekin;                     << 
116   G4LorentzVector lv(0.0,0.0,plab,e1+mass2);   << 
117   G4ThreeVector bst = lv.boostVector();        << 
118   G4double momentumCMS = plab*mass2/std::sqrt( << 
119                                                   158 
120   pLocalTmax = 4.0*momentumCMS*momentumCMS;    << 159   if(Z == 1 && A == 1)       theDef = theProton;
                                                   >> 160   else if (Z == 1 && A == 2) theDef = theDeuteron;
                                                   >> 161   else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
                                                   >> 162   else if (Z == 2 && A == 3) theDef = G4He3::He3();
                                                   >> 163   else if (Z == 2 && A == 4) theDef = theAlpha;
                                                   >> 164   else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
                                                   >> 165  
                                                   >> 166   G4double m2 = theDef->GetPDGMass();
                                                   >> 167   G4LorentzVector lv1 = aParticle->Get4Momentum();
                                                   >> 168   G4LorentzVector lv(0.0,0.0,0.0,m2);   
                                                   >> 169   lv += lv1;
121                                                   170 
122   // Sampling in CM system                     << 171   G4ThreeVector bst = lv.boostVector();
123   G4double t = SampleInvariantT(theParticle, p << 172   lv1.boost(-bst);
124                                                   173 
125   if(t < 0.0 || t > pLocalTmax) {              << 174   G4ThreeVector p1 = lv1.vect();
126     // For the very rare cases where cos(theta << 175   G4double ptot = p1.mag();
127     // print some debugging information via a  << 176   G4double tmax = 4.0*ptot*ptot;
128     // using the default algorithm             << 177   G4double t = 0.0;
129 #ifdef G4VERBOSE                               << 178 
130     if(nwarn < 2) {                            << 179   // Choose generator
131       G4ExceptionDescription ed;               << 180   G4ElasticGenerator gtype = fLElastic;
132       ed << GetModelName() << " wrong sampling << 181 
133    << " for " << aParticle->GetDefinition()->G << 182   // Q-elastic for p,n scattering on H and He
134    << " ekin=" << ekin << " MeV"               << 183   if ((theParticle == theProton || theParticle == theNeutron)
135    << " off (Z,A)=(" << Z << "," << A << ") -  << 184        && Z <= 2 && ekin >= lowEnergyLimitQ)  
136       G4Exception( "G4HadronElastic::ApplyYour << 185     gtype = fQElastic;
137       ++nwarn;                                 << 186 
                                                   >> 187   else {
                                                   >> 188     // S-wave for very low energy
                                                   >> 189     if(plab < plabLowLimit) gtype = fSWave;
                                                   >> 190     // HE-elastic for energetic projectiles
                                                   >> 191     else if(ekin >= lowEnergyLimitHE && A < 238 && Z>= 2) gtype = fHElastic;
                                                   >> 192   }
                                                   >> 193 
                                                   >> 194   //
                                                   >> 195   // Sample t
                                                   >> 196   //
                                                   >> 197   if(gtype == fQElastic) {
                                                   >> 198     if (verboseLevel >1) 
                                                   >> 199       G4cout << "G4HadronElastic: Z= " << Z << " N= " 
                                                   >> 200        << N << " pdg= " <<  projPDG
                                                   >> 201        << " mom(GeV)= " << plab/GeV << "  " << qCManager << G4endl; 
                                                   >> 202     if(Z == 1 && N == 2) N = 1;
                                                   >> 203     else if(Z == 2 && N == 1) N = 2;
                                                   >> 204     G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
                                                   >> 205 
                                                   >> 206     // check if cross section is reasonable
                                                   >> 207     if(cs > 0.0) t = qCManager->GetExchangeT(Z,N,projPDG);
                                                   >> 208     else if(plab > plabLowLimit) gtype = fLElastic;
                                                   >> 209     else gtype = fSWave;
                                                   >> 210   }
                                                   >> 211 
                                                   >> 212   if(gtype == fLElastic) {
                                                   >> 213     t = GeV*GeV*SampleT(ptot,m1,m2,aTarget);
                                                   >> 214     if(t > tmax) gtype = fSWave;
                                                   >> 215   }
                                                   >> 216 
                                                   >> 217   // use mean atomic number
                                                   >> 218   if(gtype == fHElastic) {
                                                   >> 219     G4int A0 = static_cast<G4int>(nistManager->GetAtomicMassAmu(Z)+0.5);
                                                   >> 220     t = hElastic->SampleT(theParticle,plab,Z,A0);
                                                   >> 221     if(t > tmax) gtype = fSWave;
                                                   >> 222   }
                                                   >> 223 
                                                   >> 224   // NaN finder
                                                   >> 225   if(!(t < 0.0 || t >= 0.0)) {
                                                   >> 226     if (verboseLevel > 0) {
                                                   >> 227       G4cout << "G4HadronElastic:WARNING: Z= " << Z << " N= " 
                                                   >> 228        << N << " pdg= " <<  projPDG
                                                   >> 229        << " mom(GeV)= " << plab/GeV 
                                                   >> 230        << " the model type " << gtype;
                                                   >> 231       if(gtype ==  fQElastic) G4cout << " CHIPS ";
                                                   >> 232       else if(gtype ==  fLElastic) G4cout << " LElastic ";
                                                   >> 233       else if(gtype ==  fHElastic) G4cout << " HElastic ";
                                                   >> 234       G4cout << " S-wave will be sampled" 
                                                   >> 235        << G4endl; 
138     }                                             236     }
139 #endif                                         << 237     gtype = fSWave;
140     t = G4HadronElastic::SampleInvariantT(theP << 
141   }                                               238   }
142                                                   239 
143   G4double phi  = G4UniformRand()*CLHEP::twopi << 240   if(gtype == fSWave) t = G4UniformRand()*tmax;
144   G4double cost = 1. - 2.0*t/pLocalTmax;       << 
145                                                   241 
146   if (cost > 1.0) { cost = 1.0; }              << 242   if(verboseLevel>1)
147   else if(cost < -1.0) { cost = -1.0; }        << 243     G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax 
                                                   >> 244      << " ptot= " << ptot << G4endl;
148                                                   245 
                                                   >> 246   // Sampling in CM system
                                                   >> 247   G4double phi  = G4UniformRand()*twopi;
                                                   >> 248   G4double cost = 1. - 2.0*t/tmax;
                                                   >> 249   if(std::abs(cost) > 1.0) cost = -1.0 + 2.0*G4UniformRand();
149   G4double sint = std::sqrt((1.0-cost)*(1.0+co    250   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
                                                   >> 251   
                                                   >> 252   if (verboseLevel>1) 
                                                   >> 253     G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
150                                                   254 
151   if (verboseLevel>1) {                        << 255   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
152     G4cout << " t= " << t << " tmax(GeV^2)= "  << 256   v1 *= ptot;
153      << " Pcms(GeV)= " << momentumCMS/GeV << " << 257   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                                                   258 
161   nlv1.boost(bst);                                259   nlv1.boost(bst); 
162                                                   260 
163   G4double eFinal = nlv1.e() - m1;                261   G4double eFinal = nlv1.e() - m1;
164   if (verboseLevel > 1) {                      << 262   if (verboseLevel > 1) 
165     G4cout <<"G4HadronElastic: m= " << m1 << " << 263     G4cout << "Scattered: "
166      << " 4-M Final: " << nlv1                 << 264      << nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal 
                                                   >> 265      << " Proj: 4-mom " << lv1 
                                                   >> 266      <<G4endl;
                                                   >> 267   if(eFinal < 0.0) {
                                                   >> 268     G4cout << "G4HadronElastic WARNING ekin= " << eFinal
                                                   >> 269      << " after scattering of " 
                                                   >> 270      << aParticle->GetDefinition()->GetParticleName()
                                                   >> 271      << " p(GeV/c)= " << plab
                                                   >> 272      << " on " << theDef->GetParticleName()
167      << G4endl;                                   273      << G4endl;
                                                   >> 274     eFinal = 0.0;
                                                   >> 275     nlv1.setE(m1);
168   }                                               276   }
169                                                   277 
170   if(eFinal <= 0.0) {                          << 278   theParticleChange.SetMomentumChange(nlv1.vect().unit());
171     theParticleChange.SetMomentumChange(0.0,0. << 279   theParticleChange.SetEnergyChange(eFinal);
172     theParticleChange.SetEnergyChange(0.0);    << 280   
173   } else {                                     << 281   G4LorentzVector nlv0 = lv - nlv1;
174     theParticleChange.SetMomentumChange(nlv1.v << 282   G4double erec =  nlv0.e() - m2;
175     theParticleChange.SetEnergyChange(eFinal); << 283   if (verboseLevel > 1) 
176   }                                            << 284     G4cout << "Recoil: "
177   lv -= nlv1;                                  << 285      << nlv0<<" m= " << m2 << " ekin(MeV)= " << erec 
178   G4double erec =  std::max(lv.e() - mass2, 0. << 286      <<G4endl;
179   if (verboseLevel > 1) {                      << 287 
180     G4cout << "Recoil: " <<" m= " << mass2 <<  << 288   if(erec > lowEnergyRecoilLimit) {
181      << " 4-mom: " << lv                       << 289     G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
182      << G4endl;                                << 290     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 {                                        291   } else {
                                                   >> 292     if(erec < 0.0) erec = 0.0;
200     theParticleChange.SetLocalEnergyDeposit(er    293     theParticleChange.SetLocalEnergyDeposit(erec);
201   }                                               294   }
202                                                   295 
203   return &theParticleChange;                      296   return &theParticleChange;
204 }                                                 297 }
205                                                   298 
206 // sample momentum transfer in the CMS system  << 
207 G4double                                          299 G4double 
208 G4HadronElastic::SampleInvariantT(const G4Part << 300 G4HadronElastic::SampleT(G4double, G4double, G4double, G4double atno2)
209           G4double mom, G4int, G4int A)        << 
210 {                                                 301 {
211   const G4double plabLowLimit = 400.0*CLHEP::M << 302   // G4cout << "Entering elastic scattering 2"<<G4endl;
212   const G4double GeV2 = GeV*GeV;               << 303   // Compute the direction of elastic scattering.
213   const G4double z07in13 = std::pow(0.7, 0.333 << 304   // It is planned to replace this code with a method based on
214   const G4double numLimit = 18.;               << 305   // parameterized functions and a Monte Carlo method to invert the CDF.
215                                                << 306 
216   G4int pdg = std::abs(part->GetPDGEncoding()) << 307   G4double ran = G4UniformRand();
217   G4double tmax = pLocalTmax/GeV2;             << 308   G4double aa, bb, cc, dd, rr;
218                                                << 309   if (atno2 <= 62.) {
219   G4double aa, bb, cc, dd;                     << 310     aa = std::pow(atno2, 1.63);
220   G4Pow* g4pow = G4Pow::GetInstance();         << 311     bb = 14.5*std::pow(atno2, 0.66);
221   if (A <= 62) {                               << 312     cc = 1.4*std::pow(atno2, 0.33);
222     if (pdg == 211){ //Pions                   << 313     dd = 10.;
223       if(mom >= plabLowLimit){     //High ener << 314   } else {
224   bb = 14.5*g4pow->Z23(A);/*14.5*/             << 315     aa = std::pow(atno2, 1.33);
225   dd = 10.;                                    << 316     bb = 60.*std::pow(atno2, 0.33);
226   cc = 0.075*g4pow->Z13(A)/dd;//1.4            << 317     cc = 0.4*std::pow(atno2, 0.40);
227   //aa = g4pow->powZ(A, 1.93)/bb;//1.63        << 318     dd = 10.;
228   aa = (A*A)/bb;//1.63                         << 319   }
229       } else {                       //Low ene << 320   aa = aa/bb;
230   bb = 29.*z07in13*z07in13*g4pow->Z23(A);      << 321   cc = cc/dd;
231   dd = 15.;                                    << 322   rr = (aa + cc)*ran;
232   cc = 0.04*g4pow->Z13(A)/dd;//1.4             << 323   if (verboseLevel > 1) {
233   aa = g4pow->powZ(A, 1.63)/bb;//1.63          << 324     G4cout << "DoIt: aa,bb,cc,dd,rr" << G4endl;
234       }                                        << 325     G4cout << aa << " " << bb << " " << cc << " " << dd << " " << rr << G4endl;
235     } else { //Other particles                 << 326   }
236       bb = 14.5*g4pow->Z23(A);                 << 327   G4double t1 = -std::log(ran)/bb;
237       dd = 20.;                                << 328   G4double t2 = -std::log(ran)/dd;
238       aa = (A*A)/bb;//1.63                     << 329   if (verboseLevel > 1) {
239       cc = 1.4*g4pow->Z13(A)/dd;               << 330     G4cout << "t1,Fctcos " << t1 << " " << Fctcos(t1, aa, bb, cc, dd, rr) << G4endl;
240     }                                          << 331     G4cout << "t2,Fctcos " << t2 << " " << Fctcos(t2, aa, bb, cc, dd, rr) << G4endl;
241       //===========================            << 
242   } else { //(A>62)                            << 
243     if (pdg == 211) {                          << 
244       if(mom >= plabLowLimit){ //high          << 
245   bb = 60.*z07in13*g4pow->Z13(A);//60          << 
246   dd = 30.;                                    << 
247   aa = 0.5*(A*A)/bb;//1.33                     << 
248   cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4     -- << 
249       } else { //low                           << 
250   bb = 120.*z07in13*g4pow->Z13(A);//60         << 
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   G4double eps = 0.001;
263   G4double q2 = 1.0 - G4Exp(-std::min(dd*tmax, << 334   G4int ind1 = 10;
264   G4double s1 = q1*aa;                         << 335   G4double t = 0.0;
265   G4double s2 = q2*cc;                         << 336   G4int ier1;
266   if((s1 + s2)*G4UniformRand() < s2) {         << 337   ier1 = Rtmi(&t, t1, t2, eps, ind1,
267     q1 = q2;                                   << 338         aa, bb, cc, dd, rr);
268     bb = dd;                                   << 339   if (verboseLevel > 1) {
                                                   >> 340     G4cout << "From Rtmi, ier1=" << ier1 << G4endl;
                                                   >> 341     G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << G4endl;
269   }                                               342   }
270   return -GeV2*G4Log(1.0 - G4UniformRand()*q1) << 343   if (ier1 != 0) t = 0.25*(3.*t1 + t2);
                                                   >> 344   if (verboseLevel > 1) {
                                                   >> 345       G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << 
                                                   >> 346               G4endl;
                                                   >> 347   }
                                                   >> 348   return t;
271 }                                                 349 }
272                                                   350 
273 ////////////////////////////////////////////// << 351 // The following is a "translation" of a root-finding routine
274 //                                             << 352 // from GEANT3.21/GHEISHA.  Some of the labelled block structure has
275 // Cofs for s-,c-,b-particles ds/dt slopes     << 353 // been retained for clarity.  This routine will not be needed after
276                                                << 354 // the planned revisions to DoIt().
277 G4double G4HadronElastic::GetSlopeCof(const G4 << 355 
                                                   >> 356 G4int
                                                   >> 357 G4HadronElastic::Rtmi(G4double* x, G4double xli, G4double xri, G4double eps, 
                                                   >> 358           G4int iend, 
                                                   >> 359           G4double aa, G4double bb, G4double cc, G4double dd, 
                                                   >> 360           G4double rr)
278 {                                                 361 {
279   // The input parameter "pdg" should be the a << 362    G4int ier = 0;
280   // (i.e. the same value for a particle and i << 363    G4double xl = xli;
                                                   >> 364    G4double xr = xri;
                                                   >> 365    *x = xl;
                                                   >> 366    G4double tol = *x;
                                                   >> 367    G4double f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 368    if (f == 0.) return ier;
                                                   >> 369    G4double fl, fr;
                                                   >> 370    fl = f;
                                                   >> 371    *x = xr;
                                                   >> 372    tol = *x;
                                                   >> 373    f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 374    if (f == 0.) return ier;
                                                   >> 375    fr = f;
                                                   >> 376 
                                                   >> 377 // Error return in case of wrong input data
                                                   >> 378    if (fl*fr >= 0.) {
                                                   >> 379       ier = 2;
                                                   >> 380       return ier;
                                                   >> 381    }
                                                   >> 382 
                                                   >> 383 // Basic assumption fl*fr less than 0 is satisfied.
                                                   >> 384 // Generate tolerance for function values.
                                                   >> 385    G4int i = 0;
                                                   >> 386    G4double tolf = 100.*eps;
                                                   >> 387 
                                                   >> 388 // Start iteration loop
                                                   >> 389 label4:
                                                   >> 390    i++;
                                                   >> 391 
                                                   >> 392 // Start bisection loop
                                                   >> 393    for (G4int k = 1; k <= iend; k++) {
                                                   >> 394       *x = 0.5*(xl + xr);
                                                   >> 395       tol = *x;
                                                   >> 396       f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 397       if (f == 0.) return 0;
                                                   >> 398       if (f*fr < 0.) {      // Interchange xl and xr in order to get the
                                                   >> 399          tol = xl;          // same Sign in f and fr
                                                   >> 400          xl = xr;
                                                   >> 401          xr = tol;
                                                   >> 402          tol = fl;
                                                   >> 403          fl = fr;
                                                   >> 404          fr = tol;
                                                   >> 405       }
                                                   >> 406       tol = f - fl;
                                                   >> 407       G4double a = f*tol;
                                                   >> 408       a = a + a;
                                                   >> 409       if (a < fr*(fr - fl) && i <= iend) goto label17;
                                                   >> 410       xr = *x;
                                                   >> 411       fr = f;
                                                   >> 412 
                                                   >> 413 // Test on satisfactory accuracy in bisection loop
                                                   >> 414       tol = eps;
                                                   >> 415       a = std::abs(xr);
                                                   >> 416       if (a > 1.) tol = tol*a;
                                                   >> 417       if (std::abs(xr - xl) <= tol && std::abs(fr - fl) <= tolf) goto label14;
                                                   >> 418    }
                                                   >> 419 // End of bisection loop
                                                   >> 420 
                                                   >> 421 // No convergence after iend iteration steps followed by iend
                                                   >> 422 // successive steps of bisection or steadily increasing function
                                                   >> 423 // values at right bounds.  Error return.
                                                   >> 424    ier = 1;
                                                   >> 425 
                                                   >> 426 label14:
                                                   >> 427    if (std::abs(fr) > std::abs(fl)) {
                                                   >> 428       *x = xl;
                                                   >> 429       f = fl;
                                                   >> 430    }
                                                   >> 431    return ier;
                                                   >> 432 
                                                   >> 433 // Computation of iterated x-value by inverse parabolic interp
                                                   >> 434 label17:
                                                   >> 435    G4double a = fr - f;
                                                   >> 436    G4double dx = (*x - xl)*fl*(1. + f*(a - tol)/(a*(fr - fl)))/tol;
                                                   >> 437    G4double xm = *x;
                                                   >> 438    G4double fm = f;
                                                   >> 439    *x = xl - dx;
                                                   >> 440    tol = *x;
                                                   >> 441    f = Fctcos(tol, aa, bb, cc, dd, rr);
                                                   >> 442    if (f == 0.) return ier;
                                                   >> 443 
                                                   >> 444 // Test on satisfactory accuracy in iteration loop
                                                   >> 445    tol = eps;
                                                   >> 446    a = std::abs(*x);
                                                   >> 447    if (a > 1) tol = tol*a;
                                                   >> 448    if (std::abs(dx) <= tol && std::abs(f) <= tolf) return ier;
                                                   >> 449 
                                                   >> 450 // Preparation of next bisection loop
                                                   >> 451    if (f*fl < 0.) {
                                                   >> 452       xr = *x;
                                                   >> 453       fr = f;
                                                   >> 454    }
                                                   >> 455    else {
                                                   >> 456       xl = *x;
                                                   >> 457       fl = f;
                                                   >> 458       xr = xm;
                                                   >> 459       fr = fm;
                                                   >> 460    }
                                                   >> 461    goto label4;
                                                   >> 462 }
281                                                   463 
282   G4double coeff = 1.0;                        << 
283                                                   464 
284   // heavy barions                             << 465 // Test function for root-finder
285                                                   466 
286   static const G4double  lBarCof1S  = 0.88;    << 467 G4double
287   static const G4double  lBarCof2S  = 0.76;    << 468 G4HadronElastic::Fctcos(G4double t, 
288   static const G4double  lBarCof3S  = 0.64;    << 469       G4double aa, G4double bb, G4double cc, G4double dd, 
289   static const G4double  lBarCof1C  = 0.784378 << 470       G4double rr)
290   static const G4double  lBarCofSC  = 0.664378 << 471 {
291   static const G4double  lBarCof2SC = 0.544378 << 472    const G4double expxl = -82.;
292   static const G4double  lBarCof1B  = 0.740659 << 473    const G4double expxu = 82.;
293   static const G4double  lBarCofSB  = 0.620659 << 474 
294   static const G4double  lBarCof2SB = 0.500659 << 475    G4double test1 = -bb*t;
295                                                << 476    if (test1 > expxu) test1 = expxu;
296   if( pdg == 3122 || pdg == 3222 ||  pdg == 31 << 477    if (test1 < expxl) test1 = expxl;
297   {                                            << 478 
298     coeff = lBarCof1S; // Lambda, Sigma+, Sigm << 479    G4double test2 = -dd*t;
299                                                << 480    if (test2 > expxu) test2 = expxu;
300   } else if( pdg == 3322 || pdg == 3312   )    << 481    if (test2 < expxl) test2 = expxl;
301   {                                            << 482 
302     coeff = lBarCof2S; // Xi-, Xi0             << 483    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 }                                                 484 }
386                                                   485 
387                                                   486 
                                                   >> 487 void
                                                   >> 488 G4HadronElastic::Defs1(G4double p, G4double px, G4double py, G4double pz, 
                                                   >> 489            G4double pxinc, G4double pyinc, G4double pzinc, 
                                                   >> 490            G4double* pxnew, G4double* pynew, G4double* pznew)
                                                   >> 491 {
                                                   >> 492 // Transform scattered particle to reflect direction of incident particle
                                                   >> 493    G4double pt2 = pxinc*pxinc + pyinc*pyinc;
                                                   >> 494    if (pt2 > 0.) {
                                                   >> 495       G4double cost = pzinc/p;
                                                   >> 496       G4double sint1 = std::sqrt(std::abs((1. - cost )*(1.+cost)));
                                                   >> 497       G4double sint2 = std::sqrt(pt2)/p;
                                                   >> 498       G4double sint = 0.5*(sint1 + sint2);
                                                   >> 499       G4double ph = pi*0.5;
                                                   >> 500       if (pyinc < 0.) ph = pi*1.5;
                                                   >> 501       if (std::abs(pxinc) > 1.e-6) ph = std::atan2(pyinc, pxinc);
                                                   >> 502       G4double cosp = std::cos(ph);
                                                   >> 503       G4double sinp = std::sin(ph);
                                                   >> 504       if (verboseLevel > 1) {
                                                   >> 505          G4cout << "cost sint " << cost << " " << sint << G4endl;
                                                   >> 506          G4cout << "cosp sinp " << cosp << " " << sinp << G4endl;
                                                   >> 507       }
                                                   >> 508       *pxnew = cost*cosp*px - sinp*py + sint*cosp*pz;
                                                   >> 509       *pynew = cost*sinp*px + cosp*py + sint*sinp*pz;
                                                   >> 510       *pznew =     -sint*px                 +cost*pz;
                                                   >> 511    }
                                                   >> 512    else {
                                                   >> 513        G4double cost=pzinc/p;
                                                   >> 514        *pxnew = cost*px;
                                                   >> 515        *pynew = py;
                                                   >> 516        *pznew = cost*pz;
                                                   >> 517    }
                                                   >> 518 }
388                                                   519