From e126a3451bdbdca20bd5011c685d8d2e227fbfbe Mon Sep 17 00:00:00 2001 From: Cassie Date: Thu, 21 Apr 2016 17:12:56 -0400 Subject: [PATCH 1/4] Added a fitting routine to analyzer.C, can be toggled with the old version by changing ROOFIT_ON to 1. Also added hv channels, fit parameters, errors on the fit parameters and the rate of the background as branches in the analysis tree --- calibration/analyzer.C | 411 ++++++++++++++++++++++++++++++++++++++--- calibration/analyzer.h | 53 +++++- 2 files changed, 432 insertions(+), 32 deletions(-) diff --git a/calibration/analyzer.C b/calibration/analyzer.C index 648235b..9b5c726 100644 --- a/calibration/analyzer.C +++ b/calibration/analyzer.C @@ -47,6 +47,8 @@ using namespace RooFit; #define MAX_PEAKS 5 +const string analysis_path("/Users/cassiereuter/Documents/dev/Modulation/analysis/calibration/"); + /*---------------------------------------------------------------------------------------------------*/ float source_energy[NUMBER_OF_SOURCES][MAX_PEAKS] = // @@ -57,7 +59,7 @@ float source_energy[NUMBER_OF_SOURCES][MAX_PEAKS] = {511.,1157.020,511.+1157.020,-1,-1}, // ID1: Ti44 {1173.2,1332.5,1173.2+1332.5,-1,-1}, // ID2: Co60 {661.7,-1,-1,-1,-1}, // ID3: CS137 - {-1,-1,-1,-1,-1}, // ID4: MN54 + {834.,-1,-1,-1,-1}, // ID4: MN54 {1460.,-1,-1,-1,-1} // ID5: K40 }; @@ -74,6 +76,12 @@ const float emax = 3000.; // in keV const float adc_max_volt = 2.; const float base_max_val = 2000; +// +// for cas_fitter +// +TH1 *h_bg; +Double_t slope, b; + /*----------------------------------------------------------------------------------------------------*/ @@ -90,6 +98,33 @@ Double_t fitf(Double_t *v, Double_t *par) return fitval; } + +/*----------------------------------------------------------------------------------------------------*/ + +// +// Gaussian function + 2nd order polynomial for simple rate fitting +// +// define a function with 3 parameters +Double_t fitMC(Double_t *ibin, Double_t *par) +{ + int bin = (int)round(((par[0]*ibin[0]-par[1])-b)/slope); + // get MC value at the appropriate bin (y) + Double_t bg = h_bg->GetBinContent(bin); + + // scale background appropriately + Double_t bgval = par[2]*(bg-par[3]); + + // get energy value for Gaussian + Double_t E = (Double_t)h_bg->GetBinCenter(bin); + Double_t arg = (E-par[5])/par[6]; + Double_t gaus = par[4]/(sqrt(2*TMath::Pi())*par[6])*TMath::Exp(-0.5*arg*arg); + + // BG + gauss + Double_t fitval = bgval+gaus; + + return fitval; +} + /*----------------------------------------------------------------------------------------------------*/ void analyzer::fit_spectrum(int ichannel, double *fit_range){ // @@ -130,15 +165,15 @@ void analyzer::fit_spectrum(int ichannel, double *fit_range){ // string mc_file=""; if (id == TI44){ - mc_file = "MC_ti44_modulation.root"; + mc_file = analysis_path + "MC_ti44_modulation.root"; } else if(id == CO60){ - mc_file = "MC_co60_modulation.root"; + mc_file = analysis_path + "MC_co60_modulation.root"; } else if(id == CS137){ - mc_file = "MC_cs137_modulation.root"; + mc_file = analysis_path + "MC_cs137_modulation.root"; } else if(id == MN54){ - mc_file = "MC_mn54_modulation.root"; + mc_file = analysis_path + "MC_mn54_modulation.root"; } else if(id == K40){ - mc_file = "MC_k40_modulation.root"; + mc_file = analysis_path + "MC_k40_modulation.root"; } else { cout <<"fit_spectrum:: BAD source identifier"<GetXaxis()->SetRangeUser(start,stop); + + TSpectrum *s = new TSpectrum(1); + s->Search(_pk_tmp[ipeak], 5, "new"); + Double_t *peakpos = s->GetPositionX(); + if (ipeak != 0) e_start = peakpos[ipeak]*source_energy[id][ipeak] / source_energy[id][0]; + else e_start = peakpos[0]; + + maxbin = _pk_tmp[ichannel]->GetMaximumBin(); + maxval = _pk_tmp[ichannel]->GetBinContent(maxbin); + e_start = _pk_tmp[ichannel]->GetBinCenter(maxbin); + + _pk_tmp[ichannel]->GetXaxis()->SetRangeUser(0.,3000.); + + // + // the background template for each of the sources obtained from a GEANT4 simulation + // + string mc_file=""; + if (id == TI44){ + mc_file = analysis_path + "MC_ti44_modulation.root"; + } else if(id == CO60){ + mc_file = analysis_path + "MC_co60_modulation.root"; + } else if(id == CS137){ + mc_file = analysis_path + "MC_cs137_modulation.root"; + } else if(id == MN54){ + mc_file = analysis_path + "MC_mn54_modulation.root"; + } else if(id == K40){ + mc_file = analysis_path + "MC_k40_modulation.root"; + } else { + cout <<"fit_spectrum:: BAD source identifier"<Get("h2"); + // get start and stop bins for MC + + int mc_start, mc_end; + int nentries = h_bg->GetNbinsX(); + mc_start = 1; + // get start + for (int i = 1; i <= nentries; i++) + { + if (h_bg->GetBinCenter(i) >= start) { + mc_start = i; + break; + } + } + + // get end + for (int i = mc_start; i <= nentries; i++) + { + if (h_bg->GetBinCenter(i) >= stop) { + mc_end = i; + break; + } + } + slope = (stop-start)/(mc_end-mc_start); + b = stop-slope*mc_end; + + // fit a Gauss + background to a photopeak + + Double_t res_start = 0.06/2.35*sqrt(662.)*sqrt(e_start); + Double_t scale_guess = _pk_tmp[ichannel]->GetBinWidth(1)/h_bg->GetBinWidth(1); + Double_t vert_scaleguess = _pk_tmp[ichannel]->GetBinContent(100)/h_bg->GetBinContent(100); // some random bin + + _pk_tmp[ichannel]->Draw(); + TF1 *func = new TF1("fit",fitMC,start,stop,7); + + cout << "GUESSES: [0] = " << 1 << " [1] = " << 0 << " [2] = " << vert_scaleguess << " [3] = " << 0. << " [4] = " << maxval << " [5] = " << e_start << " [6] = " << res_start << endl; + + func->SetParameters(1., 0., vert_scaleguess, 0., maxval,e_start,res_start); + + func->SetParLimits(5, e_start-100, e_start+100); + + func->SetParNames("horzscalefactor","horz_offset","bgamp", "bgvertoffset", "C","mean","sigma"); + + _pk_tmp[ichannel]->Fit("fit","R Q","",start,stop); + _pk_tmp[ichannel]->Fit("fit","R Q","",start,stop); + _pk_tmp[ichannel]->Fit("fit","R","same",start,stop); + + + Double_t peak = func->GetParameter(4); + _t_energy = func->GetParameter(5); + if(ipeak ==0) e0 = _t_energy; + Double_t sigma = func->GetParameter(6); + _t_res = 0; + if(_t_energy>0) _t_res = 2.355*sigma/_t_energy ; + + Double_t bin_width = (emax-emin)/nbin[ichannel]; + Double_t dR1 = func->GetParError(4)/TIME_INTERVAL/bin_width; + _t_rate = peak / TIME_INTERVAL / bin_width; + cout <<"get_interval_data:: ich ="<GetParameter(4) << " BG AMP: " << func->GetParameter(2) << endl; + + Double_t parameters[7]; + Double_t errors[7]; + + for (int i = 0; i < 7; i++) { + parameters[i] = func->GetParameter(i); + errors[i] = func->GetParError(i); + } + + Double_t bg_int(0); + Double_t bg_err(0); + for (int i = mc_start; i < mc_end; i++) { + bg_int += func->GetParameter(2)*(h_bg->GetBinContent((int)round(func->GetParameter(1)*(i-func->GetParameter(1))) - func->GetParameter(3))); + bg_err += func->GetParameter(2)*(h_bg->GetBinError((int)round(func->GetParameter(1)*(i-func->GetParameter(1))) - func->GetParameter(3))); + } + + bg_int /= TIME_INTERVAL; + bg_int /= bin_width; + bg_err /= TIME_INTERVAL; + bg_err /= bin_width; + + // fill the output tree..... + addTreeEntry(_t_energy,_t_rate,dR1,_t_res,ichannel,ipeak, parameters, errors, bg_int, bg_err); + + TF1 *fitresult = _pk_tmp[ichannel]->GetFunction("fit"); + + if (PLOT_ON_SCREEN == 1) { + TCanvas *c1 = new TCanvas("c1","c1",900,900); + c1->cd(); + TF1* fittedgaus = new TF1("fittedgaus", "[0]/(2*TMath::Pi()*[2])*TMath::Exp((-1*(x-[1])*(x-[1]))/([2]*[2]))", fit_range[0], fit_range[1]); + fittedgaus->SetParameters(func->GetParameter(4), func->GetParameter(5), func->GetParameter(6)); + TF1* fittedbg = new TF1("fittedbg", "[2]*(h_bg->GetBinContent((int)round((([0]*x-[1]-b)/slope)))-[3])", mc_start, mc_end); + fittedbg->SetParameters(func->GetParameter(0), func->GetParameter(1), func->GetParameter(2), func->GetParameter(3)); + fittedgaus->SetLineColor(2); + fittedbg->SetLineColor(4); + fittedgaus->Draw("same"); + fittedbg->Draw("same"); + + c1->Update(); + + + PlotResiduals(_pk_tmp[ichannel], fitresult, fit_range[0], fit_range[1]); + } + + + delete func; + f_mc->Close(); + return; + +} + /*----------------------------------------------------------------------------------------------------*/ void analyzer::fit_spectrum(int ichannel){ // @@ -309,15 +507,15 @@ void analyzer::fit_spectrum(int ichannel){ // string mc_file=""; if (id == TI44){ - mc_file = "/user/z37/Modulation/analysis/calibration/MC_ti44_modulation.root"; + mc_file = analysis_path + "MC_ti44_modulation.root"; } else if(id == CO60){ - mc_file = "/user/z37/Modulation/analysis/calibration/MC_co60_modulation.root"; + mc_file = analysis_path + "MC_co60_modulation.root"; } else if(id == CS137){ - mc_file = "/user/z37/Modulation/analysis/calibration/MC_cs137_modulation.root"; + mc_file = analysis_path + "MC_cs137_modulation.root"; } else if(id == MN54){ - mc_file = "/user/z37/Modulation/analysis/calibration/MC_mn54_modulation.root"; + mc_file = analysis_path + "MC_mn54_modulation.root"; } else if(id == K40){ - mc_file = "/user/z37/Modulation/analysis/calibration/MC_k40_modulation.root"; + mc_file = analysis_path + "MC_k40_modulation.root"; } else { cout <<"fit_spectrum:: BAD source identifier"<Fill(); } @@ -543,6 +766,35 @@ void analyzer::book_histograms(){ tree->Branch("bz", &_t_bz, "bz/D"); tree->Branch("btot", &_t_btot, "btot/D"); tree->Branch("humid", &_t_humid, "humid/D"); + + tree->Branch("hv0", &_t_hv0, "hv0/D"); + tree->Branch("hv1", &_t_hv1, "hv1/D"); + tree->Branch("hv2", &_t_hv2, "hv2/D"); + tree->Branch("hv3", &_t_hv3, "hv3/D"); + tree->Branch("hv4", &_t_hv4, "hv4/D"); + tree->Branch("hv5", &_t_hv5, "hv5/D"); + tree->Branch("hv6", &_t_hv6, "hv6/D"); + tree->Branch("hv7", &_t_hv7, "hv7/D"); + + tree->Branch("bg_int", &_t_bg_int, "bg_int/D"); + tree->Branch("bg_err", &_t_bg_err, "bg_err/D"); + + tree->Branch("par0",&_t_par0, "par0/D"); + tree->Branch("par1",&_t_par1, "par1/D"); + tree->Branch("par2",&_t_par2, "par2/D"); + tree->Branch("par3",&_t_par3, "par3/D"); + tree->Branch("par4",&_t_par4, "par4/D"); + tree->Branch("par5",&_t_par5, "par5/D"); + tree->Branch("par6",&_t_par6, "par6/D"); + + tree->Branch("err0",&_t_err0, "err0/D"); + tree->Branch("err1",&_t_err1, "err1/D"); + tree->Branch("err2",&_t_err2, "err2/D"); + tree->Branch("err3",&_t_err3, "err3/D"); + tree->Branch("err4",&_t_err4, "err4/D"); + tree->Branch("err5",&_t_err5, "err5/D"); + tree->Branch("err6",&_t_err6, "err6/D"); + } /*----------------------------------------------------------------------------------------------------*/ void analyzer::fill_histograms(){ @@ -607,25 +859,33 @@ void analyzer::get_interval_data(){ int id = source_id[ich]; if (id == TI44) { range[0] = 400; range[1] = 620; - fit_spectrum(ich, range); + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 0, range); range[0] = 1000; range[1] = 1300; - fit_spectrum(ich, range); + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 1, range); range[0] = 1500; range[1] = 2000; - fit_spectrum(ich, range); + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 2, range); } else if (id == CO60 ) { range[0] = 900; range[1] = 1600; - fit_spectrum(ich, range); + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 0, range); range[0] = 2200; range[1] = 2800; - fit_spectrum(ich, range); + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 1, range); } else if (id == CS137 ) { range[0] = 400; range[1] = 1000; - fit_spectrum(ich, range); + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 0, range);; } else if (id == MN54) { - range[0] = 0; range[1] = 2000; - fit_spectrum(ich, range); + range[0] = 700; range[1] = 900; + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 0, range); } else if (id == K40) { range[0] = 1300; range[1] = 1600; - fit_spectrum(ich, range); + if (ROOFIT_ON == 1) fit_spectrum(ich, range); + else cas_fitter(ich, 0, range); } else { // // if we deal with the background detectors we use a linear fit + gauss to fit the signal @@ -719,8 +979,16 @@ void analyzer::fit_spectrum_simple(int ichannel){ cout <<"get_interval_data:: ich ="<GetParameter(i); + errors[i] = func->GetParError(i); + } + // fille the output tree..... - addTreeEntry(_t_energy,_t_rate,1.0,_t_res,ichannel,ipeak); + addTreeEntry(_t_energy,_t_rate,1.0,_t_res,ichannel,ipeak, parameters, errors, 0.0, 0.0); // c1->Update(); delete func; @@ -783,6 +1051,15 @@ void analyzer::reset_interval_data(){ _t_btot = 0; _t_humid = 0; + _t_hv0 = 0; + _t_hv1 = 0; + _t_hv2 = 0; + _t_hv3 = 0; + _t_hv4 = 0; + _t_hv5 = 0; + _t_hv6 = 0; + _t_hv7 = 0; + } /*----------------------------------------------------------------------------------------------------*/ void analyzer::add_interval_data(){ @@ -796,18 +1073,34 @@ void analyzer::add_interval_data(){ _t_bz += bz; _t_btot += btot; _t_humid += humid; + _t_hv0 += hv0; + _t_hv1 += hv1; + _t_hv2 += hv2; + _t_hv3 += hv3; + _t_hv4 += hv4; + _t_hv5 += hv5; + _t_hv6 += hv6; + _t_hv7 += hv7; } /*----------------------------------------------------------------------------------------------------*/ void analyzer::calculate_interval_data(){ if(n_interval>0){ - _t_temp /= n_interval; - _t_pres /= n_interval; - _t_bx /= n_interval; - _t_by /= n_interval; - _t_bz /= n_interval; - _t_btot /= n_interval; - _t_humid /= n_interval; + _t_temp += temp; + _t_pres += pres; + _t_bx += bx; + _t_by += by; + _t_bz += bz; + _t_btot += btot; + _t_humid += humid; + _t_hv0 += hv0; + _t_hv1 += hv1; + _t_hv2 += hv2; + _t_hv3 += hv3; + _t_hv4 += hv4; + _t_hv5 += hv5; + _t_hv6 += hv6; + _t_hv7 += hv7; } } /*---------------------------------------------------------------------------------------------------*/ @@ -870,6 +1163,64 @@ void analyzer::get_source_id() cout <<"analyzer::get_source_id ... done"<GetNbinsX(); + + cout << "PLOT RESIDUALS" << endl; + Float_t resid_i(0); + std::vector residuals; + std::vector energies; + std::vector eerrors; + std::vector rerrors; + + int mc_start, mc_end; + mc_start = 1; + // get start + for (int i = 1; i <= nbins; i++) + { + if (inthist->GetBinCenter(i) >= Elow) { + mc_start = i; + break; + } + } + + // get end + for (int i = mc_start; i <= nbins; i++) + { + if (inthist->GetBinCenter(i) >= Ehigh) { + mc_end = i; + break; + } + } + + for (int i = mc_start; i <= mc_end; i++) { + if (inthist->GetBinCenter(i) >= Elow && inthist->GetBinCenter(i) < Ehigh) { + cout << "Bin content: " << inthist->GetBinContent(i) << endl; + cout << "Bin center: " << inthist->GetBinCenter(i) << endl; + cout << "Eval: " << fitresult->Eval(inthist->GetBinCenter(i)); + resid_i = (inthist->GetBinContent(i) - fitresult->Eval(inthist->GetBinCenter(i)))/inthist->GetBinError(i); + cout <<"i: " << i << " measured: " << inthist->GetBinContent(i) << " expected: " << fitresult->Eval(inthist->GetBinCenter(i)) << " diff : " << inthist->GetBinContent(i) - fitresult->Eval(inthist->GetBinCenter(i)) << " sigma: " << inthist->GetBinError(i) << " residual: " << resid_i << endl; + residuals.push_back(resid_i); + energies.push_back(inthist->GetBinCenter(i)); + eerrors.push_back(inthist->GetBinError(i)); + rerrors.push_back(inthist->GetBinError(i)); + } + } + + TCanvas *c2 = new TCanvas("c2", "c2"); + c2->cd(); + // draw residuals here + + TGraph *residualplot = new TGraph(energies.size(), energies.data(), residuals.data()); + residualplot->GetXaxis()->SetTitle("Energy "); + residualplot->GetYaxis()->SetTitle("Residuals"); + residualplot->Draw("AP"); + + return; +} + /*----------------------------------------------------------------------------------------------------*/ // // MAIN:: Loop routine diff --git a/calibration/analyzer.h b/calibration/analyzer.h index 1026776..6b1ec0f 100644 --- a/calibration/analyzer.h +++ b/calibration/analyzer.h @@ -14,8 +14,9 @@ /*----------------------------------------------------------------------------*/ -#define TIME_INTERVAL 3600 +#define TIME_INTERVAL 300 #define PLOT_ON_SCREEN 0 +#define ROOFIT_ON 0 /*----------------------------------------------------------------------------*/ #define NUMBER_OF_CHANNELS 8 @@ -62,6 +63,14 @@ class analyzer { Double_t bz; Double_t btot; Double_t humid; + Double_t hv0; + Double_t hv1; + Double_t hv2; + Double_t hv3; + Double_t hv4; + Double_t hv5; + Double_t hv6; + Double_t hv7; @@ -82,6 +91,14 @@ class analyzer { TBranch *b_bz; //! TBranch *b_btot; //! TBranch *b_humid; //! + TBranch *b_hv0; //! + TBranch *b_hv1; //! + TBranch *b_hv2; //! + TBranch *b_hv3; //! + TBranch *b_hv4; //! + TBranch *b_hv5; //! + TBranch *b_hv6; //! + TBranch *b_hv7; //! analyzer(string fname, string cname); virtual ~analyzer(); @@ -104,9 +121,11 @@ class analyzer { void fit_spectrum(int ichannel); void fit_spectrum_simple(int ichannel); void fit_spectrum(int ichannel, double *fit_range); - void addTreeEntry(Double_t E, Double_t R, Double_t dR, Double_t res, Int_t ich, Int_t ipk); + void addTreeEntry(Double_t E, Double_t R, Double_t dR, Double_t res, Int_t ich, Int_t ipk, Double_t* parameters, Double_t* errors, Double_t bg_int, Double_t bg_err); void processFitData(RooRealVar N, RooRealVar f, RooRealVar E, RooRealVar sig, int ichannel, int ipeak); void get_source_id(); + void cas_fitter(int ichannel, int ipeak, double *fit_range); + void PlotResiduals(TH1F *inthist, TF1* fitresult, Float_t Elow, Float_t Ehigh); double covariance(int i, int j); // histograms @@ -135,6 +154,36 @@ class analyzer { Double_t _t_btot; Double_t _t_humid; + Double_t _t_hv0; + Double_t _t_hv1; + Double_t _t_hv2; + Double_t _t_hv3; + Double_t _t_hv4; + Double_t _t_hv5; + Double_t _t_hv6; + Double_t _t_hv7; + + Double_t _t_par0; + Double_t _t_par1; + Double_t _t_par2; + Double_t _t_par3; + Double_t _t_par4; + Double_t _t_par5; + Double_t _t_par6; + Double_t _t_par7; + + Double_t _t_err0; + Double_t _t_err1; + Double_t _t_err2; + Double_t _t_err3; + Double_t _t_err4; + Double_t _t_err5; + Double_t _t_err6; + Double_t _t_err7; + + Double_t _t_bg_int; + Double_t _t_bg_err; + // time information Double_t t0,tstart,time_since_start,delta_t; From 0fdf4b0bbd97dddca8635c7b9c5097c128b9d5da Mon Sep 17 00:00:00 2001 From: Cassie Date: Mon, 25 Apr 2016 13:12:16 -0400 Subject: [PATCH 2/4] Added a fitting routine to analyzer.C, can be toggled with the old version by changing ROOFIT_ON to 1. Also added hv channels, fit parameters, errors on the fit parameters and the rate of the background as branches in the analysis tree --- calibration/analyzer.C | 53 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 10 deletions(-) diff --git a/calibration/analyzer.C b/calibration/analyzer.C index 9b5c726..0bcd889 100644 --- a/calibration/analyzer.C +++ b/calibration/analyzer.C @@ -125,6 +125,19 @@ Double_t fitMC(Double_t *ibin, Double_t *par) return fitval; } +// just the BG +Double_t BGmodel(Double_t *x, Double_t *par) +{ + int bin = (int)round(((par[0]*(x[0]-par[1]))-b)/slope); + // get MC value at the appropriate bin (y) + Double_t bg = h_bg->GetBinContent(bin); + + // scale background appropriately + Double_t bgval = par[2]*(bg-par[3]); + + return bgval; +} + /*----------------------------------------------------------------------------------------------------*/ void analyzer::fit_spectrum(int ichannel, double *fit_range){ // @@ -405,8 +418,7 @@ void analyzer::cas_fitter(int ichannel, int ipeak, double *fit_range){ Double_t res_start = 0.06/2.35*sqrt(662.)*sqrt(e_start); Double_t scale_guess = _pk_tmp[ichannel]->GetBinWidth(1)/h_bg->GetBinWidth(1); Double_t vert_scaleguess = _pk_tmp[ichannel]->GetBinContent(100)/h_bg->GetBinContent(100); // some random bin - - _pk_tmp[ichannel]->Draw(); + TF1 *func = new TF1("fit",fitMC,start,stop,7); cout << "GUESSES: [0] = " << 1 << " [1] = " << 0 << " [2] = " << vert_scaleguess << " [3] = " << 0. << " [4] = " << maxval << " [5] = " << e_start << " [6] = " << res_start << endl; @@ -416,10 +428,10 @@ void analyzer::cas_fitter(int ichannel, int ipeak, double *fit_range){ func->SetParLimits(5, e_start-100, e_start+100); func->SetParNames("horzscalefactor","horz_offset","bgamp", "bgvertoffset", "C","mean","sigma"); - + _pk_tmp[ichannel]->Fit("fit","R Q","",start,stop); _pk_tmp[ichannel]->Fit("fit","R Q","",start,stop); - _pk_tmp[ichannel]->Fit("fit","R","same",start,stop); + _pk_tmp[ichannel]->Fit("fit","R","",start,stop); Double_t peak = func->GetParameter(4); @@ -462,21 +474,39 @@ void analyzer::cas_fitter(int ichannel, int ipeak, double *fit_range){ TF1 *fitresult = _pk_tmp[ichannel]->GetFunction("fit"); if (PLOT_ON_SCREEN == 1) { - TCanvas *c1 = new TCanvas("c1","c1",900,900); + + TCanvas *c1 = new TCanvas("c1", "c1", 600, 300); c1->cd(); - TF1* fittedgaus = new TF1("fittedgaus", "[0]/(2*TMath::Pi()*[2])*TMath::Exp((-1*(x-[1])*(x-[1]))/([2]*[2]))", fit_range[0], fit_range[1]); + c1->SetLogy(); + _pk_tmp[ichannel]->GetXaxis()->SetRangeUser(fit_range[0]-100,fit_range[1]+100); + _pk_tmp[ichannel]->SetTitle("energy spectrum"); + cout << "ENTRIES: " << _pk_tmp[ichannel]->GetEntries() << endl; + _pk_tmp[ichannel]->Draw(); + + TF1* fittedgaus = new TF1("fittedgaus", "[0]/(2*TMath::Pi()*[2])*TMath::Exp((-1*(x-[1])*(x-[1]))/([2]*[2]))", fit_range[0]-100, fit_range[1]+100); fittedgaus->SetParameters(func->GetParameter(4), func->GetParameter(5), func->GetParameter(6)); - TF1* fittedbg = new TF1("fittedbg", "[2]*(h_bg->GetBinContent((int)round((([0]*x-[1]-b)/slope)))-[3])", mc_start, mc_end); + TF1* fittedbg = new TF1("fittedbg", BGmodel, fit_range[0]-100, fit_range[1]+100, 4); fittedbg->SetParameters(func->GetParameter(0), func->GetParameter(1), func->GetParameter(2), func->GetParameter(3)); fittedgaus->SetLineColor(2); fittedbg->SetLineColor(4); + fitresult->Draw("same"); fittedgaus->Draw("same"); fittedbg->Draw("same"); - + fitresult->SetRange(fit_range[0]-100, fit_range[1]+100); + fitresult->Draw("same"); + + /* + char pdfname[256]; + sprintf(pdfname,"fitted_mn54.png"); + c1->Print(pdfname); + sprintf(pdfname,"fitted_mn54.pdf"); + c1->Print(pdfname); + */ c1->Update(); - + PlotResiduals(_pk_tmp[ichannel], fitresult, fit_range[0], fit_range[1]); + } @@ -893,7 +923,7 @@ void analyzer::get_interval_data(){ fit_spectrum_simple(ich); } } - _pk_tmp[ich]->Reset(); // reset the histogram + //_pk_tmp[ich]->Reset(); // reset the histogram // TODO } // loop over channels } /*----------------------------------------------------------------------------------------------------*/ @@ -1248,6 +1278,8 @@ void analyzer::Loop() Long64_t nbytes = 0, nb = 0; Bool_t last_event = false; + for (int ichannel = 0; ichannel < NUMBER_OF_CHANNELS; ichannel++) _pk_tmp[ichannel]->Sumw2(); + for (Long64_t jentry=0; jentry Date: Thu, 23 Feb 2017 17:04:31 -0500 Subject: [PATCH 3/4] Update ecal.C --- calibration/ecal.C | 92 +++++++++++++++++++++++++++++++--------------- 1 file changed, 63 insertions(+), 29 deletions(-) diff --git a/calibration/ecal.C b/calibration/ecal.C index 7c2bcf8..662f88d 100644 --- a/calibration/ecal.C +++ b/calibration/ecal.C @@ -10,6 +10,8 @@ #include #include #include +#include +#include /*---------------------------------------------------------------------------------------------------*/ // @@ -29,8 +31,8 @@ // (2) Set range_low and range_high to find the range in which the peak with cal_energy should be // CH0 CH1 CH2 CH3 CH4 CH5 CH6 CH7 -const double range_low[NUMBER_OF_CHANNELS] ={0.16e-6 , 0.14e-6 , 0.05e-6 , 0.05e-6 , 0. , 0. , 0. , 0. }; -const double range_high[NUMBER_OF_CHANNELS]={1e-6 , 1e-6 , 1e-6 , 1e-6 , 1e-6 , 1e-6 , 1e-6 , 1e6 }; +const double range_low[NUMBER_OF_CHANNELS] ={0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. }; +const double range_high[NUMBER_OF_CHANNELS]={.1e-6 , .1e-6 , 1e-6 , 1e-6 , 2e-6 , 0.2e-6 , 0.4e-6 , 0.4e-6 }; /*---------------------------------------------------------------------------------------------------*/ float source_energy[NUMBER_OF_SOURCES][MAX_PEAKS] = @@ -39,10 +41,10 @@ float source_energy[NUMBER_OF_SOURCES][MAX_PEAKS] = // { {1460.,-1,-1,-1,-1}, // ID0: Background - {511.,1157.020,511.+1157.020,-1,-1}, // ID1: Ti44 + {511.,1157.020,(1157.020+511.), -1, -1}, // ID1: Ti44 {1173.2,1332.5,1173.2+1332.5,-1,-1}, // ID2: Co60 {661.7,-1,-1,-1,-1}, // ID3: CS137 - {-1,-1,-1,-1,-1}, // ID4: MN54 + {834.8,-1,-1,-1,-1}, // ID4: MN54 {1460.,-1,-1,-1,-1} // ID5: K40 }; @@ -110,7 +112,7 @@ void ecal::book_histograms(){ /*---------------------------------------------------------------------------------------------------*/ void ecal::fill_histograms(int ilevel){ - Long64_t nentries = fChain->GetEntriesFast(); + Long64_t nentries = fChain->GetEntries(); cout<<"Start calibration loop.... nentries ="<GetEntriesFast(); - cout<<"Start calibration loop.... nentries ="<GetEntries(); + cout<<"Start calibration loop.... nentries ="<GetEntriesFast() <GetEntry(jentry); nbytes += nb; // @@ -261,44 +267,53 @@ void ecal::ecal_continuous(){ // // fill the histogram // - if(error == 0 && !istestpulse) _integral[channel]->Fill(integral); + if (!istestpulse) _integral[channel]->Fill(integral); - if(time_since_start - t0 > TIME_INTERVAL) { + if(time_since_start - t0 > TIME_INTERVAL || last_event) { // maximum validity time of this calibration.... _cal_tmax = time; // // do the calibration // + do_calibration(); + // // file the output tree // + fill_tree(_cal_tmin,_cal_tmax); + // // reset the histograms // + reset_histograms(); + // // reset the time for the start of the next interval // t0 = time_since_start; // SET THE MINIMUM TIME OF THE NEXT INTERVAL _cal_tmin = time; + } - if(jentry%100000 == 0) cout<<"Processed "<Write(); _f->Close(); + } /*---------------------------------------------------------------------------------------------------*/ void ecal::reset_histograms(){ @@ -338,6 +353,7 @@ void ecal::ecal_single(){ sprintf(parname,"cal_ch%02d_c%i",ich,ipar); TParameter *p1 = new TParameter(parname,ccal[ich][ipar]); p1->Write(); + delete p1; } } // @@ -350,7 +366,6 @@ void ecal::ecal_single(){ // fill_histograms(AFTER_CALIBRATION); - _f->Write(); _f->Close(); @@ -410,22 +425,34 @@ void ecal::do_calibration(){ // the source identifier for this channel int id = source_id[ich]; + for (int i=0; i0) npeak++; } - /// if(npeak == 0){ // no calibration possible ..... E = 1.0*integral. - // if the calibration fails I want a simple linear realtion between integral and energy - /// } - // STEP1: find the highest peak // define the proper range for searching the calibration constant _integral[ich]->GetXaxis()->SetRangeUser(range_low[ich],range_high[ich]); - // find the bin with maximum value in the range - int ibin = _integral[ich]->GetMaximumBin(); - double val = _integral[ich]->GetBinCenter(ibin); + + // Cassie changed: step 1 + // need the primary peak, not necessarily the tallest (e.g. Ti-44, where the tallest peak is 68-78 keV, but the primary peak is 511 keV + TSpectrum *s = new TSpectrum(MAX_PEAKS); // 5 peaks just in case... we're really only looking for the biggest one + s->Search(_integral[ich], 2.5, "new"); + double *peakpos = s->GetPositionX(); + double val; + int ibin = 0; + val = peakpos[0]; + + int nentries = _integral[ich]->GetNbinsX(); + for (int i = 0; i = _integral[ich]->GetBinCenter(i)){ + ibin = i; + break; + } + } double maxval = _integral[ich]->GetBinContent(ibin); + delete s; Double_t ee[MAX_PEAKS]; Double_t dee[MAX_PEAKS]; @@ -450,13 +477,14 @@ void ecal::do_calibration(){ vlow = val - 5*sig_expected;//0.025e-6; vhigh = val + 5*sig_expected;//0.025e-6; + /* if(ich == 4 || ich ==5){ if(ipeak == 0 ) { vhigh = val+3*sig_expected;//0.015e-6; } if(ipeak == 1 ) vlow = val-3*sig_expected;//0.015e-6; } - + */ // // fit a Gauss around the desired location in the spectrum // @@ -473,12 +501,6 @@ void ecal::do_calibration(){ func->SetParameters(maxval,val,sig_expected); func->SetParNames("C","mean","sigma"); _integral[ich]->Fit("fit","Q","",vlow,vhigh); - // if(ich==2 || ich==4){ - // cout << ich<<" " < Date: Thu, 23 Feb 2017 17:05:52 -0500 Subject: [PATCH 4/4] Update ecal.h --- calibration/ecal.h | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/calibration/ecal.h b/calibration/ecal.h index 410088b..8356319 100644 --- a/calibration/ecal.h +++ b/calibration/ecal.h @@ -33,7 +33,7 @@ // CALIBRATION_MODE 1 = one calibration every TIME_INTERVAL seconds // #define CALIBRATION_MODE 1 -#define TIME_INTERVAL 1800 +#define TIME_INTERVAL 600 /*----------------------------------------------------------------------------*/ @@ -41,6 +41,8 @@ // Fixed size dimensions of array or collections stored in the TTree if any. +using namespace std; + class ecal { public : string calFile; @@ -53,6 +55,7 @@ public : Float_t height; Double_t time; UChar_t istestpulse; + UChar_t istrigger; Int_t error; Float_t baseline; Float_t rms; @@ -63,6 +66,7 @@ public : TBranch *b_height; //! TBranch *b_time; //! TBranch *b_istestpulse; //! + TBranch *b_istrigger; //! TBranch *b_error; //! TBranch *b_baseline; //! TBranch *b_baselineRMS; //! @@ -128,6 +132,7 @@ ecal::ecal(string fname,string cname) : fChain(0) char cmd[256]; TChain *tree = new TChain("T"); sprintf(cmd,"%s*.root",fname.c_str()); + cout << cmd << endl; tree->Add(cmd); calFile = cname; @@ -178,6 +183,7 @@ void ecal::Init(TChain *tree) // Set branch addresses and branch pointers if (!tree) return; + fChain = tree; fCurrent = -1; fChain->SetMakeClass(1); @@ -190,6 +196,7 @@ void ecal::Init(TChain *tree) fChain->SetBranchAddress("error", &error, &b_error); fChain->SetBranchAddress("baseline", &baseline, &b_baseline); fChain->SetBranchAddress("rms", &rms, &b_baselineRMS); + //fChain->SetBranchAddress("istrigger", &istrigger, &b_istrigger); Notify(); }