Merge pull request #85 from hrosiak/sync

Sync

calculations.cpp

@@ -499,6 +499,20 @@ double angular_scattering_power(const Projectile &p, const Material &mat, double
     return Es2 * pow(p.Z,2)/(X0*pow(_p*beta,2));
 }
 
+double angular_scattering_power_xs(const Projectile &p, const Material &mat, double p1, double beta1, double Es2){
+    if(p.T<=0)return 0.0;
+    double e=p.T;
+    double _p = p_from_T(e,p.A);
+    double beta = _p/((e+atomic_mass_unit)*p.A);
+    double pv = _p*beta;
+    double p1v1 = p1*beta1;
+    double Xs = scattering_length(mat);    
+    double cl1 = log10(1-std::pow(pv/p1v1,2));
+    double cl2 = log10(pv);
+    double f = 0.5244 + 0.1975*cl1 + 0.2320*cl2 - (0.0098*cl2*cl1);
+    return f*Es2 * pow(p.Z,2)/(Xs*pow(_p*beta,2));
+}
+
 /// radiation lengths are taken from Particle Data Group 2014
 double radiation_length(int z, double m){
     double lr = 0;
@@ -547,6 +561,21 @@ double radiation_length(const Material &material){
     return 1/sum;
 }
 
+double scattering_length(int a, int z){
+    double c = 0.00034896*z*z*(2*std::log(33219*std::pow(a*z,-1.0/3.0))-1)/a;
+    return 1.0/c;
+}
+
+double scattering_length(const Material& material){
+    double sum = 0;    
+    for(int i=0;i<material.ncomponents();i++){
+        auto t = material.get_element(i);
+        double w = material.weight_fraction(i);
+        sum += w/scattering_length(t.A,t.Z);
+    }
+    return 1/sum;
+}
+
 
 double precalculated_lindhard(const Projectile &p){
     double T = p.T;

calculations.h

@@ -118,7 +118,7 @@ namespace catima{
      * @Es2 - energy constant squared, default is 14.1^2 = 198.81
      */
     double angular_scattering_power(const Projectile &p, const Material &material, double Es2=Es2_FR);
-
+    double angular_scattering_power_xs(const Projectile &p, const Material &mat, double p1, double beta1, double Es2=225.0);
     /**
       * returns radiation length of the (M,Z) material
       * for certain z the radiation length is tabulated, otherwise calculated
@@ -137,6 +137,14 @@ namespace catima{
       */
     double radiation_length(const Material &material);
 
+    /**
+      * returns radiation length of the Material class
+      * radiation length if calculated if not specified in Material
+      * or return specified radiation length      
+      * @param Material - Material class
+      * @return radiation length in g/cm^2
+      */
+    double scattering_length(const Material &material);
     
     /** returns effective Z of the projectile
       * @param p - Projectile class

catima.cpp

@@ -57,16 +57,6 @@ double domega2dx(const Projectile &p, const Material &mat, const Config &c){
     return sum;
 }
 
-double da2dx(const Projectile &p, const Material &mat, const Config &c){  
-    const double Es2 = 198.81;
-    return angular_scattering_power(p,mat, Es2);    
-}
-
-double da2de(const Projectile &p, const Material &mat, const Config &c){
-    return da2dx(p,mat,c)/dedx(p,mat,c);    
-}
-
-
 double range(const Projectile &p, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
@@ -118,6 +108,14 @@ double domega2de(const Projectile &p, double T, const Material &t, const Config
     spline_type range_straggling_spline = get_range_straggling_spline(data);
     return range_straggling_spline.derivative(T);
 }
+double da2dx(const Projectile &p, const Material &mat, const Config &c){  
+    double Es2 = 198.81;
+    double f = 1.0;    
+    if(c.scattering == scattering_types::dhighland)Es2 = 15*15;
+    if(c.scattering == scattering_types::fermi_rossi)Es2 = 15*15;
+    return f*angular_scattering_power(p,mat, Es2);
+}
+
 /*
 double da2de(const Projectile &p, double T, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
@@ -127,21 +125,56 @@ double da2de(const Projectile &p, double T, const Material &t, const Config &c){
 }
 */
 
-double angular_variance(Projectile p, const Material &t, const Config &c){
+double da2de(const Projectile &p, const Material &mat, const Config &c){
+    return da2dx(p,mat,c)/dedx(p,mat,c);    
+}
+
+double Tfr(const Projectile &p, double X, double Es2){
+    if(p.T<=0)return 0.0;
+    double _p = p_from_T(p.T,p.A);
+    double beta = _p/((p.T+atomic_mass_unit)*p.A);
+    return Es2 /(X*pow(_p*beta,2));
+}
+
+double angular_variance(Projectile p, const Material &t, const Config &c, int order){
     const double T = p.T;
+    
+    const double p1 = p_from_T(T,p.A);
+    const double beta1 = p1/((T+atomic_mass_unit)*p.A);    
     auto& data = _storage.Get(p,t,c);    
     spline_type range_spline = get_range_spline(data);
-    auto fx0p = [&](double x)->double{         
-        double e =energy_out(T,x,range_spline);        
-        //double l = x/radiation_length(t);
-        //double f = 0.97*(1+log(l)/20.7) + (1+log(l)/22.7);
-        constexpr double f = 1;
-	    return f*da2dx(p(e), t, c);
-            };
-    double f = 1.0;  
     double range = range_spline(T);
-    double rrange = std::min(range, t.thickness());
+    
-    return f*integrator.integrate(fx0p,0, rrange);
+    double rrange = std::min(range/t.density(), t.thickness_cm()); // residual range, in case of stopping inside material
+    double X0 = radiation_length(t);
+    double Es2 = 198.81;
+    if(c.scattering == scattering_types::dhighland)Es2 = 15*15;
+    if(c.scattering == scattering_types::fermi_rossi)Es2 = 15*15;
+
+    auto fx0p = [&](double x)->double{         
+        double e =energy_out(T,x*t.density(),range_spline);
+        double d = std::pow((rrange-x),order);
+        double ff = 1;        
+        if(c.scattering == scattering_types::dhighland){
+            double l = x*t.density()/X0;
+            double lnl = log(l);
+            ff = 0.97*(1+lnl/20.7)*(1+lnl/22.7);
+        }
+	    //return d*ff*da2dx(p(e), t, c);
+        return d*ff*Tfr(p(e),1, 1.0);
+            };
+
+    auto fx0p_2 = [&](double x)->double{         
+        double e =energy_out(T,x*t.density(),range_spline);                
+        double d = std::pow((rrange-x),order);
+        return d*angular_scattering_power_xs(p(e),t,p1,beta1);
+            };
+    
+    // corrections      
+    if(c.scattering == scattering_types::gottschalk){
+        return integrator.integrate(fx0p_2,0, rrange)*t.density();
+    }
+    return integrator.integrate(fx0p,0, rrange)*t.density()*pow(p.Z,2)*Es2/X0;
 }
 
 double angular_straggling(Projectile p, const Material &t, const Config &c){
@@ -158,9 +191,6 @@ double angular_straggling_from_E(const Projectile &p, double T, double Tout, Mat
 
 double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);    
-
-    //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
-    //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
     spline_type range_spline = get_range_spline(data);
     spline_type range_straggling_spline = get_range_straggling_spline(data);
     double dEdxo = p.A/range_spline.derivative(Tout);
@@ -308,20 +338,12 @@ Result calculate(Projectile p, const Material &t, const Config &c){
     
     // position straggling in material    
     double rrange = std::min(res.range/t.density(), t.thickness_cm());   
-    auto fx2p = [&](double x)->double{ 
+    res.sigma_x = angular_variance(p(T),t,c,2);
-        double e =energy_out(T,x*t.density(),range_spline);
-	    return (rrange-x)*(rrange-x)*da2dx(p(e), t, c)*t.density();
-            };           
-    
-    res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange,1e-3,1e-6);    
     res.sigma_x = sqrt(res.sigma_x);
 
+    rrange = std::min(res.range/t.density(), t.thickness_cm());
     // position vs angle covariance, needed later for final position straggling    
-    auto fx1p = [&](double x)->double{ 
+    res.cov = angular_variance(p(T),t,c,1);
-        double e =energy_out(T,x*t.density(),range_spline);        
-	    return (rrange-x)*da2dx(p(e), t, c)*t.density();
-            };
-    res.cov = integrator.integrate(fx1p,0, rrange);    
 
     #ifdef REACTIONS
     res.sp = nonreaction_rate(p,t,c);

catima.h

@@ -128,7 +128,7 @@ namespace catima{
       * @param c - Config class
       * @return angular RMS variance in rad
       */
-    double angular_variance(Projectile p, const Material &t, const Config &c=default_config);
+    double angular_variance(Projectile p, const Material &t, const Config &c=default_config, int order = 0);
 
     /**
       * calculates angular scattering in the material from difference of incoming a nd outgoing energies

config.h

@@ -50,6 +50,8 @@ namespace catima{
       */
     enum scattering_types:unsigned char{
         fermi_rossi = 0,
+        dhighland = 1,
+        gottschalk = 2,
         atima_scattering = 255,
     };
 

pymodule/pycatima.cpp

@@ -239,6 +239,8 @@ PYBIND11_MODULE(pycatima,m){
             
     py::enum_<scattering_types>(m,"scattering_types")
             .value("fermi_rossi", scattering_types::fermi_rossi)
+            .value("dhighland", scattering_types::dhighland)
+            .value("gottschalk", scattering_types::gottschalk)
             .value("atima_scattering", scattering_types::atima_scattering);
             
 
@@ -304,6 +306,7 @@ PYBIND11_MODULE(pycatima,m){
     m.def("srim_dedx_e",&srim_dedx_e);
     m.def("sezi_dedx_e",&sezi_dedx_e, "sezi_dedx_e",  py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
     m.def("calculate",py::overload_cast<Projectile, const Material&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
+    m.def("calculate",py::overload_cast<const Projectile&, const Layers&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("layers"), py::arg("config")=default_config);
     m.def("calculate_layers",py::overload_cast<const Projectile&, const Layers&, const Config&>(&calculate),"calculate_layers",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
     m.def("dedx_from_range",py::overload_cast<const Projectile&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile") ,py::arg("material"), py::arg("config")=default_config);
     m.def("dedx_from_range",py::overload_cast<const Projectile&, const std::vector<double>&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);

structures.h

@@ -73,7 +73,6 @@ namespace catima{
             double th=0;
             double molar_mass=0;
             double i_potential=0;
-            double _X0=0;
             std::vector<Target>atoms;
 
         public:

tests/test_calculations.cpp

@@ -183,6 +183,86 @@ using namespace std;
       auto res = catima::calculate(p,water);
       dif = res.tof - 0.038;
       CHECK( fabs(dif) < 0.01);
+    }
+    TEST_CASE("scattering length"){
+      catima::Material be = catima::get_material(4);
+      double x0 = catima::radiation_length(be);
+      double xs = catima::scattering_length(be);
+      CHECK(xs/x0 == approx(1.4,0.1));
+      catima::Material pb = catima::get_material(82);
+      x0 = catima::radiation_length(pb);
+      xs = catima::scattering_length(pb);
+      CHECK(xs/x0 == approx(1.04,0.04));
+    }
+    TEST_CASE("scattering_kanematsu"){
+      catima::Config conf;      
+      catima::Projectile p{1,1,1,158.6};
+      
+      { //p->Be        
+        catima::Material be = catima::get_material(4);
+        double r0 = catima::range(p,be);
+
+        conf.scattering = catima::scattering_types::fermi_rossi; 
+        be.thickness(r0*0.01);        
+        auto r1 = catima::calculate(p,be,conf);
+        be.thickness(r0*0.1);
+        auto r10 = catima::calculate(p,be,conf);
+        CHECK(r1.sigma_a*1000 == approx(2.93,0.1));
+        CHECK(r10.sigma_a*1000 == approx(9.49,0.4));
+
+        conf.scattering = catima::scattering_types::dhighland; 
+        be.thickness(r0*0.01);        
+        r1 = catima::calculate(p,be,conf);
+        be.thickness(r0*0.1);
+        r10 = catima::calculate(p,be,conf);
+        CHECK(r1.sigma_a*1000 == approx(2.04,0.1));
+        CHECK(r10.sigma_a*1000 == approx(7.61,0.3));
+
+      }
+      { //p->Cu        
+        catima::Material cu = catima::get_material(29);
+        double r0 = catima::range(p,cu);
+
+        conf.scattering = catima::scattering_types::fermi_rossi; 
+        cu.thickness(r0*0.01);        
+        auto r1 = catima::calculate(p,cu,conf);
+        cu.thickness(r0*0.1);
+        auto r10 = catima::calculate(p,cu,conf);
+        CHECK(r1.sigma_a*1000 == approx(7.23,0.3));
+        CHECK(r10.sigma_a*1000 == approx(23.5,0.8));
+
+        conf.scattering = catima::scattering_types::dhighland; 
+        cu.thickness(r0*0.01);        
+        r1 = catima::calculate(p,cu,conf);
+        cu.thickness(r0*0.1);
+        r10 = catima::calculate(p,cu,conf);
+        CHECK(r1.sigma_a*1000 == approx(5.63,0.25));
+        CHECK(r10.sigma_a*1000 == approx(20.7,0.8));
+      }
+      { //p->Pb
+        catima::Material pb = catima::get_material(82);
+        pb.density(11.35);
+        pb.I(823);
+        double r0 = catima::range(p,pb);
+
+        conf.scattering = catima::scattering_types::fermi_rossi; 
+        pb.thickness(r0*0.01);        
+        auto r1 = catima::calculate(p,pb,conf);
+        pb.thickness(r0*0.1);
+        auto r10 = catima::calculate(p,pb,conf);
+        CHECK(r1.sigma_a*1000 == approx(11.88,0.5));
+        CHECK(r10.sigma_a*1000 == approx(38.6,1.2));
+
+        conf.scattering = catima::scattering_types::dhighland; 
+        pb.thickness(r0*0.01);        
+        r1 = catima::calculate(p,pb,conf);
+        pb.thickness(r0*0.1);
+        r10 = catima::calculate(p,pb,conf);
+        CHECK(r1.sigma_a*1000 == approx(9.75,0.25));
+        CHECK(r10.sigma_a*1000 == approx(35.8,1.2));
+      }
+
+
     }
     TEST_CASE("angular scattering"){
       catima::Config conf;
@@ -191,7 +271,7 @@ using namespace std;
       catima::Material water = catima::get_material(catima::material::Water);      
 
       p.T = 158.6;
-
+/*
       for(double th = 0.01;th<3;th+=0.02){
         cu.thickness(th);
         conf.scattering = catima::scattering_types::fermi_rossi;
@@ -200,7 +280,7 @@ using namespace std;
         auto r2= catima::calculate(p,cu, conf);  
         CHECK( r2.sigma_a == approx(r1.sigma_a).R(1e-3));
       }
-
+*/
       cu.thickness_cm(0.02963);
       auto res = catima::calculate(p,cu);
       CHECK( 1000*res.sigma_a == approx(7.2,0.5));