Go to the source code of this file.
Functions | |
| void | ge_eu152_calib () |
| void ge_eu152_calib | ( | ) |
Definition at line 1 of file ge_eu152_calib.c.
00001 { 00002 // The file to run on in the sum of calibration runs 3375, 3378 and 3379. 00003 TFile* f = new TFile("calib_ge.root", "READ"); 00004 TSpectrum* s = new TSpectrum(100, 3.); 00005 TH1* spec = (TH1*)f->Get("PlotTAP_Amplitude/hGe-S#FirstComplete#{constant_fraction=0.60}{no_time_shift=true}_Amplitude"); 00006 spec->Sumw2(); 00007 00008 // Activity 00009 unsigned int src_date[] = {1, 3, 2000}; // 1 March 2000? 00010 unsigned int meas_date[] = {19, 12, 2013}; // 19 Dec 2013 00011 Double_t age = 13.8; // Years according to WolframAlpha 00012 Double_t halflife = 13.537; 00013 Double_t activity0 = 46.4e3; // RunPSI2013:661 00014 TTimeStamp t0(src_date[2], src_date[1], src_date[0], 0, 0, 0); 00015 TTimeStamp t(meas_date[2], meas_date[1], meas_date[0], 0, 0, 0); 00016 Double_t dt = (t.GetSec() - t0.GetSec())/60./60./24./365.2425; 00017 Double_t activity = activity0*TMath::Power(0.5,dt/halflife); 00018 Double_t livetime = 3245.5; // Analysis-R13:129 00019 Double_t events = activity*livetime; 00020 00021 // Spectrum 00022 const unsigned int NPeaks = 15; 00023 const unsigned int Peak511Index = 5, Peak122Index = 0; 00024 Double_t Energy[NPeaks] = {121.7817, 244.6975, 344.2785, 411.1163, 443.965, 510.999, 778.9040, 867.378, 964.079, 1085.869, 1089.737, 1112.074, 1212.948, 1299.140, 1408.006}; 00025 Double_t RedEn[NPeaks-1] = {121.7817, 244.6975, 344.2785, 411.1163, 443.965, 778.9040, 867.378, 964.079, 1085.869, 1089.737, 1112.074, 1212.948, 1299.140, 1408.006}; 00026 Double_t Intensity[NPeaks-1] = {0.2858, 0.07583, 0.265, 0.02234, 0.02821, 0.12942, 0.04245, 0.14605, 0.10207, 0.01726, 0.13644, 0.01422, 0.01623, 0.21005}; 00027 Double_t LogEnergy[NPeaks]; for (unsigned int i = 0; i < NPeaks; ++i) LogEnergy[i] = TMath::Log(Energy[i]); 00028 Double_t ExpectedCounts[NPeaks-1]; for (unsigned int i = 0; i < NPeaks-1; ++i) ExpectedCounts[i] = events*Intensity[i]; 00029 Double_t ADC[NPeaks] = {1000.75, 2015.48, 2831.37, 3379.41, 3648.90, 4196.38, 6392.91, 7118.00, 7909.78, 8907.53, 8941.08, 9121.85, 9948.85, 10654.3, 11547.5}; // First guess 00030 Double_t ADC_Meas[NPeaks]; 00031 const char* Functions[NPeaks] = { "gaus(0)+[3]*x+[4]", 00032 "gaus(0)+gaus(3)+[6]*x+[7]", 00033 "gaus(0)+gaus(3)+[6]*x+[7]", 00034 "gaus(0)+[3]*x+[4]", 00035 "gaus(0)+[3]*x+[4]", 00036 "gaus(0)+[3]*x+[4]", 00037 "gaus(0)+[3]*x+[4]", 00038 "gaus(0)+[3]*x+[4]", 00039 "gaus(0)+[3]*x+[4]", 00040 "gaus(0)+gaus(3)+[6]*x+[7]", 00041 "gaus(0)+gaus(3)+[6]*x+[7]", 00042 "gaus(0)+gaus(3)+[6]*x+[7]", 00043 "gaus(0)+[3]*x+[4]", 00044 "gaus(0)+gaus(3)+[6]*x+[7]", 00045 "gaus(0)+[3]*x+[4]" }; 00046 const unsigned int NParam[NPeaks] = { 5, 8, 8, 5, 5, 5, 5, 5, 5, 8, 8, 8, 5, 8, 5 }; 00047 const Double_t FitParam[NPeaks][11] = { {1250., ADC[0], 8., 0., 0., 0., 0., 0., 0., 0., 0.}, 00048 {200., ADC[1], 8., 20., 1966., 8., 0., 100., 0., 0., 0.}, 00049 {70., ADC[2], 8., 10., 2894., 8., 0., 70., 0., 0., 0.}, 00050 {50., ADC[3], 8., 0., 40., 0., 0., 0., 0., 0., 0.}, 00051 {55., ADC[4], 8., 0., 35., 0., 0., 0., 0., 0., 0.}, 00052 {450., ADC[5], 10., 0., 25., 0., 0., 0., 0., 0., 0.}, 00053 {125., ADC[6], 8., 0., 15., 0., 0., 0., 0., 0., 0.}, 00054 {40., ADC[7], 8., 0., 10., 0., 0., 0., 0., 0., 0.}, 00055 {120., ADC[8], 8., 0., 10., 0., 0., 0., 0., 0., 0.}, 00056 {75., ADC[9], 8., 15., ADC[10], 8., 0., 5., 0., 0., 0.}, 00057 {15., ADC[10], 8., 75., ADC[9], 8., 0., 5., 0., 0., 0.}, 00058 {95., ADC[11], 8., 10., 9170., 8., 0., 5., 0., 0., 0.}, 00059 {10., ADC[12], 8., 0., 5., 0.. 0., 0., 0., 0., 0.}, 00060 {10., ADC[13], 8., 25., 10610., 8., 0., 2., 0., 0., 0.}, 00061 {120., ADC[14], 8., 0., 0., 0., 0., 0., 0., 0., 0.} }; 00062 Double_t FitParamRes[NPeaks][11]; 00063 Double_t FitParamErr[NPeaks][11]; 00064 Double_t Counts[NPeaks-1]; 00065 Double_t Efficiency[NPeaks-1]; 00066 Double_t LogEfficiency[NPeaks-1]; 00067 00068 // Errors 00069 Double_t Chi2Ndf[NPeaks]; 00070 Double_t Error_ADC[NPeaks]; 00071 Double_t Error_Counts[NPeaks-1]; 00072 Double_t Error_Eff[NPeaks-1]; 00073 Double_t Error_LogEff[NPeaks-1]; 00074 00075 // Count the events (integral of fit gaussian) 00076 const Double_t rt2pi = TMath::Sqrt(2.*TMath::Pi()); 00077 for (unsigned int i = 0; i < NPeaks; ++i) { 00078 TF1* fit = new TF1("fit", Functions[i]); 00079 const Double_t fit_window = 100.; 00080 spec->GetXaxis()->SetRangeUser(ADC[i] - fit_window, ADC[i] + fit_window); 00081 fit->SetParameters(FitParam[i]); 00082 TFitResultPtr res = spec->Fit(fit,"SMENQ"); 00083 Chi2Ndf[i] = res->Chi2()/(Double_t)res->Ndf(); 00084 for (unsigned int j = 0; j < NParam[i]; ++j) { 00085 FitParamRes[i][j] = res->Value(j); 00086 FitParamErr[i][j] = res->ParError(j); 00087 } 00088 } 00089 for (unsigned int i = 0; i < NPeaks; ++i) { 00090 // Get energy calibration 00091 ADC_Meas[i] = FitParamRes[i][1]; 00092 Error_ADC[i] = FitParamErr[i][1]; 00093 // Get efficiency calibration 00094 // Adjust for lack of prediction for 511 peak 00095 unsigned int j = i; 00096 if (i == Peak511Index) 00097 continue; 00098 else if (i > Peak511Index) 00099 --j; 00100 Counts[j] = FitParamRes[i][0]*FitParamRes[i][2]*rt2pi; 00101 Error_Counts[j] = Counts[j]*TMath::Sqrt(TMath::Power(FitParamErr[i][0]/FitParamRes[i][0],2.)+TMath::Power(FitParamErr[i][2]/FitParamRes[i][2],2.)); 00102 Efficiency[j] = Counts[j] / ExpectedCounts[j]; 00103 Error_Eff[j] = Error_Counts[j] / ExpectedCounts[j]; 00104 LogEfficiency[j] = TMath::Log(Efficiency[j]); 00105 Error_LogEff[j] = Error_Eff[j] / Efficiency[j]; 00106 } 00107 00108 // Draw 00109 Double_t en_par[2], en_par_err[2]; 00110 Double_t eff_par[2], eff_par_err[2], Chi2[2]; 00111 unsigned int Ndf[2]; 00112 /* TGraphErrors* e_cal = new TGraphErrors(NPeaks, Energy, ADC_Meas, 0, Error_ADC); */ 00113 TGraphErrors* e_cal = new TGraphErrors(NPeaks-1, ADC_Meas+1, Energy+1, Error_ADC+1, 0); 00114 /* TGraphErrors* eff_cal = new TGraphErrors(NPeaks-1, LogEnergy+1, LogEfficiency+1, 0, Error_LogEff+1); */ 00115 TGraphErrors* eff_cal = new TGraphErrors(NPeaks-1, Energy+1, Efficiency+1, 0, Error_Eff+1); 00116 e_cal->SetTitle("Ge Energy Calibration"); e_cal->GetXaxis()->SetTitle("ADC Value [ADC]"); e_cal->GetYaxis()->SetTitle("Energy [keV]"); 00117 eff_cal->SetTitle("Ge Efficiency Calibration"); eff_cal->GetXaxis()->SetTitle("Energy [keV]"); eff_cal->GetYaxis()->SetTitle("Counts / Expected Counts"); 00118 TCanvas* c_en = new TCanvas("c_en"); 00119 TF1* fit = new TF1("fit", "[0] + [1]*x"); 00120 fit->SetParameter(0, 0); 00121 fit->SetParameter(1, 100); 00122 fit->SetParName(0, "Offset"); 00123 fit->SetParName(1, "Gradient"); 00124 TFitResultPtr res = e_cal->Fit(fit,"SME"); e_cal->Draw("A*"); 00125 TPaveText* text_en = new TPaveText(0.2, 0.65, 0.5, 0.75, "nb NDC"); text_en->SetBorderSize(0); text_en->SetTextSize(0.04); text_en->SetFillColor(kWhite); text_en->AddText("Energy = Offset + Gradient * ADC"); text_en->Draw(); 00126 Chi2[0] = res->Chi2(); Ndf[0] = res->Ndf(); 00127 en_par[0] = res->Value(0); en_par[1] = res->Value(1); 00128 en_par_err[0] = res->ParError(0); en_par_err[1] = res->ParError(1); 00129 TCanvas* c_eff = new TCanvas("c_eff"); 00130 // gStyle->SetOptFit(111); 00131 TF1* logfit = new TF1("func", "[0]*(x**[1])"); 00132 logfit->SetParName(0, "Coefficient"); 00133 logfit->SetParName(1, "Exponent"); 00134 logfit->SetParameters(TMath::Exp(-2.6), -0.84); 00135 res = eff_cal->Fit(logfit, "SME"); eff_cal->Draw("A*"); 00136 TPaveText* text_eff = new TPaveText(0.5, 0.5, 0.65, 0.6, "nb NDC"); text_eff->SetBorderSize(0); text_eff->SetTextSize(0.05); text_eff->AddText("#varepsilon = aE^{b}"); text_eff->SetFillColor(kWhite); text_eff->Draw(); 00137 /* res = eff_cal->Fit("pol1", "SME"); eff_cal->Draw("A*"); */ 00138 Chi2[1] = res->Chi2(); Ndf[1] = res->Ndf(); 00139 eff_par[0] = res->Value(0); eff_par[1] = res->Value(1); 00140 eff_par_err[0] = res->ParError(0); eff_par_err[1] = res->ParError(1); 00141 // Print out 00142 printf("\n"); 00143 printf(" %-10s | %-8s | %-5s | %-10s | %-5s | %-7s\n", "Energy", "ADC", "+/-", "Eff", "+/-", "Chi2"); 00144 for (unsigned int i = 0; i < NPeaks; ++i) { 00145 int adc_exp = TMath::Floor(TMath::Log10(Error_ADC[i])); 00146 int eff_exp = TMath::Floor(TMath::Log10(Error_Eff[i])); 00147 printf(" %-10.10g | %-8.10g | %-5.1g | %-10.10g | %-5.1g | %-7.3g\n", 00148 Energy[i], 00149 TMath::Floor(ADC_Meas[i] / TMath::Power(10., adc_exp)) * TMath::Power(10., adc_exp), 00150 Error_ADC[i], 00151 TMath::Floor(Efficiency[i] / TMath::Power(10., eff_exp)) * TMath::Power(10., eff_exp), 00152 Error_Eff[i], 00153 Chi2Ndf[i]); 00154 } 00155 printf("\n"); 00156 printf("ADC = %g(%g) * Energy + %g(%g) (%g/%u)\n", en_par[1], en_par_err[1], en_par[0], en_par_err[0], Chi2[0], Ndf[0]); 00157 /* printf("Eff = %g(%g)*Energy^[%g(%g)] (%g/%u)\n", TMath::Exp(eff_par[0]), TMath::Exp(eff_par[0])*eff_par_err[0], eff_par[1], eff_par_err[1], Chi2[1], Ndf[1]); */ 00158 printf("Eff = %g(%g)*Energy^[%g(%g)] (%g/%u)\n", eff_par[0], eff_par_err[0], eff_par[1], eff_par_err[1], Chi2[1], Ndf[1]); 00159 00160 // c_en->SaveAs("~/plots/ThesisPlots/ge-calibration-curve.pdf"); 00161 // c_eff->SaveAs("~/plots/ThesisPlots/ge-efficiency-curve.pdf"); 00162 00163 double si_peak = 400.177; 00164 double al_peak = 346.828; 00165 std::cout << "Eff(" << si_peak << ") = " << logfit->Eval(si_peak) << std::endl; 00166 std::cout << "Eff(" << al_peak << ") = " << logfit->Eval(al_peak) << std::endl; 00167 }
1.6.1