Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.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/G4ElasticHadrNucleusHE.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.cc (Version 9.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 //
                                                   >>  27 // $Id: G4ElasticHadrNucleusHE.cc,v 1.82 2010-11-18 22:49:57 vnivanch Exp $
                                                   >>  28 // GEANT4 tag $Name: not supported by cvs2svn $
                                                   >>  29 //
                                                   >>  30 //
 26 //  The generator of high energy hadron-nucleu     31 //  The generator of high energy hadron-nucleus elastic scattering
 27 //  The hadron kinetic energy T > 1 GeV            32 //  The hadron kinetic energy T > 1 GeV
 28 //  N.Starkov 2003.                            <<  33 //  N.  Starkov 2003.
 29 //                                                 34 //
                                                   >>  35 //  19.05.04 Variant for G4 6.1: The 'ApplyYourself' was changed
 30 //  19.11.05 The HE elastic scattering on prot     36 //  19.11.05 The HE elastic scattering on proton is added (N.Starkov)
 31 //  16.11.06 The low energy boundary is shifte     37 //  16.11.06 The low energy boundary is shifted to T = 400 MeV (N.Starkov)
 32 //  23.11.06 General cleanup, ONQ0=3, use poin     38 //  23.11.06 General cleanup, ONQ0=3, use pointer instead of particle name (VI)
 33 //  02.05.07 Scale sampled t as p^2 (VI)           39 //  02.05.07 Scale sampled t as p^2 (VI)
 34 //  15.05.07 Redesign and cleanup (V.Ivanchenk     40 //  15.05.07 Redesign and cleanup (V.Ivanchenko)
 35 //  17.05.07 cleanup (V.Grichine)                  41 //  17.05.07 cleanup (V.Grichine)
 36 //  19.04.12 Fixed reproducibility violation ( << 
 37 //  12.06.12 Fixed warnings of shadowed variab << 
 38 //                                                 42 //
 39                                                    43 
 40 #include  "G4ElasticHadrNucleusHE.hh"              44 #include  "G4ElasticHadrNucleusHE.hh"
 41 #include  "G4PhysicalConstants.hh"             << 
 42 #include  "G4SystemOfUnits.hh"                 << 
 43 #include  "Randomize.hh"                           45 #include  "Randomize.hh"
 44 #include  "G4ios.hh"                               46 #include  "G4ios.hh"
 45 #include  "G4ParticleTable.hh"                     47 #include  "G4ParticleTable.hh"
 46 #include  "G4NucleiProperties.hh"              << 
 47 #include  "G4IonTable.hh"                          48 #include  "G4IonTable.hh"
 48 #include  "G4Proton.hh"                            49 #include  "G4Proton.hh"
 49 #include  "G4PionPlus.hh"                      << 
 50 #include  "G4PionMinus.hh"                     << 
 51 #include  "G4NistManager.hh"                       50 #include  "G4NistManager.hh"
 52 #include  "G4ProductionCutsTable.hh"           << 
 53 #include  "G4MaterialCutsCouple.hh"            << 
 54 #include  "G4Material.hh"                      << 
 55 #include  "G4Element.hh"                       << 
 56 #include  "G4Log.hh"                           << 
 57 #include  "G4Exp.hh"                           << 
 58                                                    51 
 59 using namespace std;                               52 using namespace std;
 60                                                    53 
 61 const G4int G4ElasticHadrNucleusHE::fHadronCod << 
 62 {211,-211,2112,2212,321,-321,130,310,311,-311, << 
 63  3122,3222,3112,3212,3312,3322,3334,           << 
 64  -2212,-2112,-3122,-3222,-3112,-3212,-3312,-33 << 
 65                                                << 
 66 const G4int G4ElasticHadrNucleusHE::fHadronTyp << 
 67 {2,3,6,0,4,5,4,4,4,5,                          << 
 68  0,0,0,0,0,0,0,                                << 
 69  1,7,1,1,1,1,1,1,1};                           << 
 70                                                << 
 71 const G4int G4ElasticHadrNucleusHE::fHadronTyp << 
 72 {3,4,1,0,5,6,5,5,5,6,                          << 
 73  0,0,0,0,0,0,0,                                << 
 74  2,2,2,2,2,2,2,2,2};                           << 
 75                                                << 
 76 G4double G4ElasticHadrNucleusHE::fLineF[]  = { << 
 77 G4double G4ElasticHadrNucleusHE::fEnergy[] = { << 
 78 G4double G4ElasticHadrNucleusHE::fLowEdgeEnerg << 
 79 G4double G4ElasticHadrNucleusHE::fBinom[240][2 << 
 80                                                << 
 81 G4ElasticData*                                 << 
 82 G4ElasticHadrNucleusHE::fElasticData[NHADRONS] << 
 83                                                << 
 84 #ifdef G4MULTITHREADED                         << 
 85   G4Mutex G4ElasticHadrNucleusHE::elasticMutex << 
 86 #endif                                         << 
 87                                                << 
 88 G4bool G4ElasticHadrNucleusHE::fStoreToFile =  << 
 89 G4bool G4ElasticHadrNucleusHE::fRetrieveFromFi << 
 90                                                << 
 91 const G4double invGeV    =  1.0/CLHEP::GeV;    << 
 92 const G4double MbToGeV2  =  2.568;             << 
 93 const G4double GeV2      =  CLHEP::GeV*CLHEP:: << 
 94 const G4double invGeV2   =  1.0/GeV2;          << 
 95 const G4double protonM   =  CLHEP::proton_mass << 
 96 const G4double protonM2  =  protonM*protonM;   << 
 97                                                << 
 98 //////////////////////////////////////////////     54 ///////////////////////////////////////////////////////////////
                                                   >>  55 //
                                                   >>  56 //
 99                                                    57 
100 G4ElasticData::G4ElasticData(const G4ParticleD     58 G4ElasticData::G4ElasticData(const G4ParticleDefinition* p, 
101            G4int Z, G4int A, const G4double* e <<  59            G4int Z, G4double AWeight, G4double* eGeV)
102 {                                                  60 { 
103   G4double massGeV  = p->GetPDGMass()*invGeV;  <<  61   hadr     = p;
104   G4double mass2GeV2= massGeV*massGeV;         <<  62   massGeV  = p->GetPDGMass()/GeV;
                                                   >>  63   mass2GeV2= massGeV*massGeV;
                                                   >>  64   AtomicWeight = G4int(AWeight + 0.5);
105                                                    65 
106   DefineNucleusParameters(A);                  <<  66   DefineNucleusParameters(AWeight);
107   G4double limitQ2 = 35./(R1*R1);     //  (GeV << 
108                                                    67 
109   massA  = G4NucleiProperties::GetNuclearMass( <<  68   limitQ2 = 35./(R1*R1);     //  (GeV/c)^2
110   massA2 = massA*massA;                        <<  69 
111   /*                                           <<  70   G4double dQ2 = limitQ2/(ONQ2 - 1.);
112   G4cout << " G4ElasticData for " << p->GetPar <<  71 
113    << " Z= " << Z << " A= " << A << " R1= " << <<  72   TableQ2[0] = 0.0;
114    << " R2= " << R2 << G4endl;                 <<  73 
115   */                                           <<  74   for(G4int ii = 1; ii < ONQ2; ii++) 
116   for(G4int kk = 0; kk<NENERGY; ++kk)          << 
117   {                                                75   {
118     G4double elab = e[kk] + massGeV;           <<  76     TableQ2[ii] = TableQ2[ii-1]+dQ2;
119     G4double plab2= e[kk]*(e[kk] + 2.0*massGeV << 
120     G4double Q2m  = 4.0*plab2*massA2/(mass2GeV << 
121                                                << 
122     if(Z == 1 && p == G4Proton::Proton()) { Q2 << 
123                                                << 
124     maxQ2[kk] = Q2m;                           << 
125     /*                                         << 
126     G4cout << " Ekin= " << e[kk] << " Q2m= " < << 
127      << " limitQ2= " << limitQ2 << G4endl;     << 
128     */                                         << 
129   }                                                77   }
130                                                    78 
131   dQ2 = limitQ2/(G4double)(ONQ2-2);            <<  79   massA  = AWeight*amu_c2/GeV;
                                                   >>  80   massA2 = massA*massA; 
                                                   >>  81 
                                                   >>  82   for(G4int kk = 0; kk < NENERGY; kk++) 
                                                   >>  83   {
                                                   >>  84     dnkE[kk] = 0;
                                                   >>  85     G4double elab = eGeV[kk] + massGeV;
                                                   >>  86     G4double plab2= eGeV[kk]*(eGeV[kk] + 2.0*massGeV);
                                                   >>  87     G4double Q2m  = 4.0*plab2*massA2/(mass2GeV2 + massA2 + 2.*massA2*elab);
                                                   >>  88 
                                                   >>  89     if(Z == 1 && p == G4Proton::Proton()) Q2m *= 0.5;
                                                   >>  90 
                                                   >>  91     maxQ2[kk] = std::min(limitQ2, Q2m);
                                                   >>  92     TableCrossSec[ONQ2*kk] = 0.0;
                                                   >>  93 
                                                   >>  94 //    G4cout<<" kk  eGeV[kk] "<<kk<<"  "<<eGeV[kk]<<G4endl;
                                                   >>  95   }
132 }                                                  96 }
133                                                    97 
134 //////////////////////////////////////////////     98 /////////////////////////////////////////////////////////////////////////
                                                   >>  99 //
                                                   >> 100 //
135                                                   101 
136 void G4ElasticData::DefineNucleusParameters(G4 << 102 void G4ElasticData::DefineNucleusParameters(G4double A)
137 {                                                 103 {
138   switch (A) {                                 << 104   switch (AtomicWeight)
                                                   >> 105     {
139     case 207:                                     106     case 207:
140     case 208:                                     107     case 208:
141       R1       = 20.5;                            108       R1       = 20.5;
                                                   >> 109 //      R1       = 17.5;
                                                   >> 110 //      R1       = 21.3;    
142       R2       = 15.74;                           111       R2       = 15.74;
                                                   >> 112 //      R2       = 10.74;
                                                   >> 113 
143       Pnucl    = 0.4;                             114       Pnucl    = 0.4;
144       Aeff     = 0.7;                             115       Aeff     = 0.7;
145       break;                                      116       break;
146     case 237:                                     117     case 237:
147     case 238:                                     118     case 238:
148       R1       = 21.7;                            119       R1       = 21.7;    
149       R2       = 16.5;                            120       R2       = 16.5;
150       Pnucl    = 0.4;                             121       Pnucl    = 0.4;
151       Aeff     = 0.7;                             122       Aeff     = 0.7;
152       break;                                      123       break;
153     case 90:                                      124     case 90:
154     case 91:                                      125     case 91:
155       R1    = 16.5;                            << 126       R1    = 16.5*1.0;
156       R2    = 11.62;                              127       R2    = 11.62;
157       Pnucl = 0.4;                                128       Pnucl = 0.4;
158       Aeff  = 0.7;                                129       Aeff  = 0.7;
159       break;                                      130       break;
160     case 58:                                      131     case 58:
161     case 59:                                      132     case 59:
162       R1    = 15.75;                           << 133       R1    = 15.0*1.05;
163       R2    = 9.9;                                134       R2    = 9.9;
164       Pnucl = 0.45;                               135       Pnucl = 0.45;
165       Aeff  = 0.85;                               136       Aeff  = 0.85;
166       break;                                      137       break;
167     case 48:                                      138     case 48:
168     case 47:                                      139     case 47:
169       R1    = 14.0;                               140       R1    = 14.0;
170       R2    = 9.26;                               141       R2    = 9.26;
171       Pnucl = 0.31;                               142       Pnucl = 0.31;
172       Aeff  = 0.75;                               143       Aeff  = 0.75;
173       break;                                      144       break;
174     case 40:                                      145     case 40:
175     case 41:                                      146     case 41:
176       R1    = 13.3;                               147       R1    = 13.3;
177       R2    = 9.26;                               148       R2    = 9.26;
178       Pnucl = 0.31;                               149       Pnucl = 0.31;
179       Aeff  = 0.75;                               150       Aeff  = 0.75;
180       break;                                      151       break;
181     case 28:                                      152     case 28:
182     case 29:                                      153     case 29:
183       R1    = 12.0;                               154       R1    = 12.0;
184       R2    = 7.64;                               155       R2    = 7.64;
185       Pnucl = 0.253;                              156       Pnucl = 0.253;
186       Aeff  = 0.8;                                157       Aeff  = 0.8;
187       break;                                      158       break;
188     case 16:                                      159     case 16:
189       R1    = 10.50;                              160       R1    = 10.50;
190       R2    = 5.5;                                161       R2    = 5.5;
191       Pnucl = 0.7;                                162       Pnucl = 0.7;
192       Aeff  = 0.98;                               163       Aeff  = 0.98;
193       break;                                      164       break;
194     case 12:                                      165     case 12:
195       R1    = 9.3936;                             166       R1    = 9.3936;
196       R2    = 4.63;                               167       R2    = 4.63;
197       Pnucl = 0.7;                                168       Pnucl = 0.7;
                                                   >> 169 //      Pnucl = 0.5397;
198       Aeff  = 1.0;                                170       Aeff  = 1.0;
199       break;                                      171       break;
200     case 11:                                      172     case 11:
201       R1    = 9.0;                                173       R1    = 9.0;
202       R2    = 5.42;                               174       R2    = 5.42;
203       Pnucl = 0.19;                               175       Pnucl = 0.19;
                                                   >> 176 //      Pnucl = 0.5397;
204       Aeff  = 0.9;                                177       Aeff  = 0.9;
205       break;                                      178       break;
206     case 9:                                       179     case 9:
207       R1    = 9.9;                                180       R1    = 9.9;
208       R2    = 6.5;                                181       R2    = 6.5;
209       Pnucl = 0.690;                              182       Pnucl = 0.690;
210       Aeff  = 0.95;                               183       Aeff  = 0.95;
211       break;                                      184       break;
212     case 4:                                       185     case 4:
213       R1    = 5.3;                                186       R1    = 5.3;   
214       R2    = 3.7;                                187       R2    = 3.7;
215       Pnucl = 0.4;                                188       Pnucl = 0.4;
216       Aeff  = 0.75;                               189       Aeff  = 0.75;
217       break;                                      190       break;
218     case 1:                                       191     case 1:
219       R1    = 4.5;                                192       R1    = 4.5;   
220       R2    = 2.3;                                193       R2    = 2.3;
221       Pnucl = 0.177;                              194       Pnucl = 0.177;
222       Aeff  = 0.9;                                195       Aeff  = 0.9;
223       break;                                      196       break;
224     default:                                      197     default:
225       R1    = 4.45*G4Exp(G4Log((G4double)(A -  << 198       R1    = 4.45*std::pow(A - 1.,0.309)*0.9;
226       R2    = 2.3 *G4Exp(G4Log((G4double)A)* 0 << 199       R2    = 2.3 *std::pow(A, 0.36);
227                                                   200 
228       if(A < 100 && A > 3) { Pnucl = 0.176 + 0 << 201       if(A < 100 && A > 3) Pnucl = 0.176 + 0.00275*A;
229       else                 { Pnucl = 0.4; }    << 202       else                 Pnucl = 0.4;
230       //G4cout<<" Deault: A= "<<A<<"  R1 R2 Ae << 203 
231       //      <<Aeff<<"  "<<Pnucl<<G4endl;     << 204 //G4cout<<" Deault: A= "<<A<<"  R1 R2 Aeff Pnucl "<<R1<<"  "<<R2<<"  "
232                                                << 205 //      <<Aeff<<"  "<<Pnucl<<G4endl;
233       if(A >= 100)               { Aeff = 0.7; << 206 
234       else if(A < 100 && A > 75) { Aeff = 1.5  << 207       if(A >= 100)               Aeff = 0.7;
235       else                       { Aeff = 0.9; << 208       else if(A < 100 && A > 75) Aeff = 1.5 - 0.008*A;
                                                   >> 209       else                       Aeff = 0.9;
236       break;                                      210       break;
237   }                                            << 211     }
238   //G4cout<<" Result: A= "<<A<<"  R1 R2 Aeff P << 212 //G4cout<<" Result: A= "<<A<<"  R1 R2 Aeff Pnucl "<<R1<<"  "<<R2<<"  "
239   //      <<Aeff<<"  "<<Pnucl<<G4endl;         << 213 //      <<Aeff<<"  "<<Pnucl<<G4endl;
240 }                                                 214 }
241                                                   215 
242 //////////////////////////////////////////////    216 ////////////////////////////////////////////////////////////////////
                                                   >> 217 //
                                                   >> 218 //  The constructor for the generating of events
                                                   >> 219 //
243                                                   220 
244 G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE    221 G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE(const G4String& name)
245   : G4HadronElastic(name), fDirectory(nullptr) << 222   : G4HadronElastic(name)
246 {                                                 223 {
247   dQ2 = hMass = hMass2 = hLabMomentum = hLabMo << 224   dQ2 = hMass = hMass2 =  hLabMomentum = hLabMomentum2 = MomentumCM = HadrEnergy 
248     = R1 = R2 = Pnucl = Aeff = HadrTot = HadrS    225     = R1 = R2 = Pnucl = Aeff = HadrTot = HadrSlope = HadrReIm = TotP = DDSect2
249     = DDSect3 = ConstU = Slope1 = Slope2 = Coe << 226     = DDSect3 = ConstU = FmaxT = Slope1 = Slope2 = Coeff1 = Coeff2 = MaxTR 
250     = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = << 227     = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = 0.0;
251   iHadrCode = iHadron = iHadron1 = 0;          << 228   NumbN = iHadrCode = iHadron = 0;
252                                                   229 
253   verboseLevel = 0;                               230   verboseLevel = 0;
254   ekinLowLimit = 400.0*CLHEP::MeV;             << 231   plabLowLimit = 20.0*MeV;
                                                   >> 232   lowestEnergyLimit = 0.0;
                                                   >> 233   //Description();
                                                   >> 234 
                                                   >> 235   MbToGeV2  =  2.568;
                                                   >> 236   sqMbToGeV =  1.602;
                                                   >> 237   Fm2ToGeV2 =  25.68;
                                                   >> 238   GeV2      =  GeV*GeV;
                                                   >> 239   protonM   =  proton_mass_c2/GeV;
                                                   >> 240   protonM2  =  protonM*protonM;
255                                                   241 
256   BoundaryP[0]=9.0; BoundaryTG[0]=5.0;Boundary << 242   BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.;
257   BoundaryP[1]=20.0;BoundaryTG[1]=1.5;Boundary    243   BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.;
258   BoundaryP[2]=5.0; BoundaryTG[2]=1.0;Boundary    244   BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5;
259   BoundaryP[3]=8.0; BoundaryTG[3]=3.0;Boundary    245   BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.;
260   BoundaryP[4]=7.0; BoundaryTG[4]=3.0;Boundary    246   BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.;
261   BoundaryP[5]=5.0; BoundaryTG[5]=2.0;Boundary    247   BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.;
262   BoundaryP[6]=5.0; BoundaryTG[6]=1.5;Boundary    248   BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0;
263                                                   249 
                                                   >> 250   Binom();
                                                   >> 251   // energy in GeV
                                                   >> 252   Energy[0] = 0.4;
                                                   >> 253   Energy[1] = 0.6;
                                                   >> 254   Energy[2] = 0.8;
                                                   >> 255   LowEdgeEnergy[0] = 0.0;
                                                   >> 256   LowEdgeEnergy[1] = 0.5;
                                                   >> 257   LowEdgeEnergy[2] = 0.7;
                                                   >> 258   G4double e = 1.0;
                                                   >> 259   G4double f = std::pow(10.,0.1);
                                                   >> 260   for(G4int i=3; i<NENERGY; i++) {
                                                   >> 261     Energy[i] = e;
                                                   >> 262     LowEdgeEnergy[i] = e/f;
                                                   >> 263     e *= f*f;
                                                   >> 264   }
264   nistManager = G4NistManager::Instance();        265   nistManager = G4NistManager::Instance();
265                                                   266 
266   if(fEnergy[0] == 0.0) {                      << 267   // PDG code for hadrons
267 #ifdef G4MULTITHREADED                         << 268   G4int ic[NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311,
268     G4MUTEXLOCK(&elasticMutex);                << 269       3122,3222,3112,3212,3312,3322,3334,
269     if(fEnergy[0] == 0.0) {                    << 270       -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334};
270 #endif                                         << 271   // internal index 
271       isMaster = true;                         << 272   G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
272       Binom();                                 << 273       0,0,0,0,0,0,0,
273       // energy in GeV                         << 274       1,7,1,1,1,1,1,1,1};
274       fEnergy[0] = 0.4;                        << 275 
275       fEnergy[1] = 0.6;                        << 276   G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
276       fEnergy[2] = 0.8;                        << 277                         0,0,0,0,0,0,0,
277       fEnergy[3] = 1.0;                        << 278                         2,2,2,2,2,2,2,2,2};
278       fLowEdgeEnergy[0] = 0.0;                 << 
279       fLowEdgeEnergy[1] = 0.5;                 << 
280       fLowEdgeEnergy[2] = 0.7;                 << 
281       fLowEdgeEnergy[3] = 0.9;                 << 
282       G4double f = G4Exp(G4Log(10.)*0.1);      << 
283       G4double e = f*f;                        << 
284       for(G4int i=4; i<NENERGY; ++i) {         << 
285   fEnergy[i] = e;                              << 
286   fLowEdgeEnergy[i] = e/f;                     << 
287   e *= f*f;                                    << 
288       }                                        << 
289       if(verboseLevel > 0) {                   << 
290   G4cout << "### G4ElasticHadrNucleusHE: energ << 
291   for(G4int i=0; i<NENERGY; ++i) {             << 
292     G4cout << "  " << i << "   " << fLowEdgeEn << 
293      << "  " << fEnergy[i] << G4endl;          << 
294   }                                            << 
295       }                                        << 
296 #ifdef G4MULTITHREADED                         << 
297     }                                          << 
298     G4MUTEXUNLOCK(&elasticMutex);              << 
299 #endif                                         << 
300   }                                            << 
301 }                                              << 
302                                                   279 
303 ////////////////////////////////////////////// << 280   for(G4int j=0; j<NHADRONS; j++) 
                                                   >> 281   {
                                                   >> 282     HadronCode[j]  = ic[j];
                                                   >> 283     HadronType[j]  = id[j];
                                                   >> 284     HadronType1[j] = id1[j];
304                                                   285 
305 void G4ElasticHadrNucleusHE::ModelDescription( << 286     for(G4int k = 0; k < 93; k++) { SetOfElasticData[j][k] = 0; }
306 {                                              << 287   } 
307   outFile << "G4ElasticHadrNucleusHE is a hadr << 
308     << "model developed by N. Starkov which us << 
309     << "parameterization to calculate the fina << 
310     << "for all hadrons with incident momentum << 
311 }                                                 288 }
312                                                   289 
313 ////////////////////////////////////////////// << 
314                                                   290 
315 G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusH << 291 void G4ElasticHadrNucleusHE::Description() const
316 {                                                 292 {
317   if(isMaster) {                               << 293   char* dirName = getenv("G4PhysListDocDir");
318     for(G4int j = 0; j < NHADRONS; ++j) {      << 294   if (dirName) {
319       for(G4int k = 0; k < ZMAX; ++k) {        << 295     std::ofstream outFile;
320   G4ElasticData* ptr = fElasticData[j][k];     << 296     G4String outFileName = GetModelName() + ".html";
321   if(ptr) {                                    << 297     G4String pathName = G4String(dirName) + "/" + outFileName;
322     delete ptr;                                << 298     outFile.open(pathName);
323     fElasticData[j][k] = nullptr;              << 299     outFile << "<html>\n";
324     for(G4int l = j+1; l < NHADRONS; ++l) {    << 300     outFile << "<head>\n";
325       if(ptr == fElasticData[l][k]) { fElastic << 301 
326     }                                          << 302     outFile << "<title>Description of G4ElasticHadrNucleusHE Model</title>\n";
327   }                                            << 303     outFile << "</head>\n";
328       }                                        << 304     outFile << "<body>\n";
329     }                                          << 305 
330     delete fDirectory;                         << 306     outFile << "G4ElasticHadrNucleusHE is a hadron-nucleus elastic scattering\n"
331     fDirectory = nullptr;                      << 307             << "model developed by N. Starkov which uses a Glauber model\n"
                                                   >> 308             << "parameterization to calculate the final state.  It is valid\n"
                                                   >> 309             << "for all hadrons with incident energies above 1 GeV.\n";
                                                   >> 310 
                                                   >> 311     outFile << "</body>\n";
                                                   >> 312     outFile << "</html>\n";
                                                   >> 313     outFile.close();
332   }                                               314   }
333 }                                                 315 }
334                                                   316 
                                                   >> 317 
335 //////////////////////////////////////////////    318 ///////////////////////////////////////////////////////////////////
                                                   >> 319 //
                                                   >> 320 //
336                                                   321 
337 void G4ElasticHadrNucleusHE::InitialiseModel() << 322 G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE()
338 {                                                 323 {
339   if(!isMaster) { return; }                    << 324   for(G4int j = 0; j < NHADRONS; j++) 
340   G4ProductionCutsTable* theCoupleTable=       << 325   {
341     G4ProductionCutsTable::GetProductionCutsTa << 326     for(G4int k = 0; k < 93; k++) 
342   G4int numOfCouples = (G4int)theCoupleTable-> << 327     {
343                                                << 328       if(SetOfElasticData[j][k]) delete SetOfElasticData[j][k];
344   for(G4int i=0; i<2; ++i) {                   << 
345     const G4ParticleDefinition* p = G4PionPlus << 
346     if(1 == i) { p = G4PionMinus::PionMinus(); << 
347     iHadrCode = fHadronCode[i];                << 
348     iHadron   = fHadronType[i];                << 
349     iHadron1  = fHadronType1[i];               << 
350     hMass     = p->GetPDGMass()*invGeV;        << 
351     hMass2    = hMass*hMass;                   << 
352     for(G4int j=0; j<numOfCouples; ++j) {      << 
353       auto mat = theCoupleTable->GetMaterialCu << 
354       auto elmVec = mat->GetElementVector();   << 
355       std::size_t numOfElem = mat->GetNumberOf << 
356       for(std::size_t k=0; k<numOfElem; ++k) { << 
357         G4int Z = std::min((*elmVec)[k]->GetZa << 
358         if(!fElasticData[i][Z]) {              << 
359           if(1 == i && Z > 1) {                << 
360             fElasticData[1][Z] = fElasticData[ << 
361           } else {                             << 
362             FillData(p, i, Z);                 << 
363           }                                    << 
364         }                                      << 
365       }                                        << 
366     }                                             329     }
367   }                                            << 330   } 
368 }                                                 331 }
369                                                   332 
370 //////////////////////////////////////////////    333 ////////////////////////////////////////////////////////////////////
                                                   >> 334 //
                                                   >> 335 //
371                                                   336 
372 G4double                                          337 G4double 
373 G4ElasticHadrNucleusHE::SampleInvariantT(const    338 G4ElasticHadrNucleusHE::SampleInvariantT(const G4ParticleDefinition* p,
374            G4double inLabMom,                     339            G4double inLabMom, 
375            G4int iZ, G4int A)                  << 340            G4int Z, G4int N)
376 {                                                 341 {
377   G4double mass = p->GetPDGMass();             << 342   G4double plab  = inLabMom/GeV;   // (GeV/c)
378   G4double kine = sqrt(inLabMom*inLabMom + mas << 343   G4double Q2 = 0;
379   if(kine <= ekinLowLimit) {                   << 344 
380     return G4HadronElastic::SampleInvariantT(p << 
381   }                                            << 
382   G4int Z = std::min(iZ,ZMAX-1);               << 
383   G4double Q2 = 0.0;                           << 
384   iHadrCode = p->GetPDGEncoding();                345   iHadrCode = p->GetPDGEncoding();
385                                                   346 
386   // below computations in GeV/c               << 347   NumbN = N;
387   hMass  = mass*invGeV;                        << 
388   hMass2 = hMass*hMass;                        << 
389   G4double plab = inLabMom*invGeV;             << 
390   G4double tmax = pLocalTmax*invGeV2;          << 
391                                                   348 
392   if(verboseLevel > 1) {                       << 349   if(verboseLevel > 1)
393     G4cout<< "G4ElasticHadrNucleusHE::SampleT: << 350   {
                                                   >> 351     G4cout<< " G4ElasticHadrNucleusHE::SampleT: " 
394     << " for " << p->GetParticleName()            352     << " for " << p->GetParticleName() 
395     << " at Z= " << Z << " A= " << A           << 353     << " at Z= " << Z << " A= " << N
396     << " plab(GeV)= " << plab                     354     << " plab(GeV)= " << plab
397     << " hadrCode= " << iHadrCode              << 
398     << G4endl;                                    355     << G4endl;
399   }                                               356   }
400   iHadron = -1;                                << 
401   G4int idx;                                      357   G4int idx;
402   for(idx=0; idx<NHADRONS; ++idx) {            << 358 
403     if(iHadrCode == fHadronCode[idx]) {        << 359   for( idx = 0 ; idx < NHADRONS; idx++) 
404       iHadron = fHadronType[idx];              << 360   {
405       iHadron1 = fHadronType1[idx];            << 361     if(iHadrCode == HadronCode[idx]) break;
406       break;                                   << 
407     }                                          << 
408   }                                               362   }
409   // Hadron is not in the list                 << 
410   if(0 > iHadron) { return 0.0; }              << 
411                                                   363 
412   if(Z==1) {                                   << 364   // Hadron is not in the list
413     Q2 = HadronProtonQ2(plab, tmax);           << 
414                                                   365 
415     if (verboseLevel>1) {                      << 366   if( idx >= NHADRONS ) return Q2;
416       G4cout<<"  Proton : Q2  "<<Q2<<G4endl;   << 
417     }                                          << 
418   } else {                                     << 
419     const G4ElasticData* ElD1 = fElasticData[i << 
420                                                   367 
421     // Construct elastic data                  << 368   iHadron = HadronType[idx];
422     if(!ElD1) {                                << 369   iHadrCode = HadronCode[idx];
423       FillData(p, idx, Z);                     << 
424       ElD1 = fElasticData[idx][Z];             << 
425       if(!ElD1) { return 0.0; }                << 
426     }                                          << 
427                                                   370 
428     // sample scattering                       << 371   if(Z==1)
429     Q2 = HadronNucleusQ2_2(ElD1, plab, tmax);  << 372     {
                                                   >> 373       hMass  = p->GetPDGMass()/GeV;
                                                   >> 374       hMass2 = hMass*hMass;
430                                                   375 
431     if(verboseLevel > 1) {                     << 376       G4double T = sqrt(plab*plab+hMass2)-hMass;
432       G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< "  << 
433             << Q2/tmax <<G4endl;               << 
434     }                                          << 
435   }                                            << 
436   return Q2*GeV2;                              << 
437 }                                              << 
438                                                   377 
439 ////////////////////////////////////////////// << 378       if(T > 0.4) Q2 = HadronProtonQ2(p, plab);
440                                                   379 
441 void G4ElasticHadrNucleusHE::FillData(const G4 << 380       if (verboseLevel>1)
442                                       G4int id << 381   G4cout<<"  Proton : Q2  "<<Q2<<G4endl;
443 {                                              << 
444 #ifdef G4MULTITHREADED                         << 
445   G4MUTEXLOCK(&elasticMutex);                  << 
446   if(!fElasticData[idx][Z]) {                  << 
447 #endif                                         << 
448     G4int A = G4lrint(nistManager->GetAtomicMa << 
449     G4ElasticData* pElD = new G4ElasticData(p, << 
450     if(fRetrieveFromFile) {                    << 
451       std::ostringstream ss;                   << 
452       InFileName(ss, p, Z);                    << 
453       std::ifstream infile(ss.str(), std::ios: << 
454       for(G4int i=0; i<NENERGY; ++i) {         << 
455   if(ReadLine(infile, pElD->fCumProb[i])) {    << 
456     continue;                                  << 
457   } else {                                     << 
458     fRetrieveFromFile = false;                 << 
459           break;                               << 
460   }                                            << 
461       }                                        << 
462       infile.close();                          << 
463     }                                          << 
464     R1     = pElD->R1;                         << 
465     R2     = pElD->R2;                         << 
466     Aeff   = pElD->Aeff;                       << 
467     Pnucl  = pElD->Pnucl;                      << 
468     dQ2    = pElD->dQ2;                        << 
469     if(verboseLevel > 0) {                     << 
470       G4cout<<"### FillData for " << p->GetPar << 
471       << " Z= " << Z << " idx= " << idx << " i << 
472       <<" iHadron1= " << iHadron1 << " iHadrCo << 
473             <<"\n   R1= " << R1 << " R2= " <<  << 
474       <<" Pnucl= " << Pnucl << G4endl;         << 
475     }                                             382     }
                                                   >> 383   else
                                                   >> 384     {
                                                   >> 385       G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
476                                                   386 
477     if(!fRetrieveFromFile) {                   << 387       // Construct elastic data
478       for(G4int i=0; i<NENERGY; ++i) {         << 388       if(!ElD1) 
479   G4double T = fEnergy[i];                     << 389   {
480   hLabMomentum2 = T*(T + 2.*hMass);            << 390     G4double AWeight = nistManager->GetAtomicMassAmu(Z);
481   hLabMomentum  = std::sqrt(hLabMomentum2);    << 391     ElD1 = new  G4ElasticData(p, Z, AWeight, Energy);
482   HadrEnergy = hMass + T;                      << 392     SetOfElasticData[idx][Z] = ElD1;
483   DefineHadronValues(Z);                       << 393     
484   Q2max = pElD->maxQ2[i];                      << 394     if(verboseLevel > 1)
485                                                << 395       {
486   G4int length  = FillFq2(A);                  << 396         G4cout<< " G4ElasticHadrNucleusHE::SampleT:  new record " << idx
487   (pElD->fCumProb[i]).reserve(length);         << 397         << " for " << p->GetParticleName() << " Z= " << Z
488   G4double norm = 1.0/fLineF[length-1];        << 398         << G4endl;
489                                                << 399       }
490   if(verboseLevel > 0) {                       << 400   }  
491     G4cout << "### i= " << i << " Z= " << Z << << 401       hMass          = ElD1->massGeV;
492      << " length= " << length << " Q2max= " << << 402       hMass2         = ElD1->mass2GeV2;
493   }                                            << 403       G4double M     = ElD1->massA;
                                                   >> 404       G4double M2    = ElD1->massA2;
                                                   >> 405       G4double plab2 = plab*plab;
                                                   >> 406       G4double Q2max = 4.*plab2*M2/
                                                   >> 407   (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
                                                   >> 408 
                                                   >> 409       // sample scattering
                                                   >> 410       G4double T = sqrt(plab2+hMass2)-hMass;
494                                                   411 
495   (pElD->fCumProb[i]).push_back(0.0);          << 412       if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
496   for(G4int ii=1; ii<length-1; ++ii) {         << 
497     (pElD->fCumProb[i]).push_back(fLineF[ii]*n << 
498     if(verboseLevel > 2) {                     << 
499       G4cout << "    ii= " << ii << " val= "   << 
500        << (pElD->fCumProb[i])[ii] << G4endl;   << 
501     }                                          << 
502   }                                            << 
503   (pElD->fCumProb[i]).push_back(1.0);          << 
504       }                                        << 
505     }                                          << 
506                                                   413 
507     if(fStoreToFile) {                         << 414       if(verboseLevel > 1)
508       std::ostringstream ss;                   << 415   G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< "  t/tmax= " << Q2/Q2max <<G4endl;
509       OutFileName(ss, p, Z);                   << 
510       std::ofstream fileout(ss.str());         << 
511       for(G4int i=0; i<NENERGY; ++i) {         << 
512   WriteLine(fileout, pElD->fCumProb[i]);       << 
513       }                                        << 
514       fileout.close();                         << 
515     }                                          << 
516                                                << 
517     if(verboseLevel > 0) {                     << 
518       G4cout << " G4ElasticHadrNucleusHE::Fill << 
519        << " for " << p->GetParticleName() << " << 
520        << " A= " << A << G4endl;               << 
521     }                                             416     }
522     fElasticData[idx][Z] = pElD;               << 417   return  Q2*GeV2;
523                                                << 
524 #ifdef G4MULTITHREADED                         << 
525   }                                            << 
526   G4MUTEXUNLOCK(&elasticMutex);                << 
527 #endif                                         << 
528 }                                                 418 }
529                                                   419 
530 ////////////////////////////////////////////// << 420 //////////////////////////////////////////////////////////////////////////
                                                   >> 421 //
                                                   >> 422 //
531                                                   423 
532 void G4ElasticHadrNucleusHE::InterpolateHN(G4i << 424 G4double 
533                    const G4double C0P[], const << 425 G4ElasticHadrNucleusHE::SampleT(const G4ParticleDefinition* p,
534                    const G4double B0P[], const << 426         G4double inLabMom, 
                                                   >> 427         G4int Z, G4int N)
535 {                                                 428 {
536   G4int i;                                     << 429   return SampleInvariantT(p, inLabMom, Z, N);
537                                                << 
538   for(i=1; i<n; ++i) { if(hLabMomentum <= EnP[ << 
539   if(i == n) { i = n - 1; }                    << 
540                                                << 
541   Coeff0 = LineInterpol(EnP[i], EnP[i-1], C0P[ << 
542   Coeff1 = LineInterpol(EnP[i], EnP[i-1], C1P[ << 
543   Slope0 = LineInterpol(EnP[i], EnP[i-1], B0P[ << 
544   Slope1 = LineInterpol(EnP[i], EnP[i-1], B1P[ << 
545                                                << 
546 //  G4cout<<"  InterpolHN:  n i "<<n<<"  "<<i< << 
547 //        <<hLabMomentum<<G4endl;              << 
548 }                                                 430 }
549                                                   431 
550 //////////////////////////////////////////////    432 //////////////////////////////////////////////////////////////////////////
                                                   >> 433 //
                                                   >> 434 //
551                                                   435 
552 G4double                                       << 436 G4double G4ElasticHadrNucleusHE::
553 G4ElasticHadrNucleusHE::HadronNucleusQ2_2(cons << 437                           HadronNucleusQ2_2(G4ElasticData* pElD, G4int Z, 
554                                           G4do << 438                                             G4double plab, G4double tmax)
555 {                                                 439 {
556   G4double ekin  = std::sqrt(hMass2 + plab*pla << 440   G4double LineFq2[ONQ2];
                                                   >> 441 
                                                   >> 442   G4double Rand = G4UniformRand();
                                                   >> 443 
                                                   >> 444   G4int      iNumbQ2 = 0;
                                                   >> 445   G4double   Q2 = 0.0;
                                                   >> 446 
                                                   >> 447   G4double ptot2 = plab*plab;
                                                   >> 448   G4double ekin  = std::sqrt(hMass2 + ptot2) - hMass;
                                                   >> 449 
                                                   >> 450   if(verboseLevel > 1)
                                                   >> 451     G4cout<<"Q2_2: ekin  plab  "<<ekin<<"    "<<plab<<"  tmax "<<tmax<<G4endl;
557                                                   452 
558   if(verboseLevel > 1) {                       << 
559     G4cout<<"Q2_2: ekin(GeV)= " << ekin << "   << 
560           <<"  tmax(GeV2)= " << tmax <<G4endl; << 
561   }                                            << 
562   // Find closest energy bin                      453   // Find closest energy bin
563   G4int idx;                                   << 454   G4int NumbOnE; 
564   for(idx=0; idx<NENERGY-1; ++idx) {           << 455   for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ ) 
565     if(ekin <= fLowEdgeEnergy[idx+1]) { break; << 456   {
                                                   >> 457     if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break;
566   }                                               458   }
567   //G4cout << "   idx= " << idx << G4endl;     << 459   G4double* dNumbQ2 = pElD->TableQ2;
                                                   >> 460 
                                                   >> 461   G4int index = NumbOnE*ONQ2;
568                                                   462 
569   // Select kinematics for node energy            463   // Select kinematics for node energy
570   R1    = pElD->R1;                            << 464   G4double T     = Energy[NumbOnE];
571   dQ2   = pElD->dQ2;                           << 465   hLabMomentum2  = T*(T + 2.*hMass);
572   Q2max = pElD->maxQ2[idx];                    << 466   G4double Q2max = pElD->maxQ2[NumbOnE];
573   G4int length = (G4int)(pElD->fCumProb[idx]). << 467   G4int length   = pElD->dnkE[NumbOnE];
574                                                   468 
575   G4double Rand = G4UniformRand();             << 469   // Build vector
                                                   >> 470   if(length == 0) 
                                                   >> 471     {
                                                   >> 472       R1    = pElD->R1;
                                                   >> 473       R2    = pElD->R2;
                                                   >> 474       Aeff  = pElD->Aeff;
                                                   >> 475       Pnucl = pElD->Pnucl;
                                                   >> 476       hLabMomentum = std::sqrt(hLabMomentum2);
                                                   >> 477  
                                                   >> 478       DefineHadronValues(Z);
576                                                   479 
577   G4int iNumbQ2 = 0;                           << 480       if(verboseLevel >0)
578   for(iNumbQ2=1; iNumbQ2<length; ++iNumbQ2) {  << 481   {
579     if(Rand <= (pElD->fCumProb[idx])[iNumbQ2]) << 482     G4cout<<"1  plab  T  "<<plab<<"  "<<T<<"  sigTot  B  ReIm  "
580   }                                            << 483     <<HadrTot<<"  "<<HadrSlope<<"  "<<HadrReIm<<G4endl;
581   iNumbQ2 = std::min(iNumbQ2, length - 1);     << 484     G4cout<<"  R1  R2  Aeff  p  "<<R1<<"  "<<R2<<"  "<<Aeff<<"  "
582   G4double Q2 = GetQ2_2(iNumbQ2, length, pElD- << 485     <<Pnucl<<G4endl;
583   Q2 = std::min(Q2, Q2max);                    << 486   }
584   Q2 *= tmax/Q2max;                            << 
585                                                   487 
586   if(verboseLevel > 1) {                       << 488       //pElD->CrossSecMaxQ2[NumbOnE] = 1.0;
                                                   >> 489 
                                                   >> 490       if(verboseLevel > 1)
                                                   >> 491   G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE 
                                                   >> 492         << " length= " << length 
                                                   >> 493         << " Q2max= " << Q2max 
                                                   >> 494         << " ekin= " << ekin <<G4endl;
                                                   >> 495     
                                                   >> 496       pElD->TableCrossSec[index] = 0;
                                                   >> 497 
                                                   >> 498 
                                                   >> 499       dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
                                                   >> 500 
                                                   >> 501       GetHeavyFq2(NumbN, LineFq2);  //  %%%%%%%%%%%%%%%%%%%%%%%%%
                                                   >> 502 
                                                   >> 503       for(G4int ii=0; ii<ONQ2; ii++)
                                                   >> 504   {
                                                   >> 505     //if(verboseLevel > 2)
                                                   >> 506     //  G4cout<<"  ii LineFq2  "<<ii<<"  "<<LineFq2[ii]/LineFq2[ONQ2-1]
                                                   >> 507     //  <<"  dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl;
                                                   >> 508 
                                                   >> 509     pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
                                                   >> 510   }
                                                   >> 511     
                                                   >> 512       pElD->dnkE[NumbOnE] = ONQ2;
                                                   >> 513       length = ONQ2;
                                                   >> 514     } 
                                                   >> 515 
                                                   >> 516   G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
                                                   >> 517 
                                                   >> 518   for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ ) 
                                                   >> 519     {
                                                   >> 520       if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
                                                   >> 521     }
                                                   >> 522   Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
                                                   >> 523 
                                                   >> 524   if(tmax < Q2max) Q2 *= tmax/Q2max;
                                                   >> 525 
                                                   >> 526   if(verboseLevel > 1)
587     G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<"    527     G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2 
588     << " rand= " << Rand << " Q2max= " << Q2ma << 528     << " rand= " << Rand << G4endl;
589           << " tmax= " << tmax << G4endl;      << 529   
590   }                                            << 
591   return Q2;                                      530   return Q2;
592 }                                                 531 }       
593                                                   532 
594 //////////////////////////////////////////////    533 ///////////////////////////////////////////////////////////////////////
595 //                                                534 //
596 //  The randomization of one dimensional array    535 //  The randomization of one dimensional array 
597 //                                                536 //
598                                                << 537 G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int kk, G4double * Q,
599 G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int << 538            G4double * F, G4double ranUni)
600            const std::vector<G4double>& F,     << 
601                                          G4dou << 
602 {                                                 539 {
603   //G4cout << "GetQ2_2 kk= " << kk << " kmax=  << 540   G4double ranQ2;
604   //   << F.size() << "  rand= " << ranUni <<  << 
605   if(kk == kmax-1) {                           << 
606     G4double X1 = dQ2*kk;                      << 
607     G4double F1 = F[kk-1];                     << 
608     G4double X2 = Q2max;                       << 
609     G4double xx = R1*(X2 - X1);                << 
610     xx = (xx > 20.) ? 0.0 : G4Exp(-xx);        << 
611     G4double Y = X1 - G4Log(1.0 - (ranUni - F1 << 
612     return Y;                                  << 
613   }                                            << 
614   G4double F1, F2, F3, X1, X2, X3;             << 
615                                                   541 
616   if(kk == 1 || kk == 0) {                     << 542   G4double F1  = F[kk-2];
617     F1 = F[0];                                 << 543   G4double F2  = F[kk-1];
618     F2 = F[1];                                 << 544   G4double F3  = F[kk];
619     F3 = F[2];                                 << 545   G4double X1  = Q[kk-2];
620     X1 = 0.0;                                  << 546   G4double X2  = Q[kk-1];
621     X2 = dQ2;                                  << 547   G4double X3  = Q[kk];
622     X3 = dQ2*2;                                << 548 
623   } else {                                     << 549   if(verboseLevel > 2) 
624     F1 = F[kk-2];                              << 
625     F2 = F[kk-1];                              << 
626     F3 = F[kk];                                << 
627     X1 = dQ2*(kk-2);                           << 
628     X2 = dQ2*(kk-1);                           << 
629     X3 = dQ2*kk;                               << 
630   }                                            << 
631   if(verboseLevel > 1) {                       << 
632     G4cout << "GetQ2_2 kk= " << kk << " X2= "     550     G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3 
633      << " F2= " << F2 << " F3= " << F3 << " Rn    551      << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl;
                                                   >> 552 
                                                   >> 553   if(kk == 1 || kk == 0)
                                                   >> 554   {
                                                   >> 555      F1  = F[0]; 
                                                   >> 556      F2  = F[1];
                                                   >> 557      F3  = F[2];
                                                   >> 558      X1  = Q[0];
                                                   >> 559      X2  = Q[1];
                                                   >> 560      X3  = Q[2];
634   }                                               561   }
635                                                   562 
636   G4double F12 = F1*F1;                           563   G4double F12 = F1*F1;
637   G4double F22 = F2*F2;                           564   G4double F22 = F2*F2;
638   G4double F32 = F3*F3;                           565   G4double F32 = F3*F3;
639                                                   566 
640   G4double D0  = F12*F2+F1*F32+F3*F22-F32*F2-F    567   G4double D0  = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
641                                                   568 
642   if(verboseLevel > 2) {                       << 569   if(verboseLevel > 2) 
643     G4cout << "       X1= " << X1 << " F1= " <    570     G4cout << "       X1= " << X1 << " F1= " << F1 << "  D0= " 
644            << D0 << G4endl;                       571            << D0 << G4endl; 
645   }                                            << 572 
646   G4double Y;                                  << 573   if(std::abs(D0) < 0.00000001)
647   if(std::abs(D0) < 1.e-9) {                   << 574     { 
648     Y = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2) << 575       ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
649   } else {                                     << 576     }
650     G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F << 577   else    
651     G4double DB = X2*F12+X1*F32+X3*F22-X2*F32- << 578     {
652     G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F2 << 579       G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
                                                   >> 580       G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22;
                                                   >> 581       G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
653              -X1*F2*F32-X2*F3*F12-X3*F1*F22;      582              -X1*F2*F32-X2*F3*F12-X3*F1*F22;
654     Y = (DA*ranUni*ranUni + DB*ranUni + DC)/D0 << 583       ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
655   }                                            << 584     }
656   return Y;                                    << 585   return ranQ2;         //  MeV^2
657 }                                                 586 }
658                                                   587 
659 //////////////////////////////////////////////    588 ////////////////////////////////////////////////////////////////////////
660                                                << 589 //
661 G4int G4ElasticHadrNucleusHE::FillFq2(G4int A) << 590 //
                                                   >> 591 G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Nucleus, G4double* LineF)
662 {                                                 592 {
                                                   >> 593   G4int ii, jj, aSimp;
663   G4double curQ2, curSec;                         594   G4double curQ2, curSec;
664   G4double curSum = 0.0;                          595   G4double curSum = 0.0;
665   G4double totSum = 0.0;                          596   G4double totSum = 0.0;
666                                                   597 
667   G4double ddQ2 = dQ2*0.1;                     << 598   G4double ddQ2 = dQ2/20;
668   G4double Q2l  = 0.0;                         << 599   G4double Q2l  = 0;
669                                                   600 
670   G4int ii = 0;                                << 601   LineF[0] = 0;
671   for(ii=1; ii<ONQ2-1; ++ii) {                 << 602   for(ii = 1; ii<ONQ2; ii++)
672     curSum = curSec = 0.0;                     << 603     {
673                                                << 604       curSum = 0;
674     for(G4int jj=0; jj<10; ++jj) {             << 605       aSimp  = 4;   
675       curQ2 = Q2l+(jj + 0.5)*ddQ2;             << 606 
676       if(curQ2 >= Q2max) { break; }            << 607       for(jj = 0; jj<20; jj++)
677       curSec = HadrNucDifferCrSec(A, curQ2);   << 608   {
678       curSum += curSec;                        << 609     curQ2 = Q2l+jj*ddQ2;
679     }                                          << 610 
680     G4double del = (curQ2 >= Q2max) ? Q2max -  << 611     curSec  = HadrNucDifferCrSec(Nucleus, curQ2);
681     Q2l    += del;                             << 612     curSum += curSec*aSimp;
682     curSum *= del*0.1;                         << 613 
683     totSum += curSum;                          << 614     if(aSimp > 3) aSimp = 2;
684     fLineF[ii] = totSum;                       << 615     else          aSimp = 4;
685     if (verboseLevel>2) {                      << 616 
686       G4cout<<ii << ". FillFq2: A= " << A << " << 617     if(jj == 0 && verboseLevel>2)
687       <<dQ2<<" Tot= "<<totSum << " dTot " <<cu << 618       G4cout<<"  Q2  "<<curQ2<<"  AIm  "<<aAIm<<"  DIm  "<<aDIm
688       <<" curSec= " <<curSec<<G4endl;          << 619       <<"  Diff  "<<curSec<<"  totSum  "<<totSum<<G4endl;
689     }                                          << 620   }
690     if(totSum*1.e-4 > curSum || Q2l >= Q2max)  << 621 
691   }                                            << 622       Q2l    += dQ2;
692   ii = std::min(ii, ONQ2-2);                   << 623       curSum *= ddQ2/2.3;   //  $$$$$$$$$$$$$$$$$$$$$$$
693   curQ2 = Q2l;                                 << 624       totSum += curSum;
694   G4double xx = R1*(Q2max - curQ2);            << 625 
695   if(xx > 0.0) {                               << 626       LineF[ii] = totSum;
696     xx = (xx > 20.) ? 0.0 : G4Exp(-xx);        << 627   
697     curSec = HadrNucDifferCrSec(A, curQ2);     << 628       if (verboseLevel>2)
698     totSum += curSec*(1.0 - xx)/R1;            << 629   G4cout<<"  GetHeavy: Q2  dQ2  totSum  "<<Q2l<<"  "<<dQ2<<"  "<<totSum
699   }                                            << 630         <<"  curSec  "
700   fLineF[ii + 1] = totSum;                     << 631         <<curSec<<"  totSum  "<< totSum<<"  DTot "
701   if (verboseLevel>1) {                        << 632         <<curSum<<G4endl;
702     G4cout << "### FillFq2 done curQ2= " << cu << 633     }      
703            << " sumG= " << fLineF[ONQ2-2] << " << 634   return totSum;
704      << " Nbins= " << ii + 1 << G4endl;        << 
705   }                                            << 
706   return ii + 2;                               << 
707 }                                                 635 }
708                                                   636 
709 //////////////////////////////////////////////    637 ////////////////////////////////////////////////////////////////////////
                                                   >> 638 //
                                                   >> 639 //
710                                                   640 
711 G4double G4ElasticHadrNucleusHE::GetLightFq2(G << 641 G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus, 
                                                   >> 642                                              G4double Q2)
712 {                                                 643 {
713   // Scattering off proton                     << 644   // Scattering of proton
714   if(Z == 1)                                      645   if(Z == 1) 
715   {                                               646   {
716     G4double SqrQ2  = std::sqrt(Q2);              647     G4double SqrQ2  = std::sqrt(Q2);
717     G4double valueConstU = 2.*(hMass2 + proton << 648     G4double ConstU = 2.*(hMass2 + protonM2) - Q2;
718                                                   649 
719     G4double y = (1.-Coeff1-Coeff0)/HadrSlope* << 650     G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-std::exp(-HadrSlope*Q2))
720       + Coeff0*(1.-G4Exp(-Slope0*Q2))          << 651       + Coeff0*(1.-std::exp(-Slope0*Q2))
721       + Coeff2/Slope2*G4Exp(Slope2*valueConstU << 652       + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1.)
722       + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1 << 653       + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
723                                                   654 
724     return y;                                     655     return y;
725   }                                               656   }
726                                                   657 
727   // The preparing of probability function        658   // The preparing of probability function  
728                                                   659 
729   G4double prec = A > 208  ?  1.0e-7 : 1.0e-6; << 660   G4double prec = Nucleus > 208  ?  1.0e-7 : 1.0e-6;
730                                                   661 
731   G4double    Stot     = HadrTot*MbToGeV2;        662   G4double    Stot     = HadrTot*MbToGeV2;     //  Gev^-2
732   G4double    Bhad     = HadrSlope;         //    663   G4double    Bhad     = HadrSlope;         //  GeV^-2
733   G4double    Asq      = 1+HadrReIm*HadrReIm;     664   G4double    Asq      = 1+HadrReIm*HadrReIm;
734   G4double    Rho2     = std::sqrt(Asq);          665   G4double    Rho2     = std::sqrt(Asq);
735                                                   666 
736   if(verboseLevel >1) {                        << 667 //  if(verboseLevel >1)
737     G4cout<<" Fq2 Before for i Tot B Im "<<Had    668     G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<"  "<<HadrSlope<<"  "
738       <<HadrReIm<<G4endl;                         669       <<HadrReIm<<G4endl;
739   }                                            << 670 
740   if(verboseLevel > 1) {                          671   if(verboseLevel > 1) {
741     G4cout << "GetFq2: Stot= " << Stot << " Bh    672     G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad 
742            <<"  Im "<<HadrReIm                    673            <<"  Im "<<HadrReIm 
743            << " Asq= " << Asq << G4endl;          674            << " Asq= " << Asq << G4endl;
744     G4cout << "R1= " << R1 << " R2= " << R2 <<    675     G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl;
745   }                                               676   }
746   G4double    R12      = R1*R1;                   677   G4double    R12      = R1*R1;
747   G4double    R22      = R2*R2;                   678   G4double    R22      = R2*R2;
748   G4double    R12B     = R12+2*Bhad;              679   G4double    R12B     = R12+2*Bhad;
749   G4double    R22B     = R22+2*Bhad;              680   G4double    R22B     = R22+2*Bhad;
750                                                   681 
751   G4double    Norm     = (R12*R1-Pnucl*R22*R2)    682   G4double    Norm     = (R12*R1-Pnucl*R22*R2); // HP->Aeff;
752                                                   683 
753   G4double    R13      = R12*R1/R12B;             684   G4double    R13      = R12*R1/R12B;
754   G4double    R23      = Pnucl*R22*R2/R22B;       685   G4double    R23      = Pnucl*R22*R2/R22B;
755   G4double    Unucl    = Stot/twopi*R13/Norm;  << 686   G4double    Unucl    = Stot/twopi/Norm*R13;
756   G4double    UnucRho2 = -Unucl*Rho2;             687   G4double    UnucRho2 = -Unucl*Rho2;
757                                                   688 
758   G4double    FiH      = std::asin(HadrReIm/Rh    689   G4double    FiH      = std::asin(HadrReIm/Rho2);
759   G4double    NN2      = R23/R13;                 690   G4double    NN2      = R23/R13;
760                                                   691 
761   if(verboseLevel > 2) {                       << 692   if(verboseLevel > 2) 
762     G4cout << "UnucRho2= " << UnucRho2 << " Fi    693     G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2 
763      << " Norm= " << Norm << G4endl;              694      << " Norm= " << Norm << G4endl;
764   }                                            << 
765   G4double    Prod0 = 0.;                      << 
766   G4double    N1    = -1.0;                    << 
767                                                   695 
768   for(G4int i1 = 1; i1<= A; ++i1) ////++++++++ << 696   G4double    dddd;
                                                   >> 697  
                                                   >> 698   G4double    Prod0    = 0;
                                                   >> 699   G4double    N1       = -1.0;
                                                   >> 700   //G4double    Tot0     = 0;
                                                   >> 701   G4double    exp1;
                                                   >> 702 
                                                   >> 703   G4double    Prod3 ;
                                                   >> 704   G4double    exp2  ;
                                                   >> 705   G4double    N4, N5, N2, Prod1, Prod2;
                                                   >> 706   G4int    i1, i2, m1, m2;
                                                   >> 707 
                                                   >> 708   for(i1 = 1; i1<= Nucleus; i1++) ////++++++++++  i1
769     {                                             709     {
770       N1 *= (-Unucl*Rho2*(A-i1+1)/(G4double)i1 << 710       N1    = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
771       G4double Prod1 = 0.;                     << 711       Prod1 = 0;
772       G4double N2    = -1.;                    << 712       //Tot0  = 0;
                                                   >> 713       N2    = -1;
773                                                   714 
774       for(G4int i2 = 1; i2<=A; ++i2) ////+++++ << 715       for(i2 = 1; i2<=Nucleus; i2++) ////+++++++++ i2
775         {                                         716         {
776           N2 *= (-Unucl*Rho2*(A-i2+1)/(G4doubl << 717           N2    = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
777           G4double Prod2 = 0;                  << 718           Prod2 = 0; 
778           G4double N5    = -1/NN2;             << 719           N5    = -1/NN2;
779     for(G4int j2=0; j2<= i2; ++j2) ////+++++++ << 720     for(m2=0; m2<= i2; m2++) ////+++++++++ m2
780             {                                     721             {
781               G4double Prod3 = 0;              << 722               Prod3 = 0;
782               G4double exp2  = 1./((G4double)j << 723               exp2  = 1/(m2/R22B+(i2-m2)/R12B);
783               N5 *= (-NN2);                    << 724               N5    = -N5*NN2;
784               G4double N4 = -1./NN2;           << 725               N4    = -1/NN2;
785         for(G4int j1=0; j1<=i1; ++j1) ////++++ << 726         for(m1=0; m1<=i1; m1++) ////++++++++ m1
786     {                                             727     {
787       G4double exp1  = 1./((G4double)j1/R22B+( << 728       exp1  = 1/(m1/R22B+(i1-m1)/R12B);
788       G4double dddd  = 0.25*(exp1+exp2);       << 729       dddd  = exp1+exp2;
789       N4    *= (-NN2);                         << 730       N4    = -N4*NN2;
790       Prod3 +=                                 << 731       Prod3 = Prod3+N4*exp1*exp2*
791                     N4*exp1*exp2*(1.-G4Exp(-Q2 << 732         (1-std::exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][m1];
792     }                                   // j1  << 733                }                                   // m1
793         Prod2 += Prod3*N5*GetBinomCof(i2,j2);  << 734         Prod2 = Prod2 +Prod3*N5*SetBinom[i2][m2];
794       }                                      / << 735       }                                      // m2
795     Prod1 += Prod2*N2*std::cos(FiH*(i1-i2));   << 736     Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
796                                                   737 
797     if (std::abs(Prod2*N2/Prod1)<prec) break;  << 738     if (std::fabs(Prod2*N2/Prod1)<prec) break;
798         }                                         739         }                                         // i2
799       Prod0 += Prod1*N1;                       << 740       Prod0   = Prod0 + Prod1*N1;
800       if(std::abs(N1*Prod1/Prod0) < prec) brea << 741       if(std::fabs(N1*Prod1/Prod0) < prec) break;
                                                   >> 742 
801     }                                             743     }                                           // i1
802                                                   744 
803   const G4double fact = 0.25*CLHEP::pi/MbToGeV << 745   Prod0 *= 0.25*pi/MbToGeV2;  //  This is in mb
804   Prod0 *= fact;  //  This is in mb            << 
805                                                   746 
806   if(verboseLevel>1) {                         << 747   if(verboseLevel>1) 
807     G4cout << "GetLightFq2 Z= " << Z << " A= " << 748     G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus 
808      <<" Q2= " << Q2 << " Res= " << Prod0 << G    749      <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl;
809   }                                            << 
810   return Prod0;                                   750   return Prod0;
811 }                                                 751 }
812                                                << 752 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
813 ////////////////////////////////////////////// << 753 G4double G4ElasticHadrNucleusHE::
814                                                << 754    HadrNucDifferCrSec(G4int Nucleus, G4double aQ2)
815 G4double                                       << 
816 G4ElasticHadrNucleusHE::HadrNucDifferCrSec(G4i << 
817 {                                                 755 {
818   //   ------ All external kinematical variabl << 756 //   ------ All external kinematical variables are in MeV -------
819   //            ------ but internal in GeV !!! << 757 //            ------ but internal in GeV !!!!  ------
                                                   >> 758 
                                                   >> 759   G4double    theQ2 = aQ2;   ///GeV/GeV;  
820                                                   760 
821   // Scattering of proton                         761   // Scattering of proton
822   if(A == 1)                                   << 762   if(Nucleus == 1) 
823   {                                               763   {
824     G4double SqrQ2  = std::sqrt(aQ2);             764     G4double SqrQ2  = std::sqrt(aQ2);
825     G4double valueConstU = hMass2 + protonM2-2 << 765     G4double ConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2;
826                                                << 766 
827     BoundaryTL[0] = Q2max;                     << 767     G4double MaxT = 4*MomentumCM*MomentumCM;
828     BoundaryTL[1] = Q2max;                     << 768 
829     BoundaryTL[3] = Q2max;                     << 769      BoundaryTL[0] = MaxT;
830     BoundaryTL[4] = Q2max;                     << 770      BoundaryTL[1] = MaxT;
831     BoundaryTL[5] = Q2max;                     << 771      BoundaryTL[3] = MaxT;
832                                                << 772      BoundaryTL[4] = MaxT;
833     G4double dSigPodT = HadrTot*HadrTot*(1+Had << 773      BoundaryTL[5] = MaxT;
834       ( Coeff1*G4Exp(-Slope1*SqrQ2)+           << 774 
835   Coeff2*G4Exp( Slope2*(valueConstU)+aQ2)+     << 775     G4double dSigPodT;
836   (1-Coeff1-Coeff0)*G4Exp(-HadrSlope*aQ2)+     << 776 
837   Coeff0*G4Exp(-Slope0*aQ2) )*2.568/(16*pi);   << 777     dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)*
                                                   >> 778                  (
                                                   >> 779                   Coeff1*std::exp(-Slope1*SqrQ2)+
                                                   >> 780                   Coeff2*std::exp( Slope2*(ConstU)+aQ2)+
                                                   >> 781                   (1-Coeff1-Coeff0)*std::exp(-HadrSlope*aQ2)+
                                                   >> 782                  +Coeff0*std::exp(-Slope0*aQ2)
                                                   >> 783 //                +0.1*(1-std::fabs(CosTh))
                                                   >> 784                   )/16/3.1416*2.568;
838                                                   785 
839     return dSigPodT;                              786     return dSigPodT;
840   }                                               787   }
841                                                   788 
842   G4double    Stot     = HadrTot*MbToGeV2;     << 789     G4double    Stot     = HadrTot*MbToGeV2; 
843   G4double    Bhad     = HadrSlope;            << 790     G4double    Bhad     = HadrSlope; 
844   G4double    Asq      = 1+HadrReIm*HadrReIm;  << 791     G4double    Asq      = 1+HadrReIm*HadrReIm;
845   G4double    Rho2     = std::sqrt(Asq);       << 792     G4double    Rho2     = std::sqrt(Asq);
846   G4double    R12      = R1*R1;                << 793     G4double    Pnuclp   = 0.001;
847   G4double    R22      = R2*R2;                << 794                 Pnuclp   = Pnucl;
848   G4double    R12B     = R12+2*Bhad;           << 795     G4double    R12      = R1*R1;
849   G4double    R22B     = R22+2*Bhad;           << 796     G4double    R22      = R2*R2;
850   G4double    R12Ap    = R12+20;               << 797     G4double    R12B     = R12+2*Bhad;
851   G4double    R22Ap    = R22+20;               << 798     G4double    R22B     = R22+2*Bhad;
852   G4double    R13Ap    = R12*R1/R12Ap;         << 799     G4double    R12Ap    = R12+20;
853   G4double    R23Ap    = R22*R2*Pnucl/R22Ap;   << 800     G4double    R22Ap    = R22+20;
854   G4double    R23dR13  = R23Ap/R13Ap;          << 801     G4double    R13Ap    = R12*R1/R12Ap;
855   G4double    R12Apd   = 2/R12Ap;              << 802     G4double    R23Ap    = R22*R2/R22Ap*Pnuclp;
856   G4double    R22Apd   = 2/R22Ap;              << 803     G4double    R23dR13  = R23Ap/R13Ap;
857   G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);  << 804     G4double    R12Apd   = 2/R12Ap;
858                                                << 805     G4double    R22Apd   = 2/R22Ap;
859   G4double DDSec1p  = (DDSect2+DDSect3*G4Log(0 << 806     G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
860   G4double DDSec2p  = (DDSect2+DDSect3*G4Log(0 << 807 
861                              std::sqrt((R12+R2 << 808     G4double DDSec1p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R1/4));
862   G4double DDSec3p  = (DDSect2+DDSect3*G4Log(0 << 809     G4double DDSec2p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/
                                                   >> 810                              std::sqrt((R12+R22)/2)/4));
                                                   >> 811     G4double DDSec3p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R2/4));
                                                   >> 812 
                                                   >> 813     G4double    Norm     = (R12*R1-Pnucl*R22*R2)*Aeff;
                                                   >> 814     G4double    Normp    = (R12*R1-Pnuclp*R22*R2)*Aeff;
                                                   >> 815     G4double    R13      = R12*R1/R12B;
                                                   >> 816     G4double    R23      = Pnucl*R22*R2/R22B;
                                                   >> 817     G4double    Unucl    = Stot/twopi/Norm*R13;
                                                   >> 818     G4double    UnuclScr = Stot/twopi/Normp*R13Ap;
                                                   >> 819     G4double    SinFi    = HadrReIm/Rho2;
                                                   >> 820     G4double    FiH      = std::asin(SinFi);
                                                   >> 821     G4double    N        = -1;
                                                   >> 822     G4double    N2       = R23/R13;
                                                   >> 823 
                                                   >> 824     G4double    ImElasticAmpl0  = 0;
                                                   >> 825     G4double    ReElasticAmpl0  = 0;
                                                   >> 826 
                                                   >> 827     G4double    exp1;
                                                   >> 828     G4double    N4;
                                                   >> 829     G4double    Prod1, Tot1=0, medTot, DTot1, DmedTot;
                                                   >> 830     G4int       i;
863                                                   831 
864   G4double    Norm     = (R12*R1-Pnucl*R22*R2) << 832     for( i=1; i<=Nucleus; i++)
865   G4double    R13      = R12*R1/R12B;          << 833     {
866   G4double    R23      = Pnucl*R22*R2/R22B;    << 834       N       = -N*Unucl*(Nucleus-i+1)/i*Rho2;
867   G4double    Unucl    = Stot/(twopi*Norm)*R13 << 835       N4      = 1;
868   G4double    UnuclScr = Stot/(twopi*Norm)*R13 << 836       Prod1   = std::exp(-theQ2/i*R12B/4)/i*R12B;
869   G4double    SinFi    = HadrReIm/Rho2;        << 837       medTot  = R12B/i;
870   G4double    FiH      = std::asin(SinFi);     << 838 
871   G4double    N        = -1;                   << 839        for(G4int l=1; l<=i; l++)
872   G4double    N2       = R23/R13;              << 840        {
873                                                << 841          exp1    = l/R22B+(i-l)/R12B;
874   G4double ImElasticAmpl0 = 0;                 << 842          N4      = -N4*(i-l+1)/l*N2;
875   G4double ReElasticAmpl0 = 0;                 << 843          Prod1   = Prod1+N4/exp1*std::exp(-theQ2/exp1/4);
876   G4double exp1;                               << 844          medTot  = medTot+N4/exp1;
877                                                << 845        }  // end l
878   for(G4int i=1; i<=A; ++i) {                  << 846 
879     N  *= (-Unucl*Rho2*(A-i+1)/(G4double)i);   << 847       ReElasticAmpl0  = ReElasticAmpl0+Prod1*N*std::sin(FiH*i);
880     G4double N4 = 1;                           << 848       ImElasticAmpl0  = ImElasticAmpl0+Prod1*N*std::cos(FiH*i);
881     G4double medTot = R12B/(G4double)i;        << 849       Tot1            = Tot1+medTot*N*std::cos(FiH*i);
882     G4double Prod1  = G4Exp(-aQ2*R12B/(G4doubl << 850       if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001) break;
883                                                << 851     }      // i
884     for(G4int l=1; l<=i; ++l) {                << 852 
885       exp1 = l/R22B+(i-l)/R12B;                << 853     ImElasticAmpl0 = ImElasticAmpl0*pi/2.568;   // The amplitude in mB
886       N4 *= (-N2*(i-l+1)/(G4double)l);         << 854     ReElasticAmpl0 = ReElasticAmpl0*pi/2.568;   // The amplitude in mB
887       G4double expn4 = N4/exp1;                << 855     Tot1           = Tot1*twopi/2.568;
888       Prod1  += expn4*G4Exp(-aQ2/(exp1*4));    << 856 
889       medTot += expn4;                         << 857     G4double C1 = R13Ap*R13Ap/2*DDSec1p;
890     }  // end l                                << 858     G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p;
891                                                << 859     G4double C3 = R23Ap*R23Ap/2*DDSec3p;
892     G4double dcos = N*std::cos(FiH*i);         << 860 
893     ReElasticAmpl0  += Prod1*N*std::sin(FiH*i) << 861     G4double N1p  = 1;
894     ImElasticAmpl0  += Prod1*dcos;             << 862 
895     if(std::abs(Prod1*N/ImElasticAmpl0) < 0.00 << 863     G4double Din1 = 0.5;     
896   }      // i                                  << 864 
897                                                << 865     Din1  = 0.5*(C1*std::exp(-theQ2/8*R12Ap)/2*R12Ap-
898   static const G4double pi25 = CLHEP::pi/2.568 << 866                  C2/R12ApdR22Ap*std::exp(-theQ2/4/R12ApdR22Ap)+
899   ImElasticAmpl0 *= pi25;   // The amplitude i << 867                  C3*R22Ap/2*std::exp(-theQ2/8*R22Ap));
900   ReElasticAmpl0 *= pi25;   // The amplitude i << 868 
901                                                << 869     DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
902   G4double C1 = R13Ap*R13Ap*0.5*DDSec1p;       << 870 
903   G4double C2 = 2*R23Ap*R13Ap*0.5*DDSec2p;     << 871     G4double exp1p;
904   G4double C3 = R23Ap*R23Ap*0.5*DDSec3p;       << 872     G4double exp2p;
905                                                << 873     G4double exp3p;
906   G4double N1p  = 1;                           << 874     G4double N2p;
907   G4double Din1 = 0.5*(C1*G4Exp(-aQ2/8*R12Ap)/ << 875     G4double Din2, BinCoeff;
908            C2/R12ApdR22Ap*G4Exp(-aQ2/(4*R12Apd << 876 
909            C3*R22Ap/2*G4Exp(-aQ2/8*R22Ap));    << 877     BinCoeff = 1;
910                                                << 878 
911   G4double DTot1 = 0.5*(C1*0.5*R12Ap-C2/R12Apd << 879     for( i = 1; i<= Nucleus-2; i++)
912                                                << 880     {
913   for(G4int i=1; i<= A-2; ++i) {               << 881       N1p     = -N1p*UnuclScr*(Nucleus-i-1)/i*Rho2;
914     N1p *= (-UnuclScr*Rho2*(A-i-1)/(G4double)i << 882       N2p     = 1;
915     G4double N2p  = 1;                         << 883       Din2    = 0;
916     G4double Din2 = 0;                         << 884       DmedTot = 0;
917     G4double DmedTot = 0;                      << 885         for(G4int l = 0; l<=i; l++) 
918     G4double BinCoeff = 1.0;                   << 886         {
919     for(G4int l=0; l<=i; ++l) {                << 887           if(l == 0)      BinCoeff = 1;
920       if(l > 0) { BinCoeff *= (i-l+1)/(G4doubl << 888           else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
921                                                << 889 
922       exp1  = l/R22B+(i-l)/R12B;               << 890           exp1  = l/R22B+(i-l)/R12B;
923       G4double exp1p = exp1+R12Apd;            << 891           exp1p = exp1+R12Apd;
924       G4double exp2p = exp1+R12ApdR22Ap;       << 892           exp2p = exp1+R12ApdR22Ap;
925       G4double exp3p = exp1+R22Apd;            << 893           exp3p = exp1+R22Apd;
926                                                << 894 
927       Din2 += N2p*BinCoeff*(C1/exp1p*G4Exp(-aQ << 895           Din2  = Din2 + N2p*BinCoeff*
928           C2/exp2p*G4Exp(-aQ2/(4*exp2p))+      << 896             (C1/exp1p*std::exp(-theQ2/4/exp1p)-
929           C3/exp3p*G4Exp(-aQ2/(4*exp3p)));     << 897                    C2/exp2p*std::exp(-theQ2/4/exp2p)+
930                                                << 898                    C3/exp3p*std::exp(-theQ2/4/exp3p));
931       DmedTot += N2p*BinCoeff*(C1/exp1p-C2/exp << 
932                                                << 
933       N2p *= -R23dR13;                         << 
934     }     // l                                 << 
935                                                << 
936     G4double dcos = N1p*std::cos(FiH*i)/(G4dou << 
937     Din1  += Din2*dcos;                        << 
938     DTot1 += DmedTot*dcos;                     << 
939                                                << 
940     if(std::abs(Din2*N1p/Din1) < 0.000001) bre << 
941   }           //  i                            << 
942   G4double gg = (G4double)(A*(A-1)*4)/(Norm*No << 
943                                                   899 
944   Din1  *= (-gg);                              << 900     DmedTot = DmedTot + N2p*BinCoeff*
945   DTot1 *= 5*gg;                               << 901       (C1/exp1p-C2/exp2p+C3/exp3p);
946                                                   902 
947   //  ----------------  dSigma/d|-t|,  mb/(GeV << 903     N2p   = -N2p*R23dR13;
                                                   >> 904   }     // l
948                                                   905 
949   G4double DiffCrSec2 = (ReElasticAmpl0*ReElas << 906   Din1  = Din1+Din2*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
950        (ImElasticAmpl0+Din1)*                  << 907   DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
951        (ImElasticAmpl0+Din1))/twopi;           << 908  
                                                   >> 909   if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
                                                   >> 910     }           //  i
                                                   >> 911 
                                                   >> 912     Din1 = -Din1*Nucleus*(Nucleus-1)
                                                   >> 913                  /2/pi/Normp/2/pi/Normp*16*pi*pi;
952                                                   914 
953   Dtot11 = DTot1;                              << 915     DTot1 = DTot1*Nucleus*(Nucleus-1)
954   aAIm   = ImElasticAmpl0;                     << 916                  /2/pi/Normp/2/pi/Normp*16*pi*pi;
955   aDIm   = Din1;                               << 
956                                                   917 
957   return DiffCrSec2;  //  dSig/d|-t|,  mb/(GeV << 918     DTot1 *= 5;   //  $$$$$$$$$$$$$$$$$$$$$$$$ 
958 }                                              << 919 //     Din1  *= 0.2;    //   %%%%%%%%%%%%%%%%%%%%%%%   proton
                                                   >> 920 //     Din1 *= 0.05;    //   %%%%%%%%%%%%%%%%%%%%%%%  pi+
                                                   >> 921 //  ----------------  dSigma/d|-t|,  mb/(GeV/c)^-2  -----------------
                                                   >> 922 
                                                   >> 923     G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
                                                   >> 924                            (ImElasticAmpl0+Din1)*
                                                   >> 925                            (ImElasticAmpl0+Din1))*2/4/pi;
                                                   >> 926 
                                                   >> 927     Tot1   = Tot1-DTot1;
                                                   >> 928      //  Tott1  = Tot1*1.0;
                                                   >> 929     Dtot11 = DTot1;
                                                   >> 930     aAIm   = ImElasticAmpl0;
                                                   >> 931     aDIm   = Din1;
                                                   >> 932 
                                                   >> 933     return DiffCrSec2*1.0;  //  dSig/d|-t|,  mb/(GeV/c)^-2
                                                   >> 934 }   // function
                                                   >> 935 //  ##############################################
959                                                   936 
960 //////////////////////////////////////////////    937 ////////////////////////////////////////////////////////////////
                                                   >> 938 //
                                                   >> 939 //
961                                                   940 
962 void G4ElasticHadrNucleusHE::DefineHadronValue << 941 void  G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
963 {                                                 942 {
                                                   >> 943   HadrEnergy = std::sqrt(hMass2 + hLabMomentum2);
                                                   >> 944 
964   G4double sHadr = 2.*HadrEnergy*protonM+proto    945   G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2;
965   G4double sqrS  = std::sqrt(sHadr);              946   G4double sqrS  = std::sqrt(sHadr);
                                                   >> 947   G4double Ecm   = 0.5*(sHadr-hMass2+protonM2)/sqrS;
                                                   >> 948   MomentumCM     = std::sqrt(Ecm*Ecm-protonM2);
966                                                   949   
967   if(verboseLevel>2) {                         << 950   if(verboseLevel>2)
968     G4cout << "GetHadrValues: Z= " << Z << " i << 951     G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron 
969      << " E(GeV)= " << HadrEnergy << " sqrS= "    952      << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS
970      << " plab= " << hLabMomentum                 953      << " plab= " << hLabMomentum   
971      <<"  E - m  "<<HadrEnergy - hMass<< G4end    954      <<"  E - m  "<<HadrEnergy - hMass<< G4endl;
972   }                                            << 955 
973   G4double TotN = 0.0;                            956   G4double TotN = 0.0;
974   G4double logE = G4Log(HadrEnergy);           << 957   G4double logE = std::log(HadrEnergy);
975   G4double logS = G4Log(sHadr);                << 958   G4double logS = std::log(sHadr);
976            TotP = 0.0;                            959            TotP = 0.0;
977                                                   960 
978   switch (iHadron) {                           << 961   switch (iHadron)
979   case 0:                  //  proton, neutron << 962     {
980   case 6:                                      << 963     case 0:                  //  proton, neutron
981                                                << 964     case 6:
982     if(hLabMomentum > 10) {                    << 965 
983       TotP = TotN = 7.5*logE - 40.12525 + 103* << 966       if(hLabMomentum > 10)
984                                                << 967   TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); //  mb
985     } else {                                   << 968 
986       // ==================  neutron  ======== << 969       else
987                                                << 970   {
988       if( hLabMomentum > 1.4 ) {               << 971 // ==================  neutron  ================
989   TotN = 33.3+15.2*(hLabMomentum2-1.35)/       << 972 
990     (G4Exp(G4Log(hLabMomentum)*2.37)+0.95);    << 973 ////    if(iHadrCode == 2112) 
                                                   >> 974 
                                                   >> 975 
                                                   >> 976     if( hLabMomentum > 1.4 )
                                                   >> 977       TotN = 33.3+15.2*(hLabMomentum2-1.35)/
                                                   >> 978         (std::pow(hLabMomentum,2.37)+0.95);
991                                                   979     
992       } else if(hLabMomentum > 0.8) {          << 980     else if(hLabMomentum > 0.8)
993   G4double A0 = logE + 0.0513;                 << 981       {
994   TotN = 33.0 + 25.5*A0*A0;                    << 982         G4double A0 = logE + 0.0513;
995       } else {                                 << 983         TotN = 33.0 + 25.5*A0*A0;  
996   G4double A0 = logE - 0.2634;  // log(1.3)    << 984       }
997   TotN = 33.0 + 30.*A0*A0*A0*A0;               << 985     else 
998       }                                        << 986       {
999       //  =================  proton  ========= << 987         G4double A0 = logE - 0.2634;  // log(1.3)
                                                   >> 988         TotN = 33.0 + 30.*A0*A0*A0*A0;
                                                   >> 989       }
                                                   >> 990 //  =================  proton  ===============
                                                   >> 991 //       else if(iHadrCode == 2212) 
                                                   >> 992     {
                                                   >> 993       if(hLabMomentum >= 1.05)
                                                   >> 994               {
                                                   >> 995     TotP = 39.0+75.*(hLabMomentum-1.2)/
                                                   >> 996       (hLabMomentum2*hLabMomentum+0.15);
                                                   >> 997               }
                                                   >> 998 
                                                   >> 999       else if(hLabMomentum >= 0.7)
                                                   >> 1000         {
                                                   >> 1001      G4double A0 = logE + 0.3147;
                                                   >> 1002      TotP = 23.0 + 40.*A0*A0;
                                                   >> 1003         }
                                                   >> 1004       else 
                                                   >> 1005               {
                                                   >> 1006     TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5);
                                                   >> 1007         }
                                                   >> 1008     }
                                                   >> 1009   }
1000                                                  1010 
1001       if(hLabMomentum >= 1.05) {              << 1011 //      HadrTot = 0.5*(82*TotP+126*TotN)/104;  //  $$$$$$$$$$$$$$$$$$
1002   TotP = 39.0+75.*(hLabMomentum-1.2)/(hLabMom << 1012       HadrTot = 0.5*(TotP+TotN);
1003       } else if(hLabMomentum >= 0.7) {        << 1013 //  ...................................................
1004   G4double A0 = logE + 0.3147;                << 1014       //  Proton slope
1005   TotP = 23.0 + 40.*A0*A0;                    << 1015       if(hLabMomentum >= 2.)       HadrSlope = 5.44 + 0.88*logS;
1006       } else {                                << 1016 
1007   TotP = 23.+50.*G4Exp(G4Log(G4Log(0.73/hLabM << 1017       else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37;
1008       }                                       << 1018 
1009     }                                         << 1019       else                         HadrSlope = 1.5;
1010     HadrTot = 0.5*(TotP+TotN);                << 1020 
1011     //  ..................................... << 1021 //  ...................................................
1012     //  Proton slope                          << 1022       if(hLabMomentum >= 1.2)
1013     if(hLabMomentum >= 2.)       { HadrSlope  << 1023    HadrReIm  = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18);
1014     else if(hLabMomentum >= 0.5) { HadrSlope  << 1024        
1015     else                         { HadrSlope  << 1025       else if(hLabMomentum >= 0.6) 
1016                                               << 1026   HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/
1017     //  ..................................... << 1027     (std::pow(3*hLabMomentum,2.2)+1);     
1018     if(hLabMomentum >= 1.2) {                 << 1028 
1019       HadrReIm  = 0.13*(logS - 5.8579332)*G4E << 1029       else 
1020     } else if(hLabMomentum >= 0.6) {          << 1030   HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2);
1021       HadrReIm = -75.5*(G4Exp(G4Log(hLabMomen << 1031 //  ...................................................
1022   (G4Exp(G4Log(3*hLabMomentum)*2.2)+1);       << 1032       DDSect2   = 2.2;                              //mb*GeV-2
1023     } else {                                  << 1033       DDSect3   = 0.6;                               //mb*GeV-2
1024       HadrReIm = 15.5*hLabMomentum/(27*hLabMo << 1034       //  ================== lambda  ==================
1025     }                                         << 1035       if( iHadrCode == 3122)
1026     //  ..................................... << 1036   {
1027     DDSect2   = 2.2;                          << 1037     HadrTot   *= 0.88;
1028     DDSect3   = 0.6;                          << 1038     HadrSlope *=0.85;
1029     //  ================== lambda  ========== << 1039   }
1030     if( iHadrCode == 3122) {                  << 
1031       HadrTot   *= 0.88;                      << 
1032       HadrSlope *=0.85;                       << 
1033       //  ================== sigma +  =======    1040       //  ================== sigma +  ==================
1034     } else if( iHadrCode == 3222) {           << 1041       else if( iHadrCode == 3222)
1035       HadrTot   *=0.81;                       << 1042   {
1036       HadrSlope *=0.85;                       << 1043     HadrTot   *=0.81;
                                                   >> 1044     HadrSlope *=0.85;
                                                   >> 1045   }
1037       //  ================== sigma 0,-  =====    1046       //  ================== sigma 0,-  ==================
1038     } else if(iHadrCode == 3112 || iHadrCode  << 1047       else if(iHadrCode == 3112 || iHadrCode == 3212 )
1039       HadrTot   *=0.88;                       << 1048   {
1040       HadrSlope *=0.85;                       << 1049     HadrTot   *=0.88;
                                                   >> 1050     HadrSlope *=0.85;
                                                   >> 1051   }
1041       //  ===================  xi  ==========    1052       //  ===================  xi  =================
1042     } else if( iHadrCode == 3312 || iHadrCode << 1053       else if( iHadrCode == 3312 || iHadrCode == 3322 )
1043       HadrTot   *=0.77;                       << 1054   {
1044       HadrSlope *=0.75;                       << 1055     HadrTot   *=0.77;
                                                   >> 1056     HadrSlope *=0.75;
                                                   >> 1057   }
1045       //  =================  omega  =========    1058       //  =================  omega  =================
1046     } else if( iHadrCode == 3334) {           << 1059       else if( iHadrCode == 3334)
1047       HadrTot   *=0.78;                       << 1060   {
1048       HadrSlope *=0.7;                        << 1061     HadrTot   *=0.78;
1049     }                                         << 1062     HadrSlope *=0.7;
1050     break;                                    << 1063   }
1051     //  ===================================== << 1064 
1052   case 1:              //   antiproton        << 1065       break;
1053   case 7:              //   antineutron       << 1066 //  ===========================================================
1054                                               << 1067     case 1:              //   antiproton
1055     HadrTot   = 5.2+5.2*logE + 123.2/sqrS;    << 1068     case 7:              //   antineutron
1056     HadrSlope = 8.32+0.57*logS;               << 1069 
1057                                               << 1070       HadrTot   = 5.2+5.2*logE + 123.2/sqrS;     //  mb
1058     if( HadrEnergy < 1000 ) {                 << 1071       HadrSlope = 8.32+0.57*logS;                //(GeV/c)^-2
1059       HadrReIm  = 0.06*(sqrS-2.236)*(sqrS-14. << 1072 
1060     } else {                                  << 1073       if( HadrEnergy < 1000 )
1061       HadrReIm  = 0.6*(logS - 5.8579332)*G4Ex << 1074   HadrReIm  = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
1062     }                                         << 1075       else
1063     DDSect2   = 11;                           << 1076   HadrReIm  = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25);
1064     DDSect3   = 3;                            << 1077 
1065     //  ================== lambda  ========== << 1078       DDSect2   = 11;                            //mb*(GeV/c)^-2
1066     if( iHadrCode == -3122) {                 << 1079       DDSect3   = 3;                             //mb*(GeV/c)^-2
1067       HadrTot   *= 0.88;                      << 1080       //  ================== lambda  ==================
1068       HadrSlope *=0.85;                       << 1081       if( iHadrCode == -3122)
                                                   >> 1082   {
                                                   >> 1083     HadrTot   *= 0.88;
                                                   >> 1084     HadrSlope *=0.85;
                                                   >> 1085   }
1069       //  ================== sigma +  =======    1086       //  ================== sigma +  ==================
1070     } else if( iHadrCode == -3222) {          << 1087       else if( iHadrCode == -3222)
1071       HadrTot   *=0.81;                       << 1088   {
1072       HadrSlope *=0.85;                       << 1089     HadrTot   *=0.81;
                                                   >> 1090     HadrSlope *=0.85;
                                                   >> 1091   }
1073       //  ================== sigma 0,-  =====    1092       //  ================== sigma 0,-  ==================
1074     } else if(iHadrCode == -3112 || iHadrCode << 1093       else if(iHadrCode == -3112 || iHadrCode == -3212 )
1075       HadrTot   *=0.88;                       << 1094   {
1076       HadrSlope *=0.85;                       << 1095     HadrTot   *=0.88;
1077     //  ===================  xi  ============ << 1096     HadrSlope *=0.85;
1078     } else if( iHadrCode == -3312 || iHadrCod << 1097   }
1079       HadrTot   *=0.77;                       << 1098       //  ===================  xi  =================
1080       HadrSlope *=0.75;                       << 1099       else if( iHadrCode == -3312 || iHadrCode == -3322 )
                                                   >> 1100   {
                                                   >> 1101     HadrTot   *=0.77;
                                                   >> 1102     HadrSlope *=0.75;
                                                   >> 1103   }
1081       //  =================  omega  =========    1104       //  =================  omega  =================
1082     } else if( iHadrCode == -3334) {          << 1105       else if( iHadrCode == -3334)
1083       HadrTot   *=0.78;                       << 1106   {
1084       HadrSlope *=0.7;                        << 1107     HadrTot   *=0.78;
1085     }                                         << 1108           HadrSlope *=0.7;
1086     break;                                    << 1109   }
1087     //  ------------------------------------- << 
1088   case 2:             //   pi plus, pi minus  << 
1089   case 3:                                     << 
1090                                               << 
1091     if(hLabMomentum >= 3.5) {                 << 
1092       TotP = 10.6+2.*logE + 25.*G4Exp(-logE*0 << 
1093       //  =================================== << 
1094     } else if(hLabMomentum >= 1.15) {         << 
1095       G4double x = (hLabMomentum - 2.55)/0.55 << 
1096       G4double y = (hLabMomentum - 1.47)/0.22 << 
1097       TotP = 3.2*G4Exp(-x*x) + 12.*G4Exp(-y*y << 
1098       //  =================================== << 
1099     } else if(hLabMomentum >= 0.4) {          << 
1100       TotP  = 88*(logE+0.2877)*(logE+0.2877)+ << 
1101     //  ===================================== << 
1102     } else {                                  << 
1103       G4double x = (hLabMomentum - 0.29)/0.08 << 
1104       TotP = 20. + 180.*G4Exp(-x*x);          << 
1105     }                                         << 
1106     //  ------------------------------------- << 
1107                                                  1110 
1108     if(hLabMomentum >= 3.0 ) {                << 1111       break;
1109       TotN = 10.6 + 2.*logE + 30.*G4Exp(-logE << 1112 //  -------------------------------------------
1110     } else if(hLabMomentum >= 1.3) {          << 1113     case 2:             //   pi plus, pi minus
1111       G4double x = (hLabMomentum - 2.1)/0.4;  << 1114     case 3:
1112       G4double y = (hLabMomentum - 1.4)/0.12; << 1115 
1113       TotN = 36.1+0.079 - 4.313*logE + 3.*G4E << 1116       if(hLabMomentum >= 3.5)
1114     } else if(hLabMomentum >= 0.65) {         << 1117   TotP = 10.6+2.*logE + 25.*std::pow(HadrEnergy,-0.43); // mb
1115       G4double x = (hLabMomentum - 0.72)/0.06 << 1118 //  =========================================
1116       G4double y = (hLabMomentum - 1.015)/0.0 << 1119       else if(hLabMomentum >= 1.15)
1117       TotN = 36.1 + 10.*G4Exp(-x*x) + 24*G4Ex << 1120   {
1118     } else if(hLabMomentum >= 0.37) {         << 1121           G4double x = (hLabMomentum - 2.55)/0.55; 
1119       G4double x = G4Log(hLabMomentum/0.48);  << 1122     G4double y = (hLabMomentum - 1.47)/0.225;
1120       TotN = 26. + 110.*x*x;                  << 1123     TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5;
1121     } else {                                  << 1124   }
1122       G4double x = (hLabMomentum - 0.29)/0.07 << 1125 //  =========================================
1123       TotN = 28.0 + 40.*G4Exp(-x*x);          << 1126       else if(hLabMomentum >= 0.4)
1124     }                                         << 1127   {
1125     HadrTot = (TotP+TotN)*0.5;                << 1128   TotP  = 88*(logE+0.2877)*(logE+0.2877)+14.0;
1126     //  ..................................... << 1129         }
1127     HadrSlope = 7.28+0.245*logS;        // Ge << 1130 //  =========================================
1128     HadrReIm  = 0.2*(logS - 4.6051702)*G4Exp( << 1131       else 
1129                                               << 1132   {
1130     DDSect2   = 0.7;                          << 1133     G4double x = (hLabMomentum - 0.29)/0.085;
1131     DDSect3   = 0.27;                         << 1134     TotP = 20. + 180.*std::exp(-x*x);
1132                                               << 1135   }
1133     break;                                    << 1136 //  -------------------------------------------
1134     //  ===================================== << 1137 
1135   case 4:            //  K plus               << 1138       if(hLabMomentum >= 3.0 )
1136                                               << 1139   TotN = 10.6 + 2.*logE + 30.*std::pow(HadrEnergy,-0.43); // mb
1137     HadrTot   = 10.6+1.8*logE + 9.0*G4Exp(-lo << 1140 
1138     if(HadrEnergy>100) { HadrSlope = 15.0; }  << 1141       else if(hLabMomentum >= 1.3) 
1139     else { HadrSlope = 1.0+1.76*logS - 2.84/s << 1142   {
1140                                               << 1143           G4double x = (hLabMomentum - 2.1)/0.4;
1141     HadrReIm  = 0.4*(sHadr-20)*(sHadr-150)*G4 << 1144           G4double y = (hLabMomentum - 1.4)/0.12;
1142     DDSect2   = 0.7;                          << 1145     TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) + 
1143     DDSect3   = 0.21;                         << 1146                                               1.5*std::exp(-y*y);
1144     break;                                    << 1147   }
1145     //  ===================================== << 1148       else if(hLabMomentum >= 0.65)
1146   case 5:              //   K minus           << 1149   {
1147                                               << 1150           G4double x = (hLabMomentum - 0.72)/0.06;
1148     HadrTot   = 10+1.8*logE + 25./sqrS; // mb << 1151           G4double y = (hLabMomentum - 1.015)/0.075;
1149     HadrSlope = 6.98+0.127*logS;         // G << 1152     TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y);
1150     HadrReIm  = 0.4*(sHadr-20)*(sHadr-20)*G4E << 1153   }
1151     DDSect2   = 0.7;                          << 1154       else if(hLabMomentum >= 0.37)
1152     DDSect3   = 0.27;                         << 1155   {
1153     break;                                    << 1156     G4double x = std::log(hLabMomentum/0.48);
                                                   >> 1157     TotN = 26. + 110.*x*x;
                                                   >> 1158   }
                                                   >> 1159       else 
                                                   >> 1160   {
                                                   >> 1161           G4double x = (hLabMomentum - 0.29)/0.07;
                                                   >> 1162     TotN = 28.0 + 40.*std::exp(-x*x);
                                                   >> 1163   }
                                                   >> 1164       HadrTot = (TotP+TotN)/2;
                                                   >> 1165 //  ........................................
                                                   >> 1166       HadrSlope = 7.28+0.245*logS;        // GeV-2
                                                   >> 1167       HadrReIm  = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15);
                                                   >> 1168 
                                                   >> 1169       DDSect2   = 0.7;                               //mb*GeV-2
                                                   >> 1170       DDSect3   = 0.27;                              //mb*GeV-2
                                                   >> 1171 
                                                   >> 1172       break;
                                                   >> 1173 //  ==========================================================
                                                   >> 1174     case 4:            //  K plus
                                                   >> 1175 
                                                   >> 1176       HadrTot   = 10.6+1.8*logE + 9.0*std::pow(HadrEnergy,-0.55);  // mb
                                                   >> 1177       if(HadrEnergy>100) HadrSlope = 15.0;
                                                   >> 1178       else HadrSlope = 1.0+1.76*logS - 2.84/sqrS;   // GeV-2
                                                   >> 1179 
                                                   >> 1180       HadrReIm  = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
                                                   >> 1181       DDSect2   = 0.7;                             //mb*GeV-2
                                                   >> 1182       DDSect3   = 0.21;                            //mb*GeV-2
                                                   >> 1183       break;
                                                   >> 1184 //  =========================================================
                                                   >> 1185      case 5:              //   K minus
                                                   >> 1186 
                                                   >> 1187        HadrTot   = 10+1.8*logE + 25./sqrS; // mb
                                                   >> 1188        HadrSlope = 6.98+0.127*logS;         // GeV-2
                                                   >> 1189        HadrReIm  = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
                                                   >> 1190        DDSect2   = 0.7;                             //mb*GeV-2
                                                   >> 1191        DDSect3   = 0.27;                            //mb*GeV-2
                                                   >> 1192        break;
1154   }                                              1193   }   
1155   //  ======================================= << 1194 //  =========================================================
1156   if(verboseLevel>2) {                        << 1195   if(verboseLevel>2)
1157     G4cout << "HadrTot= " << HadrTot << " Had    1196     G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope
1158      << " HadrReIm= " << HadrReIm << " DDSect    1197      << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2
1159      << " DDSect3= " << DDSect3 << G4endl;       1198      << " DDSect3= " << DDSect3 << G4endl;
1160   }                                           << 1199 
1161   if(Z != 1) return;                             1200   if(Z != 1) return;
1162                                                  1201 
1163   // Scattering of protons                       1202   // Scattering of protons
1164                                                  1203 
1165   Coeff0 = Coeff1 = Coeff2 = 0.0;                1204   Coeff0 = Coeff1 = Coeff2 = 0.0;
1166   Slope0 = Slope1 = 1.0;                         1205   Slope0 = Slope1 = 1.0;
1167   Slope2 = 5.0;                                  1206   Slope2 = 5.0;
1168                                                  1207 
1169   // data for iHadron=0                          1208   // data for iHadron=0
1170   static const G4double EnP0[6]={1.5,3.0,5.0,    1209   static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0};
1171   static const G4double C0P0[6]={0.15,0.02,0.    1210   static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002};
1172   static const G4double C1P0[6]={0.05,0.02,0.    1211   static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0};
1173   static const G4double B0P0[6]={1.5,2.5,3.0,    1212   static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25};
1174   static const G4double B1P0[6]={5.0,1.0,3.5,    1213   static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8};
1175                                                  1214       
1176   // data for iHadron=6,7                        1215   // data for iHadron=6,7
1177   static const G4double EnN[5]={1.5,5.0,10.0,    1216   static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0};
1178   static const G4double C0N[5]={0.0,0.0,0.02,    1217   static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01};
1179   static const G4double C1N[5]={0.06,0.008,0.    1218   static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003};
1180   static const G4double B0N[5]={1.5,2.5,3.8,3    1219   static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5};
1181   static const G4double B1N[5]={1.5,2.2,3.6,4    1220   static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8};
1182                                                  1221 
1183   // data for iHadron=1                          1222   // data for iHadron=1
1184   static const G4double EnP[2]={1.5,4.0};        1223   static const G4double EnP[2]={1.5,4.0};
1185   static const G4double C0P[2]={0.001,0.0005}    1224   static const G4double C0P[2]={0.001,0.0005};
1186   static const G4double C1P[2]={0.003,0.001};    1225   static const G4double C1P[2]={0.003,0.001};
1187   static const G4double B0P[2]={2.5,4.5};        1226   static const G4double B0P[2]={2.5,4.5};
1188   static const G4double B1P[2]={1.0,4.0};        1227   static const G4double B1P[2]={1.0,4.0};
1189                                                  1228 
1190   // data for iHadron=2                          1229   // data for iHadron=2
1191   static const G4double EnPP[4]={1.0,2.0,3.0,    1230   static const G4double EnPP[4]={1.0,2.0,3.0,4.0};
1192   static const G4double C0PP[4]={0.0,0.0,0.0,    1231   static const G4double C0PP[4]={0.0,0.0,0.0,0.0};
1193   static const G4double C1PP[4]={0.15,0.08,0.    1232   static const G4double C1PP[4]={0.15,0.08,0.02,0.01};
1194   static const G4double B0PP[4]={1.5,2.8,3.8,    1233   static const G4double B0PP[4]={1.5,2.8,3.8,3.8};
1195   static const G4double B1PP[4]={0.8,1.6,3.6,    1234   static const G4double B1PP[4]={0.8,1.6,3.6,4.6};
1196                                                  1235 
1197   // data for iHadron=3                          1236   // data for iHadron=3
1198   static const G4double EnPPN[4]={1.0,2.0,3.0    1237   static const G4double EnPPN[4]={1.0,2.0,3.0,4.0};
1199   static const G4double C0PPN[4]={0.0,0.0,0.0    1238   static const G4double C0PPN[4]={0.0,0.0,0.0,0.0};
1200   static const G4double C1PPN[4]={0.0,0.0,0.0    1239   static const G4double C1PPN[4]={0.0,0.0,0.0,0.0};
1201   static const G4double B0PPN[4]={1.5,2.8,3.8    1240   static const G4double B0PPN[4]={1.5,2.8,3.8,3.8};
1202   static const G4double B1PPN[4]={0.8,1.6,3.6    1241   static const G4double B1PPN[4]={0.8,1.6,3.6,4.6};
1203                                                  1242 
1204   // data for iHadron=4                          1243   // data for iHadron=4
1205   static const G4double EnK[4]={1.4,2.33,3.0,    1244   static const G4double EnK[4]={1.4,2.33,3.0,5.0};
1206   static const G4double C0K[4]={0.0,0.0,0.0,0    1245   static const G4double C0K[4]={0.0,0.0,0.0,0.0};
1207   static const G4double C1K[4]={0.01,0.007,0.    1246   static const G4double C1K[4]={0.01,0.007,0.005,0.003};
1208   static const G4double B0K[4]={1.5,2.0,3.8,3    1247   static const G4double B0K[4]={1.5,2.0,3.8,3.8};
1209   static const G4double B1K[4]={1.6,1.6,1.6,1    1248   static const G4double B1K[4]={1.6,1.6,1.6,1.6};
1210                                                  1249 
1211   // data for iHadron=5                          1250   // data for iHadron=5
1212   static const G4double EnKM[2]={1.4,4.0};       1251   static const G4double EnKM[2]={1.4,4.0};
1213   static const G4double C0KM[2]={0.006,0.002}    1252   static const G4double C0KM[2]={0.006,0.002};
1214   static const G4double C1KM[2]={0.00,0.00};     1253   static const G4double C1KM[2]={0.00,0.00};
1215   static const G4double B0KM[2]={2.5,3.5};       1254   static const G4double B0KM[2]={2.5,3.5};
1216   static const G4double B1KM[2]={1.6,1.6};       1255   static const G4double B1KM[2]={1.6,1.6};
1217                                                  1256 
1218   switch(iHadron) {                           << 1257   switch(iHadron)
1219   case 0:                                     << 1258     {
                                                   >> 1259     case 0 :
1220                                                  1260 
1221     if(hLabMomentum <BoundaryP[0]) {          << 1261       if(hLabMomentum <BoundaryP[0])
1222       InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P << 1262   InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0);
1223     }                                         << 
1224     Coeff2 = 0.8/hLabMomentum2;               << 
1225     break;                                    << 
1226                                                  1263 
1227   case 6:                                     << 1264       Coeff2 = 0.8/hLabMomentum/hLabMomentum;
                                                   >> 1265       break; 
1228                                                  1266 
1229     if(hLabMomentum < BoundaryP[1]) {         << 1267     case  6 :
1230       InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);   << 
1231     }                                         << 
1232     Coeff2 = 0.8/hLabMomentum2;               << 
1233     break;                                    << 
1234                                                  1268 
1235   case 1:                                     << 1269       if(hLabMomentum < BoundaryP[1])
1236   case 7:                                     << 1270   InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);
1237     if(hLabMomentum <  BoundaryP[2]) {        << 
1238       InterpolateHN(2,EnP,C0P,C1P,B0P,B1P);   << 
1239     }                                         << 
1240     break;                                    << 
1241                                                  1271 
1242   case 2:                                     << 1272       Coeff2 = 0.8/hLabMomentum/hLabMomentum;
                                                   >> 1273       break; 
1243                                                  1274 
1244     if(hLabMomentum < BoundaryP[3]) {         << 1275     case 1 :
1245       InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1P << 1276     case 7 :
1246     }                                         << 1277       if(hLabMomentum <  BoundaryP[2])
1247     Coeff2 = 0.02/hLabMomentum;               << 1278   InterpolateHN(2,EnP,C0P,C1P,B0P,B1P);
1248     break;                                    << 1279       break; 
1249                                                  1280 
1250   case 3:                                     << 1281     case 2 :
1251                                                  1282 
1252     if(hLabMomentum < BoundaryP[4]) {         << 1283       if(hLabMomentum < BoundaryP[3])
1253       InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN << 1284   InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP);
1254     }                                         << 1285 
1255     Coeff2 = 0.02/hLabMomentum;               << 1286       Coeff2 = 0.02/hLabMomentum;
1256     break;                                    << 1287       break; 
                                                   >> 1288 
                                                   >> 1289     case 3 :
                                                   >> 1290 
                                                   >> 1291       if(hLabMomentum < BoundaryP[4])
                                                   >> 1292   InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN);
                                                   >> 1293 
                                                   >> 1294       Coeff2 = 0.02/hLabMomentum;
                                                   >> 1295       break;
1257                                                  1296  
1258   case 4:                                     << 1297     case 4 :
1259                                                  1298 
1260     if(hLabMomentum < BoundaryP[5]) {         << 1299       if(hLabMomentum < BoundaryP[5])
1261       InterpolateHN(4,EnK,C0K,C1K,B0K,B1K);   << 1300   InterpolateHN(4,EnK,C0K,C1K,B0K,B1K);
1262     }                                         << 1301 
1263     if(hLabMomentum < 1) { Coeff2 = 0.34; }   << 1302       if(hLabMomentum < 1) Coeff2 = 0.34;
1264     else  { Coeff2 = 0.34/(hLabMomentum2*hLab << 1303       else  Coeff2 = 0.34/hLabMomentum2/hLabMomentum;
1265     break;                                    << 1304       break; 
1266                                               << 1305 
1267   case 5:                                     << 1306     case 5 :
1268     if(hLabMomentum < BoundaryP[6]) {         << 1307       if(hLabMomentum < BoundaryP[6])
1269       InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1K << 1308   InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM);
                                                   >> 1309 
                                                   >> 1310       if(hLabMomentum < 1) Coeff2 = 0.01;
                                                   >> 1311       else  Coeff2 = 0.01/hLabMomentum2/hLabMomentum;
                                                   >> 1312       break; 
1270     }                                            1313     }
1271     if(hLabMomentum < 1) { Coeff2 = 0.01; }   << 
1272     else  { Coeff2 = 0.01/(hLabMomentum2*hLab << 
1273     break;                                    << 
1274   }                                           << 
1275                                                  1314 
1276   if(verboseLevel > 2) {                      << 1315   if(verboseLevel > 2) 
1277     G4cout<<"  HadrVal : Plasb  "<<hLabMoment    1316     G4cout<<"  HadrVal : Plasb  "<<hLabMomentum
1278     <<"  iHadron  "<<iHadron<<"  HadrTot  "<<    1317     <<"  iHadron  "<<iHadron<<"  HadrTot  "<<HadrTot<<G4endl;
1279   }                                           << 
1280 }                                                1318 }
1281                                                  1319 
1282 ///////////////////////////////////////////// << 1320 //  =====================================================
                                                   >> 1321 void  G4ElasticHadrNucleusHE::
                                                   >> 1322        GetKinematics(const G4ParticleDefinition * aHadron,
                                                   >> 1323                            G4double MomentumH)
                                                   >> 1324 {
                                                   >> 1325   if (verboseLevel>1)
                                                   >> 1326     G4cout<<"1  GetKin.: HadronName MomentumH "
                                                   >> 1327     <<aHadron->GetParticleName()<<"  "<<MomentumH<<G4endl;
                                                   >> 1328 
                                                   >> 1329   DefineHadronValues(1);
1283                                                  1330 
                                                   >> 1331   G4double Sh     = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV
                                                   >> 1332 
                                                   >> 1333   ConstU = 2*protonM2+2*hMass2-Sh;
                                                   >> 1334 
                                                   >> 1335   G4double MaxT = 4*MomentumCM*MomentumCM;
                                                   >> 1336 
                                                   >> 1337   BoundaryTL[0] = MaxT; //2.0;
                                                   >> 1338   BoundaryTL[1] = MaxT;
                                                   >> 1339   BoundaryTL[3] = MaxT;
                                                   >> 1340   BoundaryTL[4] = MaxT;
                                                   >> 1341   BoundaryTL[5] = MaxT;
                                                   >> 1342 
                                                   >> 1343   G4int NumberH=0;
                                                   >> 1344 
                                                   >> 1345   while(iHadrCode!=HadronCode[NumberH]) NumberH++;
                                                   >> 1346 
                                                   >> 1347   NumberH = HadronType1[NumberH];   
                                                   >> 1348 
                                                   >> 1349   if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
                                                   >> 1350   else MaxTR = BoundaryTG[NumberH];
                                                   >> 1351 
                                                   >> 1352   if (verboseLevel>1)
                                                   >> 1353     G4cout<<"3  GetKin. : NumberH  "<<NumberH
                                                   >> 1354     <<"  Bound.P[NumberH] "<<BoundaryP[NumberH]
                                                   >> 1355     <<"  Bound.TL[NumberH] "<<BoundaryTL[NumberH]
                                                   >> 1356     <<"  Bound.TG[NumberH] "<<BoundaryTG[NumberH]
                                                   >> 1357     <<"  MaxT MaxTR "<<MaxT<<"  "<<MaxTR<<G4endl;
                                                   >> 1358 
                                                   >> 1359 //     GetParametersHP(aHadron, MomentumH);
                                                   >> 1360 }
                                                   >> 1361 //  ============================================================
1284 G4double G4ElasticHadrNucleusHE::GetFt(G4doub    1362 G4double G4ElasticHadrNucleusHE::GetFt(G4double Q2)
1285 {                                                1363 {
1286   G4double Fdistr=0;                             1364   G4double Fdistr=0;
1287   G4double SqrQ2 = std::sqrt(Q2);                1365   G4double SqrQ2 = std::sqrt(Q2);
1288                                                  1366  
1289   Fdistr = (1-Coeff1-Coeff0) / HadrSlope*(1-G << 1367   Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU))
1290     + Coeff0*(1-G4Exp(-Slope0*Q2))            << 1368     /HadrSlope*(1-std::exp(-HadrSlope*Q2))
1291     + Coeff2/Slope2*G4Exp(Slope2*ConstU)*(G4E << 1369     + Coeff0*(1-std::exp(-Slope0*Q2))
1292     + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+Sqr << 1370     + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1)
                                                   >> 1371     + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
1293                                                  1372 
1294   if (verboseLevel>1) {                       << 1373   if (verboseLevel>1)
1295     G4cout<<"Old:  Coeff0 Coeff1 Coeff2 "<<Co    1374     G4cout<<"Old:  Coeff0 Coeff1 Coeff2 "<<Coeff0<<"  "
1296           <<Coeff1<<"  "<<Coeff2<<"  Slope Sl    1375           <<Coeff1<<"  "<<Coeff2<<"  Slope Slope0 Slope1 Slope2 "
1297           <<HadrSlope<<"  "<<Slope0<<"  "<<Sl    1376           <<HadrSlope<<"  "<<Slope0<<"  "<<Slope1<<"  "<<Slope2
1298           <<"  Fdistr "<<Fdistr<<G4endl;         1377           <<"  Fdistr "<<Fdistr<<G4endl;
1299   }                                           << 
1300   return Fdistr;                                 1378   return Fdistr;
1301 }                                                1379 }
1302                                               << 1380 //  +++++++++++++++++++++++++++++++++++++++
1303 ///////////////////////////////////////////// << 1381 G4double G4ElasticHadrNucleusHE::GetQ2(G4double Ran)
1304                                               << 
1305 G4double                                      << 
1306 G4ElasticHadrNucleusHE::HadronProtonQ2(G4doub << 
1307 {                                                1382 {
1308   hLabMomentum  = plab;                       << 1383   G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
1309   hLabMomentum2 = hLabMomentum*hLabMomentum;  << 1384   G4double Q2=0;
1310   HadrEnergy    = std::sqrt(hMass2 + hLabMome << 
1311   DefineHadronValues(1);                      << 
1312                                               << 
1313   G4double Sh = 2.0*protonM*HadrEnergy+proton << 
1314   ConstU = 2*protonM2+2*hMass2-Sh;            << 
1315                                               << 
1316   BoundaryTL[0] = tmax;                       << 
1317   BoundaryTL[1] = tmax;                       << 
1318   BoundaryTL[3] = tmax;                       << 
1319   BoundaryTL[4] = tmax;                       << 
1320   BoundaryTL[5] = tmax;                       << 
1321                                               << 
1322   G4double MaxTR = (plab < BoundaryP[iHadron1 << 
1323     BoundaryTL[iHadron1] : BoundaryTG[iHadron << 
1324                                               << 
1325   if (verboseLevel>1) {                       << 
1326     G4cout<<"3  GetKin. : iHadron1  "<<iHadro << 
1327     <<"  Bound.P[iHadron1] "<<BoundaryP[iHadr << 
1328     <<"  Bound.TL[iHadron1] "<<BoundaryTL[iHa << 
1329     <<"  Bound.TG[iHadron1] "<<BoundaryTG[iHa << 
1330     <<"  MaxT MaxTR "<<tmax<<"  "<<MaxTR<<G4e << 
1331   }                                           << 
1332                                               << 
1333   G4double rand = G4UniformRand();            << 
1334                                               << 
1335   G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=Max << 
1336                                                  1385 
1337   G4double norm  = 1.0/GetFt(MaxTR);          << 1386   FmaxT = GetFt(MaxTR);
1338   G4double delta = GetFt(DDD0)*norm - rand;   << 1387   delta = GetDistrFun(DDD0)-Ran;
1339                                                  1388 
1340   static const G4int maxNumberOfLoops = 10000 << 1389   while(std::fabs(delta) > 0.0001)
1341   G4int loopCounter = -1;                     << 
1342   while ( (std::abs(delta) > 0.0001) &&       << 
1343           ++loopCounter < maxNumberOfLoops )  << 
1344     {                                            1390     {
1345       if(delta>0)                                1391       if(delta>0) 
1346       {                                          1392       {
1347         DDD2 = DDD0;                             1393         DDD2 = DDD0;
1348         DDD0 = (DDD0+DDD1)*0.5;                  1394         DDD0 = (DDD0+DDD1)*0.5;
1349       }                                          1395       }
1350       else if(delta<0.0)                      << 1396       else if(delta<0)
1351       {                                          1397       {
1352         DDD1 = DDD0;                             1398         DDD1 = DDD0; 
1353         DDD0 = (DDD0+DDD2)*0.5;                  1399         DDD0 = (DDD0+DDD2)*0.5;
1354       }                                          1400       }
1355       delta = GetFt(DDD0)*norm - rand;        << 1401       delta = GetDistrFun(DDD0)-Ran;
1356     }                                            1402     }
1357   return (loopCounter >= maxNumberOfLoops) ?  << 
1358 }                                             << 
1359                                                  1403 
1360 ///////////////////////////////////////////// << 1404   Q2 = DDD0;
1361                                                  1405 
1362 void G4ElasticHadrNucleusHE::Binom()          << 1406   return Q2;
1363 {                                             << 
1364   for(G4int N = 0; N < 240; ++N) {            << 
1365     G4double J = 1.0;                         << 
1366     for(G4int M = 0; M <= N; ++M) {           << 
1367       G4double Fact1 = 1.0;                   << 
1368       if (N > 0 && N > M && M > 0 ) {         << 
1369   J *= (G4double)(N-M+1)/(G4double)M;         << 
1370   Fact1 = J;                                  << 
1371       }                                       << 
1372       fBinom[N][M] = Fact1;                   << 
1373     }                                         << 
1374   }                                           << 
1375 }                                                1407 }
                                                   >> 1408 //  ++++++++++++++++++++++++++++++++++++++++++
                                                   >> 1409 G4double G4ElasticHadrNucleusHE::
                                                   >> 1410        HadronProtonQ2(const G4ParticleDefinition * p,
                                                   >> 1411                             G4double inLabMom)
                                                   >> 1412 {
1376                                                  1413 
1377 ///////////////////////////////////////////// << 1414   hMass         = p->GetPDGMass()/GeV;
                                                   >> 1415   hMass2        = hMass*hMass;
                                                   >> 1416   hLabMomentum  = inLabMom;
                                                   >> 1417   hLabMomentum2 = hLabMomentum*hLabMomentum;
                                                   >> 1418   HadrEnergy    = sqrt(hLabMomentum2+hMass2); 
1378                                                  1419 
1379 void                                          << 1420   G4double Rand = G4UniformRand();
1380 G4ElasticHadrNucleusHE::InFileName(std::ostri << 
1381            const G4ParticleDefinition* p, G4i << 
1382 {                                             << 
1383   if(!fDirectory) {                           << 
1384     fDirectory = G4FindDataDir("G4LEDATA");   << 
1385     if (fDirectory) {                         << 
1386       ss << fDirectory << "/";                << 
1387     }                                         << 
1388   }                                           << 
1389   OutFileName(ss, p, Z);                      << 
1390 }                                             << 
1391                                                  1421 
1392 ///////////////////////////////////////////// << 1422   GetKinematics(p, inLabMom);
1393                                                  1423 
1394 void                                          << 1424   G4double Q2 = GetQ2(Rand);
1395 G4ElasticHadrNucleusHE::OutFileName(std::ostr << 1425 
1396             const G4ParticleDefinition* p, G4 << 1426   return Q2;
1397 {                                             << 
1398   ss << "hedata/" << p->GetParticleName() <<  << 
1399 }                                                1427 }
1400                                                  1428 
1401 ///////////////////////////////////////////// << 1429 //  ===========================================
                                                   >> 1430 
                                                   >> 1431 ///////////////////////////////////////////////////////////////////
                                                   >> 1432 //
                                                   >> 1433 //  
1402                                                  1434 
1403 G4bool G4ElasticHadrNucleusHE::ReadLine(std:: << 1435 void G4ElasticHadrNucleusHE::Binom()
1404           std::vector<G4double>& v)           << 
1405 {                                                1436 {
1406   G4int n(0);                                 << 1437   G4int N, M;
1407   infile >> n;                                << 1438   G4double  Fact1, J;
1408   if (infile.fail()) { return false; }        << 
1409   if(n > 0) {                                 << 
1410     v.reserve(n);                             << 
1411     G4double x(0.0);                          << 
1412     for(G4int i=0; i<n; ++i) {                << 
1413       infile >> x;                            << 
1414       if (infile.fail()) { return false; }    << 
1415       v.emplace_back(x);                      << 
1416     }                                         << 
1417   }                                           << 
1418   return true;                                << 
1419 }                                             << 
1420                                                  1439 
1421 ///////////////////////////////////////////// << 1440   for(N = 0; N < 240; N++)
                                                   >> 1441   {
                                                   >> 1442     J = 1;
1422                                                  1443 
1423 void G4ElasticHadrNucleusHE::WriteLine(std::o << 1444       for( M = 0; M <= N; M++ )
1424                std::vector<G4double>& v)      << 1445       {
1425 {                                             << 1446     Fact1 = 1;
1426   std::size_t n = v.size();                   << 1447 
1427   outfile << n << G4endl;                     << 1448     if ( ( N > 0 ) && ( N > M ) && M > 0 )
1428   if(n > 0) {                                 << 1449           {
1429     for(std::size_t i=0; i<n; ++i) {          << 1450               J *= G4double(N-M+1)/G4double(M);
1430       outfile << v[i] << " ";                 << 1451               Fact1 = J;
1431     }                                         << 1452           }
1432     outfile << G4endl;                        << 1453     SetBinom[N][M] = Fact1;
                                                   >> 1454       }
1433   }                                              1455   }
                                                   >> 1456   return;
1434 }                                                1457 }
1435                                                  1458 
                                                   >> 1459 
                                                   >> 1460 //
                                                   >> 1461 //
1436 /////////////////////////////////////////////    1462 ///////////////////////////////////////////////////////////
                                                   >> 1463 
1437                                                  1464