source: trunk/source/processes/hadronic/cross_sections/src/G4GGNuclNuclCrossSection.cc @ 990

Last change on this file since 990 was 968, checked in by garnier, 15 years ago

fichier ajoutes
File size: 26.6 KB
Line 
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// 24.11.08 V. Grichine - first implementation
28//
29
30#include "G4GGNuclNuclCrossSection.hh"
31
32#include "G4ParticleTable.hh"
33#include "G4IonTable.hh"
34#include "G4ParticleDefinition.hh"
35
36
37
38////////////////////////////////////////////////////////////////////////////////
39//
40//
41
42G4GGNuclNuclCrossSection::G4GGNuclNuclCrossSection() 
43: fUpperLimit( 100000 * GeV ),
44  fLowerLimit( 0.1 * GeV ),
45  fRadiusConst( 1.08*fermi )  // 1.1, 1.3 ?
46{
47  theProton   = G4Proton::Proton();
48  theNeutron  = G4Neutron::Neutron();
49}
50
51///////////////////////////////////////////////////////////////////////////////////////
52//
53//
54
55G4GGNuclNuclCrossSection::~G4GGNuclNuclCrossSection()
56{
57}
58
59
60////////////////////////////////////////////////////////////////////////////////////////
61//
62//
63
64
65G4bool
66G4GGNuclNuclCrossSection::IsApplicable(const G4DynamicParticle* aDP, 
67                                          const G4Element*  anElement)
68{
69  return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN());
70} 
71
72////////////////////////////////////////////////////////////////////////////////////////
73//
74//
75
76G4bool
77G4GGNuclNuclCrossSection::IsZAApplicable(const G4DynamicParticle* aDP, 
78                                            G4double Z, G4double)
79{
80  G4bool applicable      = false;
81  // G4int baryonNumber     = aDP->GetDefinition()->GetBaryonNumber();
82  G4double kineticEnergy = aDP->GetKineticEnergy();
83
84  //  const G4ParticleDefinition* theParticle = aDP->GetDefinition();
85 
86  if ( kineticEnergy  >= fLowerLimit && Z > 1.5 ) applicable = true;
87
88  return applicable;
89}
90
91////////////////////////////////////////////////////////////////////////////////////////
92//
93// Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
94// Glauber model with Gribov correction calculated in the dipole approximation on
95// light cone. Gaussian density helps to calculate rest integrals of the model.
96// [1] B.Z. Kopeliovich, nucl-th/0306044
97
98
99G4double G4GGNuclNuclCrossSection::
100GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement, G4double T)
101{
102  return GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(), T);
103}
104
105////////////////////////////////////////////////////////////////////////////////////////
106//
107// Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
108// Glauber model with Gribov correction calculated in the dipole approximation on
109// light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model.
110// [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
111
112
113
114G4double G4GGNuclNuclCrossSection::
115GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double tZ, G4double tA, G4double)
116{
117  G4double xsection, sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
118
119  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
120  G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
121
122  G4double pTkin = aParticle->GetKineticEnergy(); 
123  pTkin /= pA;
124
125  G4double pN = pA - pZ;
126  if( pN < 0. ) pN = 0.;
127
128  G4double tN = tA - tZ;
129  if( tN < 0. ) tN = 0.;
130
131  G4double tR = GetNucleusRadius(tA); 
132  G4double pR = GetNucleusRadius(pA); 
133
134  sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
135          (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
136
137  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );   // basically 2piRR
138
139  ratio = sigma/nucleusSquare;
140
141  xsection =  nucleusSquare*std::log( 1. + ratio );
142
143  fTotalXsc = xsection;
144
145  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
146
147  fElasticXsc   = fTotalXsc - fInelasticXsc;
148   
149  G4double difratio = ratio/(1.+ratio);
150
151  fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
152
153  // production to be checked !!! edit MK xsc
154
155  sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscMK(theProton, pTkin, theProton) +
156          (pZ*tN+pN*tZ)*GetHadronNucleonXscMK(theProton, pTkin, theNeutron);
157 
158  ratio = sigma/nucleusSquare;
159
160  fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
161
162  if (fElasticXsc < 0.) fElasticXsc = 0.;
163
164  return xsection; 
165}
166
167//////////////////////////////////////////////////////////////////////////
168//
169// Return single-diffraction/inelastic cross-section ratio
170
171G4double G4GGNuclNuclCrossSection::
172GetRatioSD(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
173{
174  G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
175
176  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
177  G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
178
179  G4double pTkin = aParticle->GetKineticEnergy(); 
180  pTkin /= pA;
181
182  G4double pN = pA - pZ;
183  if( pN < 0. ) pN = 0.;
184
185  G4double tN = tA - tZ;
186  if( tN < 0. ) tN = 0.;
187
188  G4double tR = GetNucleusRadius(tA); 
189  G4double pR = GetNucleusRadius(pA); 
190
191  sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
192          (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
193
194  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );   // basically 2piRR
195
196  ratio = sigma/nucleusSquare;
197
198
199  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
200   
201  G4double difratio = ratio/(1.+ratio);
202
203  fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
204
205  if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
206  else                    ratio = 0.;
207
208  return ratio; 
209}
210
211//////////////////////////////////////////////////////////////////////////
212//
213// Return suasi-elastic/inelastic cross-section ratio
214
215G4double G4GGNuclNuclCrossSection::
216GetRatioQE(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
217{
218  G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
219
220  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
221  G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
222
223  G4double pTkin = aParticle->GetKineticEnergy(); 
224  pTkin /= pA;
225
226  G4double pN = pA - pZ;
227  if( pN < 0. ) pN = 0.;
228
229  G4double tN = tA - tZ;
230  if( tN < 0. ) tN = 0.;
231
232  G4double tR = GetNucleusRadius(tA); 
233  G4double pR = GetNucleusRadius(pA); 
234
235  sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
236          (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
237
238  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );   // basically 2piRR
239
240  ratio = sigma/nucleusSquare;
241
242  fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
243
244  //  sigma = GetHNinelasticXsc(aParticle, tA, tZ);
245  ratio = sigma/nucleusSquare;
246
247  fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
248
249  if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
250  else                                ratio = 0.;
251  if ( ratio < 0. )                   ratio = 0.;
252
253  return ratio; 
254}
255
256/////////////////////////////////////////////////////////////////////////////////////
257//
258// Returns hadron-nucleon Xsc according to differnt parametrisations:
259// [2] E. Levin, hep-ph/9710546
260// [3] U. Dersch, et al, hep-ex/9910052
261// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
262
263G4double
264G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle, 
265                                                  const G4Element* anElement          )
266{
267  G4double At = anElement->GetN();  // number of nucleons
268  G4double Zt = anElement->GetZ();  // number of protons
269
270
271  return GetHadronNucleonXsc( aParticle, At, Zt );
272}
273
274
275
276
277/////////////////////////////////////////////////////////////////////////////////////
278//
279// Returns hadron-nucleon Xsc according to differnt parametrisations:
280// [2] E. Levin, hep-ph/9710546
281// [3] U. Dersch, et al, hep-ex/9910052
282// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
283
284G4double
285G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle, 
286                                                   G4double At,  G4double Zt       )
287{
288  G4double xsection = 0.;
289
290
291  G4double targ_mass = G4ParticleTable::GetParticleTable()->
292  GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
293
294  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
295
296  G4double proj_mass     = aParticle->GetMass();
297  G4double proj_momentum = aParticle->GetMomentum().mag();
298  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
299
300  sMand /= GeV*GeV;  // in GeV for parametrisation
301  proj_momentum /= GeV;
302
303  const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
304 
305
306  if(pParticle == theNeutron) // as proton ???
307  {
308    xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
309  } 
310  else if(pParticle == theProton) 
311  {
312    xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
313    // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
314    // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
315  } 
316 
317  xsection *= millibarn;
318  return xsection;
319}
320
321
322/////////////////////////////////////////////////////////////////////////////////////
323//
324// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
325// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
326
327G4double
328G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle, 
329                                                  const G4Element* anElement          )
330{
331  G4double At = anElement->GetN();  // number of nucleons
332  G4double Zt = anElement->GetZ();  // number of protons
333
334
335  return GetHadronNucleonXscPDG( aParticle, At, Zt );
336}
337
338
339
340
341/////////////////////////////////////////////////////////////////////////////////////
342//
343// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
344// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
345//  At = number of nucleons,  Zt = number of protons
346
347G4double
348G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle, 
349                                                     G4double At,  G4double Zt )
350{
351  G4double xsection = 0.;
352
353  G4double Nt = At-Zt;              // number of neutrons
354  if (Nt < 0.) Nt = 0.; 
355
356
357  G4double targ_mass = G4ParticleTable::GetParticleTable()->
358  GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
359
360  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
361
362  G4double proj_mass     = aParticle->GetMass(); 
363  G4double proj_momentum = aParticle->GetMomentum().mag();
364
365  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
366
367  sMand         /= GeV*GeV;  // in GeV for parametrisation
368
369  // General PDG fit constants
370
371  G4double s0   = 5.38*5.38; // in Gev^2
372  G4double eta1 = 0.458;
373  G4double eta2 = 0.458;
374  G4double B    = 0.308;
375
376
377  const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
378 
379
380  if(pParticle == theNeutron) // proton-neutron fit
381  {
382    xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
383                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
384    xsection  += Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
385                      + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
386  } 
387  else if(pParticle == theProton) 
388  {
389     
390      xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.) 
391                          + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
392
393      xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.) 
394                          + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
395  } 
396  xsection *= millibarn; // parametrised in mb
397  return xsection;
398}
399
400
401
402
403/////////////////////////////////////////////////////////////////////////////////////
404//
405// Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
406// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
407// projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
408
409G4double
410G4GGNuclNuclCrossSection::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle, 
411                                                 G4double pTkin, 
412                                                 G4ParticleDefinition* tParticle)
413{
414  G4double xsection(0), Delta, A0, B0;
415  G4double hpXsc(0);
416  G4double hnXsc(0);
417
418
419  G4double targ_mass     = tParticle->GetPDGMass();
420  G4double proj_mass     = pParticle->GetPDGMass(); 
421
422  G4double proj_energy   = proj_mass + pTkin; 
423  G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
424
425  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
426
427  sMand         /= GeV*GeV;  // in GeV for parametrisation
428  proj_momentum /= GeV;
429  proj_energy   /= GeV;
430  proj_mass     /= GeV;
431
432  // General PDG fit constants
433
434  //  G4double s0   = 5.38*5.38; // in Gev^2
435  //  G4double eta1 = 0.458;
436  //  G4double eta2 = 0.458;
437  //  G4double B    = 0.308;
438
439 
440
441
442 
443  if( proj_momentum >= 10. ) // high energy: pp = nn = np
444    // if( proj_momentum >= 2.)
445  {
446    Delta = 1.;
447
448    if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
449
450    if( proj_momentum >= 10.)
451    {
452        B0 = 7.5;
453        A0 = 100. - B0*std::log(3.0e7);
454
455        xsection = A0 + B0*std::log(proj_energy) - 11
456                  + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
457                     0.93827*0.93827,-0.165);        //  mb
458    }
459  }
460  else // low energy pp = nn != np
461  {
462      if(pParticle == tParticle) // pp or nn      // nn to be pp
463      {
464        if( proj_momentum < 0.73 )
465        {
466          hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
467        }
468        else if( proj_momentum < 1.05  )
469        {
470          hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
471                         (std::log(proj_momentum/0.73));
472        }
473        else  // if( proj_momentum < 10.  )
474        {
475          hnXsc = 39.0 + 
476              75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
477        }
478        xsection = hnXsc;
479      }
480      else  // pn to be np
481      {
482        if( proj_momentum < 0.8 )
483        {
484          hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
485        }     
486        else if( proj_momentum < 1.4 )
487        {
488          hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
489        }
490        else    // if( proj_momentum < 10.  )
491        {
492          hpXsc = 33.3+
493              20.8*(std::pow(proj_momentum,2.0)-1.35)/
494                 (std::pow(proj_momentum,2.50)+0.95);
495        }
496        xsection = hpXsc;
497      }
498  }
499  xsection *= millibarn; // parametrised in mb
500  return xsection;
501}
502
503/*
504/////////////////////////////////////////////////////////////////////////////////////
505//
506// Returns hadron-nucleon inelastic cross-section based on proper parametrisation
507
508G4double
509G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
510                                                  const G4Element* anElement          )
511{
512  G4double At = anElement->GetN();  // number of nucleons
513  G4double Zt = anElement->GetZ();  // number of protons
514
515
516  return GetHNinelasticXsc( aParticle, At, Zt );
517}
518
519/////////////////////////////////////////////////////////////////////////////////////
520//
521// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
522
523G4double
524G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
525                                                     G4double At,  G4double Zt )
526{
527  //  G4ParticleDefinition* hadron = aParticle->GetDefinition();
528  G4double sumInelastic, Nt = At - Zt;
529
530  if(Nt < 0.) Nt = 0.;
531 
532  sumInelastic  = Zt*GetHadronNucleonXscMK(aParticle, theProton);
533  sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);   
534 
535  return sumInelastic;
536}
537*/
538
539/////////////////////////////////////////////////////////////////////////////////////
540//
541// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
542
543G4double
544G4GGNuclNuclCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle, 
545                                                     G4double At,  G4double Zt )
546{
547  G4int PDGcode    = aParticle->GetDefinition()->GetPDGEncoding();
548  G4int absPDGcode = std::abs(PDGcode);
549
550  G4double Elab = aParticle->GetTotalEnergy();             
551                          // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
552  G4double Plab = aParticle->GetMomentum().mag();           
553                          // std::sqrt(Elab * Elab - 0.88);
554
555  Elab /= GeV;
556  Plab /= GeV;
557
558  G4double LogPlab    = std::log( Plab );
559  G4double sqrLogPlab = LogPlab * LogPlab;
560
561  //G4cout<<"Plab = "<<Plab<<G4endl;
562
563  G4double NumberOfTargetProtons  = Zt; 
564  G4double NumberOfTargetNucleons = At;
565  G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
566
567  if(NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;
568
569  G4double Xtotal = 0., Xelastic = 0., Xinelastic =0.;
570
571  if( absPDGcode > 1000 )  //------Projectile is baryon --------
572  {
573       G4double XtotPP = 48.0 +  0. *std::pow(Plab, 0.  ) +
574                         0.522*sqrLogPlab - 4.51*LogPlab;
575
576       G4double XtotPN = 47.3 +  0. *std::pow(Plab, 0.  ) +
577                         0.513*sqrLogPlab - 4.27*LogPlab;
578
579       G4double XelPP  = 11.9 + 26.9*std::pow(Plab,-1.21) +
580                         0.169*sqrLogPlab - 1.85*LogPlab;
581
582       G4double XelPN  = 11.9 + 26.9*std::pow(Plab,-1.21) +
583                         0.169*sqrLogPlab - 1.85*LogPlab;
584
585       Xtotal          = ( NumberOfTargetProtons  * XtotPP +
586                           NumberOfTargetNeutrons * XtotPN  );
587
588       Xelastic        = ( NumberOfTargetProtons  * XelPP  +
589                           NumberOfTargetNeutrons * XelPN   );
590  }
591  Xinelastic = Xtotal - Xelastic;
592
593  if(Xinelastic < 0.) Xinelastic = 0.;
594
595  return Xinelastic*= millibarn;
596}
597
598/////////////////////////////////////////////////////////////////////////////////////
599//
600// Returns hadron-nucleon cross-section based on Mikhail Kossov CHIPS parametrisation of
601// data from G4QuasiFreeRatios class
602
603G4double
604G4GGNuclNuclCrossSection::GetHadronNucleonXscMK(G4ParticleDefinition* pParticle, G4double pTkin,
605                                                G4ParticleDefinition* nucleon  )
606{
607  G4int I = -1;
608  G4int PDG = pParticle->GetPDGEncoding();
609  G4double totalXsc = 0;
610  G4double elasticXsc = 0;
611  G4double inelasticXsc;
612  // G4int absPDG = std::abs(PDG);
613
614  G4double pM = pParticle->GetPDGMass();
615  G4double p  = std::sqrt(pTkin*(pTkin+2*pM))/GeV;
616
617  G4bool F = false;           
618  if(nucleon == theProton)       F = true;
619  else if(nucleon == theNeutron) F = false;
620  else
621  {
622    G4cout << "nucleon is not proton or neutron, return xsc for proton" << G4endl;
623    F = true;
624  }
625
626  G4bool kfl = true;                             // Flag of K0/aK0 oscillation
627  G4bool kf  = false;
628
629  if( PDG == 130 || PDG == 310 )
630  {
631    kf = true;
632    if( G4UniformRand() > .5 ) kfl = false;
633  }
634  if     ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) I = 0; // pp/nn
635  else if( (PDG == 2112 && F) || (PDG == 2212 && !F) ) I = 1; // np/pn
636  else
637  {
638    G4cout<<"MK PDG = "<<PDG
639          <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
640    G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash");
641  }
642
643  // Each parameter set can have not more than nPoints = 128 parameters
644
645  static const G4double lmi = 3.5;       // min of (lnP-lmi)^2 parabola
646  static const G4double pbe = .0557;     // elastic (lnP-lmi)^2 parabola coefficient
647  static const G4double pbt = .3;        // total (lnP-lmi)^2 parabola coefficient
648  static const G4double pmi = .1;        // Below that fast LE calculation is made
649  static const G4double pma = 1000.;     // Above that fast HE calculation is made
650                 
651  if( p <= 0.)
652  {
653    G4cout<<" p = "<<p<<" is zero or negative"<<G4endl;
654
655    elasticXsc   = 0.;
656    inelasticXsc = 0.;
657    totalXsc     = 0.;
658
659    return totalXsc;
660  }
661  if (!I)                          // pp/nn
662  {
663    if( p < pmi )
664    {
665      G4double p2 = p*p;
666      elasticXsc          = 1./(.00012 + p2*.2);
667      totalXsc          = elasticXsc;
668    }
669    else if(p>pma)
670    {
671      G4double lp  = std::log(p)-lmi;
672      G4double lp2 = lp*lp;
673      elasticXsc  = pbe*lp2 + 6.72;
674      totalXsc    = pbt*lp2 + 38.2;
675    }
676    else
677    {
678      G4double p2  = p*p;
679      G4double LE  = 1./( .00012 + p2*.2);
680      G4double lp  = std::log(p) - lmi;
681      G4double lp2 = lp*lp;
682      G4double rp2 = 1./p2;
683      elasticXsc  = LE + ( pbe*lp2 + 6.72+32.6/p)/( 1. + rp2/p);
684      totalXsc    = LE + ( pbt*lp2 + 38.2+52.7*rp2)/( 1. + 2.72*rp2*rp2);
685    }
686  }
687  else if( I==1 )                        // np/pn
688  {
689    if( p < pmi )
690    {
691      G4double p2 = p*p;
692      elasticXsc = 1./( .00012 + p2*( .051 + .1*p2));
693      totalXsc   = elasticXsc;
694    }
695    else if( p > pma )
696    {
697      G4double lp  = std::log(p) - lmi;
698      G4double lp2 = lp*lp;
699      elasticXsc  = pbe*lp2 + 6.72;
700      totalXsc    = pbt*lp2 + 38.2;
701    }
702    else
703    {
704      G4double p2  = p*p;
705      G4double LE  = 1./( .00012 + p2*( .051 + .1*p2 ) );
706      G4double lp  = std::log(p) - lmi;
707      G4double lp2 = lp*lp;
708      G4double rp2 = 1./p2;
709      elasticXsc  = LE + (pbe*lp2 + 6.72 + 30./p)/( 1. + .49*rp2/p);
710      totalXsc    = LE + (pbt*lp2 + 38.2)/( 1. + .54*rp2*rp2);
711    }
712  }
713  else
714  {
715    G4cout<<"PDG incoding = "<<I<<" is not defined (0-1)"<<G4endl;
716 
717  }
718  if( elasticXsc > totalXsc ) elasticXsc = totalXsc;
719
720  totalXsc   *= millibarn;
721  elasticXsc *= millibarn;
722  inelasticXsc   = totalXsc - elasticXsc;
723  if (inelasticXsc < 0.) inelasticXsc = 0.;
724
725  return inelasticXsc;
726}
727
728////////////////////////////////////////////////////////////////////////////////////
729//
730//
731
732G4double
733G4GGNuclNuclCrossSection::GetNucleusRadius( const G4DynamicParticle* , 
734                                               const G4Element* anElement)
735{
736  G4double At       = anElement->GetN();
737  G4double oneThird = 1.0/3.0;
738  G4double cubicrAt = std::pow (At, oneThird); 
739
740
741  G4double R;  // = fRadiusConst*cubicrAt;
742  /* 
743  G4double tmp = std::pow( cubicrAt-1., 3.);
744  tmp         += At;
745  tmp         *= 0.5;
746
747  if (At > 20.)   // 20.
748  {
749    R = fRadiusConst*std::pow (tmp, oneThird);
750  }
751  else
752  {
753    R = fRadiusConst*cubicrAt;
754  }
755  */
756 
757  R = fRadiusConst*cubicrAt;
758
759  // return R;  // !!!!
760
761
762 
763  G4double meanA  = 21.;
764
765  G4double tauA1  = 40.; 
766  G4double tauA2  = 10.; 
767  G4double tauA3  = 5.; 
768
769  G4double a1 = 0.85;
770  G4double b1 = 1. - a1;
771
772  G4double b2 = 0.3;
773  G4double b3 = 4.;
774
775  if (At > 20.)   // 20.
776  {
777    R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) ); 
778  }
779  else if (At > 3.5)
780  {
781    R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) ); 
782  }
783  else 
784  {
785    R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) ); 
786  } 
787  return R;
788 
789}
790
791////////////////////////////////////////////////////////////////////////////////////
792//
793//
794
795G4double
796G4GGNuclNuclCrossSection::GetNucleusRadius(G4double At)
797{
798  G4double R;
799
800  // R = GetNucleusRadiusGG(At);
801
802  R = GetNucleusRadiusDE(At);
803
804  return R;
805}
806
807///////////////////////////////////////////////////////////////////
808
809G4double
810G4GGNuclNuclCrossSection::GetNucleusRadiusGG(G4double At)
811{
812 
813  G4double oneThird = 1.0/3.0;
814  G4double cubicrAt = std::pow (At, oneThird); 
815
816
817  G4double R;  // = fRadiusConst*cubicrAt;
818 
819  /*
820  G4double tmp = std::pow( cubicrAt-1., 3.);
821  tmp         += At;
822  tmp         *= 0.5;
823
824  if (At > 20.)
825  {
826    R = fRadiusConst*std::pow (tmp, oneThird);
827  }
828  else
829  {
830    R = fRadiusConst*cubicrAt;
831  }
832  */
833 
834  R = fRadiusConst*cubicrAt;
835
836  G4double meanA = 20.;
837  G4double tauA  = 20.; 
838
839  if ( At > 20.)   // 20.
840  {
841    R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) ); 
842  }
843  else
844  {
845    R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) ); 
846  }
847
848  return R;
849
850}
851
852
853G4double
854G4GGNuclNuclCrossSection::GetNucleusRadiusDE(G4double A)
855{
856
857  // algorithm from diffuse-elastic
858
859  G4double R, r0, a11, a12, a13, a2, a3;
860
861  a11 = 1.26;  // 1.08, 1.16
862  a12 = 1.;  // 1.08, 1.16
863  a13 = 1.12;  // 1.08, 1.16
864  a2 = 1.1;
865  a3 = 1.;
866
867
868  if( A < 50. )
869  {
870    if( 10  < A && A <= 15. )     r0  = a11*( 1 - std::pow(A, -2./3.) )*fermi;   // 1.08*fermi;
871    else if( 15  < A && A <= 20 ) r0  = a12*( 1 - std::pow(A, -2./3.) )*fermi;
872    else if( 20  < A && A <= 30 ) r0  = a13*( 1 - std::pow(A, -2./3.) )*fermi;
873    else                          r0  = a2*fermi;
874
875    R = r0*std::pow( A, 1./3. );
876  }
877  else
878  {
879    r0 = a3*fermi;
880
881    R  = r0*std::pow(A, 0.27);
882  }
883  return R;
884
885
886
887}
888
889
890
891
892
893
894////////////////////////////////////////////////////////////////////////////////////
895//
896//
897
898G4double G4GGNuclNuclCrossSection::CalculateEcmValue( const G4double mp , 
899                                                         const G4double mt , 
900                                                         const G4double Plab )
901{
902  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
903  G4double Ecm  = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
904  // G4double Pcm  = Plab * mt / Ecm;
905  // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
906
907  return Ecm ; // KEcm;
908}
909
910
911////////////////////////////////////////////////////////////////////////////////////
912//
913//
914
915G4double G4GGNuclNuclCrossSection::CalcMandelstamS( const G4double mp , 
916                                                       const G4double mt , 
917                                                       const G4double Plab )
918{
919  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
920  G4double sMand  = mp*mp + mt*mt + 2*Elab*mt ;
921
922  return sMand;
923}
924
925
926//
927//
928///////////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.