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 11.0.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 //                                                 26 //
 27 // Geant4 Header : G4HadronElastic                 27 // Geant4 Header : G4HadronElastic
 28 //                                                 28 //
 29 // Author : V.Ivanchenko 29 June 2009 (redesig     29 // Author : V.Ivanchenko 29 June 2009 (redesign old elastic model)
 30 //                                                 30 //  
 31                                                    31 
 32 #include "G4HadronElastic.hh"                      32 #include "G4HadronElastic.hh"
 33 #include "G4SystemOfUnits.hh"                      33 #include "G4SystemOfUnits.hh"
 34 #include "G4ParticleTable.hh"                      34 #include "G4ParticleTable.hh"
 35 #include "G4ParticleDefinition.hh"                 35 #include "G4ParticleDefinition.hh"
 36 #include "G4IonTable.hh"                           36 #include "G4IonTable.hh"
 37 #include "Randomize.hh"                            37 #include "Randomize.hh"
 38 #include "G4Proton.hh"                             38 #include "G4Proton.hh"
 39 #include "G4Neutron.hh"                            39 #include "G4Neutron.hh"
 40 #include "G4Deuteron.hh"                           40 #include "G4Deuteron.hh"
 41 #include "G4Alpha.hh"                              41 #include "G4Alpha.hh"
 42 #include "G4Pow.hh"                                42 #include "G4Pow.hh"
 43 #include "G4Exp.hh"                                43 #include "G4Exp.hh"
 44 #include "G4Log.hh"                                44 #include "G4Log.hh"
 45 #include "G4HadronicParameters.hh"                 45 #include "G4HadronicParameters.hh"
 46 #include "G4PhysicsModelCatalog.hh"                46 #include "G4PhysicsModelCatalog.hh"
 47                                                    47 
 48                                                    48 
 49 G4HadronElastic::G4HadronElastic(const G4Strin     49 G4HadronElastic::G4HadronElastic(const G4String& name) 
 50   : G4HadronicInteraction(name), secID(-1)         50   : G4HadronicInteraction(name), secID(-1)
 51 {                                                  51 {
 52   SetMinEnergy( 0.0*GeV );                         52   SetMinEnergy( 0.0*GeV );
 53   SetMaxEnergy( G4HadronicParameters::Instance     53   SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
 54   lowestEnergyLimit= 1.e-6*eV;                     54   lowestEnergyLimit= 1.e-6*eV;
 55   pLocalTmax  = 0.0;                               55   pLocalTmax  = 0.0;
 56   nwarn = 0;                                       56   nwarn = 0;
 57                                                    57 
 58   theProton   = G4Proton::Proton();                58   theProton   = G4Proton::Proton();
 59   theNeutron  = G4Neutron::Neutron();              59   theNeutron  = G4Neutron::Neutron();
 60   theDeuteron = G4Deuteron::Deuteron();            60   theDeuteron = G4Deuteron::Deuteron();
 61   theAlpha    = G4Alpha::Alpha();                  61   theAlpha    = G4Alpha::Alpha();
 62                                                    62 
 63   secID = G4PhysicsModelCatalog::GetModelID( "     63   secID = G4PhysicsModelCatalog::GetModelID( "model_" + name );
 64 }                                                  64 }
 65                                                    65 
 66 G4HadronElastic::~G4HadronElastic()                66 G4HadronElastic::~G4HadronElastic()
 67 {}                                                 67 {}
 68                                                    68 
 69                                                    69 
 70 void G4HadronElastic::ModelDescription(std::os     70 void G4HadronElastic::ModelDescription(std::ostream& outFile) const
 71 {                                                  71 {
 72   outFile << "G4HadronElastic is the base clas     72   outFile << "G4HadronElastic is the base class for all hadron-nucleus\n" 
 73           << "elastic scattering models except     73           << "elastic scattering models except HP.\n" 
 74           << "By default it uses the Gheisha t     74           << "By default it uses the Gheisha two-exponential momentum\n"
 75     << "transfer parameterization.  The model      75     << "transfer parameterization.  The model is fully relativistic\n"
 76     << "as opposed to the original Gheisha mod     76     << "as opposed to the original Gheisha model which was not.\n"
 77     << "This model may be used for all long-li     77     << "This model may be used for all long-lived hadrons at all\n"
 78     << "incident energies but fit the data onl     78     << "incident energies but fit the data only for relativistic scattering.\n";
 79 }                                                  79 }
 80                                                    80 
 81 G4HadFinalState* G4HadronElastic::ApplyYoursel     81 G4HadFinalState* G4HadronElastic::ApplyYourself(
 82      const G4HadProjectile& aTrack, G4Nucleus&     82      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
 83 {                                                  83 {
 84   theParticleChange.Clear();                       84   theParticleChange.Clear();
 85                                                    85 
 86   const G4HadProjectile* aParticle = &aTrack;      86   const G4HadProjectile* aParticle = &aTrack;
 87   G4double ekin = aParticle->GetKineticEnergy(     87   G4double ekin = aParticle->GetKineticEnergy();
 88                                                    88 
 89   // no scattering below the limit                 89   // no scattering below the limit
 90   if(ekin <= lowestEnergyLimit) {                  90   if(ekin <= lowestEnergyLimit) {
 91     theParticleChange.SetEnergyChange(ekin);       91     theParticleChange.SetEnergyChange(ekin);
 92     theParticleChange.SetMomentumChange(0.,0.,     92     theParticleChange.SetMomentumChange(0.,0.,1.);
 93     return &theParticleChange;                     93     return &theParticleChange;
 94   }                                                94   }
 95                                                    95 
 96   G4int A = targetNucleus.GetA_asInt();            96   G4int A = targetNucleus.GetA_asInt();
 97   G4int Z = targetNucleus.GetZ_asInt();            97   G4int Z = targetNucleus.GetZ_asInt();
 98                                                    98 
 99   // Scattered particle referred to axis of in     99   // Scattered particle referred to axis of incident particle
100   const G4ParticleDefinition* theParticle = aP    100   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
101   G4double m1 = theParticle->GetPDGMass();        101   G4double m1 = theParticle->GetPDGMass();
102   G4double plab = std::sqrt(ekin*(ekin + 2.0*m    102   G4double plab = std::sqrt(ekin*(ekin + 2.0*m1));
103                                                   103 
104   if (verboseLevel>1) {                           104   if (verboseLevel>1) {
105     G4cout << "G4HadronElastic: "                 105     G4cout << "G4HadronElastic: " 
106      << aParticle->GetDefinition()->GetParticl    106      << aParticle->GetDefinition()->GetParticleName() 
107      << " Plab(GeV/c)= " << plab/GeV              107      << " Plab(GeV/c)= " << plab/GeV  
108      << " Ekin(MeV) = " << ekin/MeV               108      << " Ekin(MeV) = " << ekin/MeV 
109      << " scattered off Z= " << Z                 109      << " scattered off Z= " << Z 
110      << " A= " << A                               110      << " A= " << A 
111      << G4endl;                                   111      << G4endl;
112   }                                               112   }
113                                                   113 
114   G4double mass2 = G4NucleiProperties::GetNucl    114   G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
115   G4double e1 = m1 + ekin;                        115   G4double e1 = m1 + ekin;
116   G4LorentzVector lv(0.0,0.0,plab,e1+mass2);      116   G4LorentzVector lv(0.0,0.0,plab,e1+mass2);
117   G4ThreeVector bst = lv.boostVector();           117   G4ThreeVector bst = lv.boostVector();
118   G4double momentumCMS = plab*mass2/std::sqrt(    118   G4double momentumCMS = plab*mass2/std::sqrt(m1*m1 + mass2*mass2 + 2.*mass2*e1);
119                                                   119 
120   pLocalTmax = 4.0*momentumCMS*momentumCMS;       120   pLocalTmax = 4.0*momentumCMS*momentumCMS;
121                                                   121 
122   // Sampling in CM system                        122   // Sampling in CM system
123   G4double t = SampleInvariantT(theParticle, p    123   G4double t = SampleInvariantT(theParticle, plab, Z, A);
124                                                   124 
125   if(t < 0.0 || t > pLocalTmax) {                 125   if(t < 0.0 || t > pLocalTmax) {
126     // For the very rare cases where cos(theta    126     // For the very rare cases where cos(theta) is greater than 1 or smaller than -1,
127     // print some debugging information via a     127     // print some debugging information via a "JustWarning" exception, and resample
128     // using the default algorithm                128     // using the default algorithm
129 #ifdef G4VERBOSE                                  129 #ifdef G4VERBOSE
130     if(nwarn < 2) {                               130     if(nwarn < 2) {
131       G4ExceptionDescription ed;                  131       G4ExceptionDescription ed;
132       ed << GetModelName() << " wrong sampling    132       ed << GetModelName() << " wrong sampling t= " << t << " tmax= " << pLocalTmax
133    << " for " << aParticle->GetDefinition()->G    133    << " for " << aParticle->GetDefinition()->GetParticleName() 
134    << " ekin=" << ekin << " MeV"                  134    << " ekin=" << ekin << " MeV" 
135    << " off (Z,A)=(" << Z << "," << A << ") -     135    << " off (Z,A)=(" << Z << "," << A << ") - will be resampled" << G4endl;
136       G4Exception( "G4HadronElastic::ApplyYour    136       G4Exception( "G4HadronElastic::ApplyYourself", "hadEla001", JustWarning, ed);
137       ++nwarn;                                    137       ++nwarn;
138     }                                             138     }
139 #endif                                            139 #endif
140     t = G4HadronElastic::SampleInvariantT(theP    140     t = G4HadronElastic::SampleInvariantT(theParticle, plab, Z, A);
141   }                                               141   }
142                                                   142 
143   G4double phi  = G4UniformRand()*CLHEP::twopi    143   G4double phi  = G4UniformRand()*CLHEP::twopi;
144   G4double cost = 1. - 2.0*t/pLocalTmax;          144   G4double cost = 1. - 2.0*t/pLocalTmax;
145                                                   145 
146   if (cost > 1.0) { cost = 1.0; }                 146   if (cost > 1.0) { cost = 1.0; }
147   else if(cost < -1.0) { cost = -1.0; }           147   else if(cost < -1.0) { cost = -1.0; } 
148                                                   148 
149   G4double sint = std::sqrt((1.0-cost)*(1.0+co    149   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
150                                                   150 
151   if (verboseLevel>1) {                           151   if (verboseLevel>1) {
152     G4cout << " t= " << t << " tmax(GeV^2)= "     152     G4cout << " t= " << t << " tmax(GeV^2)= " << pLocalTmax/(GeV*GeV) 
153      << " Pcms(GeV)= " << momentumCMS/GeV << "    153      << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost 
154      << " sin(t)=" << sint << G4endl;             154      << " sin(t)=" << sint << G4endl;
155   }                                               155   }
156   G4LorentzVector nlv1(momentumCMS*sint*std::c    156   G4LorentzVector nlv1(momentumCMS*sint*std::cos(phi),
157            momentumCMS*sint*std::sin(phi),        157            momentumCMS*sint*std::sin(phi),
158                        momentumCMS*cost,          158                        momentumCMS*cost,
159            std::sqrt(momentumCMS*momentumCMS +    159            std::sqrt(momentumCMS*momentumCMS + m1*m1));
160                                                   160 
161   nlv1.boost(bst);                                161   nlv1.boost(bst); 
162                                                   162 
163   G4double eFinal = nlv1.e() - m1;                163   G4double eFinal = nlv1.e() - m1;
164   if (verboseLevel > 1) {                         164   if (verboseLevel > 1) {
165     G4cout <<"G4HadronElastic: m= " << m1 << "    165     G4cout <<"G4HadronElastic: m= " << m1 << " Efin(MeV)= " << eFinal 
166      << " 4-M Final: " << nlv1                    166      << " 4-M Final: " << nlv1 
167      << G4endl;                                   167      << G4endl;
168   }                                               168   }
169                                                   169 
170   if(eFinal <= 0.0) {                             170   if(eFinal <= 0.0) { 
171     theParticleChange.SetMomentumChange(0.0,0.    171     theParticleChange.SetMomentumChange(0.0,0.0,1.0);
172     theParticleChange.SetEnergyChange(0.0);       172     theParticleChange.SetEnergyChange(0.0);
173   } else {                                        173   } else {
174     theParticleChange.SetMomentumChange(nlv1.v    174     theParticleChange.SetMomentumChange(nlv1.vect().unit());
175     theParticleChange.SetEnergyChange(eFinal);    175     theParticleChange.SetEnergyChange(eFinal);
176   }                                               176   }
177   lv -= nlv1;                                     177   lv -= nlv1;
178   G4double erec =  std::max(lv.e() - mass2, 0.    178   G4double erec =  std::max(lv.e() - mass2, 0.0);
179   if (verboseLevel > 1) {                         179   if (verboseLevel > 1) {
180     G4cout << "Recoil: " <<" m= " << mass2 <<     180     G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
181      << " 4-mom: " << lv                          181      << " 4-mom: " << lv 
182      << G4endl;                                   182      << G4endl;
183   }                                               183   }
184                                                   184  
185   // the recoil is created if kinetic energy a    185   // the recoil is created if kinetic energy above the threshold
186   if(erec > GetRecoilEnergyThreshold()) {         186   if(erec > GetRecoilEnergyThreshold()) {
187     G4ParticleDefinition * theDef = nullptr;      187     G4ParticleDefinition * theDef = nullptr;
188     if(Z == 1 && A == 1)       { theDef = theP    188     if(Z == 1 && A == 1)       { theDef = theProton; }
189     else if (Z == 1 && A == 2) { theDef = theD    189     else if (Z == 1 && A == 2) { theDef = theDeuteron; }
190     else if (Z == 1 && A == 3) { theDef = G4Tr    190     else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
191     else if (Z == 2 && A == 3) { theDef = G4He    191     else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
192     else if (Z == 2 && A == 4) { theDef = theA    192     else if (Z == 2 && A == 4) { theDef = theAlpha; }
193     else {                                        193     else {
194       theDef =                                    194       theDef = 
195   G4ParticleTable::GetParticleTable()->GetIonT    195   G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
196     }                                             196     }
197     G4DynamicParticle * aSec = new G4DynamicPa    197     G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv.vect().unit(), erec);
198     theParticleChange.AddSecondary(aSec, secID    198     theParticleChange.AddSecondary(aSec, secID);
199   } else {                                        199   } else {
200     theParticleChange.SetLocalEnergyDeposit(er    200     theParticleChange.SetLocalEnergyDeposit(erec);
201   }                                               201   }
202                                                   202 
203   return &theParticleChange;                      203   return &theParticleChange;
204 }                                                 204 }
205                                                   205 
206 // sample momentum transfer in the CMS system     206 // sample momentum transfer in the CMS system 
207 G4double                                          207 G4double 
208 G4HadronElastic::SampleInvariantT(const G4Part    208 G4HadronElastic::SampleInvariantT(const G4ParticleDefinition* part,
209           G4double mom, G4int, G4int A)           209           G4double mom, G4int, G4int A)
210 {                                                 210 {
211   const G4double plabLowLimit = 400.0*CLHEP::M    211   const G4double plabLowLimit = 400.0*CLHEP::MeV;
212   const G4double GeV2 = GeV*GeV;                  212   const G4double GeV2 = GeV*GeV;
213   const G4double z07in13 = std::pow(0.7, 0.333    213   const G4double z07in13 = std::pow(0.7, 0.3333333333);
214   const G4double numLimit = 18.;                  214   const G4double numLimit = 18.;
215                                                   215 
216   G4int pdg = std::abs(part->GetPDGEncoding())    216   G4int pdg = std::abs(part->GetPDGEncoding());
217   G4double tmax = pLocalTmax/GeV2;                217   G4double tmax = pLocalTmax/GeV2;
218                                                   218 
219   G4double aa, bb, cc, dd;                        219   G4double aa, bb, cc, dd;
220   G4Pow* g4pow = G4Pow::GetInstance();            220   G4Pow* g4pow = G4Pow::GetInstance();
221   if (A <= 62) {                                  221   if (A <= 62) {
222     if (pdg == 211){ //Pions                      222     if (pdg == 211){ //Pions
223       if(mom >= plabLowLimit){     //High ener    223       if(mom >= plabLowLimit){     //High energy
224   bb = 14.5*g4pow->Z23(A);/*14.5*/                224   bb = 14.5*g4pow->Z23(A);/*14.5*/
225   dd = 10.;                                       225   dd = 10.;
226   cc = 0.075*g4pow->Z13(A)/dd;//1.4               226   cc = 0.075*g4pow->Z13(A)/dd;//1.4
227   //aa = g4pow->powZ(A, 1.93)/bb;//1.63           227   //aa = g4pow->powZ(A, 1.93)/bb;//1.63
228   aa = (A*A)/bb;//1.63                            228   aa = (A*A)/bb;//1.63
229       } else {                       //Low ene    229       } else {                       //Low energy
230   bb = 29.*z07in13*z07in13*g4pow->Z23(A);         230   bb = 29.*z07in13*z07in13*g4pow->Z23(A);
231   dd = 15.;                                       231   dd = 15.;
232   cc = 0.04*g4pow->Z13(A)/dd;//1.4                232   cc = 0.04*g4pow->Z13(A)/dd;//1.4
233   aa = g4pow->powZ(A, 1.63)/bb;//1.63             233   aa = g4pow->powZ(A, 1.63)/bb;//1.63
234       }                                           234       }
235     } else { //Other particles                    235     } else { //Other particles
236       bb = 14.5*g4pow->Z23(A);                    236       bb = 14.5*g4pow->Z23(A);
237       dd = 20.;                                   237       dd = 20.;
238       aa = (A*A)/bb;//1.63                        238       aa = (A*A)/bb;//1.63
239       cc = 1.4*g4pow->Z13(A)/dd;                  239       cc = 1.4*g4pow->Z13(A)/dd;          
240     }                                             240     }
241       //===========================               241       //===========================
242   } else { //(A>62)                               242   } else { //(A>62)
243     if (pdg == 211) {                             243     if (pdg == 211) {
244       if(mom >= plabLowLimit){ //high             244       if(mom >= plabLowLimit){ //high
245   bb = 60.*z07in13*g4pow->Z13(A);//60             245   bb = 60.*z07in13*g4pow->Z13(A);//60
246   dd = 30.;                                       246   dd = 30.;
247   aa = 0.5*(A*A)/bb;//1.33                        247   aa = 0.5*(A*A)/bb;//1.33
248   cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4     --    248   cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4     ---    2: 0.4
249       } else { //low                              249       } else { //low
250   bb = 120.*z07in13*g4pow->Z13(A);//60            250   bb = 120.*z07in13*g4pow->Z13(A);//60
251   dd = 30.;                                       251   dd = 30.;
252   aa = 2.*g4pow->powZ(A,1.33)/bb;                 252   aa = 2.*g4pow->powZ(A,1.33)/bb;
253   cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4     --    253   cc = 4.*g4pow->powZ(A,0.4)/dd;//1:0.4     ---    2: 0.4
254       }                                           254       }
255     } else {                                      255     } else {
256       bb = 60.*g4pow->Z13(A);                     256       bb = 60.*g4pow->Z13(A);
257       dd = 25.;                                   257       dd = 25.;
258       aa = g4pow->powZ(A,1.33)/bb;//1.33          258       aa = g4pow->powZ(A,1.33)/bb;//1.33
259       cc = 0.2*g4pow->powZ(A,0.4)/dd;//1:0.4      259       cc = 0.2*g4pow->powZ(A,0.4)/dd;//1:0.4     ---    2: 0.4
260     }                                             260     }
261   }                                               261   }
262   G4double q1 = 1.0 - G4Exp(-std::min(bb*tmax,    262   G4double q1 = 1.0 - G4Exp(-std::min(bb*tmax, numLimit));
263   G4double q2 = 1.0 - G4Exp(-std::min(dd*tmax,    263   G4double q2 = 1.0 - G4Exp(-std::min(dd*tmax, numLimit));
264   G4double s1 = q1*aa;                            264   G4double s1 = q1*aa;
265   G4double s2 = q2*cc;                            265   G4double s2 = q2*cc;
266   if((s1 + s2)*G4UniformRand() < s2) {            266   if((s1 + s2)*G4UniformRand() < s2) {
267     q1 = q2;                                      267     q1 = q2;
268     bb = dd;                                      268     bb = dd;
269   }                                               269   }
270   return -GeV2*G4Log(1.0 - G4UniformRand()*q1)    270   return -GeV2*G4Log(1.0 - G4UniformRand()*q1)/bb;
271 }                                                 271 }
272                                                   272 
273 //////////////////////////////////////////////    273 //////////////////////////////////////////////
274 //                                                274 //
275 // Cofs for s-,c-,b-particles ds/dt slopes        275 // Cofs for s-,c-,b-particles ds/dt slopes
276                                                   276 
277 G4double G4HadronElastic::GetSlopeCof(const G4    277 G4double G4HadronElastic::GetSlopeCof(const G4int pdg )
278 {                                                 278 {
279   // The input parameter "pdg" should be the a    279   // The input parameter "pdg" should be the absolute value of the PDG code
280   // (i.e. the same value for a particle and i    280   // (i.e. the same value for a particle and its antiparticle).
281                                                   281 
282   G4double coeff = 1.0;                           282   G4double coeff = 1.0;
283                                                   283 
284   // heavy barions                                284   // heavy barions
285                                                   285 
286   static const G4double  lBarCof1S  = 0.88;       286   static const G4double  lBarCof1S  = 0.88;
287   static const G4double  lBarCof2S  = 0.76;       287   static const G4double  lBarCof2S  = 0.76;
288   static const G4double  lBarCof3S  = 0.64;       288   static const G4double  lBarCof3S  = 0.64;
289   static const G4double  lBarCof1C  = 0.784378    289   static const G4double  lBarCof1C  = 0.784378;
290   static const G4double  lBarCofSC  = 0.664378    290   static const G4double  lBarCofSC  = 0.664378;
291   static const G4double  lBarCof2SC = 0.544378    291   static const G4double  lBarCof2SC = 0.544378;
292   static const G4double  lBarCof1B  = 0.740659    292   static const G4double  lBarCof1B  = 0.740659;
293   static const G4double  lBarCofSB  = 0.620659    293   static const G4double  lBarCofSB  = 0.620659;
294   static const G4double  lBarCof2SB = 0.500659    294   static const G4double  lBarCof2SB = 0.500659;
295                                                   295   
296   if( pdg == 3122 || pdg == 3222 ||  pdg == 31    296   if( pdg == 3122 || pdg == 3222 ||  pdg == 3112 || pdg == 3212  )
297   {                                               297   {
298     coeff = lBarCof1S; // Lambda, Sigma+, Sigm    298     coeff = lBarCof1S; // Lambda, Sigma+, Sigma-, Sigma0
299                                                   299 
300   } else if( pdg == 3322 || pdg == 3312   )       300   } else if( pdg == 3322 || pdg == 3312   )
301   {                                               301   {
302     coeff = lBarCof2S; // Xi-, Xi0                302     coeff = lBarCof2S; // Xi-, Xi0
303   }                                               303   }
304   else if( pdg == 3324)                           304   else if( pdg == 3324)
305   {                                               305   {
306     coeff = lBarCof3S; // Omega                   306     coeff = lBarCof3S; // Omega
307   }                                               307   }
308   else if( pdg == 4122 ||  pdg == 4212 ||   pd    308   else if( pdg == 4122 ||  pdg == 4212 ||   pdg == 4222 ||   pdg == 4112   )
309   {                                               309   {
310     coeff = lBarCof1C; // LambdaC+, SigmaC+, S    310     coeff = lBarCof1C; // LambdaC+, SigmaC+, SigmaC++, SigmaC0
311   }                                               311   }
312   else if( pdg == 4332 )                          312   else if( pdg == 4332 )
313   {                                               313   {
314     coeff = lBarCof2SC; // OmegaC                 314     coeff = lBarCof2SC; // OmegaC
315   }                                               315   }
316   else if( pdg == 4232 || pdg == 4132 )           316   else if( pdg == 4232 || pdg == 4132 )
317   {                                               317   {
318     coeff = lBarCofSC; // XiC+, XiC0              318     coeff = lBarCofSC; // XiC+, XiC0
319   }                                               319   }
320   else if( pdg == 5122 || pdg == 5222 || pdg =    320   else if( pdg == 5122 || pdg == 5222 || pdg == 5112 || pdg == 5212    )
321   {                                               321   {
322     coeff = lBarCof1B; // LambdaB, SigmaB+, Si    322     coeff = lBarCof1B; // LambdaB, SigmaB+, SigmaB-, SigmaB0
323   }                                               323   }
324   else if( pdg == 5332 )                          324   else if( pdg == 5332 )
325   {                                               325   {
326     coeff = lBarCof2SB; // OmegaB-                326     coeff = lBarCof2SB; // OmegaB-
327   }                                               327   }
328   else if( pdg == 5132 || pdg == 5232 ) // XiB    328   else if( pdg == 5132 || pdg == 5232 ) // XiB-, XiB0
329   {                                               329   {
330     coeff = lBarCofSB;                            330     coeff = lBarCofSB;
331   }                                               331   }
332   // heavy mesons Kaons?                          332   // heavy mesons Kaons?
333   static const G4double lMesCof1S = 0.82; // K    333   static const G4double lMesCof1S = 0.82; // Kp/piP kaons?
334   static const G4double llMesCof1C = 0.676568;    334   static const G4double llMesCof1C = 0.676568;
335   static const G4double llMesCof1B = 0.610989;    335   static const G4double llMesCof1B = 0.610989;
336   static const G4double llMesCof2C = 0.353135;    336   static const G4double llMesCof2C = 0.353135;
337   static const G4double llMesCof2B = 0.221978;    337   static const G4double llMesCof2B = 0.221978;
338   static const G4double llMesCofSC = 0.496568;    338   static const G4double llMesCofSC = 0.496568;
339   static const G4double llMesCofSB = 0.430989;    339   static const G4double llMesCofSB = 0.430989;
340   static const G4double llMesCofCB = 0.287557;    340   static const G4double llMesCofCB = 0.287557;
341   static const G4double llMesCofEtaP = 0.88;      341   static const G4double llMesCofEtaP = 0.88;
342   static const G4double llMesCofEta = 0.76;       342   static const G4double llMesCofEta = 0.76;
343                                                   343 
344   if( pdg == 321 || pdg == 311 || pdg == 310 )    344   if( pdg == 321 || pdg == 311 || pdg == 310 )
345   {                                               345   {
346     coeff = lMesCof1S; //K+-0                     346     coeff = lMesCof1S; //K+-0
347   }                                               347   }
348   else if( pdg == 511 ||  pdg == 521  )           348   else if( pdg == 511 ||  pdg == 521  )
349   {                                               349   {
350     coeff = llMesCof1B; // BMeson0, BMeson+       350     coeff = llMesCof1B; // BMeson0, BMeson+
351   }                                               351   }
352   else if(pdg == 421 ||  pdg == 411 )             352   else if(pdg == 421 ||  pdg == 411 )
353   {                                               353   {
354     coeff = llMesCof1C; // DMeson+, DMeson0       354     coeff = llMesCof1C; // DMeson+, DMeson0
355   }                                               355   }
356   else if( pdg == 531  )                          356   else if( pdg == 531  )
357   {                                               357   {
358     coeff = llMesCofSB; // BSMeson0               358     coeff = llMesCofSB; // BSMeson0
359   }                                               359   }
360   else if( pdg == 541 )                           360   else if( pdg == 541 )
361   {                                               361   {
362     coeff = llMesCofCB; // BCMeson+-              362     coeff = llMesCofCB; // BCMeson+-
363   }                                               363   }
364   else if(pdg == 431 )                            364   else if(pdg == 431 ) 
365   {                                               365   {
366     coeff = llMesCofSC; // DSMeson+-              366     coeff = llMesCofSC; // DSMeson+-
367   }                                               367   }
368   else if(pdg == 441 || pdg == 443 )              368   else if(pdg == 441 || pdg == 443 )
369   {                                               369   {
370     coeff = llMesCof2C; // Etac, JPsi             370     coeff = llMesCof2C; // Etac, JPsi
371   }                                               371   }
372   else if(pdg == 553 )                            372   else if(pdg == 553 )
373   {                                               373   {
374     coeff = llMesCof2B; // Upsilon                374     coeff = llMesCof2B; // Upsilon
375   }                                               375   }
376   else if(pdg == 221 )                            376   else if(pdg == 221 )
377   {                                               377   {
378     coeff = llMesCofEta; // Eta                   378     coeff = llMesCofEta; // Eta
379   }                                               379   }
380   else if(pdg == 331 )                            380   else if(pdg == 331 )
381   {                                               381   {
382     coeff = llMesCofEtaP; // Eta'                 382     coeff = llMesCofEtaP; // Eta'
383   }                                               383   } 
384   return coeff;                                   384   return coeff;
385 }                                                 385 }
386                                                   386 
387                                                   387 
388                                                   388