#include #include #include #include #include using namespace std; #include #include #include #include #include #include #include #include double rest_mass = 0.511 ; // rest mass of electron double deta_by_ds(double* s, double* p){ //S as function of rho double ess = s[0]; double E0 = p[0]; double R = p[1]; double B = p[2]; double d50 = p[3] ; double a = p[4] ; double A = p[5] ; return -3*B/sqrt((E0*1.08/(1+A*exp((ess-d50)/a)))*(E0*1.08/(1+A*exp((ess-d50)/a))+2*rest_mass)); } double dz_by_ds(double* s, double* p){ //S as function of rho double ess = s[0] ; double E0 = p[0]; double R = p[1]; double B = p[2]; double d50 = p[3] ; double a = p[4] ; double A = p[5] ; TF1 deta_ds("deta_by_ds", deta_by_ds, 0.0, R, 6); // denominator function deta_ds.SetParameters(E0, R, B, d50, a, A); // Set A0, mu, R, E0 deta_ds.FixParameter(0, E0); deta_ds.FixParameter(1, R); deta_ds.FixParameter(2, B); return cos(deta_ds.Integral(0, ess)); } double zee(double* s, double *p){ double ess = s[0] ; double E0 = p[0]; double R = p[1]; double B = p[2]; double d50 = p[3] ; double a = p[4] ; double A = p[5] ; TF1 dzed_ds("dz_by_ds", dz_by_ds,0.0, R, 6); // denominator function dzed_ds.SetParameters(E0, R, B, d50, a, A); // Set A0, mu, R, E0 dzed_ds.FixParameter(0, E0); dzed_ds.FixParameter(1, R); dzed_ds.FixParameter(2, B); return dzed_ds.Integral(0.000, ess); // integrate up rg = E0/qB } double zee_max(double* b, double *p){ double B = b[0] ; double E0 = p[0]; double R= p[1]; double d50 = p[2] ; double a = p[3] ; double A = p[4] ; TF1 z_as_func_B("zee", zee, 0, R, 6); // numerator function z_as_func_B.SetParameters(E0, R, B, d50, a, A); // Set A0, mu, R, E z_as_func_B.FixParameter(0, E0); z_as_func_B.FixParameter(1, R); z_as_func_B.FixParameter(2, B); return z_as_func_B.GetMaximum() ; } double dalpha_by_ds_shab(double* s, double* p){ double ess = s[0]; double E0 = p[0]; double R = p[1]; double B = p[2]; double A0 = p[3] ; //depth dose intercept double mu = p[4] ; //attenuation parameter double d50 = p[5] ; double a = p[6] ; double A = p[7] ; TF1 dzed_ds("dz_by_ds", dz_by_ds,0.0, R, 6); // denominator function dzed_ds.SetParameters(E0, R, B, d50, a, A); // Set A0, mu, R, E dzed_ds.FixParameter(0, E0); dzed_ds.FixParameter(1, R); dzed_ds.FixParameter(2, B); TF1 deta_ds("deta_by_ds", deta_by_ds, 0.0, R, 6); // denominator function deta_ds.SetParameters(E0, R, B, d50, a, A); // Set A0, mu, R, E0 deta_ds.FixParameter(0, E0); deta_ds.FixParameter(1, R); deta_ds.FixParameter(2, B); return (-1+A0*exp(mu*dzed_ds.Integral(0.000, ess)))*cos(deta_ds.Integral(0, ess))*E0*1.08/(1+A*exp((ess-d50)/a)); } double dalpha_by_dz_shab(double* z, double* p){ double zed = z[0]; double E0 = p[0]; double R = p[1]; double B = p[2]; double A0 = p[3] ; //depth dose intercept double mu = p[4] ; //attenuation parameter double d50 = p[5] ; double a = p[6] ; double A = p[7] ; return (-1+A0*exp(mu*zed))*E0*1.08/(1+A*exp((zed-d50)/a)); } double alpha_shab(double* b, double *p) { double B = b[0] ; double E0 = p[0]; double R = p[1]; double A0 = p[2] ; //depth dose intercept double mu = p[3] ; //attenuation parameter double d50 = p[4] ; double a = p[5] ; double A = p[6] ; TF1 dalpha_ds_shab("dalpha_by_ds_shab", dalpha_by_ds_shab ,0.0, 7, 8); // denominator function dalpha_ds_shab.SetParameters(d50, a, A); // Set A0, mu, R, E dalpha_ds_shab.FixParameter(0, E0); dalpha_ds_shab.FixParameter(1, R); dalpha_ds_shab.FixParameter(2, B); dalpha_ds_shab.FixParameter(3, A0); dalpha_ds_shab.FixParameter(4, mu); TF1 z_for_B("zee", zee, 0, R, 6); // numerator function z_for_B.SetParameters(E0, R, B, d50, a, A); // Set A0, mu, R, E TF1 dalpha_dz_shab("dalpha_dz_shab", dalpha_by_dz_shab ,0.0, 7, 8); // denominator function dalpha_dz_shab.SetParameters(d50, a, A); // Set A0, mu, R, E dalpha_dz_shab.FixParameter(0, E0); dalpha_dz_shab.FixParameter(1, R); dalpha_dz_shab.FixParameter(2, B); dalpha_dz_shab.FixParameter(3, A0); dalpha_dz_shab.FixParameter(4, mu); return dalpha_ds_shab.Integral(0.000, R)/dalpha_dz_shab.Integral(0.000, z_for_B.GetMaximum()); } int main() { const int points=23; double field[points]={0,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,2,3,4,5,6,7,8,9}; double aleph[points]={1,0.9999437,0.999734507,0.999385603,0.999000074,0.998535291,0.998014297,0.997441248,0.996780681,0.996168181,0.99546186,0.99478182,0.994043266,0.993393902,0.992664322,0.98966929,0.985912217,0.984121614,0.983207472,0.982636313,0.982315577,0.982197799,0.982103447}; TCanvas *c = new TCanvas("c","aleph", 700, 500); TGraph * graph=new TGraph(points, field, aleph); graph->GetXaxis()->SetRangeUser(0, 1); graph->GetYaxis()->SetRangeUser(0.996, 1.001); graph->GetXaxis()->SetTitle("Field / T"); graph->GetYaxis()->SetTitle("alpha "); graph->GetXaxis()->CenterTitle(); graph->GetYaxis()->CenterTitle(); graph->SetMarkerSize(1); graph->SetMarkerColor(kBlue); graph->SetMarkerStyle(20); graph->SetLineColor(2); graph->SetLineWidth(0); graph->Draw(); graph->SetTitle("Alpha 6 MV Water"); double E0 = 6; // particle energy double R = 3.052 ; // csda range double A0 = 2.52831; //depth dose intercept double mu = 0.0092 ; //attenuation parameter double d50 = 0.8*R ; double a = 100*0.4*R/4.4 ; double A = 0.2; TF1* alpha= new TF1("alpha", alpha_shab, 0.0, 0.7, 7); // numerator function alpha -> FixParameter(0, E0); alpha -> FixParameter(1, R); alpha -> FixParameter(2, A0); alpha -> FixParameter(3, mu); alpha -> SetParameter(4, d50); // E/rest_mass = 10/rest_mass, St, R //alpha -> SetParLimits(4, 0, R); // E/rest_mass = 10/rest_mass, St, R alpha -> SetParameter(5, a); // E/rest_mass = 10/rest_mass, St, R // alpha -> SetParLimits(5, 0.4*R/4.4, 100*0.4*R/4.4); // E/rest_mass = 10/rest_mass, St, R alpha ->SetParameter(6, A); // E/rest_mass = 10/rest_mass, St, R // TCanvas *c2 = new TCanvas("c2","alpha", 700, 500); graph->Fit("alpha", "R"); alpha-> Draw("AL same"); return 0; }