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