diff --git a/calibration/analyzer.C b/calibration/analyzer.C index 648235b..0bcd889 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,46 @@ 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; +} + +// 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){ // @@ -130,15 +178,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 + + 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","",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", 600, 300); + c1->cd(); + 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", 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]); + + } + + + delete func; + f_mc->Close(); + return; + +} + /*----------------------------------------------------------------------------------------------------*/ void analyzer::fit_spectrum(int ichannel){ // @@ -309,15 +537,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 +796,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 +889,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 @@ -633,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 } /*----------------------------------------------------------------------------------------------------*/ @@ -719,8 +1009,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 +1081,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 +1103,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 +1193,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 @@ -897,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 #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<<" " <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(); }