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 11.2.1)


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