Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.cc (Version 9.6.p4)


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