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