X0 Material

calculations.cpp

@@ -486,6 +486,15 @@ double angular_scattering_variance(const Projectile &p, const Target &t){
     return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
 }
 
+double angular_scattering_variance(const Projectile &p, const Material &mat){
+    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 X0 = radiation_length(mat);
+    return 198.81 * pow(p.Z,2)/(X0*pow(_p*beta,2));
+}
+
 /// radiation lengths are taken from Particle Data Group 2014
 double radiation_length(int z, double m){
     double lr = 0;
@@ -524,6 +533,16 @@ double radiation_length(int z, double m){
     return lr;
 }
 
+double radiation_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/radiation_length(t.Z,t.A);
+    }
+    return 1/sum;
+}
+
 
 double precalculated_lindhard(const Projectile &p){
     double T = p.T;

calculations.h

@@ -102,6 +102,7 @@ namespace catima{
 
 
     double angular_scattering_variance(const Projectile &p, const Target &t);
+    double angular_scattering_variance(const Projectile &p, const Material &material);
 
     /**
       * returns radiation length of the (M,Z) material
@@ -111,8 +112,20 @@ namespace catima{
       * @return radiation length in g/cm^2
       */
     double radiation_length(int z, double m);
+
+    /**
+      * 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 radiation_length(const Material &material);
+
     
     /** returns effective Z of the projectile
+      * @param p - Projectile class
+      * @param t - Target class
       * @param c - Configuration, the z effective will be calculated according to c.z_effective value
       * @return - z effective
       */ 

catima.cpp

@@ -58,14 +58,7 @@ double domega2dx(const Projectile &p, const Material &mat, const Config &c){
 }
 
 double da2dx(const Projectile &p, const Material &mat, const Config &c){
-    double sum = 0;
+    return angular_scattering_variance(p,mat);    
-    
-    for(int i=0;i<mat.ncomponents();i++){
-        auto t = mat.get_element(i);
-        double w = mat.weight_fraction(i);
-        sum += w*angular_scattering_variance(p,t);
-    }
-    return sum;
 }
 
 
@@ -272,7 +265,6 @@ Result calculate(Projectile p, const Material &t, const Config &c){
     
     double rrange = std::min(res.range/t.density(), t.thickness_cm());
     auto fx2p = [&](double x)->double{ 
-	    double range = range_spline(T);
         double e =energy_out(T,x*t.density(),range_spline);
 	    return (rrange-x)*(rrange-x)*da2dx(p(e), t, c)*t.density();
             };           

structures.h

@@ -73,6 +73,7 @@ namespace catima{
             double th=0;
             double molar_mass=0;
             double i_potential=0;
+            double _X0=0;
             std::vector<Target>atoms;
 
         public:
@@ -187,6 +188,17 @@ namespace catima{
               */
             double I() const {return i_potential;};
 
+            /**
+             * set radiation length in g/cm2 unit
+             * @param value - radiation length in g/cm2 unit, if 0 default radiation length will be calculated  
+             */ 
+            Material& X0(double value){_X0 = value;return *this;}
+
+            /**
+             * return radiation length in g/cm2 unit             
+             */ 
+            double X0(){return _X0;}
+
             /**
               * return number density of atoms/molecules per cm3 in 10^23 units
               */
@@ -215,7 +227,7 @@ namespace catima{
             double number_density_cm2(int i)const{
               if(i>=atoms.size())return 0.0;
               return number_density_cm2()*molar_fraction(i);
-            }
+            }            
 
             friend bool operator==(const Material &a, const Material&b);
     };

tests/test_calculations.cpp

@@ -404,6 +404,7 @@ using namespace std;
         auto water = catima::get_material(catima::material::Water);
         auto res2 = catima::calculate(p(600),water,600);
         CHECK(res2.dEdxi == approx(92.5).epsilon(2));
+        CHECK(catima::radiation_length(water)==approx(36.1,0.2));
     }
     TEST_CASE("z_eff"){
         using namespace catima;