Geant4_10
G4ComponentGGNuclNuclXsc.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 // 24.11.08 V. Grichine - first implementation
27 //
28 
30 
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4ParticleTable.hh"
34 #include "G4IonTable.hh"
35 #include "G4ParticleDefinition.hh"
36 #include "G4HadTmpUtil.hh"
37 #include "G4HadronNucleonXsc.hh"
38 
39 
41  : G4VComponentCrossSection("Glauber-Gribov nucleus nucleus"),
42 // fUpperLimit(100000*GeV),
43  fLowerLimit(0.1*MeV),
44  fRadiusConst(1.08*fermi), // 1.1, 1.3 ?
45  fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fProductionXsc(0.0),
46  fDiffractionXsc(0.0),
47  cacheDP(G4Proton::Proton(),G4ParticleMomentum(1.,0,0),0),
48  dProton(G4Proton::Proton(),G4ParticleMomentum(1.,0,0),0),
49  dNeutron(G4Neutron::Neutron(),G4ParticleMomentum(1.,0,0),0)
50 // , fHadronNucleonXsc(0.0)
51 {
52  theProton = G4Proton::Proton();
53  theNeutron = G4Neutron::Neutron();
54  hnXsc = new G4HadronNucleonXsc();
55 }
56 
57 
59 {
60  delete hnXsc;
61 }
62 
64 
66  G4double kinEnergy,
67  G4int Z, G4int A)
68 {
69  cacheDP.SetDefinition(aParticle);
70  cacheDP.SetKineticEnergy(kinEnergy);
71  fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
72  return fTotalXsc;
73 }
74 
76 
78  G4double kinEnergy,
79  G4int Z, G4double A)
80 {
81  cacheDP.SetDefinition(aParticle);
82  cacheDP.SetKineticEnergy(kinEnergy);
83  fInelasticXsc = GetZandACrossSection(&cacheDP, Z, G4int(A));
84  return fTotalXsc;
85 }
86 
88 
90  G4double kinEnergy,
91  G4int Z, G4int A)
92 {
93  cacheDP.SetDefinition(aParticle);
94  cacheDP.SetKineticEnergy(kinEnergy);
95  fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
96  return fInelasticXsc;
97 }
98 
100 
102  G4double kinEnergy,
103  G4int Z, G4double A)
104 {
105  cacheDP.SetDefinition(aParticle);
106  cacheDP.SetKineticEnergy(kinEnergy);
107  fInelasticXsc = GetZandACrossSection(&cacheDP, Z, G4int(A));
108  return fInelasticXsc;
109 }
110 
112 
114  G4double kinEnergy,
115  G4int Z, G4double A)
116 {
117  cacheDP.SetDefinition(aParticle);
118  cacheDP.SetKineticEnergy(kinEnergy);
119  fInelasticXsc = GetZandACrossSection(&cacheDP, Z, G4int(A));
120  return fElasticXsc;
121 }
122 
124 
126  G4double kinEnergy,
127  G4int Z, G4int A)
128 {
129  cacheDP.SetDefinition(aParticle);
130  cacheDP.SetKineticEnergy(kinEnergy);
131  fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
132  return fElasticXsc;
133 }
134 
136 
138  G4double kinEnergy,
139  G4int Z, G4int A)
140 {
141  cacheDP.SetDefinition(aParticle);
142  cacheDP.SetKineticEnergy(kinEnergy);
143  fInelasticXsc = GetZandACrossSection(&cacheDP, Z, A);
144  G4double ratio = 0.;
145 
146  if(fInelasticXsc > 0.)
147  {
148  ratio = (fInelasticXsc - fProductionXsc)/fInelasticXsc;
149  if(ratio < 0.) ratio = 0.;
150  }
151  return ratio;
152 }
153 
155 
156 void
158 {
159  outFile << "G4ComponentGGNuclNuclXsc calculates total, inelastic and\n"
160  << "elastic cross sections for nucleus-nucleus collisions using\n"
161  << "the Glauber model with Gribov corrections. It is valid for\n"
162  << "all incident energies above 100 keV./n";
163 }
164 
166 
167 G4bool
169  G4int Z, const G4Material*)
170 {
171  G4bool applicable = false;
172  G4double kineticEnergy = aDP->GetKineticEnergy();
173 
174  if (kineticEnergy >= fLowerLimit && Z > 1) applicable = true;
175  return applicable;
176 }
177 
179 //
180 // Calculates total and inelastic Xsc, derives elastic as total - inelastic
181 // accordong to Glauber model with Gribov correction calculated in the dipole
182 // approximation on light cone. Gaussian density helps to calculate rest
183 // integrals of the model. [1] B.Z. Kopeliovich, nucl-th/0306044
184 
185 
188  const G4Material*)
189 {
190  G4int A = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
191  return GetZandACrossSection(aParticle, Z, A);
192 }
193 
195 //
196 // Calculates total and inelastic Xsc, derives elastic as total - inelastic
197 // accordong to Glauber model with Gribov correction calculated in the dipole
198 // approximation on light cone. Gaussian density of point-like nucleons helps
199 // to calculate rest integrals of the model. [1] B.Z. Kopeliovich,
200 // nucl-th/0306044 + simplification above
201 
202 
205  G4int tZ, G4int tA)
206 {
207  G4double xsection;
208  G4double sigma;
209  G4double cofInelastic = 2.4;
210  G4double cofTotal = 2.0;
211  G4double nucleusSquare;
212  G4double cB;
213  G4double ratio;
214 
215  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
216  G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
217 
218  G4double pTkin = aParticle->GetKineticEnergy();
219  pTkin /= pA;
220 
221  G4double pN = pA - pZ;
222  if( pN < 0. ) pN = 0.;
223 
224  G4double tN = tA - tZ;
225  if( tN < 0. ) tN = 0.;
226 
227  G4double tR = GetNucleusRadius( G4double(tZ),G4double(tA) );
228  G4double pR = GetNucleusRadius(pZ,pA);
229 
230  cB = GetCoulombBarier(aParticle, G4double(tZ), G4double(tA), pR, tR);
231 
232  if ( cB > 0. )
233  {
234  dProton.SetKineticEnergy(pTkin);
235  dNeutron.SetKineticEnergy(pTkin);
236 
237  sigma = (pZ*tZ+pN*tN)*hnXsc->GetHadronNucleonXscNS(&dProton, theProton);
238 
239  G4double ppInXsc = hnXsc->GetInelasticHadronNucleonXsc();
240 
241  sigma += (pZ*tN+pN*tZ)*hnXsc->GetHadronNucleonXscNS(&dNeutron, theProton);
242 
243  G4double npInXsc = hnXsc->GetInelasticHadronNucleonXsc();
244 
245  // G4cout<<"ppInXsc = "<<ppInXsc/millibarn<<"; npInXsc = "<<npInXsc/millibarn<<G4endl;
246  // G4cout<<"npTotXsc = "<<hnXsc->GetTotalHadronNucleonXsc()/millibarn<<"; npElXsc = "
247  // <<hnXsc->GetElasticHadronNucleonXsc()/millibarn<<G4endl;
248 
249  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
250 
251  ratio = sigma/nucleusSquare;
252  xsection = nucleusSquare*std::log( 1. + ratio );
253  fTotalXsc = xsection;
254  fTotalXsc *= cB;
255 
256  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
257 
258  fInelasticXsc *= cB;
259  fElasticXsc = fTotalXsc - fInelasticXsc;
260 
261  // if (fElasticXsc < DBL_MIN) fElasticXsc = DBL_MIN;
262  /*
263  G4double difratio = ratio/(1.+ratio);
264 
265  fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
266  */
267  // production to be checked !!! edit MK xsc
268 
269  //sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscMK(theProton, pTkin, theProton) +
270  // (pZ*tN+pN*tZ)*GetHadronNucleonXscMK(theProton, pTkin, theNeutron);
271 
272  sigma = (pZ*tZ+pN*tN)*ppInXsc + (pZ*tN+pN*tZ)*npInXsc;
273 
274  ratio = sigma/nucleusSquare;
275  fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
276 
277  if (fElasticXsc < 0.) fElasticXsc = 0.;
278  }
279  else
280  {
281  fInelasticXsc = 0.;
282  fTotalXsc = 0.;
283  fElasticXsc = 0.;
284  fProductionXsc = 0.;
285  }
286  return fInelasticXsc; // xsection;
287 }
288 
290 //
291 //
292 
295  G4double pR, G4double tR)
296 {
297  G4double ratio;
298  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
299 
300  G4double pTkin = aParticle->GetKineticEnergy();
301  // G4double pPlab = aParticle->GetTotalMomentum();
302  G4double pM = aParticle->GetDefinition()->GetPDGMass();
303  // G4double tM = tZ*proton_mass_c2 + (tA-tZ)*neutron_mass_c2; // ~ 1% accuracy
305  G4double pElab = pTkin + pM;
306  G4double totEcm = std::sqrt(pM*pM + tM*tM + 2.*pElab*tM);
307  // G4double pPcm = pPlab*tM/totEcm;
308  // G4double pTcm = std::sqrt(pM*pM + pPcm*pPcm) - pM;
309  G4double totTcm = totEcm - pM -tM;
310 
312  bC /= pR + tR;
313  bC /= 2.; // 4., 2. parametrisation cof ??? vmg
314 
315  // G4cout<<"pTkin = "<<pTkin/GeV<<"; pPlab = "
316  // <<pPlab/GeV<<"; bC = "<<bC/GeV<<"; pTcm = "<<pTcm/GeV<<G4endl;
317 
318  if( totTcm <= bC ) ratio = 0.;
319  else ratio = 1. - bC/totTcm;
320 
321  // if(ratio < DBL_MIN) ratio = DBL_MIN;
322  if( ratio < 0.) ratio = 0.;
323 
324  // G4cout <<"ratio = "<<ratio<<G4endl;
325  return ratio;
326 }
327 
328 
330 //
331 // Return single-diffraction/inelastic cross-section ratio
332 
335 {
336  G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
337 
338  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
339  G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
340 
341  G4double pTkin = aParticle->GetKineticEnergy();
342  pTkin /= pA;
343 
344  G4double pN = pA - pZ;
345  if( pN < 0. ) pN = 0.;
346 
347  G4double tN = tA - tZ;
348  if( tN < 0. ) tN = 0.;
349 
350  G4double tR = GetNucleusRadius(tZ,tA);
351  G4double pR = GetNucleusRadius(pZ,pA);
352 
353  sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
354  (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
355 
356  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
357  ratio = sigma/nucleusSquare;
358  fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
359  G4double difratio = ratio/(1.+ratio);
360 
361  fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
362 
363  if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
364  else ratio = 0.;
365 
366  return ratio;
367 }
368 
370 //
371 // Return quasi-elastic/inelastic cross-section ratio
372 
375 {
376  G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
377 
378  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
379  G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
380 
381  G4double pTkin = aParticle->GetKineticEnergy();
382  pTkin /= pA;
383 
384  G4double pN = pA - pZ;
385  if( pN < 0. ) pN = 0.;
386 
387  G4double tN = tA - tZ;
388  if( tN < 0. ) tN = 0.;
389 
390  G4double tR = GetNucleusRadius(tZ,tA);
391  G4double pR = GetNucleusRadius(pZ,pA);
392 
393  sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
394  (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
395 
396  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
397  ratio = sigma/nucleusSquare;
398  fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
399 
400  // sigma = GetHNinelasticXsc(aParticle, tA, tZ);
401  ratio = sigma/nucleusSquare;
402  fProductionXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
403 
404  if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
405  else ratio = 0.;
406  if ( ratio < 0. ) ratio = 0.;
407 
408  return ratio;
409 }
410 
412 //
413 // Returns hadron-nucleon Xsc according to differnt parametrisations:
414 // [2] E. Levin, hep-ph/9710546
415 // [3] U. Dersch, et al, hep-ex/9910052
416 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
417 
418 G4double
420  const G4Element* anElement)
421 {
422  G4int At = G4lrint(anElement->GetN()); // number of nucleons
423  G4int Zt = G4lrint(anElement->GetZ()); // number of protons
424  return GetHadronNucleonXsc(aParticle, At, Zt);
425 }
426 
427 
428 
429 
431 //
432 // Returns hadron-nucleon Xsc according to differnt parametrisations:
433 // [2] E. Levin, hep-ph/9710546
434 // [3] U. Dersch, et al, hep-ex/9910052
435 // [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
436 
437 G4double
439  G4int At, G4int Zt)
440 {
441  G4double xsection = 0.;
442 
444  GetIonTable()->GetIonMass(Zt, At);
445  targ_mass = 0.939*GeV; // ~mean neutron and proton ???
446 
447  G4double proj_mass = aParticle->GetMass();
448  G4double proj_momentum = aParticle->GetMomentum().mag();
449  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
450 
451  sMand /= GeV*GeV; // in GeV for parametrisation
452  proj_momentum /= GeV;
453  const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
454 
455  if(pParticle == theNeutron) // as proton ???
456  {
457  xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
458  }
459  else if(pParticle == theProton)
460  {
461  xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
462  }
463 
464  xsection *= millibarn;
465  return xsection;
466 }
467 
468 
469 
471 //
472 // Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
473 // http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
474 // At = number of nucleons, Zt = number of protons
475 
476 G4double
478  G4double sMand,
479  G4ParticleDefinition* tParticle)
480 {
481  G4double xsection = 0.;
482  // G4bool pORn = (tParticle == theProton || nucleon == theNeutron );
483  G4bool proton = (tParticle == theProton);
484  G4bool neutron = (tParticle == theNeutron);
485 
486  // General PDG fit constants
487 
488  G4double s0 = 5.38*5.38; // in Gev^2
489  G4double eta1 = 0.458;
490  G4double eta2 = 0.458;
491  G4double B = 0.308;
492 
493  // const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
494 
495  if(pParticle == theNeutron) // proton-neutron fit
496  {
497  if ( proton )
498  {
499  xsection = ( 35.80 + B*std::pow(std::log(sMand/s0),2.)
500  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
501  }
502  if ( neutron )
503  {
504  xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
505  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
506  }
507  }
508  else if(pParticle == theProton)
509  {
510  if ( proton )
511  {
512  xsection = (35.45 + B*std::pow(std::log(sMand/s0),2.)
513  + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
514 
515  }
516  if ( neutron )
517  {
518  xsection = (35.80 + B*std::pow(std::log(sMand/s0),2.)
519  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
520  }
521  }
522  xsection *= millibarn; // parametrised in mb
523  return xsection;
524 }
525 
526 
528 //
529 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
530 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
531 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
532 
533 G4double
535  G4double pTkin,
536  G4ParticleDefinition* tParticle)
537 {
538  G4double xsection(0);
539  // G4double Delta; DHW 19 May 2011: variable set but not used
540  G4double A0, B0;
541  G4double hpXscv(0);
542  G4double hnXscv(0);
543 
544  G4double targ_mass = tParticle->GetPDGMass();
545  G4double proj_mass = pParticle->GetPDGMass();
546 
547  G4double proj_energy = proj_mass + pTkin;
548  G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
549 
550  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
551 
552  sMand /= GeV*GeV; // in GeV for parametrisation
553  proj_momentum /= GeV;
554  proj_energy /= GeV;
555  proj_mass /= GeV;
556 
557  // General PDG fit constants
558 
559  // G4double s0 = 5.38*5.38; // in Gev^2
560  // G4double eta1 = 0.458;
561  // G4double eta2 = 0.458;
562  // G4double B = 0.308;
563 
564  if( proj_momentum >= 373.)
565  {
566  return GetHadronNucleonXscPDG(pParticle,sMand,tParticle);
567  }
568  else if( proj_momentum >= 10. ) // high energy: pp = nn = np
569  // if( proj_momentum >= 2.)
570  {
571  // Delta = 1.; // DHW 19 May 2011: variable set but not used
572  // if (proj_energy < 40.) Delta = 0.916+0.0021*proj_energy;
573 
574  if (proj_momentum >= 10.) {
575  B0 = 7.5;
576  A0 = 100. - B0*std::log(3.0e7);
577 
578  xsection = A0 + B0*std::log(proj_energy) - 11
579  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
580  0.93827*0.93827,-0.165); // mb
581  }
582  }
583  else // low energy pp = nn != np
584  {
585  if(pParticle == tParticle) // pp or nn // nn to be pp
586  {
587  if( proj_momentum < 0.73 )
588  {
589  hnXscv = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
590  }
591  else if( proj_momentum < 1.05 )
592  {
593  hnXscv = 23 + 40*(std::log(proj_momentum/0.73))*
594  (std::log(proj_momentum/0.73));
595  }
596  else // if( proj_momentum < 10. )
597  {
598  hnXscv = 39.0 +
599  75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
600  }
601  xsection = hnXscv;
602  }
603  else // pn to be np
604  {
605  if( proj_momentum < 0.8 )
606  {
607  hpXscv = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
608  }
609  else if( proj_momentum < 1.4 )
610  {
611  hpXscv = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
612  }
613  else // if( proj_momentum < 10. )
614  {
615  hpXscv = 33.3+
616  20.8*(std::pow(proj_momentum,2.0)-1.35)/
617  (std::pow(proj_momentum,2.50)+0.95);
618  }
619  xsection = hpXscv;
620  }
621  }
622  xsection *= millibarn; // parametrised in mb
623  return xsection;
624 }
625 
627 //
628 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
629 
630 G4double
632  G4int At, G4int Zt)
633 {
634  G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
635  G4int absPDGcode = std::abs(PDGcode);
636  G4double Elab = aParticle->GetTotalEnergy();
637  // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
638  G4double Plab = aParticle->GetMomentum().mag();
639  // std::sqrt(Elab * Elab - 0.88);
640 
641  Elab /= GeV;
642  Plab /= GeV;
643 
644  G4double LogPlab = std::log( Plab );
645  G4double sqrLogPlab = LogPlab * LogPlab;
646 
647  //G4cout<<"Plab = "<<Plab<<G4endl;
648 
649  G4double NumberOfTargetProtons = Zt;
650  G4double NumberOfTargetNucleons = At;
651  G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
652 
653  if(NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;
654 
655  G4double Xtotal = 0., Xelastic = 0., Xinelastic =0.;
656 
657  if( absPDGcode > 1000 ) //------Projectile is baryon --------
658  {
659  G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
660  0.522*sqrLogPlab - 4.51*LogPlab;
661 
662  G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
663  0.513*sqrLogPlab - 4.27*LogPlab;
664 
665  G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
666  0.169*sqrLogPlab - 1.85*LogPlab;
667 
668  G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
669  0.169*sqrLogPlab - 1.85*LogPlab;
670 
671  Xtotal = ( NumberOfTargetProtons * XtotPP +
672  NumberOfTargetNeutrons * XtotPN );
673 
674  Xelastic = ( NumberOfTargetProtons * XelPP +
675  NumberOfTargetNeutrons * XelPN );
676  }
677 
678  Xinelastic = Xtotal - Xelastic;
679  if(Xinelastic < 0.) Xinelastic = 0.;
680 
681  return Xinelastic*= millibarn;
682 }
683 
685 //
686 //
687 
688 G4double
690  const G4Element* anElement)
691 {
692  G4double At = anElement->GetN();
693  G4double oneThird = 1.0/3.0;
694  G4double cubicrAt = std::pow (At, oneThird);
695 
696  G4double R; // = fRadiusConst*cubicrAt;
697  R = fRadiusConst*cubicrAt;
698 
699  G4double meanA = 21.;
700  G4double tauA1 = 40.;
701  G4double tauA2 = 10.;
702  G4double tauA3 = 5.;
703 
704  G4double a1 = 0.85;
705  G4double b1 = 1. - a1;
706 
707  G4double b2 = 0.3;
708  G4double b3 = 4.;
709 
710  if (At > 20.) // 20.
711  {
712  R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
713  }
714  else if (At > 3.5)
715  {
716  R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
717  }
718  else
719  {
720  R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
721  }
722 
723  return R;
724 }
725 
727 //
728 //
729 
730 G4double
732 {
733  G4double R;
734  R = GetNucleusRadiusDE(Zt,At);
735  // R = GetNucleusRadiusRMS(Zt,At);
736 
737  return R;
738 }
739 
741 
742 G4double
744 {
745  G4double oneThird = 1.0/3.0;
746  G4double cubicrAt = std::pow (At, oneThird);
747 
748  G4double R; // = fRadiusConst*cubicrAt;
749  R = fRadiusConst*cubicrAt;
750 
751  G4double meanA = 20.;
752  G4double tauA = 20.;
753 
754  if ( At > 20.) // 20.
755  {
756  R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
757  }
758  else
759  {
760  R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
761  }
762 
763  return R;
764 }
765 
767 //
768 //
769 
770 G4double
772 {
773  // algorithm from diffuse-elastic
774 
775  G4double R, r0, a11, a12, a13, a2, a3;
776 
777  a11 = 1.26; // 1.08, 1.16
778  a12 = 1.; // 1.08, 1.16
779  a13 = 1.12; // 1.08, 1.16
780  a2 = 1.1;
781  a3 = 1.;
782 
783  // Special rms radii for light nucleii
784 
785  if (A < 50.)
786  {
787  if (std::abs(A-1.) < 0.5) return 0.89*fermi; // p
788  else if(std::abs(A-2.) < 0.5) return 2.13*fermi; // d
789  else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5) return 1.80*fermi; // t
790 
791  else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96*fermi; // He3
792  else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5) return 1.68*fermi; // He4
793 
794  else if(std::abs(Z-3.) < 0.5) return 2.40*fermi; // Li7
795  else if(std::abs(Z-4.) < 0.5) return 2.51*fermi; // Be9
796 
797  else if( 10. < A && A <= 16. ) r0 = a11*( 1 - std::pow(A, -2./3.) )*fermi; // 1.08*fermi;
798  else if( 15. < A && A <= 20. ) r0 = a12*( 1 - std::pow(A, -2./3.) )*fermi;
799  else if( 20. < A && A <= 30. ) r0 = a13*( 1 - std::pow(A, -2./3.) )*fermi;
800  else r0 = a2*fermi;
801 
802  R = r0*std::pow( A, 1./3. );
803  }
804  else
805  {
806  r0 = a3*fermi;
807 
808  R = r0*std::pow(A, 0.27);
809  }
810  return R;
811 }
812 
813 
815 //
816 // RMS radii from e-A scattering data
817 
818 G4double
820 {
821 
822  if (std::abs(A-1.) < 0.5) return 0.89*fermi; // p
823  else if(std::abs(A-2.) < 0.5) return 2.13*fermi; // d
824  else if(std::abs(Z-1.) < 0.5 && std::abs(A-3.) < 0.5) return 1.80*fermi; // t
825 
826  else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96*fermi; // He3
827  else if(std::abs(Z-2.) < 0.5 && std::abs(A-4.) < 0.5) return 1.68*fermi; // He4
828 
829  else if(std::abs(Z-3.) < 0.5) return 2.40*fermi; // Li7
830  else if(std::abs(Z-4.) < 0.5) return 2.51*fermi; // Be9
831 
832  else return 1.24*std::pow(A, 0.28 )*fermi; // A > 9
833 }
834 
835 
837 //
838 //
839 
841  const G4double mt,
842  const G4double Plab)
843 {
844  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
845  G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
846  // G4double Pcm = Plab * mt / Ecm;
847  // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
848 
849  return Ecm ; // KEcm;
850 }
851 
852 
854 //
855 //
856 
858  const G4double mt,
859  const G4double Plab)
860 {
861  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
862  G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
863 
864  return sMand;
865 }
866 
867 //
868 //
G4double GetNucleusRadiusGG(G4double At)
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4double GetN() const
Definition: G4Element.hh:134
G4double GetZ() const
Definition: G4Element.hh:131
G4double GetNucleusRadius(const G4DynamicParticle *, const G4Element *)
virtual G4double GetInelasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
G4ParticleDefinition * GetDefinition() const
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
int G4int
Definition: G4Types.hh:78
int millibarn
Definition: hepunit.py:40
static G4NistManager * Instance()
G4double GetHadronNucleonXscNS(const G4DynamicParticle *, const G4ParticleDefinition *)
G4IonTable * GetIonTable() const
G4double GetHadronNucleonXscPDG(G4ParticleDefinition *, G4double sMand, G4ParticleDefinition *)
virtual void CrossSectionDescription(std::ostream &) const
G4double CalcMandelstamS(const G4double, const G4double, const G4double)
Float_t Z
Definition: plot.C:39
G4double GetMass() const
bool G4bool
Definition: G4Types.hh:79
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1232
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double GetHadronNucleonXsc(const G4DynamicParticle *, const G4Element *)
G4double GetNucleusRadiusRMS(G4double Zt, G4double At)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetRatioQE(const G4DynamicParticle *, G4double At, G4double Zt)
void SetKineticEnergy(G4double aEnergy)
G4double CalculateEcmValue(const G4double, const G4double, const G4double)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
int G4lrint(double ad)
Definition: templates.hh:163
G4double GetRatioSD(const G4DynamicParticle *, G4double At, G4double Zt)
virtual G4double GetTotalIsotopeCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
G4double GetHadronNucleonXscNS(G4ParticleDefinition *, G4double pTkin, G4ParticleDefinition *)
G4double GetHNinelasticXscVU(const G4DynamicParticle *, G4int At, G4int Zt)
G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
virtual G4double GetInelasticIsotopeCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
G4double GetPDGCharge() const
G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *)
G4double GetNucleusRadiusDE(G4double Zt, G4double At)
double mag() const
const G4double oneThird
G4double GetCoulombBarier(const G4DynamicParticle *, G4double Z, G4double A, G4double pR, G4double tR)
virtual G4double GetElasticIsotopeCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
G4double GetInelasticHadronNucleonXsc()
G4ThreeVector GetMomentum() const
virtual G4double ComputeQuasiElasticRatio(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4int A)
G4double GetZandACrossSection(const G4DynamicParticle *, G4int Z, G4int A)