00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <math.h>
00013 #include <iostream>
00014 #include "TrkReco/RkFitMaterial.h"
00015
00016 using namespace std;
00017
00018 RkFitMaterial::RkFitMaterial(double z, double a, double i,
00019 double rho, double x0)
00020 : x0_(x0), z_(z)
00021
00022 {
00023 rza_ = rho * z / a;
00024 isq_ = i * i * 1e-18;
00025 }
00026
00027 RkFitMaterial::RkFitMaterial(const RkFitMaterial& mat)
00028 : rza_(mat.rza_), isq_(mat.isq_),
00029 x0_(mat.x0_), z_(mat.z_)
00030 {
00031
00032 }
00033
00034
00035 double RkFitMaterial::dE(double mass, double path,
00036 double p) const
00037 {
00038 if (!(p>0))
00039 return 0;
00040
00041
00042
00043
00044 const double Me = 0.000510999;
00045 double psq = p * p;
00046 double bsq = psq / (psq + mass * mass);
00047 double esq = psq / (mass * mass);
00048
00049 double s = Me / mass;
00050 double w = (4 * Me * esq
00051 / (1 + 2 * s * sqrt(1 + esq)
00052 + s * s));
00053
00054 double cc, x0;
00055 cc = 1+2*log(sqrt(isq_)/(28.8E-09*sqrt(rza_)));
00056 if (cc < 5.215)
00057 x0 = 0.2;
00058 else
00059 x0 = 0.326*cc-1.5;
00060 double x1(3), xa(cc/4.606), aa;
00061 aa = 4.606*(xa-x0)/((x1-x0)*(x1-x0)*(x1-x0));
00062 double delta(0);
00063 double x(log10(sqrt(esq)));
00064 if (x > x0){
00065 delta = 4.606*x-cc;
00066 if (x < x1) delta=delta+aa*(x1-x)*(x1-x)*(x1-x);
00067 }
00068
00069 float f1, f2, f3, f4, f5, ce;
00070 f1 = 1/esq;
00071 f2 = f1*f1;
00072 f3 = f1*f2;
00073 f4 = (f1*0.42237+f2*0.0304-f3*0.00038)*1E12;
00074 f5 = (f1*3.858-f2*0.1668+f3*0.00158)*1E18;
00075 ce = f4*isq_+f5*isq_*sqrt(isq_);
00076 return (0.0001535 * rza_ / bsq
00077 * (log(Me * esq * w / isq_)
00078 - 2 * bsq-delta-2.0*ce/z_)) * path;
00079 }
00080
00081
00082
00083 double RkFitMaterial::mcs_angle(double mass, double path,
00084 double p) const
00085 {
00086
00087
00088 double t = path / x0_;
00089 double psq = p*p;
00090 return 0.0136 * sqrt(t * (mass*mass + psq)) / psq
00091 * (1 + 0.038 * log(t));
00092 }
00093
00094
00095 double RkFitMaterial::del_E(double mass, double path,
00096 double p) const
00097 {
00098 double sigma0_2 = 0.1569*rza_*path;
00099
00100 if (sigma0_2<0) return 0;
00101
00102 double psq = p * p;
00103 double bsq = psq / (psq + mass * mass);
00104
00105
00106 double sigma_2 = sigma0_2*(1-0.5*bsq)/(1-bsq);
00107
00108 if (sigma_2<0) return 0;
00109
00110
00111 return sqrt(sigma_2)*0.001;
00112
00113 }