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 10.5.p1)


  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"            << 
 47                                                    46 
 48                                                    47 
 49 G4HadronElastic::G4HadronElastic(const G4Strin     48 G4HadronElastic::G4HadronElastic(const G4String& name) 
 50   : G4HadronicInteraction(name), secID(-1)     <<  49   : G4HadronicInteraction(name)
 51 {                                                  50 {
 52   SetMinEnergy( 0.0*GeV );                         51   SetMinEnergy( 0.0*GeV );
 53   SetMaxEnergy( G4HadronicParameters::Instance     52   SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
 54   lowestEnergyLimit= 1.e-6*eV;                 <<  53   lowestEnergyLimit= 1.e-6*eV;  
 55   pLocalTmax  = 0.0;                           << 
 56   nwarn = 0;                                   << 
 57                                                    54 
 58   theProton   = G4Proton::Proton();                55   theProton   = G4Proton::Proton();
 59   theNeutron  = G4Neutron::Neutron();              56   theNeutron  = G4Neutron::Neutron();
 60   theDeuteron = G4Deuteron::Deuteron();            57   theDeuteron = G4Deuteron::Deuteron();
 61   theAlpha    = G4Alpha::Alpha();                  58   theAlpha    = G4Alpha::Alpha();
 62                                                << 
 63   secID = G4PhysicsModelCatalog::GetModelID( " << 
 64 }                                                  59 }
 65                                                    60 
 66 G4HadronElastic::~G4HadronElastic()                61 G4HadronElastic::~G4HadronElastic()
 67 {}                                                 62 {}
 68                                                    63 
 69                                                    64 
 70 void G4HadronElastic::ModelDescription(std::os     65 void G4HadronElastic::ModelDescription(std::ostream& outFile) const
 71 {                                                  66 {
 72   outFile << "G4HadronElastic is the base clas <<  67 
 73           << "elastic scattering models except <<  68     outFile << "G4HadronElastic is a hadron-nucleus elastic scattering\n"
 74           << "By default it uses the Gheisha t <<  69             << "model which uses the Gheisha two-exponential momentum\n"
 75     << "transfer parameterization.  The model  <<  70             << "transfer parameterization.  The model is fully relativistic\n"
 76     << "as opposed to the original Gheisha mod <<  71             << "as opposed to the original Gheisha model which was not.\n"
 77     << "This model may be used for all long-li <<  72             << "This model may be used for all long-lived hadrons at all\n"
 78     << "incident energies but fit the data onl <<  73             << "incident energies.\n";
                                                   >>  74 
 79 }                                                  75 }
 80                                                    76 
 81 G4HadFinalState* G4HadronElastic::ApplyYoursel     77 G4HadFinalState* G4HadronElastic::ApplyYourself(
 82      const G4HadProjectile& aTrack, G4Nucleus&     78      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
 83 {                                                  79 {
 84   theParticleChange.Clear();                       80   theParticleChange.Clear();
 85                                                    81 
 86   const G4HadProjectile* aParticle = &aTrack;      82   const G4HadProjectile* aParticle = &aTrack;
 87   G4double ekin = aParticle->GetKineticEnergy(     83   G4double ekin = aParticle->GetKineticEnergy();
 88                                                << 
 89   // no scattering below the limit             << 
 90   if(ekin <= lowestEnergyLimit) {                  84   if(ekin <= lowestEnergyLimit) {
 91     theParticleChange.SetEnergyChange(ekin);       85     theParticleChange.SetEnergyChange(ekin);
 92     theParticleChange.SetMomentumChange(0.,0., <<  86     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
 93     return &theParticleChange;                     87     return &theParticleChange;
 94   }                                                88   }
 95                                                    89 
 96   G4int A = targetNucleus.GetA_asInt();            90   G4int A = targetNucleus.GetA_asInt();
 97   G4int Z = targetNucleus.GetZ_asInt();            91   G4int Z = targetNucleus.GetZ_asInt();
 98                                                    92 
                                                   >>  93   G4double plab = aParticle->GetTotalMomentum();
                                                   >>  94 
 99   // Scattered particle referred to axis of in     95   // Scattered particle referred to axis of incident particle
100   const G4ParticleDefinition* theParticle = aP     96   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
101   G4double m1 = theParticle->GetPDGMass();         97   G4double m1 = theParticle->GetPDGMass();
102   G4double plab = std::sqrt(ekin*(ekin + 2.0*m << 
103                                                    98 
104   if (verboseLevel>1) {                            99   if (verboseLevel>1) {
105     G4cout << "G4HadronElastic: "                 100     G4cout << "G4HadronElastic: " 
106      << aParticle->GetDefinition()->GetParticl    101      << aParticle->GetDefinition()->GetParticleName() 
107      << " Plab(GeV/c)= " << plab/GeV              102      << " Plab(GeV/c)= " << plab/GeV  
108      << " Ekin(MeV) = " << ekin/MeV               103      << " Ekin(MeV) = " << ekin/MeV 
109      << " scattered off Z= " << Z                 104      << " scattered off Z= " << Z 
110      << " A= " << A                               105      << " A= " << A 
111      << G4endl;                                   106      << G4endl;
112   }                                               107   }
113                                                   108 
114   G4double mass2 = G4NucleiProperties::GetNucl    109   G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
115   G4double e1 = m1 + ekin;                     << 110   G4LorentzVector lv1 = aParticle->Get4Momentum();
116   G4LorentzVector lv(0.0,0.0,plab,e1+mass2);   << 111   G4LorentzVector lv(0.0,0.0,0.0,mass2);   
                                                   >> 112   lv += lv1;
                                                   >> 113 
117   G4ThreeVector bst = lv.boostVector();           114   G4ThreeVector bst = lv.boostVector();
118   G4double momentumCMS = plab*mass2/std::sqrt( << 115   lv1.boost(-bst);
119                                                   116 
120   pLocalTmax = 4.0*momentumCMS*momentumCMS;    << 117   G4ThreeVector p1 = lv1.vect();
                                                   >> 118   G4double momentumCMS = p1.mag();
                                                   >> 119   G4double tmax = 4.0*momentumCMS*momentumCMS;
121                                                   120 
122   // Sampling in CM system                        121   // Sampling in CM system
123   G4double t = SampleInvariantT(theParticle, p << 122   G4double t    = SampleInvariantT(theParticle, plab, Z, A);
124                                                << 
125   if(t < 0.0 || t > pLocalTmax) {              << 
126     // For the very rare cases where cos(theta << 
127     // print some debugging information via a  << 
128     // using the default algorithm             << 
129 #ifdef G4VERBOSE                               << 
130     if(nwarn < 2) {                            << 
131       G4ExceptionDescription ed;               << 
132       ed << GetModelName() << " wrong sampling << 
133    << " for " << aParticle->GetDefinition()->G << 
134    << " ekin=" << ekin << " MeV"               << 
135    << " off (Z,A)=(" << Z << "," << A << ") -  << 
136       G4Exception( "G4HadronElastic::ApplyYour << 
137       ++nwarn;                                 << 
138     }                                          << 
139 #endif                                         << 
140     t = G4HadronElastic::SampleInvariantT(theP << 
141   }                                            << 
142                                                << 
143   G4double phi  = G4UniformRand()*CLHEP::twopi    123   G4double phi  = G4UniformRand()*CLHEP::twopi;
144   G4double cost = 1. - 2.0*t/pLocalTmax;       << 124   G4double cost = 1. - 2.0*t/tmax;
145                                                << 125   if(cost > 1.0) { cost = 1.0; }
146   if (cost > 1.0) { cost = 1.0; }              << 126   else if(cost < -1.0) { cost = -1.0; }
147   else if(cost < -1.0) { cost = -1.0; }        << 
148                                                << 
149   G4double sint = std::sqrt((1.0-cost)*(1.0+co    127   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
150                                                   128 
151   if (verboseLevel>1) {                           129   if (verboseLevel>1) {
152     G4cout << " t= " << t << " tmax(GeV^2)= "  << 130     G4cout << " t= " << t << " tmax(GeV^2)= " << tmax/(GeV*GeV) 
153      << " Pcms(GeV)= " << momentumCMS/GeV << "    131      << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost 
154      << " sin(t)=" << sint << G4endl;             132      << " sin(t)=" << sint << G4endl;
155   }                                               133   }
156   G4LorentzVector nlv1(momentumCMS*sint*std::c << 134   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
157            momentumCMS*sint*std::sin(phi),     << 135   v1 *= momentumCMS;
158                        momentumCMS*cost,       << 136   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),
159            std::sqrt(momentumCMS*momentumCMS +    137            std::sqrt(momentumCMS*momentumCMS + m1*m1));
160                                                   138 
161   nlv1.boost(bst);                                139   nlv1.boost(bst); 
162                                                   140 
163   G4double eFinal = nlv1.e() - m1;                141   G4double eFinal = nlv1.e() - m1;
164   if (verboseLevel > 1) {                         142   if (verboseLevel > 1) {
165     G4cout <<"G4HadronElastic: m= " << m1 << "    143     G4cout <<"G4HadronElastic: m= " << m1 << " Efin(MeV)= " << eFinal 
166      << " 4-M Final: " << nlv1                 << 144      << " Proj: 4-mom " << lv1 << " Final: " << nlv1 
167      << G4endl;                                   145      << G4endl;
168   }                                               146   }
169                                                   147 
170   if(eFinal <= 0.0) {                          << 148   // precision lost in kinematics, only energy is changed
171     theParticleChange.SetMomentumChange(0.0,0. << 149   if (eFinal <= 0.0) {
172     theParticleChange.SetEnergyChange(0.0);    << 150     G4double mom = nlv1.mag();
                                                   >> 151     if(mom == 0.0) {
                                                   >> 152       nlv1.set(0.0,0.0,0.0,m1);
                                                   >> 153       theParticleChange.SetEnergyChange(0.0);
                                                   >> 154     } else {
                                                   >> 155       eFinal = mom*mom/(std::sqrt(m1*m1 + mom*mom) + m1);
                                                   >> 156       theParticleChange.SetEnergyChange(eFinal);
                                                   >> 157       theParticleChange.SetMomentumChange(nlv1.vect().unit());
                                                   >> 158     }
173   } else {                                        159   } else {
174     theParticleChange.SetMomentumChange(nlv1.v    160     theParticleChange.SetMomentumChange(nlv1.vect().unit());
175     theParticleChange.SetEnergyChange(eFinal);    161     theParticleChange.SetEnergyChange(eFinal);
176   }                                               162   }
177   lv -= nlv1;                                     163   lv -= nlv1;
178   G4double erec =  std::max(lv.e() - mass2, 0. << 164   G4double erec =  lv.e() - mass2;
179   if (verboseLevel > 1) {                         165   if (verboseLevel > 1) {
180     G4cout << "Recoil: " <<" m= " << mass2 <<     166     G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
181      << " 4-mom: " << lv                          167      << " 4-mom: " << lv 
182      << G4endl;                                   168      << G4endl;
183   }                                               169   }
184                                                   170  
185   // the recoil is created if kinetic energy a << 
186   if(erec > GetRecoilEnergyThreshold()) {         171   if(erec > GetRecoilEnergyThreshold()) {
187     G4ParticleDefinition * theDef = nullptr;      172     G4ParticleDefinition * theDef = nullptr;
188     if(Z == 1 && A == 1)       { theDef = theP    173     if(Z == 1 && A == 1)       { theDef = theProton; }
189     else if (Z == 1 && A == 2) { theDef = theD    174     else if (Z == 1 && A == 2) { theDef = theDeuteron; }
190     else if (Z == 1 && A == 3) { theDef = G4Tr    175     else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
191     else if (Z == 2 && A == 3) { theDef = G4He    176     else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
192     else if (Z == 2 && A == 4) { theDef = theA    177     else if (Z == 2 && A == 4) { theDef = theAlpha; }
193     else {                                        178     else {
194       theDef =                                    179       theDef = 
195   G4ParticleTable::GetParticleTable()->GetIonT    180   G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
196     }                                             181     }
197     G4DynamicParticle * aSec = new G4DynamicPa << 182     G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv);
198     theParticleChange.AddSecondary(aSec, secID << 183     theParticleChange.AddSecondary(aSec);
199   } else {                                     << 184   } else if(erec > 0.0) {
200     theParticleChange.SetLocalEnergyDeposit(er    185     theParticleChange.SetLocalEnergyDeposit(erec);
201   }                                               186   }
202                                                   187 
203   return &theParticleChange;                      188   return &theParticleChange;
204 }                                                 189 }
205                                                   190 
206 // sample momentum transfer in the CMS system     191 // sample momentum transfer in the CMS system 
207 G4double                                          192 G4double 
208 G4HadronElastic::SampleInvariantT(const G4Part << 193 G4HadronElastic::SampleInvariantT(const G4ParticleDefinition* p, 
209           G4double mom, G4int, G4int A)        << 194           G4double plab,
                                                   >> 195           G4int Z, G4int A)
210 {                                                 196 {
211   const G4double plabLowLimit = 400.0*CLHEP::M << 197   static const G4double GeV2 = GeV*GeV;
212   const G4double GeV2 = GeV*GeV;               << 198   G4double momentumCMS = ComputeMomentumCMS(p,plab,Z,A);
213   const G4double z07in13 = std::pow(0.7, 0.333 << 199   G4double tmax = 4.0*momentumCMS*momentumCMS/GeV2;
214   const G4double numLimit = 18.;               << 200   G4double aa, bb, cc;
215                                                << 201   G4double dd = 10.;
216   G4int pdg = std::abs(part->GetPDGEncoding()) << 
217   G4double tmax = pLocalTmax/GeV2;             << 
218                                                << 
219   G4double aa, bb, cc, dd;                     << 
220   G4Pow* g4pow = G4Pow::GetInstance();            202   G4Pow* g4pow = G4Pow::GetInstance();
221   if (A <= 62) {                                  203   if (A <= 62) {
222     if (pdg == 211){ //Pions                   << 204     bb = 14.5*g4pow->Z23(A);
223       if(mom >= plabLowLimit){     //High ener << 205     aa = g4pow->powZ(A, 1.63)/bb;
224   bb = 14.5*g4pow->Z23(A);/*14.5*/             << 206     cc = 1.4*g4pow->Z13(A)/dd;
225   dd = 10.;                                    << 207   } else {
226   cc = 0.075*g4pow->Z13(A)/dd;//1.4            << 208     bb = 60.*g4pow->Z13(A);
227   //aa = g4pow->powZ(A, 1.93)/bb;//1.63        << 209     aa = g4pow->powZ(A, 1.33)/bb;
228   aa = (A*A)/bb;//1.63                         << 210     cc = 0.4*g4pow->powZ(A, 0.4)/dd;
229       } else {                       //Low ene << 
230   bb = 29.*z07in13*z07in13*g4pow->Z23(A);      << 
231   dd = 15.;                                    << 
232   cc = 0.04*g4pow->Z13(A)/dd;//1.4             << 
233   aa = g4pow->powZ(A, 1.63)/bb;//1.63          << 
234       }                                        << 
235     } else { //Other particles                 << 
236       bb = 14.5*g4pow->Z23(A);                 << 
237       dd = 20.;                                << 
238       aa = (A*A)/bb;//1.63                     << 
239       cc = 1.4*g4pow->Z13(A)/dd;               << 
240     }                                          << 
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   }                                               211   }
262   G4double q1 = 1.0 - G4Exp(-std::min(bb*tmax, << 212   G4double q1 = 1.0 - G4Exp(-bb*tmax);
263   G4double q2 = 1.0 - G4Exp(-std::min(dd*tmax, << 213   G4double q2 = 1.0 - G4Exp(-dd*tmax);
264   G4double s1 = q1*aa;                            214   G4double s1 = q1*aa;
265   G4double s2 = q2*cc;                            215   G4double s2 = q2*cc;
266   if((s1 + s2)*G4UniformRand() < s2) {            216   if((s1 + s2)*G4UniformRand() < s2) {
267     q1 = q2;                                      217     q1 = q2;
268     bb = dd;                                      218     bb = dd;
269   }                                               219   }
270   return -GeV2*G4Log(1.0 - G4UniformRand()*q1)    220   return -GeV2*G4Log(1.0 - G4UniformRand()*q1)/bb;
271 }                                                 221 }
272                                                << 
273 ////////////////////////////////////////////// << 
274 //                                             << 
275 // Cofs for s-,c-,b-particles ds/dt slopes     << 
276                                                << 
277 G4double G4HadronElastic::GetSlopeCof(const G4 << 
278 {                                              << 
279   // The input parameter "pdg" should be the a << 
280   // (i.e. the same value for a particle and i << 
281                                                << 
282   G4double coeff = 1.0;                        << 
283                                                << 
284   // heavy barions                             << 
285                                                << 
286   static const G4double  lBarCof1S  = 0.88;    << 
287   static const G4double  lBarCof2S  = 0.76;    << 
288   static const G4double  lBarCof3S  = 0.64;    << 
289   static const G4double  lBarCof1C  = 0.784378 << 
290   static const G4double  lBarCofSC  = 0.664378 << 
291   static const G4double  lBarCof2SC = 0.544378 << 
292   static const G4double  lBarCof1B  = 0.740659 << 
293   static const G4double  lBarCofSB  = 0.620659 << 
294   static const G4double  lBarCof2SB = 0.500659 << 
295                                                << 
296   if( pdg == 3122 || pdg == 3222 ||  pdg == 31 << 
297   {                                            << 
298     coeff = lBarCof1S; // Lambda, Sigma+, Sigm << 
299                                                << 
300   } else if( pdg == 3322 || pdg == 3312   )    << 
301   {                                            << 
302     coeff = lBarCof2S; // Xi-, Xi0             << 
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 }                                              << 
386                                                << 
387                                                   222 
388                                                   223