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