#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 Z_A = 3.82; // atomic number double rho_m = 1.19; // mass density kg/m3 double M_A = 7.15; // Molar Mass kg/Mol double N_V = 0.96; // Effective number of valence electrons double E_gap = 3.7; // Band gap energy double N_A = 6.02214076E23; // Avagadro's number double q_e = 1.60217662E-19; // Charge on electron double m_e = 9.10938356E-31; // mass on electron double eps_0 = 8.85418782E-12; // permititivity of free space double h_bar = 1.054571817E-34; // reduced Planck constant double E_LO = 1E3 ; // upper limit for ESTAR data double E_HI = 20E3 ; // upper limit for ESTAR data double k = 0.8 ; // fixed constant double mec2 = 0.511E6; // rest mass of electron double d1 = 3.600; // tabata constant d1 double d2 = 0.9882; // tabata constant d2 double d3 = 1.191E-3; // tabata constant d3 double d4 = 0.8622; // tabata constant d4 double d5 = 1.02501; // tabata constant d5 double d6 = 1.0803E-4; // tabata constant d6 double d7 = 0.99628; // tabata constant d7 double d8 = 1.303E-4; // tabata constant d8 double d9 = 1.02441 ; // tabata constant d9 double d10 = 1.2986E-4; // tabata constant d10 double d11 = 1.030; // tabata constant d11 double d12 = 1.11E-2; // tabata constant d12 double d13 = 1.1E-6; // tabata constant d13 double d14 = 0.959; // tabata constant d14 double b1 = 0.3879; // tabata constant d1 double b2 = 0.2178 ; // tabata constant d2 double b3 = 0.4541 ; // tabata constant d3 double b4 = 0.03068 ; // tabata constant d4 double b5 = 3.326E-16; // tabata constant d5 double b6 = 13.24 ; // tabata constant d6 double b7 = 1.316 ; // tabata constant d7 double b8 = 14.03 ; // tabata constant d8 double b9 = 0.7406 ; // tabata constant d9 double b10 = 4.294E-3; // tabata constant d10 double b11 = 1.684 ; // tabata constant d11 double b12 = 0.2264 ; // tabata constant d12 double b13 = 0.6127 ; // tabata constant d13 double b14 = 0.1207 ; // tabata constant d14 double a1 = b1*pow(Z_A,b2) ; // tabata constant a1 double a2 = b3+b4*Z_A ; // tabata constant a2 double a3 = b5*pow(Z_A, (b6-b7*log(Z_A))) ; // tabata constant a3 double a4 = b8/pow(Z_A, b9) ; // tabata constant c1 double a5 = b10*pow(Z_A, (b11-b12*log(Z_A))) ; // tabata constant a3 double a6 = b13*pow(Z_A,b14) ; // tabata constant c1 double c1 = d1*M_A/pow(Z_A,d2) ; // tabata constant c1 double c2 = d3*pow(Z_A,d4) ; // tabata constant c2 double c3 = d5-d6*Z_A ; // tabata constant c1 double c4 = d7-d8*Z_A ; // tabata constant c1 double c5 = d9-d10*Z_A ; // tabata constant c1 double c6 = d11/pow(Z_A,d12) ; // tabata constant c1 double c7 = d13/pow(Z_A,d14) ; // tabata constant c1 double E_e = 3.75; // cut off energy for linac energy profile of double C1 = 3.824 ; double C2 = 3.522 ; double C3 = -1.222 ; double C4 = 0.308/100 ; double C5 = 1.10786E-06; double I = 74E-6/rest_mass; // PMMA mean Ionization energy as fraction of rest mass double wa_0 = 6.575; // constant for parametrization of mu_W double wa_1 = -3.623E1; // constant for parametrization of mu_W double wa_2 = -1.578; // constant for parametrization of mu_W double wb_0 = 1.0; // constant for parametrization of mu_W double wb_1 = 9.667; // constant for parametrization of mu_W double wb_2 = 7.132E-1; // constant for parametrization of mu_W double wb_3 = -3.778E-4; // constant for parametrization of mu_W double ala_0 = 6.107E-2; // constant for parametrization of mu_W double ala_1 = 1.694E-2; // constant for parametrization of mu_W double ala_2 = -2.390E-3; // constant for parametrization of mu_W double ala_3 = 3.116E-4; // constant for parametrization of mu_W double ala_4 = 5.286E-4; // constant for parametrization of mu_W double ala_5 = 2.507E-5; // constant for parametrization of mu_W double alb_0 = 1.0; // constant for parametrization of mu_W double alb_1 = 7.769E-1; // constant for parametrization of mu_W double alb_2 = 2.434E-1; // constant for parametrization of mu_W double alb_3 = 3.838E-2; // constant for parametrization of mu_W double alb_4 = 3.042E-3; // constant for parametrization of mu_W double alb_5 = 9.686E-5; // constant for parametrization of mu_W double lambda(double Energy, double Energy_m, double Energy_p_eff){ //numerator function double E = Energy; double E_m = Energy_m; double E_p_eff = Energy_p_eff; double beta = -0.1+ 2.783/E_m + 0.069*pow(rho_m,0.1); double gamma = 0.191/sqrt(rho_m); double U = E_p_eff*E_p_eff*m_e*eps_0/(N_A*h_bar*h_bar*1E6); //E_p_eff*E_p_eff/(28.8*28.8) ; double C = 1.97-0.91*U; double D = 53.4-20.8*U; // return E/(E_p_eff*E_p_eff*(beta*log(gamma*E)-C/E+D/(E*E))); // Angstroms return E/(E_p_eff*E_p_eff*(beta*log(gamma*E)-C/E+D/(E*E)))*1E-10; // meters } double R_GWJD1(double Energy){ //CSDA Range according to GWJD double E_0 = Energy*1E6; // energy converted from MeV to eV double E_p_eff = h_bar*sqrt((N_V*N_A*rho_m*1E6)/(m_e*eps_0*M_A)); // effective plasmon energy electronvolts double E_m = 2.8*sqrt(E_p_eff*E_p_eff+E_gap*E_gap) ; // mean energy lost per collision double n = 1+log(log((k+E_HI/(E_m*N_V))*sqrt(0.5*exp(1)))/log((k+E_LO/(E_m*N_V))*sqrt(0.5*exp(1))))/log(E_LO/E_HI); double b = pow(E_LO,(1-n))*lambda(E_LO,E_m,E_p_eff)*pow((1-exp(-E_LO/E_m)),-1)/(E_m*(1-pow((1+E_LO/(N_V*mec2)),-2))); if(E_0 < E_m){ return (E_0/E_m)*lambda(E_m,E_m,E_p_eff)*(1-exp(-1))/pow((1-exp(-E_0/E_m)),2)*100; //return range in cm } else if(E_m <= E_0 && E_0<= E_HI){ return (E_0/E_m)*lambda(E_0,E_m,E_p_eff)/(1-exp(-E_0/E_m))*100; //return range in cm } else return b*pow(E_0, n)*(1-pow((1+E_0/(N_V*mec2)),-2))*100; //return range in cm } double R_GWJD2(double Energy){ //Extrapolated range, according to GWJD double E_0 = Energy; double t_0 = E_0/rest_mass; double f_d = 1/(a1+a2/(1+a3/pow(t_0,a4)+a5*pow(t_0,a6))); // stopping number return f_d*R_GWJD1(E_0); // return in CM } 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] ; //return 3*B/sqrt((E0*(R-ess)/R)*(E0*(R-ess)/R+2*rest_mass)); return -3*B/sqrt((E0/(1+exp((ess-d50)/a)))*(E0/(1+exp((ess-d50)/a))+2*rest_mass)); } double dz_by_ds(double* s, double* p){ //S as function of rho double E0 = p[0]; double R = p[1]; double B = p[2]; double d50 = p[3] ; double a = p[4] ; double ess = s[0] ; TF1 deta_ds("deta_by_ds", deta_by_ds, 0.0, R, 5); // denominator function deta_ds.SetParameters(E0, R, B, d50, a); // Set A0, mu, R, E0 return cos(deta_ds.Integral(0, ess)); } double zee(double* s, double *p){ double E0 = p[0]; double R = p[1]; double B = p[2]; double d50 = p[3] ; double a = p[4] ; double ess = s[0] ; TF1 dzed_ds("dz_by_ds", dz_by_ds,0.0, R, 5); // denominator function dzed_ds.SetParameters(E0, R, B, d50, a); // Set A0, mu, R, E0 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] ; TF1 z_as_func_B("zee", zee, 0, R, 5); // numerator function z_as_func_B.SetParameters(E0, R, B, d50, a); // Set A0, mu, R, E 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] ; TF1 dzed_ds("dz_by_ds", dz_by_ds,0.0, R, 5); // denominator function dzed_ds.SetParameters(E0, R, B, d50, a); // Set A0, mu, R, E TF1 deta_ds("deta_by_ds", deta_by_ds, 0.0, R, 5); // denominator function deta_ds.SetParameters(E0, R, B, d50, a); // Set A0, mu, R, E0 return (-1+A0*exp(mu*dzed_ds.Integral(0.000, ess)))*cos(deta_ds.Integral(0, ess))*E0/(1+exp((ess-d50)/a)); } double dalpha_by_dz_shab(double* z, double* p){ 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 zed = z[0]; return (-1+A0*exp(mu*zed))*E0/(1+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]; TF1 dalpha_ds_shab("dalpha_by_ds_shab", dalpha_by_ds_shab ,0.0, 7, 7); // denominator function dalpha_ds_shab.SetParameters(E0, R, B, A0, mu, d50, a); // Set A0, mu, R, E TF1 z_for_B("zee", zee, 0, R, 5); // numerator function z_for_B.SetParameters(E0, R, B, d50, a); // Set A0, mu, R, E TF1 dalpha_dz_shab("dalpha_dz_shab", dalpha_by_dz_shab ,0.0, 7,7); // denominator function dalpha_dz_shab.SetParameters(E0, R, B, A0, mu, d50, a); // Set A0, mu, R, E return dalpha_ds_shab.Integral(0.000, R)/dalpha_dz_shab.Integral(0.000, z_for_B.GetMaximum()); //return dalpha_ds.Integral(0.000, R)/dalpha_dz.Integral(0.000, R); } int main() { gROOT->SetStyle("Plain"); //set plain TStyle gStyle->SetOptStat(0); //draw statistics on plots, 0 for no output //gStyle->SetOptStat(111111); //draw statistics on plots, 0 for no output gStyle->SetOptFit(1111); //draw fit results on plot, 0 for no output gStyle->SetPalette(57); //set color map gStyle->SetOptTitle(0); //suppress title box // ------------------------------------------ 4 MVV --------------------------------------------------------------- const int points_4MVV=8; double field_4MVV[points_4MVV]={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7}; double aleph_4MVV[points_4MVV]={1, 0.999976, 0.999951, 0.999911, 0.999849, 0.999782, 0.99969, 0.999565}; /* const int points_4MVV=22; double field_4MVV[points_4MVV]={0,0.1,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}; double aleph_4MVV[points_4MVV]={1, 0.999976, 0.999951, 0.999911, 0.999849, 0.999782, 0.99969, 0.999565, 0.999447, 0.999309, 0.999156, 0.999004, 0.998829, 0.998659, 0.998471, 0.998281, 0.997506, 0.995356, 0.993983, 0.993127, 0.992565, 0.992219}; double alpha_err[points]={0, 0.0566732, 0.0566725, 0.0566714, 0.0566696, 0.0566677, 0.0566651, 0.0566616, 0.0566582, 0.0566543, 0.05665, 0.0566457, 0.0566407, 0.0566359, 0.0566305, 0.0566252, 0.0566032, 0.0565422, 0.0565031, 0.0564788, 0.0564628, 0.056453};*/ TCanvas *c = new TCanvas("c","aleph", 700, 500); TGraph * graph_4MVV = new TGraph(points_4MVV, field_4MVV, aleph_4MVV); graph_4MVV->GetXaxis()->SetRangeUser(0, 0.72); graph_4MVV->GetYaxis()->SetRangeUser(0.9988, 1.0001); graph_4MVV->GetYaxis()->SetDecimals(3); graph_4MVV->GetXaxis()->SetTitle("B (T)"); graph_4MVV->GetYaxis()->SetTitle("#alpha = D(B)/D(B = 0 T)"); graph_4MVV->GetXaxis()->CenterTitle(); graph_4MVV->GetYaxis()->CenterTitle(); graph_4MVV->SetMarkerSize(1); graph_4MVV->SetMarkerColor(kBlack); graph_4MVV->SetMarkerStyle(20); graph_4MVV->SetLineColor(kBlack); graph_4MVV->SetLineWidth(1); graph_4MVV->Draw(); graph_4MVV->SetName("4 MVV Fit"); graph_4MVV->GetYaxis()->SetNdivisions(5,5,0); graph_4MVV->GetYaxis()->SetTitleOffset(1.8); //graph_4MVV->SetTitle("Alpha 4 MV Water"); double A0_4MVV = 3.42845/2.36941 ; //depth dose intercept double mu_4MVV = 0.000298483*40 ; //attenuation parameter double E0_4MVV = 3.75; // particle energy double R_4MVV = R_GWJD1(E0_4MVV) ; // csda range double d50_4MVV = 1.4*R_4MVV ; // original initial parameters double a_4MVV = 0.4*R_4MVV/4.4 ; ; cout<< 1.4*R_4MVV << endl; // d50 original initial parameters cout << 0.4*R_4MVV/4.4 << endl; // a original initial parameters cout<< R_4MVV << endl; // d50 original initial parameters TF1* alpha_4MVV= new TF1("alpha_4MVV", alpha_shab, 0., 0.7, 6); // numerator function alpha_4MVV -> SetParameters(E0_4MVV, R_4MVV, A0_4MVV, mu_4MVV, d50_4MVV, a_4MVV); // E/rest_mass = 10/rest_mass, St, R alpha_4MVV -> FixParameter(0, E0_4MVV); alpha_4MVV -> FixParameter(1, R_4MVV); alpha_4MVV -> FixParameter(2, A0_4MVV); alpha_4MVV -> FixParameter(3, mu_4MVV); alpha_4MVV -> SetParameter(4, 2.03629 ); alpha_4MVV -> SetParameter(5, 0.129622 ); alpha_4MVV -> SetParName(4, "x_{d}"); alpha_4MVV -> SetParName(5, "a"); // TCanvas *c2 = new TCanvas("c2","alpha", 700, 500); graph_4MVV->Fit("alpha_4MVV", "R"); alpha_4MVV-> Draw("SAME"); return 0; }