Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4HadronNucleonXsc.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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // 14.03.07 V. Grichine - first implementation
 27 //
 28 // 04.09.18 V. Ivantchenko Major revision of interfaces and implementation
 29 // 30.09.18 V. Grichine hyperon-nucleon xsc first implementation
 30 
 31 #include "G4HadronNucleonXsc.hh"
 32 #include "G4Element.hh"
 33 #include "G4Proton.hh"
 34 #include "G4Nucleus.hh"
 35 #include "G4PhysicalConstants.hh"
 36 #include "G4SystemOfUnits.hh"
 37 #include "G4ParticleTable.hh"
 38 #include "G4IonTable.hh"
 39 #include "G4HadTmpUtil.hh"
 40 #include "G4Log.hh"
 41 #include "G4Exp.hh"
 42 #include "G4Pow.hh"
 43 #include "G4NuclearRadii.hh"
 44 
 45 #include "G4Neutron.hh"
 46 #include "G4PionPlus.hh"
 47 #include "G4KaonPlus.hh"
 48 #include "G4KaonMinus.hh"
 49 #include "G4KaonZeroShort.hh"
 50 #include "G4KaonZeroLong.hh"
 51 
 52 namespace
 53 {
 54   const G4double invGeV  = 1.0/CLHEP::GeV;
 55   const G4double invGeV2 = 1.0/(CLHEP::GeV*CLHEP::GeV);
 56   // PDG fit constants
 57   const G4double minLogP = 3.5;    // min of (lnP-minLogP)^2 
 58   const G4double cofLogE = .0557;  // elastic (lnP-minLogP)^2 
 59   const G4double cofLogT = .3;     // total (lnP-minLogP)^2 
 60   const G4double pMin = .1;        // fast LE calculation 
 61   const G4double pMax = 1000.;     // fast HE calculation 
 62   const G4double ekinmin = 0.1*CLHEP::MeV;   // protection against zero ekin 
 63   const G4double ekinmaxQB = 100*CLHEP::MeV; // max kinetic energy for Coulomb barrier  
 64 }
 65 
 66 G4HadronNucleonXsc::G4HadronNucleonXsc() 
 67 {
 68   // basic hadrons
 69   theProton   = G4Proton::Proton();
 70   theNeutron  = G4Neutron::Neutron();
 71   thePiPlus   = G4PionPlus::PionPlus();
 72 
 73   // basic strange mesons
 74   theKPlus    = G4KaonPlus::KaonPlus();
 75   theKMinus   = G4KaonMinus::KaonMinus();
 76   theK0S      = G4KaonZeroShort::KaonZeroShort();
 77   theK0L      = G4KaonZeroLong::KaonZeroLong();
 78 
 79   g4calc = G4Pow::GetInstance();
 80 }
 81 
 82 void G4HadronNucleonXsc::CrossSectionDescription(std::ostream& outFile) const
 83 {
 84   outFile << "G4HadronNucleonXsc calculates the total, inelastic and elastic\n"
 85           << "hadron-nucleon cross sections using several different\n"
 86           << "parameterizations within the Glauber-Gribov approach. It is\n"
 87           << "valid for all incident gammas and long-lived hadrons at\n"
 88           << "energies above 30 keV.  This is a cross section component which\n"
 89           << "is to be used to build a cross section data set.\n"; 
 90 }
 91 
 92 G4double G4HadronNucleonXsc::HadronNucleonXsc(  const G4ParticleDefinition* theParticle,
 93                                const G4ParticleDefinition* nucleon, G4double ekin)
 94 {
 95   G4double xsc(0.);
 96   G4int pdg = std::abs(theParticle->GetPDGEncoding());
 97   
 98   // p, n, pi+-, pbar, nbar
 99   if ( pdg == 2212 || pdg == 2112 || pdg == 211 ) { 
100     xsc = HadronNucleonXscNS(theParticle, nucleon, ekin);
101   }
102   else if ( pdg == 22 ) // gamma
103   {
104     xsc = HadronNucleonXscPDG(theParticle, nucleon, ekin);
105   }
106   else if ( pdg == 321 || pdg == 310 || pdg == 130 ) // K+-, K0L, K0S
107   {
108     xsc = KaonNucleonXscNS(theParticle, nucleon, ekin);
109   }
110   else if (pdg > 3000) 
111   {
112     if (pdg == 3122 || pdg == 3222 || pdg == 3112 || pdg == 3212 || pdg == 3322 || pdg == 3312 || pdg == 3324 ||
113   pdg == 4122 || pdg == 4332 || pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 4112 || pdg == 4232 || pdg == 4132 ||
114         pdg == 5122 || pdg == 5332 || pdg == 5122 || pdg == 5112 || pdg == 5222 || pdg == 5212 || pdg == 5132 || pdg == 5232
115        ) // heavy s-,c-,b-hyperons
116     {
117       xsc = HyperonNucleonXscNS(theParticle, nucleon, ekin);
118     }
119     else 
120     {
121       xsc = HadronNucleonXscPDG(theParticle, nucleon, ekin);
122     }
123   } else if (pdg > 220) {
124     if (pdg == 511 || pdg == 421 || pdg == 531 || pdg == 541 || pdg == 431 || pdg == 411 || pdg == 521 ||
125   pdg == 221 || pdg == 331 || pdg == 441 || pdg == 443 || pdg == 543) // s-,c-,b-mesons
126     {
127       xsc = SCBMesonNucleonXscNS(theParticle, nucleon, ekin);
128     }
129     else 
130     {
131       xsc = HadronNucleonXscPDG(theParticle, nucleon, ekin);
132     }
133   } else {
134     xsc = HadronNucleonXscPDG(theParticle, nucleon, ekin);
135   }
136   return xsc;
137 }
138 
139 //////////////////////////////////////////////////////////////////////////////
140 //
141 // Returns hadron-nucleon Xsc according to PDG parametrisation (2017):
142 // http://pdg.lbl.gov/2017/reviews/hadronicrpp.pdf
143 
144 G4double G4HadronNucleonXsc::HadronNucleonXscPDG(
145                              const G4ParticleDefinition* theParticle, 
146            const G4ParticleDefinition* nucleon, G4double ekin)
147 {
148   static const G4double M    = 2.1206; // in Gev
149   static const G4double eta1 = 0.4473;
150   static const G4double eta2 = 0.5486;
151   static const G4double H    = 0.272;
152 
153   G4int pdg = theParticle->GetPDGEncoding();
154 
155   G4double mass1 = (pdg == 22) ? 770. : theParticle->GetPDGMass();
156   G4double mass2 = nucleon->GetPDGMass();
157 
158   G4double sMand = CalcMandelstamS(ekin, mass1, mass2)*invGeV2;
159   G4double x = (mass1 + mass2)*invGeV + M;
160   G4double blog = G4Log(sMand/(x*x));
161 
162   G4double P(0.0), R1(0.0), R2(0.0), del(1.0);
163 
164   G4bool proton  = (nucleon == theProton);
165   G4bool neutron = (nucleon == theNeutron);
166   
167   if(theParticle == theNeutron)
168   {
169     if ( proton )
170     {
171       P  = 34.71;
172       R1 = 12.52;
173       R2 = -6.66;
174     }
175     else
176     {
177       P  = 34.41;
178       R1 = 13.07;
179       R2 = -7.394;
180     }
181   } 
182   else if(theParticle == theProton) 
183   {
184     if ( neutron )
185     {
186       P  = 34.71;
187       R1 = 12.52;
188       R2 = -6.66;
189     }
190     else
191     {
192       P  = 34.41;
193       R1 = 13.07;
194       R2 = -7.394;
195     }
196   } 
197   // pbar
198   else if(pdg == -2212) 
199   {
200     if ( neutron )
201     {
202       P  = 34.71;
203       R1 = 12.52;
204       R2 = 6.66;
205     }
206     else
207     {
208       P  = 34.41;
209       R1 = 13.07;
210       R2 = 7.394;
211     }
212   } 
213   // nbar
214   else if(pdg == -2112) 
215   {
216     if ( proton )
217     {
218       P  = 34.71;
219       R1 = 12.52;
220       R2 = 6.66;
221     }
222     else
223     {
224       P  = 34.41;
225       R1 = 13.07;
226       R2 = 7.394;
227     }
228   } 
229   // pi+
230   else if(pdg == 211) 
231   {
232     P  = 18.75;
233     R1 = 9.56;
234     R2 = -1.767;
235   } 
236   // pi-
237   else if(pdg == -211) 
238   {
239     P  = 18.75;
240     R1 = 9.56;
241     R2 = 1.767;
242   } 
243   else if(theParticle == theKPlus) 
244   {
245     if ( proton )
246     {
247       P  = 16.36;
248       R1 = 4.29;
249       R2 = -3.408;
250     }
251     else
252     {
253       P  = 16.31;
254       R1 = 3.7;
255       R2 = -1.826;
256     }
257   } 
258   else if(theParticle == theKMinus) 
259   {
260     if ( proton )
261     {
262       P  = 16.36;
263       R1 = 4.29;
264       R2 = 3.408;
265     }
266     else
267     {
268       P  = 16.31;
269       R1 = 3.7;
270       R2 = 1.826;
271     }
272   }
273   else if(theParticle == theK0S || theParticle == theK0L) 
274   {
275     P  = 16.36;
276     R1 = 2.5;
277     R2 = 0.;
278   }
279   // sigma-
280   else if(pdg == 3112) 
281   {
282     P  = 34.7;
283     R1 = -46.;
284     R2 = 48.;
285   } 
286   // gamma
287   else if(pdg == 22) // modify later on
288   {
289     del= 0.003063;
290     P  = 34.71*del;
291     R1 = (neutron) ? 0.0231 : 0.0139;
292     R2 = 0.;
293   } 
294   else  // as proton ??? 
295   {
296     if ( neutron )
297     {
298       P  = 34.71;
299       R1 = 12.52;
300       R2 = -6.66;
301     }
302     else
303     {
304       P  = 34.41;
305       R1 = 13.07;
306       R2 = -7.394;
307     }
308   } 
309   fTotalXsc = CLHEP::millibarn* 
310     (del*(H*blog*blog + P) + R1*G4Exp(-eta1*blog) + R2*G4Exp(-eta2*blog));
311   fInelasticXsc = 0.75*fTotalXsc;
312   fElasticXsc   = fTotalXsc - fInelasticXsc;
313 
314   if( proton && theParticle->GetPDGCharge() > 0. && ekin < ekinmaxQB )
315   {
316     G4double cB = CoulombBarrier(theParticle, nucleon, ekin);
317     fTotalXsc   *= cB;
318     fElasticXsc *= cB; 
319     fInelasticXsc *= cB; 
320   }
321   /*
322   G4cout << "## HadronNucleonXscPDG: ekin(MeV)= " << ekin 
323    << " tot= " << fTotalXsc/CLHEP::millibarn 
324    << " inel= " << fInelasticXsc/CLHEP::millibarn
325    << " el= " << fElasticXsc/CLHEP::millibarn
326          << G4endl;
327   */
328   return fTotalXsc;
329 }
330 
331 //////////////////////////////////////////////////////////////////////////////
332 //
333 // Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
334 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
335 
336 G4double G4HadronNucleonXsc::HadronNucleonXscNS(
337                              const G4ParticleDefinition* theParticle, 
338            const G4ParticleDefinition* nucleon, G4double ekin0)
339 {  
340   const G4double ekin = std::max(ekin0, ekinmin);
341   G4int pdg = theParticle->GetPDGEncoding();
342   /*  
343   G4cout<< "HadronNucleonXscNS: Ekin(GeV)= " << ekin/GeV << "  "
344         << theParticle->GetParticleName() << " + " 
345         << nucleon->GetParticleName() << G4endl;
346   */
347   if(pdg == -2212 || pdg == -2112) {
348     return HadronNucleonXscPDG(theParticle, nucleon, ekin);
349   }
350 
351   G4double pM   = theParticle->GetPDGMass();
352   G4double tM   = nucleon->GetPDGMass();
353   G4double pE   = ekin + pM; 
354   G4double pLab = std::sqrt(ekin*(ekin + 2*pM));
355 
356   G4double sMand = CalcMandelstamS(ekin, pM, tM)*invGeV2;
357 
358   pLab *= invGeV;
359   pE   *= invGeV;
360 
361   if(pLab >= 10.) {
362     fTotalXsc = HadronNucleonXscPDG(theParticle, nucleon, ekin)/CLHEP::millibarn;
363   } else { fTotalXsc = 0.0; }
364   fElasticXsc = 0.0;
365   //G4cout << "Stot(mb)= " << fTotalXsc << " pLab= " << pLab 
366   //   << " Smand= " << sMand <<G4endl;
367   G4double logP = G4Log(pLab);
368 
369   G4bool proton = (nucleon == theProton);
370   G4bool neutron = (nucleon == theNeutron);
371 
372   if( theParticle == theNeutron) 
373   {
374     if( pLab >= 373.)
375     {
376       fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4Log(sMand/400.)*1.65)) 
377   + 9.19*G4Exp(-G4Log(sMand)*0.458);
378     }
379     else if( pLab >= 100.)
380     {
381       fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G4Log(sMand/28.9))*1.1) 
382   + 9.19*G4Exp(-G4Log(sMand)*0.458);
383     }
384     else if( pLab >= 10.)
385     {
386       fElasticXsc =  6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
387     }
388     else  // pLab < 10 GeV/c
389     {
390       if( neutron )      // nn to be pp
391       {
392         G4double x = G4Log(pLab/0.73);
393         if( pLab < 0.4 )
394         {
395           fTotalXsc = 23 + 50*std::sqrt(g4calc->powN(-x, 7));
396           fElasticXsc = fTotalXsc;
397         }
398         else if( pLab < 0.73 )
399         {
400           fTotalXsc = 23 + 50*std::sqrt(g4calc->powN(-x, 7));
401           fElasticXsc = fTotalXsc; 
402         }
403         else if( pLab < 1.05  )
404         {
405           fTotalXsc = 23 + 40*x*x;
406           fElasticXsc = 23 + 20*x*x;
407         }
408         else    // 1.05 - 10 GeV/c
409         {
410           fTotalXsc = 39.0+75*(pLab - 1.2)/(g4calc->powN(pLab,3) + 0.15);
411           fElasticXsc =  6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
412         }
413       }
414       if( proton )   // pn to be np
415       {
416         if( pLab < 0.02 )
417         {
418           fTotalXsc = 4100+30*G4Exp(G4Log(G4Log(1.3/pLab))*3.6); 
419     fElasticXsc = fTotalXsc;
420         }      
421         else if( pLab < 0.8 )
422         {
423           fTotalXsc = 33+30*g4calc->powN(G4Log(pLab/1.3),4);
424     fElasticXsc = fTotalXsc;
425         }      
426         else if( pLab < 1.4 )
427         {
428           fTotalXsc = 33+30*g4calc->powN(G4Log(pLab/0.95),2);
429           G4double x = G4Log(0.511/pLab);
430           fElasticXsc =  6 + 52/( x*x + 1.6 );
431         }
432         else    // 1.4 < pLab < 10.  )
433         {
434           fTotalXsc = 33.3 + 20.8*(pLab*pLab - 1.35)
435       /(std::sqrt(g4calc->powN(pLab,5)) + 0.95);
436           fElasticXsc =  6 + 20/( (logP-0.182)*(logP-0.182) + 1.0 );
437         }
438       }
439     }
440   }
441   ////// proton //////////////////////////////////////////////
442   else if(theParticle == theProton) 
443   {
444     if( pLab >= 373.) // pdg due to TOTEM data
445     {
446       fElasticXsc = 6.5 + 0.308*G4Exp(G4Log(G4Log(sMand/400.))*1.65) 
447   + 9.19*G4Exp(-G4Log(sMand)*0.458);     
448     }
449     else if( pLab >= 100.)
450     {
451       fElasticXsc = 5.53 + 0.308*G4Exp(G4Log(G4Log(sMand/28.9))*1.1) 
452   + 9.19*G4Exp(-G4Log(sMand)*0.458);     
453     }
454     else if( pLab >= 10.)
455     {
456       fElasticXsc = 6. + 20./( (logP-0.182)*(logP-0.182) + 1.0 );
457     }
458     else
459     {
460       // pp
461       if( proton )
462       {
463         if( pLab < 0.73 )
464         {
465           fTotalXsc = 23 + 50*std::sqrt(g4calc->powN(G4Log(0.73/pLab),7));
466           fElasticXsc = fTotalXsc;
467         }
468         else if( pLab < 1.05  )
469         {
470     G4double x = G4Log(pLab/0.73);
471           fTotalXsc = 23 + 40*x*x;
472           fElasticXsc = 23 + 20*x*x;
473         }
474         else    // 1.05 - 10 GeV/c
475         {
476           fTotalXsc = 39.0+75*(pLab - 1.2)/(g4calc->powN(pLab,3) + 0.15);
477           fElasticXsc = 6. + 20./( (logP-0.182)*(logP-0.182) + 1.0 );
478         }
479       }
480       else if( neutron )     // pn to be np
481       {
482         if( pLab < 0.02 )
483         {
484           fTotalXsc = 4100+30*G4Exp(G4Log(G4Log(1.3/pLab))*3.6); 
485     fElasticXsc = fTotalXsc;
486         }      
487         else if( pLab < 0.8 )
488         {
489           fTotalXsc = 33+30*g4calc->powN(G4Log(pLab/1.3),4);
490     fElasticXsc = fTotalXsc;
491         }      
492         else if( pLab < 1.4 )
493         {
494           G4double x1 = G4Log(pLab/0.95);
495           G4double x2 = G4Log(0.511/pLab);
496           fTotalXsc = 33+30*x1*x1;
497           fElasticXsc =  6 + 52/(x2*x2 + 1.6);
498         }
499         else    // 1.4 < pLab < 10.  )
500         {
501           fTotalXsc = 33.3 + 20.8*(pLab*pLab - 1.35)
502       /(std::sqrt(g4calc->powN(pLab,5)) + 0.95);
503           fElasticXsc = 6. + 20./( (logP-0.182)*(logP-0.182) + 1.0 );
504         }
505       }
506     }    
507   }
508   // pi+ p; pi- n 
509   else if((pdg == 211 && proton) || (pdg == -211 && neutron))  
510   {
511     if( pLab < 0.28 )
512     {
513       fTotalXsc = 10./((logP + 1.273)*(logP + 1.273) + 0.05);
514       fElasticXsc = fTotalXsc;
515     }
516     else if( pLab < 0.68 )
517     {
518       fTotalXsc = 14./((logP + 1.273)*(logP + 1.273) + 0.07);
519       fElasticXsc = fTotalXsc;
520     }
521     else if( pLab < 0.85 )
522     {
523       G4double x = G4Log(pLab/0.77);
524       fTotalXsc = 88.*x*x + 14.9;
525       fElasticXsc = fTotalXsc*G4Exp(-3.*(pLab - 0.68));  
526     }
527     else if( pLab < 1.15 )
528     {
529       G4double x = G4Log(pLab/0.77);
530       fTotalXsc = 88.*x*x + 14.9;
531       fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
532     }
533     else if( pLab < 1.4) // ns original
534     {
535       G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
536       G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
537       fTotalXsc       = Ex1 + Ex2 + 27.5;
538       fElasticXsc = 6.0 + 1.4/(( pLab - 1.4)*( pLab - 1.4) + 0.1);
539     }
540     else if( pLab < 2.0 ) // ns original
541     {
542       G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
543       G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
544       fTotalXsc     = Ex1 + Ex2 + 27.5;
545       fElasticXsc = 3.0 + 1.36/( (logP - 0.336)*(logP - 0.336) + 0.08);    
546     }
547     else if( pLab < 3.5 ) // ns original
548     {
549       G4double Ex1 = 3.2*G4Exp(-(pLab-2.55)*(pLab-2.55)/0.55/0.55);
550       G4double Ex2 = 12*G4Exp(-(pLab-1.47)*(pLab-1.47)/0.225/0.225);
551       fTotalXsc       = Ex1 + Ex2 + 27.5;
552       fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);    
553     }
554     else if( pLab < 10. ) // my
555     {
556       fTotalXsc = 10.6 + 2.*G4Log(pE) + 25*G4Exp(-G4Log(pE)*0.43 ); 
557       fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);    
558     }
559     else //  pLab > 10 // my
560     {
561       fElasticXsc = 3.0 + 6.20/( (logP - 0.336)*(logP - 0.336) + 0.8);    
562     }
563   }
564   // pi+ n; pi- p
565   else if((pdg == 211 && neutron) || (pdg == -211 && proton))  
566   {
567     if( pLab < 0.28 ) 
568     {
569       fTotalXsc = 0.288/((pLab - 0.28)*(pLab - 0.28) + 0.004);
570       fElasticXsc = 1.8/((logP + 1.273)*(logP + 1.273) + 0.07);
571     }
572     else if( pLab < 0.395676 ) // first peak
573     {
574       fTotalXsc = 0.648/((pLab - 0.28)*(pLab - 0.28) + 0.009);
575       fElasticXsc = 0.257/((pLab - 0.28)*(pLab - 0.28) + 0.01);
576     }
577     else if( pLab < 0.5 )
578     {
579       G4double y = G4Log(pLab/0.48);
580       fTotalXsc = 26 + 110*y*y;
581       fElasticXsc = 0.37*fTotalXsc;
582     }
583     else if( pLab < 0.65 )
584     {
585       G4double x = G4Log(pLab/0.48);
586       fTotalXsc = 26. + 110.*x*x;
587       fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
588     }
589     else if( pLab < 0.72 )
590     {
591       fTotalXsc = 36.1 + 10*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
592   24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
593       fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
594     }
595     else if( pLab < 0.88 )
596     {
597       fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
598   24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
599       fElasticXsc = 0.95/((pLab - 0.72)*(pLab - 0.72) + 0.049);
600     }
601     else if( pLab < 1.03 )
602     {
603       fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
604   24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
605       fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
606     }
607     else if( pLab < 1.15 )
608     {
609       fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
610   24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
611       fElasticXsc = 2.0 + 0.4/((pLab - 1.03)*(pLab - 1.03) + 0.016);
612     }
613     else if( pLab < 1.3 )
614     {
615       fTotalXsc = 36.1 + 10.*G4Exp(-(pLab-0.72)*(pLab-0.72)/0.06/0.06)+
616   24*G4Exp(-(pLab-1.015)*(pLab-1.015)/0.075/0.075);
617       fElasticXsc = 3. + 13./pLab;
618     }
619     else if( pLab < 10.) // < 3.0) // ns original
620     {
621       fTotalXsc = 36.1 + 0.079-4.313*logP+
622   3*G4Exp(-(pLab-2.1)*(pLab-2.1)/0.4/0.4)+
623   1.5*G4Exp(-(pLab-1.4)*(pLab-1.4)/0.12/0.12);
624       fElasticXsc = 3. + 13./pLab; 
625     }
626     else   // mb 
627     {
628       fElasticXsc = 3. + 13./pLab;
629     }
630   }
631   else if( (theParticle == theKMinus) && proton )   // K-p
632   {
633     if( pLab < pMin)
634     {
635       G4double psp = pLab*std::sqrt(pLab);
636       fElasticXsc  = 5.2/psp;
637       fTotalXsc    = 14./psp;
638     }
639     else if( pLab > pMax )
640     {
641       G4double ld  = logP - minLogP;
642       G4double ld2 = ld*ld;
643       fElasticXsc  = cofLogE*ld2 + 2.23;
644       fTotalXsc    = 1.1*cofLogT*ld2 + 19.7;
645     }
646     else
647     {
648       G4double ld  = logP - minLogP;
649       G4double ld2 = ld*ld;
650       G4double sp  = std::sqrt(pLab);
651       G4double psp = pLab*sp;
652       G4double p2  = pLab*pLab;
653       G4double p4  = p2*p2;
654 
655       G4double lh  = pLab - 1.01;
656       G4double hd  = lh*lh + .011;
657 
658       G4double lm  = pLab - .39;
659       G4double md  = lm*lm + .000356;
660       G4double lh1  = pLab - 0.78;
661       G4double hd1  = lh1*lh1 + .00166;
662       G4double lh2  = pLab - 1.63;
663       G4double hd2  = lh2*lh2 + .007;
664 
665       // small peaks were added but commented out now
666       fElasticXsc  = 5.2/psp + (1.1*cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) 
667   + .004/md + 0.005/hd1+ 0.01/hd2 +.15/hd; 
668 
669       fTotalXsc   = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4) 
670         + .006/md  + 0.01/hd1+ 0.02/hd2 + .20/hd ;
671     }
672   }
673   else if( (theParticle == theKMinus) && neutron )  // K-n
674   {
675     if( pLab > pMax )
676     {
677       G4double ld  = logP - minLogP;
678       G4double ld2 = ld*ld;
679       fElasticXsc  = cofLogE*ld2 + 2.23;
680       fTotalXsc    = 1.1*cofLogT*ld2 + 19.7;
681     }
682     else
683     { 
684       G4double lh  = pLab - 0.98;
685       G4double hd  = lh*lh + .021; 
686       G4double sqrLogPlab = logP*logP;
687 
688       fElasticXsc = 5.0 +  8.1*G4Exp(-logP*1.8) + 0.16*sqrLogPlab 
689   - 1.3*logP + .15/hd;
690       fTotalXsc = 25.2 +  0.38*sqrLogPlab - 2.9*logP + 0.30/hd;
691     }
692   }
693   else if( (theParticle == theKPlus) && proton  )  // K+p
694   {
695     // VI: modified low-energy part
696     if( pLab < 0.631 )
697     {
698       fElasticXsc = fTotalXsc = 12.03;
699     }
700     else if( pLab > pMax )
701     {
702       G4double ld  = logP - minLogP;
703       G4double ld2 = ld*ld;
704       fElasticXsc  = cofLogE*ld2 + 2.23;
705       fTotalXsc    = cofLogT*ld2 + 19.2;
706     }
707     else
708     {
709       G4double ld  = logP - minLogP;
710       G4double ld2 = ld*ld;
711       G4double lr  = pLab - .38;
712       G4double LE  = .7/(lr*lr + .076);
713       G4double sp  = std::sqrt(pLab);
714       G4double p2  = pLab*pLab;
715       G4double p4  = p2*p2;
716       // VI: tuned elastic
717       fElasticXsc  = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) 
718   + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652);
719       fTotalXsc    = LE + (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) 
720   + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392);
721     }
722   }
723   else if(  (theParticle == theKPlus) && neutron) // K+n  
724   {
725     if( pLab < pMin )
726     {
727       G4double lm = pLab - 0.94;
728       G4double md = lm*lm + .392;   
729       fElasticXsc = 2./md;
730       fTotalXsc   = 4.6/md;
731     }
732     else if( pLab > pMax )
733     {
734       G4double ld  = logP - minLogP;
735       G4double ld2 = ld*ld;
736       fElasticXsc  = cofLogE*ld2 + 2.23;
737       fTotalXsc    = cofLogT*ld2 + 19.2;
738     }
739     else
740     {
741       G4double ld  = logP - minLogP;
742       G4double ld2 = ld*ld;
743       G4double sp  = std::sqrt(pLab);
744       G4double p2  = pLab*pLab;
745       G4double p4  = p2*p2;
746       G4double lm  = pLab - 0.94;
747       G4double md  = lm*lm + .392; 
748       fElasticXsc  = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
749       fTotalXsc    = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 4.6/md;
750     }
751   }
752   fTotalXsc   *= CLHEP::millibarn; 
753   fElasticXsc *= CLHEP::millibarn;
754   fElasticXsc  = std::min(fElasticXsc, fTotalXsc);
755 
756   if( proton && theParticle->GetPDGCharge() > 0. && ekin < ekinmaxQB )
757   {
758     G4double cB = G4NuclearRadii::CoulombFactor(theParticle, nucleon, ekin);
759     fTotalXsc   *= cB;
760     fElasticXsc *= cB; 
761   }
762   fInelasticXsc = std::max(fTotalXsc - fElasticXsc,0.0);
763   /*
764   G4cout<< "HNXsc: Ekin(GeV)= " << ekin/GeV << "; tot(mb)= " << fTotalXsc/millibarn
765         <<"; el(mb)= " <<fElasticXsc/millibarn
766         <<"; inel(mb)= " <<fInelasticXsc/millibarn<<"  "
767         << theParticle->GetParticleName() << " + " 
768         << nucleon->GetParticleName() << G4endl;
769   */
770   return fTotalXsc;
771 }
772 
773 //////////////////////////////////////////////////////////////////////////////
774 //
775 // Returns kaon-nucleon cross-section based on smoothed NS 
776 // tuned for the Glauber-Gribov hadron model for Z>1
777 
778 G4double G4HadronNucleonXsc::KaonNucleonXscGG(
779                              const G4ParticleDefinition* theParticle, 
780            const G4ParticleDefinition* nucleon, G4double ekin)
781 {
782   fTotalXsc = fElasticXsc = fInelasticXsc = 0.0;
783   if(theParticle == theKMinus || theParticle == theKPlus) {
784     KaonNucleonXscVG(theParticle, nucleon, ekin);
785 
786   } else if(theParticle == theK0S || theParticle == theK0L) {
787     G4double stot  = KaonNucleonXscVG(theKMinus, nucleon, ekin);
788     G4double sel   = fElasticXsc;
789     G4double sinel = fInelasticXsc;
790     stot  += KaonNucleonXscVG(theKPlus, nucleon, ekin);
791     sel   += fElasticXsc;
792     sinel += fInelasticXsc;
793     fTotalXsc = stot*0.5; 
794     fElasticXsc = sel*0.5; 
795     fInelasticXsc = sinel*0.5;
796   }
797   return fTotalXsc;
798 }
799 
800 //////////////////////////////////////////////////////////////////////////////
801 //
802 // Returns kaon-nucleon cross-section using NS
803 
804 G4double G4HadronNucleonXsc::KaonNucleonXscNS(
805                              const G4ParticleDefinition* theParticle, 
806            const G4ParticleDefinition* nucleon, G4double ekin)
807 {
808   fTotalXsc = fElasticXsc = fInelasticXsc = 0.0;
809   if(theParticle == theKMinus || theParticle == theKPlus) {
810     HadronNucleonXscNS(theParticle, nucleon, ekin);
811 
812   } else if(theParticle == theK0S || theParticle == theK0L) {
813     G4double fact = 0.5;
814     G4double stot = 0.0;
815     G4double sel  = 0.0;
816     G4double sinel= 0.0;
817     if(ekin > ekinmaxQB) {
818       stot  = HadronNucleonXscNS(theKMinus, nucleon, ekin);
819       sel   = fElasticXsc;
820       sinel = fInelasticXsc;
821       stot  += HadronNucleonXscNS(theKPlus, nucleon, ekin);
822       sel   += fElasticXsc;
823       sinel += fInelasticXsc;
824     } else {
825       fact *= std::sqrt(ekinmaxQB/std::max(ekin, ekinmin));
826       stot  = HadronNucleonXscNS(theKMinus, nucleon, ekinmaxQB);
827       sel   = fElasticXsc;
828       sinel = fInelasticXsc;
829       stot  += HadronNucleonXscNS(theKPlus, nucleon, ekinmaxQB);
830       sel   += fElasticXsc;
831       sinel += fInelasticXsc;
832     }
833     fTotalXsc = stot*fact; 
834     fElasticXsc = sel*fact; 
835     fInelasticXsc = sinel*fact;
836   }
837   return fTotalXsc;
838 }
839 
840 //////////////////////////////////////////////////////////////////////////////
841 //
842 // Returns kaon-nucleon cross-section with smoothed NS parameterisation
843 
844 G4double G4HadronNucleonXsc::KaonNucleonXscVG(
845                  const G4ParticleDefinition* theParticle, 
846      const G4ParticleDefinition* nucleon, G4double ekin)
847 {
848   G4double pM   = theParticle->GetPDGMass();
849   G4double pLab = std::sqrt(ekin*(ekin + 2*pM));
850 
851   pLab *= invGeV;
852   G4double logP = G4Log(pLab);
853 
854   fTotalXsc = 0.0;
855 
856   G4bool proton = (nucleon == theProton);
857   G4bool neutron = (nucleon == theNeutron);
858 
859   if( (theParticle == theKMinus) && proton )   // K-p
860   {
861     if( pLab < pMin)
862     {
863       G4double psp = pLab*std::sqrt(pLab);
864       fElasticXsc  = 5.2/psp;
865       fTotalXsc    = 14./psp;
866     }
867     else if( pLab > pMax )
868     {
869       G4double ld  = logP - minLogP;
870       G4double ld2 = ld*ld;
871       fElasticXsc  = cofLogE*ld2 + 2.23;
872       fTotalXsc    = 1.1*cofLogT*ld2 + 19.7;
873     }
874     else
875     {
876       G4double ld  = logP - minLogP;
877       G4double ld2 = ld*ld;
878       G4double sp  = std::sqrt(pLab);
879       G4double psp = pLab*sp;
880       G4double p2  = pLab*pLab;
881       G4double p4  = p2*p2;
882 
883       G4double lh  = pLab - 1.01;
884       G4double hd  = lh*lh + .011;
885       fElasticXsc  = 5.2/psp + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .15/hd;
886       fTotalXsc    = 14./psp + (1.1*cofLogT*ld2 + 19.5)/(1. - .21/sp + .52/p4) + .60/hd;
887     }
888   }
889   else if( (theParticle == theKMinus) && neutron )  // K-n
890   {
891     if( pLab > pMax )
892     {
893       G4double ld  = logP - minLogP;
894       G4double ld2 = ld*ld;
895       fElasticXsc  = cofLogE*ld2 + 2.23;
896       fTotalXsc    = 1.1*cofLogT*ld2 + 19.7;
897     }
898     else
899     { 
900       G4double lh  = pLab - 0.98;
901       G4double hd  = lh*lh + .045;    // vg version
902       G4double sqrLogPlab = logP*logP;
903 
904       fElasticXsc = 5.0 +  8.1*G4Exp(-logP*1.8) + 0.16*sqrLogPlab 
905   - 1.3*logP + .15/hd;
906       fTotalXsc = 25.2 +  0.38*sqrLogPlab - 2.9*logP + 0.60/hd;  // vg version
907     }
908   }
909   else if( (theParticle == theKPlus) && proton  )  // K+p
910   {
911     // VI: modified low-energy part
912     if( pLab < 0.631 )
913     {
914       fElasticXsc = fTotalXsc = 12.03;
915     }
916     else if( pLab > pMax )
917     {
918       G4double ld  = logP - minLogP;
919       G4double ld2 = ld*ld;
920       fElasticXsc  = cofLogE*ld2 + 2.23;
921       fTotalXsc    = cofLogT*ld2 + 19.2;
922     }
923     else
924     {
925       G4double ld  = logP - minLogP;
926       G4double ld2 = ld*ld;
927       G4double lr  = pLab - .38;
928       G4double LE  = .7/(lr*lr + .076);
929       G4double sp  = std::sqrt(pLab);
930       G4double p2  = pLab*pLab;
931       G4double p4  = p2*p2;
932       // VI: tuned elastic
933       fElasticXsc  = LE + (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) 
934   + 2./((pLab - 0.8)*(pLab - 0.8) + 0.652);
935       fTotalXsc    = LE + (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) 
936   + 2.6/((pLab - 1.)*(pLab - 1.) + 0.392);
937     }
938   }
939   else if(  (theParticle == theKPlus) && neutron) // K+n  
940   {
941     if( pLab < pMin )
942     {
943       G4double lm = pLab - 0.94;
944       G4double md = lm*lm + .392;   
945       fElasticXsc = 2./md;
946       fTotalXsc   = 4.6/md;
947     }
948     else if( pLab > pMax )
949     {
950       G4double ld  = logP - minLogP;
951       G4double ld2 = ld*ld;
952       fElasticXsc  = cofLogE*ld2 + 2.23;
953       fTotalXsc    = cofLogT*ld2 + 19.2;
954     }
955     else
956     {
957       G4double ld  = logP - minLogP;
958       G4double ld2 = ld*ld;
959       G4double sp  = std::sqrt(pLab);
960       G4double p2  = pLab*pLab;
961       G4double p4  = p2*p2;
962       G4double lm  = pLab - 0.94;
963       G4double md  = lm*lm + .392; 
964       fElasticXsc  = (cofLogE*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
965       fTotalXsc    = (cofLogT*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 4.6/md;
966     }
967   }
968 
969   fTotalXsc   *= CLHEP::millibarn; 
970   fElasticXsc *= CLHEP::millibarn;
971 
972   if( proton && theParticle->GetPDGCharge() > 0. )
973   {
974     G4double cB = G4NuclearRadii::CoulombFactor(theParticle, nucleon, ekin);
975     fTotalXsc   *= cB;
976     fElasticXsc *= cB; 
977   }
978   fElasticXsc = std::min(fElasticXsc, fTotalXsc);
979   fInelasticXsc = std::max(fTotalXsc - fElasticXsc,0.0);
980   /*
981   G4cout << "HNXscVG: E= " << ekin << " " << theParticle->GetParticleName() 
982    << " P: " << proton << " xtot(b)= " << fTotalXsc/barn
983    << " xel(b)= " <<  fElasticXsc/barn << " xinel(b)= " << fInelasticXsc/barn
984    << G4endl;
985   */
986   return fTotalXsc;
987 }
988 
989 //////////////////////////////////////////////////////////////////////////////
990 //
991 // Returns hyperon-nucleon cross-section using NS x-section for protons
992 
993 G4double G4HadronNucleonXsc::HyperonNucleonXscNS(
994          const G4ParticleDefinition* theParticle, 
995    const G4ParticleDefinition* nucleon, G4double ekin)
996 {
997   G4double coeff = 1.0;
998   G4int pdg = std::abs(theParticle->GetPDGEncoding());
999     
1000   // lambda, sigma+-0 and anti-hyperons
1001   if( pdg == 3122 || pdg == 3112 || pdg == 3212 || pdg == 3222 )
1002   {
1003     coeff = 0.88;
1004   } 
1005   // Xi hyperons and anti-hyperons
1006   else if( pdg == 3312 || pdg == 3322 )
1007   {
1008     coeff = 0.76;
1009   }
1010   // omega, anti_omega
1011   else if( pdg == 3334 )
1012   {
1013     coeff = 0.64;
1014   }
1015   // lambdaC, sigmaC+-0 and anti-hyperonsC
1016   else if( pdg == 4122 || pdg == 4112 || pdg == 4212 || pdg == 4222 ) 
1017   {
1018     coeff = 0.784378;
1019   }
1020   // omegaC0, anti_omegaC0
1021   else if( pdg == 4332 )
1022   {
1023     coeff = 0.544378;
1024   }
1025   // XiC+0 and anti-hyperonC
1026   else if( pdg == 4132 || pdg == 4232 )
1027   {
1028     coeff = 0.664378;
1029   }
1030   // lambdaB, sigmaB+-0 and anti-hyperonsB
1031   else if( pdg == 5122 || pdg == 5112 || pdg == 5212 || pdg == 5222 )
1032   {
1033     coeff = 0.740659;
1034   }
1035   // omegaB0, anti_omegaB0
1036   else if( pdg == 5332 )
1037   {
1038     coeff = 0.500659;
1039   }
1040   // XiB+0 and anti-hyperonB
1041   else if( pdg == 5132 || pdg == 5232 )
1042   {
1043     coeff = 0.620659;
1044   } 
1045   fTotalXsc = coeff*HadronNucleonXscNS( theProton, nucleon, ekin);
1046   fInelasticXsc *= coeff;
1047   fElasticXsc *= coeff;
1048   
1049   return fTotalXsc;
1050 }
1051 
1052 //////////////////////////////////////////////////////////////////////////////
1053 //
1054 // Returns hyperon-nucleon cross-section using NS x-section for protons
1055 
1056 G4double G4HadronNucleonXsc::SCBMesonNucleonXscNS( 
1057          const G4ParticleDefinition* theParticle, 
1058          const G4ParticleDefinition* nucleon, G4double ekin )
1059 {
1060   G4double coeff(1.0);
1061   G4int pdg = std::abs(theParticle->GetPDGEncoding());
1062 
1063   // B+-0 anti
1064   if( pdg == 511 || pdg == 521 )
1065   {
1066     coeff = 0.610989;
1067   }
1068   // D+-0 anti
1069   else if( pdg == 411 || pdg == 421 )
1070   {
1071     coeff = 0.676568;
1072   }
1073   // Bs, antiBs
1074   else if( pdg == 531 )
1075   {
1076     coeff = 0.430989;
1077   }
1078   // Bc+-
1079   else if( pdg == 541 )
1080   {
1081     coeff = 0.287557;
1082   }
1083   // Ds+-
1084   else if( pdg == 431 )
1085   {
1086     coeff = 0.496568;
1087   }
1088   // etaC, J/Psi
1089   else if( pdg == 441 || pdg == 443 )
1090   {
1091     coeff = 0.353135;
1092   }
1093   // Upsilon
1094   else if( pdg == 553 )
1095   {
1096     coeff = 0.221978;
1097   }
1098   // eta
1099   else if( pdg == 221 )
1100   {
1101     coeff = 0.76;
1102   }
1103   // eta'
1104   else if( pdg == 331 )
1105   {
1106     coeff = 0.88;
1107   }
1108   fTotalXsc = coeff*HadronNucleonXscNS(thePiPlus, nucleon, ekin);
1109   fElasticXsc *= coeff;
1110   fInelasticXsc *= coeff;
1111   return fTotalXsc;
1112 }
1113 ////////////////////////////////////////////////////////////////////////////////
1114 //
1115 // Returns hadron-nucleon cross-section based on V. Uzjinsky parametrisation of
1116 // data from G4FTFCrossSection class
1117 
1118 G4double G4HadronNucleonXsc::HadronNucleonXscVU(
1119                              const G4ParticleDefinition* theParticle, 
1120            const G4ParticleDefinition* nucleon, G4double ekin)
1121 {
1122   G4int PDGcode = theParticle->GetPDGEncoding();
1123   G4int absPDGcode = std::abs(PDGcode);
1124   G4double mass = theParticle->GetPDGMass();
1125   G4double Plab = std::sqrt(ekin*(ekin + 2.*mass))*invGeV;
1126 
1127   G4double logPlab = G4Log( Plab );
1128   G4double sqrLogPlab = logPlab * logPlab;
1129 
1130   G4bool proton = (nucleon == theProton);
1131   G4bool neutron = (nucleon == theNeutron);
1132    
1133   if( absPDGcode > 1000)  //------Projectile is baryon -
1134   {
1135     if(proton)
1136     {
1137       fTotalXsc   = 48.0 + 0.522*sqrLogPlab - 4.51*logPlab;
1138       fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1139     }
1140     if(neutron)
1141     {    
1142       fTotalXsc   = 47.3  + 0.513*sqrLogPlab - 4.27*logPlab;
1143       fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1144     }
1145   }
1146   else if( PDGcode ==  211)  //------Projectile is PionPlus ----
1147   {
1148     if(proton)
1149     {
1150       fTotalXsc  = 16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab;
1151       fElasticXsc =  0.0 + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab;
1152     }
1153     if(neutron)
1154     {    
1155       fTotalXsc   =  33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab;
1156       fElasticXsc = 1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab;
1157     }
1158   }
1159   else if( PDGcode == -211)  //------Projectile is PionMinus ----
1160   {
1161     if(proton)
1162     {
1163       fTotalXsc   = 33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab;
1164       fElasticXsc = 1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab;
1165     }
1166     if(neutron)
1167     {    
1168       fTotalXsc   = 16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab;
1169       fElasticXsc =  0.0 + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab;
1170     }
1171   }
1172   else if( PDGcode ==  111)  //------Projectile is PionZero  --
1173   {
1174     if(proton)
1175     {
1176       fTotalXsc = (16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab +   //Pi+
1177        33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab)/2; //Pi-
1178 
1179       fElasticXsc = (0.0 + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab +   //Pi+
1180          1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab)/2;//Pi-
1181 
1182     }
1183     if(neutron)
1184     {    
1185       fTotalXsc = (33.0 + 14.0 *G4Exp(-logPlab*1.36) + 0.456*sqrLogPlab - 4.03*logPlab +   //Pi+
1186        16.4 + 19.3 *G4Exp(-logPlab*0.42) + 0.19 *sqrLogPlab - 0.0 *logPlab)/2; //Pi-
1187       fElasticXsc = (1.76 + 11.2*G4Exp(-logPlab*0.64) + 0.043*sqrLogPlab - 0.0 *logPlab +  //Pi+
1188          0.0  + 11.4*G4Exp(-logPlab*0.40) + 0.079*sqrLogPlab - 0.0 *logPlab)/2;//Pi-
1189     }
1190   }
1191   else if( PDGcode == 321 )    //------Projectile is KaonPlus --
1192   {
1193     if(proton)
1194     {
1195       fTotalXsc   = 18.1 +  0.26 *sqrLogPlab - 1.0 *logPlab;
1196       fElasticXsc =  5.0 +  8.1*G4Exp(-logPlab*1.8 ) + 0.16 *sqrLogPlab - 1.3 *logPlab;
1197     }
1198     if(neutron)
1199     {    
1200       fTotalXsc   = 18.7  + 0.21 *sqrLogPlab - 0.89*logPlab;
1201       fElasticXsc =  7.3  + 0.29 *sqrLogPlab - 2.4 *logPlab;
1202     }
1203   }
1204   else if( PDGcode ==-321 )  //------Projectile is KaonMinus ----
1205   {
1206     if(proton)
1207     {
1208       fTotalXsc   = 32.1 + 0.66*sqrLogPlab - 5.6*logPlab;
1209       fElasticXsc =  7.3 + 0.29*sqrLogPlab - 2.4*logPlab;
1210     }
1211     if(neutron)
1212     {    
1213       fTotalXsc   = 25.2 + 0.38*sqrLogPlab - 2.9*logPlab;
1214       fElasticXsc =  5.0 +  8.1*G4Exp(-logPlab*1.8 ) + 0.16*sqrLogPlab - 1.3*logPlab;
1215     }
1216   }
1217   else if( PDGcode == 311 )  //------Projectile is KaonZero -----
1218   {
1219     if(proton)
1220     {
1221       fTotalXsc   = ( 18.1 + 0.26 *sqrLogPlab - 1.0 *logPlab +   //K+
1222                       32.1 + 0.66 *sqrLogPlab - 5.6 *logPlab)/2; //K-
1223       fElasticXsc = (  5.0 +  8.1*G4Exp(-logPlab*1.8 ) + 0.16 *sqrLogPlab - 1.3 *logPlab + //K+
1224                          7.3 + 0.29 *sqrLogPlab - 2.4 *logPlab)/2; //K-
1225     }
1226     if(neutron)
1227     {    
1228       fTotalXsc   = ( 18.7 + 0.21 *sqrLogPlab - 0.89*logPlab +   //K+
1229                       25.2 + 0.38 *sqrLogPlab - 2.9 *logPlab)/2; //K-
1230       fElasticXsc = (  7.3 + 0.29 *sqrLogPlab - 2.4 *logPlab +   //K+
1231                        5.0 + 8.1*G4Exp(-logPlab*1.8 ) + 0.16 *sqrLogPlab - 1.3 *logPlab)/2; //K-
1232     }
1233   }
1234   else  //------Projectile is undefined, Nucleon assumed
1235   {
1236     if(proton)
1237     {
1238       fTotalXsc   = 48.0 + 0.522*sqrLogPlab - 4.51*logPlab;
1239       fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1240     }
1241     if(neutron)
1242     {    
1243       fTotalXsc   = 47.3 + 0.513*sqrLogPlab - 4.27*logPlab;
1244       fElasticXsc = 11.9 + 26.9*G4Exp(-logPlab*1.21) + 0.169*sqrLogPlab - 1.85*logPlab;
1245     }
1246   }
1247 
1248   fTotalXsc   *= CLHEP::millibarn;
1249   fElasticXsc *= CLHEP::millibarn;
1250   fElasticXsc = std::min(fElasticXsc, fTotalXsc);
1251   fInelasticXsc = fTotalXsc - fElasticXsc;
1252 
1253   return fTotalXsc;    
1254 }
1255 
1256 //////////////////////////////////////////////////////////////////////////////
1257 //
1258 // Returns hadron-nucleon Xsc according to different parametrisations:
1259 // [2] E. Levin, hep-ph/9710546
1260 // [3] U. Dersch, et al, hep-ex/9910052
1261 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725 
1262 
1263 G4double G4HadronNucleonXsc::HadronNucleonXscEL(
1264                              const G4ParticleDefinition* theParticle, 
1265            const G4ParticleDefinition*, G4double ekin)
1266 {
1267   G4int pdg = theParticle->GetPDGEncoding();
1268   G4double xsection(0.);
1269   static const G4double targ_mass = 
1270     0.5*(CLHEP::proton_mass_c2 + CLHEP::neutron_mass_c2);
1271   G4double sMand = 
1272     CalcMandelstamS(ekin, theParticle->GetPDGMass(), targ_mass)*invGeV2;
1273 
1274   G4double x1 = G4Exp(G4Log(sMand)*0.0808);
1275   G4double x2 = G4Exp(G4Log(-sMand)*0.4525);
1276 
1277   if(pdg == 22) 
1278   {
1279     xsection = 0.0677*x1 + 0.129*x2;
1280   } 
1281   else if(theParticle == theNeutron)  
1282   {
1283     xsection = 21.70*x1 + 56.08*x2;
1284   } 
1285   else if(theParticle == theProton) 
1286   {
1287     xsection = 21.70*x1 + 56.08*x2;
1288   } 
1289   // pbar
1290   else if(pdg == -2212) 
1291   {
1292     xsection = 21.70*x1 + 98.39*x2;
1293   } 
1294   else if(theParticle == thePiPlus) 
1295   {
1296     xsection = 13.63*x1 + 27.56*x2;
1297   } 
1298   // pi-
1299   else if(pdg == -211) 
1300   {
1301     xsection = 13.63*x1 + 36.02*x2;
1302   } 
1303   else if(theParticle == theKPlus) 
1304   {
1305     xsection = 11.82*x1 + 8.15*x2;
1306   } 
1307   else if(theParticle == theKMinus) 
1308   {
1309     xsection = 11.82*x1 + 26.36*x2;
1310   }
1311   else if(theParticle == theK0S || theParticle == theK0L) 
1312   {
1313     xsection = 11.82*x1 + 17.25*x2;
1314   }
1315   else  
1316   {
1317     xsection = 21.70*x1 + 56.08*x2;
1318   }
1319   fTotalXsc     = xsection*CLHEP::millibarn;
1320   fInelasticXsc = 0.83*fTotalXsc;
1321   fElasticXsc   = fTotalXsc - fInelasticXsc;
1322   return fTotalXsc;
1323 }
1324 
1325 ///////////////////////////////////////////////////////////////////////////////
1326 
1327 G4double 
1328 G4HadronNucleonXsc::CoulombBarrier(const G4ParticleDefinition* theParticle, 
1329            const G4ParticleDefinition* nucleon, 
1330            G4double ekin)
1331 {
1332   G4double tR = 0.895*CLHEP::fermi;
1333   G4double pR = 0.5*CLHEP::fermi;
1334 
1335   if     ( theParticle == theProton ) pR = 0.895*CLHEP::fermi;
1336   else if( theParticle == thePiPlus ) pR = 0.663*CLHEP::fermi;
1337   else if( theParticle == theKPlus )  pR = 0.340*CLHEP::fermi;
1338 
1339   G4double pZ = theParticle->GetPDGCharge();
1340   G4double tZ = nucleon->GetPDGCharge();
1341 
1342   G4double pM = theParticle->GetPDGMass(); 
1343   G4double tM = nucleon->GetPDGMass();
1344 
1345   G4double pElab = ekin + pM;
1346 
1347   G4double totEcm  = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
1348 
1349   G4double totTcm  = totEcm - pM -tM;
1350 
1351   G4double bC = fine_structure_const*hbarc*pZ*tZ/(2.*(pR + tR));
1352 
1353   //G4cout<<"##CB: ekin = "<<ekin<<"; pElab= " << pElab
1354   //  <<"; bC = "<<bC<<"; Tcm = "<< totTcm <<G4endl;
1355 
1356   G4double ratio = (totTcm > bC) ? 1. - bC/totTcm : 0.0;
1357 
1358   // G4cout <<"ratio = "<<ratio<<G4endl;
1359   return ratio;
1360 }
1361 
1362 //////////////////////////////////////////////////////////////////////////////
1363