作业3.1 HPGe刻度¶
1.能量刻度¶
1.1 从.root文件中读取能谱¶
In [1]:
%jsroot on
TCanvas *c1 = new TCanvas("c1");
//从.root文件中读取能谱
TFile * filer = new TFile("gamma.root","READ");//以只读模式打开root文件，READ可替换为RECREARE等；
if(!filer->IsOpen())
{
std::cout<<"Can't open root file!"<<std::endl;
}
TH1F * h_gammaraw = (TH1F*)filer->Get("h0");//通过指针读取原始gamma谱
//此处通过指针h_gammaraw可对直方图进行操作
h_gammaraw->Draw();
c1->Draw();
1.2 刻度文件能谱与$^{152}Eu$（黑色字体标注）、$^{133}Ba$源对应图（红色字体标注）¶
- 发现李老师让找出的h0中的几个峰并非能量，而是道址！
1.3 Gauss拟合，抽取峰位信息（需要误差！！！）¶
注意：mean error 不等于sigma：高斯拟合中，sigma通常与分辨相关。
In [2]:
// 对某个峰进行高斯拟合
h_gammaraw->Fit("gaus","","",1215,1219);// 这里 1215 为拟合的左边界，1219 为拟合的右边界
h_gammaraw->Draw();
h_gammaraw->GetXaxis()->SetRangeUser(1212,1222);//设置显示的区间
c1->Draw();
// 拟合需要选择合适的拟合区间。
// 由拟合输出信息可以获得峰位及误差信息。
// 这里需要注意，选取合适的拟合区间对于峰位的精确指认至关重要。评判拟合好坏的一个标准是拟合曲线与实际曲线越贴近越好。
In [3]:
h_gammaraw->GetXaxis()->SetRangeUser(0,2500);//返回全区间，在多次进行拟合时，必须要有这个操作，否则指定的拟合区间将存在问题
h_gammaraw->Fit("gaus","","",966.,970);//指定拟合的左右边界， 在Fit函数的第二个参数设为"L"时进行最大似然法拟合
h_gammaraw->Draw();
h_gammaraw->GetXaxis()->SetRangeUser(965,971);//设置显示的区间
c1->Draw();
峰位拟合结果与能量对应关系如下：
| 编号 | channel | error | energy | 占比 |
|---|---|---|---|---|
| 1 | 1.01646E+02 | 3.75822E-04 | 81 | 0.3331 |
| 2 | 1.35970E+02 | 3.76597E-04 | 121.78 | 0.256 |
| 3 | 2.39196E+02 | 7.60790E-04 | 244.7 | 0.076 |
| 4 | 2.65806E+02 | 9.04701E-04 | 276.4 | 0.0713 |
| 5 | 2.88047E+02 | 6.29902E-04 | 302.9 | 0.1831 |
| 6 | 3.22884E+02 | 3.68731E-04 | 344.28 | 0.265 |
| 7 | 3.32758E+02 | 2.86916E-04 | 356 | 0.625 |
| 8 | 3.56179E+02 | 8.59586E-04 | 383.8 | 0.0894 |
| 9 | 3.79109E+02 | 1.81346E-03 | 411.12 | 0.022 |
| 10 | 4.06696E+02 | 1.40930E-03 | 443.96 | 0.031 |
| 11 | 6.88208E+02 | 9.17790E-04 | 778.9 | 0.129 |
| 12 | 7.62608E+02 | 1.61396E-03 | 867.37 | 0.042 |
| 13 | 8.43845E+02 | 7.40902E-04 | 964.8 | 0.146 |
| 14 | 9.46206E+02 | 9.85851E-04 | 1085.9 | 0.102 |
| 15 | 9.68231E+02 | 8.68893E-04 | 1112.1 | 0.136 |
| 16 | 1.21705E+03 | 7.78963E-04 | 1405 | 0.21 |
此处也验证了，李老师建议找的几个峰为道址而非能量。
1.4 将拟合所得峰位（带误差）与能量存到TGraphErrors图中¶
In [4]:
TCanvas *c_cali = new TCanvas("c_cali");
const Int_t number = 16;
Float_t cali_ch[number]={1.01646E+02,1.35970E+02,2.39196E+02,2.65806E+02,2.88047E+02,3.22884E+02,3.32758E+02,3.56179E+02,3.79109E+02,4.06696E+02,6.88208E+02,7.62608E+02,8.43845E+02,9.46206E+02,9.68231E+02,1.21705E+03};
Float_t cali_ch_error[number]={3.75822E-04,3.76597E-04,7.60790E-04,9.04701E-04,6.29902E-04,3.68731E-04,2.86916E-04,8.59586E-04,1.81346E-03,1.40930E-03,9.17790E-04,1.61396E-03,7.40902E-04,9.85851E-04,8.68893E-04,7.78963E-04};
Float_t cali_E[number]={81.,121.78,244.7,276.4,302.9,344.28,356.,383.8,411.12,443.96,778.9,867.37,964.8,1085.9,1112.1,1405.};
Float_t cali_E_error[number]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
TGraphErrors * gamma_cali = new TGraphErrors(number,cali_ch,cali_E,cali_ch_error,cali_E_error);//声明TGraphError，指定点的数目，以及(x,y),(sigma_x,sigma_y)
gamma_cali->SetTitle("gamma calibration of 152Eu and 133Ba");
gamma_cali->SetMarkerColor(4);
gamma_cali->SetMarkerStyle(21);
gamma_cali->GetXaxis()->SetRangeUser(0,1500);
gamma_cali->GetYaxis()->SetRangeUser(0,1500);
gamma_cali->Draw("ALP");
c_cali->Draw();
1.5 分别用pol1与pol2拟合gamma_cali¶
In [5]:
gamma_cali->Fit("pol1","","",0,1500);
gamma_cali->Draw("ALP");
c_cali->Draw();
In [6]:
/*
gamma_cali->Fit("pol2","Q","",200,500);
gamma_cali->Fit("pol2","Q","",200,1000);
gamma_cali->Fit("pol2","Q","",200,1300);
gamma_cali->Fit("pol2","","",0,1300);
gamma_cali->Draw("ALP");
c_cali->Draw();
*/
//测试发现，在执行pol1拟合后，直接用pol2拟合（不加任何条件）拟合结果明显有偏差；
//解决办法：先集中在数据点密集处小范围拟合，逐次放大拟合范围，可得到收敛结果，拟合正确；
//此处建议自定义多项式拟合函数，拟合前给定合理的初值，才能正确收敛到正确的初值。
TF1 * fitfun_pol2 = new TF1("fitfun_pol2","[0]+[1]*x+[2]*x**2",0,1500);
fitfun_pol2->SetParameters(-40.0,1.0,0);
gamma_cali->Fit("fitfun_pol2","","",0,1300);
gamma_cali->Draw("ALP");
c_cali->Draw();
1.6 利用步骤1.5得到的刻度系数，得到刻度后的gamma能谱¶
In [7]:
void cali_gamma(TH1F * h_raw,TH1D *h_cali,Double_t p0,Double_t p1,Double_t p2)
{
// h_cali = new TH1D("","",10000,0,2500);//h1为刻度后的gamma谱
TRandom3 *r = new TRandom3(0);
int Nbins = h_raw->GetXaxis()->GetNbins(); //h0为未刻度的gamma谱。
for(int i=0; i<Nbins; i++)
{
Long64_t eN = h_raw->GetBinContent(i);// GetBinContent()提取每个bin的计数
Double_t e = h_raw->GetBinLowEdge(i);// GetBinLowEdge()提取每个bin左边界的横坐标x。
for(Long64_t j=0; j<eN; j++)
{
Double_t ea = e+r->Rndm()*0.2;// 连续化能谱
ea=p2*ea*ea+p1*ea+p0;//p0,p1,p2为上面获得的二次项刻度系数。
h_cali->Fill(ea);
}
}
//h_cali->Draw();
}
In [8]:
TCanvas* c_gamma_cali_p1=new TCanvas("c_gamma_cali_p1");
TH1D * gamma_cali_p1 = new TH1D("gamma_cali_p1","",10000,0,2500);//h1为刻度后的gamma谱
cali_gamma(h_gammaraw,gamma_cali_p1,-3.96360e+01, 1.18877e+00,0);
gamma_cali_p1->SetTitle("calibrated by Pol1");
gamma_cali_p1->GetXaxis()->SetTitle("E_{#gamma}/keV");
gamma_cali_p1->GetYaxis()->SetTitle("Counts");
gamma_cali_p1->Draw();
c_gamma_cali_p1->Draw();
In [9]:
TCanvas* c_gamma_cali_p2=new TCanvas("c_gamma_cali_p2");
TH1D * gamma_cali_p2 = new TH1D("gamma_cali_p2","",10000,0,2500);//h1为刻度后的gamma谱
cali_gamma(h_gammaraw,gamma_cali_p2,-4.05649e+01, 1.19378e+00 ,-4.25799e-06);
gamma_cali_p2->SetTitle("calibrated by Pol2");
gamma_cali_p2->GetXaxis()->SetTitle("E_{#gamma}/keV");
gamma_cali_p2->GetYaxis()->SetTitle("Counts");
gamma_cali_p2->Draw();
c_gamma_cali_p2->Draw();
1.7 残差检验¶
In [10]:
/*
//查看高斯拟合是否可行:Yes!
gaus->SetParameter(0,1e6);
gaus->SetParameter(1,81.);
gaus->SetParameter(2,1.);
gamma_cali_p1->GetXaxis()->SetRangeUser(0,2500);
gamma_cali_p1->Fit("gaus","","",80,82);
gamma_cali_p1->Draw();
c_gamma_cali_p1->Draw();
Double_t get_mean=gaus->GetParameter(1);
Double_t get_meanerr=gaus->GetParError(1);
cout<<"get_peak="<<get_mean<<"+-"<<get_meanerr<<endl;
*/
//为了得到残差图，需要批量拟合，故构建高斯拟合函数
void FitGaus(TH1D *in_h,Double_t set_con, Double_t set_mean, Double_t set_sigma, Double_t fitrangel, Double_t fitranger, Double_t &get_mean, Double_t &get_meanerr)
{
gaus->SetParameter(0,set_con);
gaus->SetParameter(1,set_mean);
gaus->SetParameter(2,set_sigma);
in_h->GetXaxis()->SetRangeUser(0,2500);
in_h->Fit("gaus","Q","",set_mean-fitrangel,set_mean+fitranger);
get_mean = gaus->GetParameter(1);
get_meanerr = gaus->GetParError(1);
in_h->Draw();
}
/*
//查看自己构造的函数拟合是否可行:无法拟合，why?????????????；
Double_t func_gaus(Double_t *x, Double_t * par)
{
return par[0]*TMath::Exp(-0.5*((x[0]-par[1])/par[2])*((x[0]-par[1])/par[2]))+par[3]+x[0]*par[4];
}
TF1 *fit_gaus = new TF1("fit_gaus",func_gaus, 0 , 2500); //使用一个高斯峰
fit_gaus->SetParameters(0,7.0e5);
fit_gaus->SetParameters(1,81.0);
fit_gaus->SetParameters(2,1.0);
fit_gaus->SetParameters(3,1.0e4);
fit_gaus->SetParameters(4,0.01);
gamma_cali_p1->GetXaxis()->SetRangeUser(0,2500);
gamma_cali_p1->Fit("fit_gaus","","",80,82);
gamma_cali_p1->Draw();
c_gamma_cali_p1->Draw();
*/
In [11]:
//拟合由pol1刻度的gamma能谱。
//初始化拟合参量
Double_t fit_in_mean[16] = {81.,121.78,244.7,276.4,302.9,344.28,356.,383.8,411.12,443.96,778.9,867.37,964.8,1085.9,1112.1,1405.};//gamma_energy[16]
Double_t fit_in_constant[16] = {1e6,1e6,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5,1e5};
Double_t fit_in_rangel[16] = {2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,1.0};
Double_t fit_in_ranger[16] = {2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,4.0};
Double_t fitpol1_get_mean[16] ;
Double_t fitpol1_get_meanerr[16];
Double_t pol1_res[16];
/*
//单个检查
FitGaus(gamma_cali_p1,fit_in_constant[15], fit_in_mean[15], 1.0, fit_in_rangel[15], fit_in_ranger[15], fit_get_mean[15], fit_get_meanerr[15]);
cout<<"get_peak="<<fit_get_mean[15]<<"+-"<<fit_get_meanerr[15]<<endl;
c_gamma_cali_p1->Draw();
*/
//批量拟合
for (Int_t i=0; i<16; i++)
{
//cout<<"E="<<fit_in_mean[i]<<" keV,";
FitGaus(gamma_cali_p1,fit_in_constant[i], fit_in_mean[i], 1.0, fit_in_rangel[i], fit_in_ranger[i], fitpol1_get_mean[i], fitpol1_get_meanerr[i]);
// cout<<"get_peak="<<fit_get_mean[i]<<"+-"<<fit_get_meanerr[i]<<endl;
pol1_res[i]=fitpol1_get_mean[i]-fit_in_mean[i];
}
TCanvas *residual_pol1 = new TCanvas("residual_pol1");
TGraphErrors *g_pol1_res = new TGraphErrors(16,fit_in_mean,pol1_res);
g_pol1_res->SetTitle("gamma energy Residual fiited by pol1 or pol2");
g_pol1_res->GetXaxis()->SetTitle("E_{#gamma}/keV");
g_pol1_res->GetYaxis()->SetTitle("#DeltaE_{#gamma}/keV");
g_pol1_res->SetMarkerStyle(20);
g_pol1_res->SetLineColor(kRed);
//g_pol1_res->Draw("ALP");
//residual_pol1->Draw();
In [12]:
//拟合由pol2刻度的gamma能谱。
Double_t fitpol2_get_mean[16] ;
Double_t fitpol2_get_meanerr[16];
Double_t pol2_res[16];
//批量拟合
for (Int_t i=0; i<16; i++)
{
//cout<<"E="<<fit_in_mean[i]<<" keV,";
FitGaus(gamma_cali_p2,fit_in_constant[i], fit_in_mean[i], 1.0, fit_in_rangel[i], fit_in_ranger[i], fitpol2_get_mean[i], fitpol2_get_meanerr[i]);
//cout<<"get_peak="<<fitpol2_get_mean[i]<<"+-"<<fitpol2_get_meanerr[i]<<endl;
pol2_res[i]=fitpol2_get_mean[i]-fit_in_mean[i];
}
//TCanvas *residual_pol2 = new TCanvas("residual_pol2");
TGraphErrors *g_pol2_res = new TGraphErrors(16,fit_in_mean,pol2_res);
//g_pol2_res->SetTitle("gamma energy Residual fiited by pol2");
//g_pol2_res->GetXaxis()->SetTitle("E_{#gamma}/keV");
//g_pol2_res->GetYaxis()->SetTitle("#DeltaE_{#gamma}/keV");
g_pol2_res->SetMarkerStyle(0);
g_pol2_res->SetLineColor(kGreen);
g_pol1_res->Draw("");
g_pol1_res->GetYaxis()->SetRangeUser(-2.5,2.5);
g_pol1_res->GetXaxis()->SetRangeUser(0,1500);
g_pol2_res->Draw("same");
auto legend = new TLegend(0.2, 0.6, 0.4, 0.8);
legend->AddEntry(g_pol1_res, "pol1", "L");
legend->AddEntry(g_pol2_res, "pol2", "L");
legend->SetTextSize(0.03);
legend->Draw("same");
residual_pol1->Draw();
//residual_pol2->Draw();
1.7 小结¶
- 一次拟合低能区残差较小，高能区残差普遍很大，甚至存在残差大于1的情况；
- 二次拟合与一次拟合类似（低能区残差小，高能区残差大）。但与一次拟合相比，二次在低能区残差更大，而高能区残差相对小一些，即二次拟合的残差相对更集中；
2. 探测器峰宽刻度¶
2.1 构造高斯+一次多项式拟合函数¶
- 由于计算峰位面积、提取半宽时需要考虑本底的影响，因此仅仅用高斯拟合不再适用。为了简单，我们仅仅考虑线性本地，故需要构造高斯+一次多项式拟合函数。
- 在第1节能量刻度结果中发现，二次多项式刻度后的残差相对集中，故本节探测器峰宽刻度使用gamma_cali_p2图。
In [13]:
double fungauspol1(double *x, double *par)//构造函数：高斯函数par[0-2]+线性函数par[3-4]
{
return par[0]*TMath::Exp(-0.5*((x[0]-par[1])/par[2]) * ((x[0]-par[1])/par[2])) + par[3]+x[0]*par[4];
}
2.2 拟合gamma_cali_p2抽取半高宽¶
In [14]:
//gamma_energy of source,参考用，为了拟合时确定初始花参量用；
//Double_t gamma_energy[16] = {81.,121.78,244.7,276.4,302.9,344.28,356.,383.8,411.12,443.96,778.9,867.37,964.8,1085.9,1112.1,1405.};
TCanvas *c_fit_gauspol1 = new TCanvas("c_fit_gauspol1");
TF1 * fitgauspol1 = new TF1("fitgauspol1",fungauspol1,0,2500,5);//构造拟合曲线
//初始化拟合参量
fitgauspol1->SetParameter(0,5.0e4);//par[0]：高斯函数中峰值（Cons）;
fitgauspol1->SetParameter(1,1405.);//par[1]：高斯函数中峰位（Mean）;
fitgauspol1->SetParameter(2,1.0);//par[2]：高斯函数中半宽（sigma）;
fitgauspol1->SetParameter(3,6.0e3);//par[3]：一次多项式函数中零次项（a）;
fitgauspol1->SetParameter(4,13.1);//par[3]：一次多项式函数中一次项（b）;
//拟合过程
gamma_cali_p2->GetXaxis()->SetRangeUser(0,2500);//画出全区间；
gamma_cali_p2->Fit("fitgauspol1","","",1398,1411);
gamma_cali_p2->Draw();
gamma_cali_p2->GetXaxis()->SetRangeUser(1395,1415);//放大拟合区间，查看拟合效果；
c_fit_gauspol1->Draw();
高斯+一次多项式拟合抽取半高宽结果：
| 序号 | mean | meanerror | sigma | sigmaerror | 半宽 | 半宽error |
|---|---|---|---|---|---|---|
| （keV） | （keV） | （keV） | （keV） | （keV） | （keV） | |
| 1 | 8.07328E+01 | 4.87834E-04 | 7.71879E-01 | 4.89804E-04 | 1.81778E+00 | 1.15349E-03 |
| 2 | 1.21674E+02 | 3.92216E-04 | 7.51012E-01 | 3.67888E-04 | 1.76863E+00 | 8.66376E-04 |
| 3 | 2.44734E+02 | 9.56016E-04 | 7.83486E-01 | 9.46195E-04 | 1.84511E+00 | 2.22829E-03 |
| 4 | 2.76466E+02 | 1.13022E-03 | 7.97680E-01 | 1.14891E-03 | 1.87854E+00 | 2.70568E-03 |
| 5 | 3.02949E+02 | 6.13366E-04 | 8.11896E-01 | 5.74021E-04 | 1.91202E+00 | 1.35182E-03 |
| 6 | 3.44444E+02 | 4.40927E-04 | 8.31793E-01 | 3.88668E-04 | 1.95887E+00 | 9.15313E-04 |
| 7 | 3.56204E+02 | 3.08243E-04 | 8.29596E-01 | 2.64285E-04 | 1.95370E+00 | 6.22391E-04 |
| 8 | 3.84092E+02 | 9.39398E-04 | 8.43247E-01 | 8.78283E-04 | 1.98585E+00 | 2.06836E-03 |
| 9 | 4.11380E+02 | 2.07285E-03 | 8.44293E-01 | 2.10468E-03 | 1.98831E+00 | 4.95652E-03 |
| 10 | 4.44235E+02 | 1.68215E-03 | 8.61726E-01 | 1.65940E-03 | 2.02936E+00 | 3.90789E-03 |
| 11 | 7.78979E+02 | 8.88147E-04 | 9.58035E-01 | 8.20380E-04 | 2.25617E+00 | 1.93199E-03 |
| 12 | 8.67345E+02 | 1.88897E-03 | 9.84850E-01 | 1.90356E-03 | 2.31932E+00 | 4.48288E-03 |
| 13 | 9.63777E+02 | 8.82106E-04 | 9.85453E-01 | 7.78298E-04 | 2.32074E+00 | 1.83289E-03 |
| 14 | 1.08516E+03 | 1.21073E-03 | 9.84573E-01 | 1.54390E-03 | 2.31867E+00 | 3.63588E-03 |
| 15 | 1.11131E+03 | 1.07066E-03 | 1.02714E+00 | 1.04491E-03 | 2.41891E+00 | 2.46076E-03 |
| 16 | 1.40600e+03 | 8.61315e-04 | 1.15341e+00 | 7.12661e-04 | 2.71628E+00 | 1.67832E-03 |
2.3 画出半宽-能量的二维散点图(TGraph)¶
In [15]:
TCanvas *c_energyvsfwhm=new TCanvas("c_energyvsfwhm");
Double_t gamma_energy_calibrated[12]={1.21674E+02,2.44734E+02,2.76466E+02,3.02949E+02,3.44444E+02,3.56204E+02,3.84092E+02,4.11380E+02,4.44235E+02,7.78979E+02,8.67345E+02,1.40600e+03};
Double_t gamma_energy_fwhm[12]={1.76863E+00,1.84511E+00,1.87854E+00,1.91202E+00,1.95887E+00,1.95370E+00,1.98585E+00,1.98831E+00,2.02936E+00,2.25617E+00,2.31932E+00,2.71628E+00};
TGraph * g_evsfwhm = new TGraph(12,gamma_energy_calibrated,gamma_energy_fwhm);
g_evsfwhm->SetMarkerStyle(50);
g_evsfwhm->SetLineColor(kGreen);
g_evsfwhm->GetYaxis()->SetRangeUser(1.75,2.8);
g_evsfwhm->GetXaxis()->SetRangeUser(0,1500);
g_evsfwhm->SetTitle("the FWHM fitted by POL2");
g_evsfwhm->Draw();
g_evsfwhm->Fit("pol2","","",100,1410);
c_energyvsfwhm->Draw();
小节：拟合发现，部分数据点偏离较大，对整体拟合影响较大。上述结果去除了的序号为1、13、14、15的数据点的拟合结果。
2.4 画出半宽拟合残差图¶
In [16]:
//2.2 中高斯拟合得到
//Double_t gamma_energy_calibrated[12]={1.21674E+02,2.44734E+02,2.76466E+02,3.02949E+02,3.44444E+02,3.56204E+02,3.84092E+02,4.11380E+02,4.44235E+02,7.78979E+02,8.67345E+02,1.40600e+03};
//Double_t gamma_energy_fwhm[12]={1.76863E+00,1.84511E+00,1.87854E+00,1.91202E+00,1.95887E+00,1.95370E+00,1.98585E+00,1.98831E+00,2.02936E+00,2.25617E+00,2.31932E+00,2.71628E+00};
Double_t fwhmvsenergy_res[12];//残差变量声明，用于分析二次拟合优略；
Double_t pol2_0,pol2_1,pol2_2;//二次多项式拟合参数声明，用于计算拟合结果；
pol2_0=pol2->GetParameter(0);//二次多项式拟合参数赋值；
pol2_1=pol2->GetParameter(1);//二次多项式拟合参数赋值；
pol2_2=pol2->GetParameter(2);//二次多项式拟合参数赋值；
//计算残差
for(Int_t i=0;i<12;i++)
{
fwhmvsenergy_res[i]=pol2_0+pol2_1*gamma_energy_calibrated[i]+pol2_2*gamma_energy_calibrated[i]*gamma_energy_calibrated[i]-gamma_energy_fwhm[i];
}
//画出残差图
TCanvas *c_fwhm_res=new TCanvas("c_fwhm_res");
TGraph *g_fwhm_res=new TGraph(12,gamma_energy_calibrated,fwhmvsenergy_res);
g_fwhm_res->SetMarkerStyle(50);
g_fwhm_res->SetLineColor(kBlack);
g_fwhm_res->SetTitle("the residual of FWHM fitted by POL2");
g_fwhm_res->GetXaxis()->SetTitle("E_{#gamma}/keV");
g_fwhm_res->GetYaxis()->SetTitle("Residual/keV");
g_fwhm_res->GetXaxis()->SetRangeUser(0,1500);
//g_fwhm_res->GetYaxis()->SetRangeUser(-2.5,2.5);
g_fwhm_res->Draw("ALP");
c_fwhm_res->Draw();
小节：残差结果整体看，结果比较集中，但个别数据点偏离相对较大；
3. gamma峰面积¶
3.1 效率的绝对刻度¶
- 思想：测量到的某个能量gamma的计数与测量时间内源发射的该能量gamma的比值。
- 要求：源性能要好，且相关数据准确，如，源种类（半衰期）、出厂日期、出厂活度、测量日期、测量时长等等；
- 效率刻度方法：
$\varepsilon(E_{\gamma})=\frac{N(E_{\gamma})_m}{N(E_{\gamma})_s}$;
$\varepsilon(E_{\gamma})$:能量为$E_{\gamma}$时的探测器效率；
$N(E_{\gamma})_m)$:测量到能量为$E_{\gamma}$的gamma计数，由3.1步卡出；
$N(E_{\gamma})_s)$:源在测试时间内发射的能量为$E_{\gamma}$的gamma计数,有$N(E_{\gamma})_s)=A_0e^{-\lambda}P_{\gamma}\Delta t$;
$A_0$:源活度；
$\lambda$:衰变常数；
$t$:出厂日期到测试日期的时间差；
$\Delta t$:刻度用时；
$P_{\gamma}$:绝对$\gamma$射线发射概率；
3.2 借由2.2的结果抽取峰面积（gamma能量峰计数）¶
- 应用高斯加一次多项式拟合得到的Mean与Sigma参数，选取Mean$\pm$1.5Sigma范围积分，得到峰位总计数S1；
- 在峰附近左右两侧各选择1.5*sigma范围进行积分，获得S2,S3，用于估计峰的本底；
- gamma峰面积S=S1-S2-S3； -使用gamma_cali_p2图;
In [17]:
//2.2 中高斯拟合得到
//Double_t gamma_energy_calibrated[12]={1.21674E+02,2.44734E+02,2.76466E+02,3.02949E+02,3.44444E+02,3.56204E+02,3.84092E+02,4.11380E+02,4.44235E+02,7.78979E+02,8.67345E+02,1.40600e+03};
//Double_t gamma_energy_fwhm[12]={1.76863E+00,1.84511E+00,1.87854E+00,1.91202E+00,1.95887E+00,1.95370E+00,1.98585E+00,1.98831E+00,2.02936E+00,2.25617E+00,2.31932E+00,2.71628E+00};//由2.2移殖而来，需要用到2.1构造的函数；
Int_t S[12],S1[12],S2[12],S3[12];//声明峰面积变量；
Int_t binL,binR;//声明左右bin变量；
for (Int_t i=0;i<12;i++)
{
binL=gamma_cali_p2->FindBin(gamma_energy_calibrated[i]-1.5*gamma_energy_fwhm[i]/2.355);
binR=gamma_cali_p2->FindBin(gamma_energy_calibrated[i]+1.5*gamma_energy_fwhm[i]/2.355);
S1[i]=gamma_cali_p2->Integral(binL,binR);
cout<<"S1="<<S1[i]<<"; ";
binL=gamma_cali_p2->FindBin(gamma_energy_calibrated[i]-4.0*gamma_energy_fwhm[i]/2.355);
binR=gamma_cali_p2->FindBin(gamma_energy_calibrated[i]-2.5*gamma_energy_fwhm[i]/2.355);
S2[i]=gamma_cali_p2->Integral(binL,binR);
cout<<"S2="<<S2[i]<<"; ";
binL=gamma_cali_p2->FindBin(gamma_energy_calibrated[i]+2.5*gamma_energy_fwhm[i]/2.355);
binR=gamma_cali_p2->FindBin(gamma_energy_calibrated[i]+4.0*gamma_energy_fwhm[i]/2.355);
S3[i]=gamma_cali_p2->Integral(binL,binR);
cout<<"S3="<<S3[i]<<"; ";
S[i]=S1[i]-S2[i]-S3[i];
cout<<"S="<<S[i]<<endl;
}
3.3 源发射粒子数¶
In [18]:
Double_t t_diff=5522.0*24.0*60.0*60.0;//单位：s；start=1998-1-1；measure time:2013-2-13;
Double_t A_152=40900.0;//152Eu活度:40.9 kBq, 1 Bq=1 次/s;
Double_t d_A152=0.05;//152Eu活度不确定度；
Double_t A_133=42200.0;//133Eu活度:42.2 kBq, 1 Bq=1 次/s;
Double_t d_A133=0.03;//133Ba活度不确定度；
Double_t t_mea=7442.0;//单位：s;
Double_t T_152=13.2*365*24*60*60;//单位：s；152Er半衰期：13.2～0.3年；
Double_t T_133=10.59*365*24*60*60;//单位：s；133Ba半衰期：10.59～0.11年；
Double_t P[12]={0.256,0.076,0.0713,0.1831,0.265,0.625,0.0894,0.022,0.031,0.129,0.042,0.21};//选取各gamma峰占比；
Double_t eff[12];
Double_t N0_152=A_152*TMath::Exp(-log(2)*t_diff/T_152)*t_mea;
//cout<<N0_152<<endl;
Double_t N0_133=A_133*TMath::Exp(-log(2)*t_diff/T_133)*t_mea;
//cout<<N0_133<<endl;
for (Int_t i=0;i<12;i++)
{
if(i==3||i==5||i==6)
{
eff[i]=S[i]/(N0_133*P[i]);
}
else
{
eff[i]=S[i]/(N0_152*P[i]);
}
//cout<<eff[i]<<endl;
}
TCanvas *c_effvsenergy=new TCanvas("c_effvsenergy");
TGraph * g_eff=new TGraph(12,gamma_energy_calibrated,eff);
g_eff->SetMarkerStyle(50);
g_eff->SetLineColor(kBlack);
g_eff->SetTitle("the efficiency of gamma with 152Eu and 133Ba");
g_eff->GetXaxis()->SetTitle("E_{#gamma}/keV");
g_eff->GetYaxis()->SetTitle("Efficiency");
g_eff->GetXaxis()->SetRangeUser(0,1500);
//g_fwhm_res->GetYaxis()->SetRangeUser(-2.5,2.5);
g_eff->Draw("ALP");
c_effvsenergy->Draw();
In [19]:
!jupyter nbconvert HPGecalibration.ipynb --to html
In [ ]: