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


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