Geant4: G4PreCompoundProton Class Reference

G4PreCompoundProton::G4PreCompoundProton ( )

Definition at line 50 of file G4PreCompoundProton.cc.

References G4VPreCompoundFragment::GetA(), G4VPreCompoundFragment::GetRestA(), G4VPreCompoundFragment::GetRestZ(), G4VPreCompoundFragment::GetZ(), and G4VPreCompoundFragment::ResidualA13().

00051   : G4PreCompoundNucleon(G4Proton::Proton(), &theProtonCoulombBarrier)
00052 {
00053   ResidualA = GetRestA();
00054   ResidualZ = GetRestZ(); 
00055   theA = GetA();
00056   theZ = GetZ();
00057   ResidualAthrd = ResidualA13();
00058   FragmentAthrd = ResidualAthrd;
00059   FragmentA = theA + ResidualA;
00060 }

G4PreCompoundProton::~G4PreCompoundProton ( ) [virtual]

Definition at line 62 of file G4PreCompoundProton.cc.

00063 {}

G4double G4PreCompoundProton::CrossSection ( G4double ekin ) [protected, virtual]

Implements G4PreCompoundNucleon.

Definition at line 81 of file G4PreCompoundProton.cc.

References G4endl, G4VPreCompoundFragment::g4pow, G4VPreCompoundFragment::GetA(), G4PreCompoundNucleon::GetOpt0(), GetOpt1(), GetOpt2(), GetOpt3(), G4VPreCompoundFragment::GetRestA(), G4VPreCompoundFragment::GetRestZ(), G4VPreCompoundFragment::GetZ(), G4VPreCompoundFragment::OPTxs, G4VPreCompoundFragment::ResidualA13(), and G4Pow::Z13().

00082 {
00083   ResidualA = GetRestA();
00084   ResidualZ = GetRestZ(); 
00085   theA = GetA();
00086   theZ = GetZ();
00087   ResidualAthrd = ResidualA13();
00088   FragmentA = theA + ResidualA;
00089   FragmentAthrd = g4pow->Z13(FragmentA);
00090 
00091   if (OPTxs==0) { return GetOpt0(K); }
00092   else if( OPTxs==1) { return GetOpt1(K); }
00093   else if( OPTxs==2|| OPTxs==4) { return GetOpt2(K); }
00094   else if (OPTxs==3)  { return GetOpt3(K); }
00095   else{
00096     std::ostringstream errOs;
00097     errOs << "BAD PROTON CROSS SECTION OPTION !!"  <<G4endl;
00098     throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00099     return 0.;
00100   }
00101 }

G4double G4PreCompoundProton::GetAlpha ( ) [protected, virtual]

Implements G4PreCompoundNucleon.

Definition at line 103 of file G4PreCompoundProton.cc.

00104 {
00105   G4int aZ = ResidualZ;
00106   G4double C = 0.0;
00107   if (aZ >= 70) 
00108     {
00109       C = 0.10;
00110     } 
00111   else 
00112     {
00113       C = ((((0.15417e-06*aZ) - 0.29875e-04)*aZ + 0.21071e-02)*aZ - 0.66612e-01)*aZ + 0.98375;
00114     }
00115   return 1.0 + C;
00116 }

G4double G4PreCompoundProton::GetBeta ( ) [protected, virtual]

Implements G4PreCompoundNucleon.

Definition at line 118 of file G4PreCompoundProton.cc.

References G4VPreCompoundFragment::GetCoulombBarrier().

00119 {
00120   return -GetCoulombBarrier();
00121 }

G4double G4PreCompoundProton::GetOpt1 ( G4double K ) [protected]

Definition at line 126 of file G4PreCompoundProton.cc.

References G4VPreCompoundFragment::g4pow, and G4Pow::powZ().

Referenced by CrossSection().

00127 {
00128   G4double Kc=K; 
00129 
00130   // JMQ  xsec is set constat above limit of validity
00131   if (K > 50*MeV) { Kc = 50*MeV; }
00132 
00133   G4double landa, landa0, landa1, mu, mm0, mu1,nu, nu0, nu1, nu2,xs;
00134   G4double p, p0, p1, p2,Ec,delta,q,r,ji;
00135   
00136   p0 = 15.72;
00137   p1 = 9.65;
00138   p2 = -449.0;
00139   landa0 = 0.00437;
00140   landa1 = -16.58;
00141   mm0 = 244.7;
00142   mu1 = 0.503;
00143   nu0 = 273.1;
00144   nu1 = -182.4;
00145   nu2 = -1.872;  
00146   delta=0.;  
00147 
00148   Ec = 1.44*theZ*ResidualZ/(1.5*ResidualAthrd+delta);
00149   p = p0 + p1/Ec + p2/(Ec*Ec);
00150   landa = landa0*ResidualA + landa1;
00151 
00152   G4double resmu1 = g4pow->powZ(ResidualA,mu1); 
00153   mu = mm0*resmu1;
00154   nu = resmu1*(nu0 + nu1*Ec + nu2*(Ec*Ec));
00155   q = landa - nu/(Ec*Ec) - 2*p*Ec;
00156   r = mu + 2*nu/Ec + p*(Ec*Ec);
00157 
00158   ji=std::max(Kc,Ec);
00159   if(Kc < Ec) { xs = p*Kc*Kc + q*Kc + r;}
00160   else {xs = p*(Kc - ji)*(Kc - ji) + landa*Kc + mu + nu*(2 - Kc/ji)/ji ;}
00161   if (xs <0.0) {xs=0.0;}
00162 
00163   return xs; 
00164 }

G4double G4PreCompoundProton::GetOpt2 ( G4double K ) [protected]

Definition at line 168 of file G4PreCompoundProton.cc.

References G4cout, G4endl, G4VPreCompoundFragment::g4pow, G4Pow::logZ(), G4INCL::Math::pi, G4VPreCompoundFragment::theCoulombBarrier, and G4VPreCompoundFragment::useSICB.

Referenced by CrossSection().

00169 {
00170 
00171   G4double eekin,ekin,ff1,ff2,ff3,r0,fac,fac1,fac2,b0,xine_th(0);
00172  
00173   // This is redundant when the Coulomb  barrier is overimposed to all 
00174   // cross sections 
00175   // It should be kept when Coulomb barrier only imposed at OPTxs=2
00176 
00177   if(!useSICB && K<=theCoulombBarrier) { return 0.0; }
00178 
00179   eekin=K;
00180   G4int rnneu=ResidualA-ResidualZ;
00181   ekin=eekin/1000;
00182   r0=1.36*1.e-15;
00183   fac=pi*r0*r0;
00184   b0=2.247-0.915*(1.-1./ResidualAthrd);
00185   fac1=b0*(1.-1./ResidualAthrd);
00186   fac2=1.;
00187   if(rnneu > 1.5) { fac2 = g4pow->logZ(rnneu); }
00188   xine_th= 1.e+31*fac*fac2*(1.+ResidualAthrd-fac1);
00189   xine_th=(1.-0.15*std::exp(-ekin))*xine_th/(1.00-0.0007*ResidualA);    
00190   ff1=0.70-0.0020*ResidualA;
00191   ff2=1.00+1/G4double(ResidualA);
00192   ff3=0.8+18/G4double(ResidualA)-0.002*ResidualA;
00193   fac=1.-(1./(1.+std::exp(-8.*ff1*(std::log10(ekin)+1.37*ff2))));
00194   xine_th=xine_th*(1.+ff3*fac);
00195   ff1=1.-1/G4double(ResidualA)-0.001*ResidualA;
00196   ff2=1.17-2.7/G4double(ResidualA)-0.0014*ResidualA;
00197   fac=-8.*ff1*(std::log10(ekin)+2.0*ff2);
00198   fac=1./(1.+std::exp(fac));
00199   xine_th=xine_th*fac;            
00200   if (xine_th < 0.0){
00201     std::ostringstream errOs;
00202     G4cout<<"WARNING:  negative Wellisch cross section "<<G4endl; 
00203     errOs << "RESIDUAL: A=" << ResidualA << " Z=" << ResidualZ <<G4endl;
00204     errOs <<"  xsec("<<ekin<<" MeV) ="<<xine_th <<G4endl;
00205     throw G4HadronicException(__FILE__, __LINE__, errOs.str());
00206   }
00207   return xine_th;
00208 }

G4double G4PreCompoundProton::GetOpt3 ( G4double K ) [protected]

Definition at line 211 of file G4PreCompoundProton.cc.

References G4VPreCompoundFragment::g4pow, and G4Pow::powZ().

Referenced by CrossSection().

00212 {
00213   //     ** p from  becchetti and greenlees (but modified with sub-barrier
00214   //     ** correction function and xp2 changed from -449)
00215 
00216   G4double landa, landa0, landa1, mu, mm0, mu1,nu, nu0, nu1, nu2;
00217   G4double p, p0, p1, p2;
00218   p0 = 15.72;
00219   p1 = 9.65;
00220   p2 = -300.;
00221   landa0 = 0.00437;
00222   landa1 = -16.58;
00223   mm0 = 244.7;
00224   mu1 = 0.503;
00225   nu0 = 273.1;
00226   nu1 = -182.4;
00227   nu2 = -1.872;
00228   
00229   // parameters for  proton cross section refinement 
00230   /*
00231   G4double afit,bfit,a2,b2;
00232   afit=-0.0785656;
00233   bfit=5.10789;
00234   a2= -0.00089076;
00235   b2= 0.0231597;  
00236   */
00237 
00238   G4double ec,ecsq,xnulam,etest(0.),ra(0.),a,w,c,signor(1.),signor2,sig; 
00239   G4double b,ecut,cut,ecut2,geom,elab;
00240     
00241   G4double      flow = 1.e-18;
00242   G4double       spill= 1.e+18; 
00243    
00244   if (ResidualA <= 60)      { signor = 0.92; }
00245   else if (ResidualA < 100) { signor = 0.8 + ResidualA*0.002; }
00246   
00247   ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
00248   ecsq = ec * ec;
00249   p = p0 + p1/ec + p2/ecsq;
00250   landa = landa0*ResidualA + landa1;
00251   a = g4pow->powZ(ResidualA,mu1);
00252   mu = mm0 * a;
00253   nu = a* (nu0+nu1*ec+nu2*ecsq);
00254   
00255   c =std::min(3.15,ec*0.5);
00256   w = 0.7 * c / 3.15; 
00257   
00258   xnulam = nu / landa;
00259   if (xnulam > spill) { xnulam=0.; }
00260   if (xnulam >= flow) { etest =std::sqrt(xnulam) + 7.; }
00261   
00262   a = -2.*p*ec + landa - nu/ecsq;
00263   b = p*ecsq + mu + 2.*nu/ec;
00264   ecut = 0.;
00265   cut = a*a - 4.*p*b;
00266   if (cut > 0.) { ecut = std::sqrt(cut); }
00267   ecut = (ecut-a) / (p+p);
00268   ecut2 = ecut;
00269   //JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
00270   // ecut<0 means that there is no cut with energy axis, i.e. xs is set 
00271   // to 0 bellow minimum
00272   //  if (cut < 0.) ecut2 = ecut - 2.;
00273   if (cut < 0.) { ecut2 = ecut; }
00274   elab = K * FragmentA /G4double(ResidualA);
00275   sig = 0.;
00276   if (elab <= ec) { //start for E<Ec 
00277     if (elab > ecut2) { sig = (p*elab*elab+a*elab+b) * signor; }
00278     
00279     signor2 = (ec-elab-c) / w;
00280     signor2 = 1. + std::exp(signor2);
00281     sig = sig / signor2;
00282   }              //end for E<=Ec
00283   else{           //start for  E>Ec
00284     sig = (landa*elab+mu+nu/elab) * signor;
00285     geom = 0.;
00286     
00287     if (xnulam < flow || elab < etest) 
00288       {
00289         if (sig <0.0) {sig=0.0;}
00290         return sig;
00291       }
00292     geom = std::sqrt(theA*K);
00293     geom = 1.23*ResidualAthrd + ra + 4.573/geom;
00294     geom = 31.416 * geom * geom;
00295     sig = std::max(geom,sig);
00296     
00297   }   //end for E>Ec
00298   return sig;
00299 }

G4double G4PreCompoundProton::GetRj	(	G4int	NumberParticles,
		G4int	NumberCharged
	)			`[protected, virtual]`

Implements G4PreCompoundNucleon.

Definition at line 65 of file G4PreCompoundProton.cc.

00066 {
00067   G4double rj = 0.0;
00068   if(nParticles > 0) { 
00069     rj = static_cast<G4double>(nCharged)/static_cast<G4double>(nParticles);
00070   }
00071   return rj;
00072 }


Public Member Functions
	G4PreCompoundProton ()
virtual	~G4PreCompoundProton ()
Protected Member Functions
virtual G4double	GetRj (G4int NumberParticles, G4int NumberCharged)
virtual G4double	CrossSection (G4double ekin)
virtual G4double	GetAlpha ()
virtual G4double	GetBeta ()
G4double	GetOpt1 (G4double K)
G4double	GetOpt2 (G4double K)
G4double	GetOpt3 (G4double K)

G4PreCompoundProton Class Reference

Public Member Functions

Protected Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation