offline/script/function_fitting.C

0001 #include <iostream>
0002 #include <string>
0003 #include <fstream>
0004 #include <TAxis.h>
0005 #include <TGraph.h>
0006 #include <TSpectrum.h>
0007 #include <TH1F.h>
0008 #include <TF1.h>
0009 #include <TMath.h>
0010 #include <TROOT.h>
0011 #include <TFile.h>
0012 #include "TError.h"
0013
0014 double gamma_function(double *x, double *par);
0015 float get_temp_shift(TH1F* hist, float y);
0016 float get_temp_scale(TH1F* hist, float x);
0017 float get_temp_N(float x, float y, float shift, float scale);
0018 float get_fit_min(TH1F* hist, float y);
0019 float get_fit_max(TH1F* hist, float y);
0020 void do_peak_finder(TH1F* hist, float& x, float& y, float& shift, float& scale, float& N, float& min, float& max);
0021 void do_peak_fit(TH1F* hist, TF1* func);
0022
0023 void function_fitting(int min_run, int max_run, std::string calo_type) {
0024   static const int n_etabin = 24;
0025   static const int n_phibin = 64;
0026
0027   // Output file for each run.

0028   std::string input_filename, output_filename, check_filename, mpv_filename;
0029   if (calo_type == "ihcal") {
0030     input_filename = "hist/ihcal_hist_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".root";
0031     output_filename = "output/fitresult_ihcal_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".root";
0032     check_filename = "output/fitcheck_ihcal_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".root";
0033     mpv_filename = "output/ihcal_mpv_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".txt";
0034   } else if (calo_type == "ohcal") {
0035     input_filename = "hist/ohcal_hist_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".root";
0036     output_filename = "output/fitresult_ohcal_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".root";
0037     check_filename = "output/fitcheck_ohcal_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".root";
0038     mpv_filename = "output/ohcal_mpv_" + std::to_string(min_run) + "_" + std::to_string(max_run) + ".txt";
0039   } else {
0040     std::cout << "Invalid calo_type specified. Use 'ihcal' or 'ohcal'." << std::endl;
0041     return;
0042   }
0043
0044   TFile* hist_input = new TFile(input_filename.c_str(), "READ");
0045   TFile* output_file = new TFile(output_filename.c_str(), "RECREATE");
0046   TFile* check_file = new TFile(check_filename.c_str(), "RECREATE");
0047   ofstream mpv_file(mpv_filename.c_str(), std::ios::trunc);
0048
0049   if (!hist_input || hist_input->IsZombie()) {
0050     std::cout << "Input file not found: " << input_filename << std::endl;
0051     return;
0052   }
0053
0054   TH1F* h_event = dynamic_cast<TH1F*>(hist_input->Get("h_event"));
0055   TH1F* h_channel[n_etabin][n_phibin];
0056   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0057     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0058       std::string hist_name = Form("h_channel_%d_%d", ieta, iphi);
0059       h_channel[ieta][iphi] = (TH1F*)hist_input->Get(hist_name.c_str());
0060     }
0061   }
0062
0063   // Setup check hist.

0064   TH1F* h_mip; h_mip = (TH1F*)hist_input->Get("h_mip");
0065   TH2F* h_2D_hit = new TH2F("h_2D_hit", ";ieta;iphi", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0066   TH2F* h_2D_mpv = new TH2F("h_2D_mpv", ";ieta;iphi", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0067   TH2F* h_2D_mean = new TH2F("h_2D_mean", ";ieta;iphi", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0068   TH2F* h_2D_mpvuncertainty = new TH2F("h_2D_mpvuncertainty", ";ieta;iphi", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0069   TH2F* h_2D_meanuncertainty = new TH2F("h_2D_meanuncertainty", ";ieta;iphi", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0070   TH2F* h_2D_fitchi2 = new TH2F("h_2D_fitchi2", ";ieta;iphi", n_etabin, 0, n_etabin, n_phibin, 0, n_phibin);
0071
0072   // Fit channel hist and fill check hist & check file.

0073   TF1 *f_gamma = new TF1("f_gamma", gamma_function, 0, 10000, 4);
0074   f_gamma->SetParName(0,"Peak(ADC)");
0075   f_gamma->SetParName(1,"Shift");
0076   f_gamma->SetParName(2,"Scale");
0077   f_gamma->SetParName(3,"N");
0078   do_peak_fit(h_mip, f_gamma);
0079   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0080     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0081       //std::cout << "Fitting ieta: " << ieta << " iphi: " << iphi << std::endl;

0082       do_peak_fit(h_channel[ieta][iphi], f_gamma);
0083       h_2D_hit->SetBinContent(ieta+1, iphi+1, h_channel[ieta][iphi]->GetEntries());
0084       h_2D_mpv->SetBinContent(ieta+1, iphi+1, f_gamma->GetParameter(0));
0085       h_2D_mean->SetBinContent(ieta+1, iphi+1, h_channel[ieta][iphi]->GetMean());
0086       h_2D_mpvuncertainty->SetBinContent(ieta+1, iphi+1, f_gamma->GetParError(0)/(float)f_gamma->GetParameter(0));
0087       h_2D_meanuncertainty->SetBinContent(ieta+1, iphi+1, h_channel[ieta][iphi]->GetMeanError()/(float)h_channel[ieta][iphi]->GetMean());
0088       h_2D_fitchi2->SetBinContent(ieta+1, iphi+1, f_gamma->GetChisquare()/(float)f_gamma->GetNDF());
0089
0090       if (h_channel[ieta][iphi]->GetEntries() == 0) {
0091         mpv_file << ieta << " " << iphi << " " << 0 << std::endl;
0092       } else {
0093         mpv_file << ieta << " " << iphi << " " << f_gamma->GetParameter(0) << std::endl;
0094       }
0095     }
0096   }
0097
0098   // Save hist.

0099   output_file->cd();
0100   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0101     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0102       h_channel[ieta][iphi]->Write();
0103     }
0104   }
0105   check_file->cd();
0106   h_mip->Write();
0107   h_2D_hit->Write();
0108   h_2D_mpv->Write();
0109   h_2D_mean->Write();
0110   h_2D_mpvuncertainty->Write();
0111   h_2D_meanuncertainty->Write();
0112   h_2D_fitchi2->Write();
0113   h_event->Write();
0114 }
0115
0116 double gamma_function(double *x, double *par) {
0117   double peak = par[0];
0118   double shift = par[1];
0119   double scale = par[2];
0120   double N = par[3];
0121
0122   if (scale==0) return 0;
0123   double arg_para = (x[0] - shift) / scale;
0124   if (arg_para < 0) arg_para=0;
0125   double peak_para = (peak - shift) / scale;
0126   double numerator = N * pow(arg_para, peak_para) * TMath::Exp(-arg_para);
0127   double denominator = ROOT::Math::tgamma(peak_para + 1) * scale;
0128   if (denominator==0) return 1e8;
0129   double val = numerator / denominator;
0130   if (isnan(val)) return 0;
0131   return val;
0132 }
0133
0134 float get_temp_shift(TH1F* hist, float y) {
0135   float shift = 0.;
0136   for (int i = 1; i <= hist->GetNbinsX(); i++) {
0137     if (hist->GetBinContent(i) > y/10.) {
0138       shift = hist->GetBinCenter(i);
0139       break;
0140     }
0141   }
0142   return shift;
0143 }
0144
0145 float get_temp_scale(TH1F* hist, float x) {
0146   float scale = (hist->GetMean() - x);
0147   return scale;
0148 }
0149
0150 float get_temp_N(float x, float y, float shift, float scale) {
0151   float temp_par0 = (x - shift) / scale;
0152   float temp_par1 = pow(temp_par0, temp_par0) * TMath::Exp(-temp_par0);
0153   float temp_par2 = ROOT::Math::tgamma(temp_par0 + 1) * scale * y;
0154   float temp_N = temp_par2 / temp_par1;
0155   return temp_N;
0156 }
0157
0158 float get_fit_min(TH1F* hist, float y) {
0159   float min = 0.;
0160   for (int i = 1; i <= hist->GetNbinsX(); i++) {
0161     if (hist->GetBinContent(i) > 2*y/3.) {
0162       min = hist->GetBinCenter(i-2);
0163       if (hist->GetBinContent(i) > 3*y/4.) {
0164         min = hist->GetBinCenter(i-3);
0165       }
0166       break;
0167     }
0168   }
0169   return min;
0170 }
0171
0172 float get_fit_max(TH1F* hist, float y) {
0173   float max = 0.;
0174   for (int i = hist->GetNbinsX(); i >= 1; i--) {
0175     if (hist->GetBinContent(i) > y/5.) {
0176       max = hist->GetBinCenter(i-3);
0177       break;
0178     }
0179   }
0180   return max;
0181 }
0182
0183 void do_peak_finder(TH1F* hist, float& x, float& y, float& shift, float& scale, float& N, float& min, float& max) {
0184   TH1F *h_peakfinder = (TH1F*)hist->Clone("h_peakfinder");
0185   h_peakfinder->Smooth(2);
0186   TSpectrum *s = new TSpectrum(10);
0187   int nfound = 0;
0188   nfound = s->Search(h_peakfinder, 2, "nobackground new", 0.6);
0189   if (nfound == 0) {
0190     int index = hist->GetMaximumBin();
0191     x = hist->GetBinCenter(index);
0192     y = hist->GetBinContent(index);
0193   } else {
0194     double *xtemp = s->GetPositionX();
0195     double *ytemp = s->GetPositionY();
0196     x = xtemp[0];
0197     y = ytemp[0];
0198   }
0199   shift = get_temp_shift(h_peakfinder, y);
0200   scale = get_temp_scale(hist, x);
0201   N = get_temp_N(x, y, shift, scale);
0202   min = get_fit_min(h_peakfinder, y);
0203   max = get_fit_max(h_peakfinder, y);
0204   delete h_peakfinder;
0205   delete s;
0206 }
0207
0208 void set_parameters(TF1* func, float peak_position, float shift, float scale, float N) {
0209   func->SetParameter(0, peak_position);
0210   func->SetParameter(1, shift);
0211   func->SetParameter(2, scale);
0212   func->SetParameter(3, N);
0213 }
0214
0215 void do_peak_fit(TH1F* hist, TF1* func) {
0216   int n_rebin = 5;
0217   float bin_width = hist->GetBinWidth(1);
0218
0219   // *** First fitting try ***//

0220   int min_nloop = 9; // loop through min ?bins

0221   int min_lowloop = 4; // loop from finder_min - ?bins

0222   int max_nloop = 13; // loop through max ?bins

0223   int max_lowloop = 3; // loop from finder_max - ?bins

0224
0225   int min_lowlimit = 2; // min limit: shift + ?bins

0226   int min_highlimit = 4; // min limit: peak - ?bins

0227   int max_lowlimit = 4; // max limit: peak + ?bins

0228
0229   gErrorIgnoreLevel = kError;
0230   hist->Rebin(n_rebin);
0231   gErrorIgnoreLevel = kWarning;
0232   float peak_position_finder, peak_value_finder, shift_finder, scale_finder, N_finder, min_finder, max_finder;
0233   do_peak_finder(hist, peak_position_finder, peak_value_finder, shift_finder, scale_finder, N_finder, min_finder, max_finder);
0234   if (std::isnan(N_finder)) N_finder = hist->GetEntries();
0235   float fitting_best_chi2 = 1e8;
0236   float fitting_min_recorder, fitting_max_recorder;
0237   float fitting_min, fitting_max;
0238   for (int imin = 0; imin < min_nloop; ++imin) {
0239     bool reach_max = false;
0240     if (min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin <= shift_finder + min_lowlimit*bin_width*n_rebin) continue;
0241     if (reach_max) continue;
0242     if (min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin >= peak_position_finder - min_highlimit*bin_width*n_rebin) {
0243       fitting_min = peak_position_finder - min_highlimit*bin_width*n_rebin;
0244       reach_max = true;
0245     } else {
0246       fitting_min = min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin;
0247     }
0248     for (int imax = 0; imax < max_nloop; ++imax) {
0249       if (max_finder + imax*bin_width*n_rebin - max_lowloop*bin_width*n_rebin <= peak_position_finder + max_lowlimit*bin_width*n_rebin) continue;
0250       set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0251       fitting_max = max_finder + imax*bin_width*n_rebin - max_lowloop*bin_width*n_rebin;
0252       hist->Fit(func, "Q", "", fitting_min, fitting_max);
0253       float chi2ndf = func->GetChisquare()/(float)func->GetNDF();
0254       if (chi2ndf < fitting_best_chi2 && !std::isnan(func->GetParError(0)) && !std::isnan(func->GetParError(1)) && !std::isnan(func->GetParError(2)) && !std::isnan(func->GetParError(3))) {
0255         fitting_best_chi2 = chi2ndf;
0256         fitting_min_recorder = fitting_min;
0257         fitting_max_recorder = fitting_max;
0258       }
0259     }
0260   }
0261   set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0262   hist->Fit(func, "Q", "", fitting_min_recorder, fitting_max_recorder);
0263
0264   // Try to fit again if the error is large. (Usually caused by low fitting range.)

0265   if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0266     float fitting_min_recorder_test1 = fitting_min_recorder - 2*min_lowloop*bin_width;
0267     float fitting_max_recorder_test1 = fitting_max_recorder + 2*min_lowloop*bin_width;
0268     set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0269     hist->Fit(func, "Q", "", fitting_min_recorder_test1, fitting_max_recorder_test1);
0270     if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0271       fitting_min_recorder_test1 = fitting_min_recorder;
0272       fitting_max_recorder_test1 = fitting_max_recorder;
0273       set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0274       hist->Fit(func, "Q", "", fitting_min_recorder_test1, fitting_max_recorder_test1);
0275     }
0276     if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0277       float fitting_min_recorder_test2 = fitting_min_recorder - 4*min_lowloop*bin_width;
0278       float fitting_max_recorder_test2 = fitting_max_recorder + 4*min_lowloop*bin_width;
0279       set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0280       hist->Fit(func, "Q", "", fitting_min_recorder_test2, fitting_max_recorder_test2);
0281       if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0282         fitting_min_recorder_test2 = fitting_min_recorder_test1;
0283         fitting_max_recorder_test2 = fitting_max_recorder_test1;
0284         set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0285         hist->Fit(func, "Q", "", fitting_min_recorder_test2, fitting_max_recorder_test2);
0286       }
0287       if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0288         float fitting_min_recorder_test3 = fitting_min_recorder - 6*min_lowloop*bin_width;
0289         float fitting_max_recorder_test3 = fitting_max_recorder + 6*min_lowloop*bin_width;
0290         set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0291         hist->Fit(func, "Q", "", fitting_min_recorder_test3, fitting_max_recorder_test3);
0292         if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0293           fitting_min_recorder_test3 = fitting_min_recorder_test2;
0294           fitting_max_recorder_test3 = fitting_max_recorder_test2;
0295           set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0296           hist->Fit(func, "Q", "", fitting_min_recorder_test3, fitting_max_recorder_test3);
0297         }
0298       }
0299     }
0300   }
0301
0302   // A quick try to fix fitting failure. (Usually caused by parameter preset.)

0303   if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0304     set_parameters(func, peak_position_finder, 2, 1.5, N_finder);
0305     hist->Fit(func, "Q", "", fitting_min_recorder - min_lowloop*bin_width, fitting_max_recorder + 15*min_lowloop*bin_width);
0306   }
0307
0308
0309   // *** Second fitting try (with larger range) ***//

0310   if (func->GetChisquare()/(float)func->GetNDF() > 2 || std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0311     min_nloop = 15; // loop through min ?bins

0312     min_lowloop = 8; // loop from finder_min - ?bins

0313     max_nloop = 25; // loop through max ?bins

0314     max_lowloop = 5; // loop from finder_max - ?bins

0315
0316     min_lowlimit = 1; // min limit: shift + ?bins

0317     min_highlimit = 3; // min limit: peak - ?bins

0318     max_lowlimit = 4; // max limit: peak + ?bins

0319
0320     fitting_best_chi2 = 1e8;
0321     for (int imin = 0; imin < min_nloop; ++imin) {
0322       bool reach_max = false;
0323       if (min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin <= shift_finder + min_lowlimit*bin_width*n_rebin) continue;
0324       if (reach_max) continue;
0325       if (min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin >= peak_position_finder - min_highlimit*bin_width*n_rebin) {
0326         fitting_min = peak_position_finder - min_highlimit*bin_width*n_rebin;
0327         reach_max = true;
0328       } else {
0329         fitting_min = min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin;
0330       }
0331       for (int imax = 0; imax < max_nloop; ++imax) {
0332       if (max_finder + imax*bin_width*n_rebin - max_lowloop*bin_width*n_rebin <= peak_position_finder + max_lowlimit*bin_width*n_rebin) continue;
0333         set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0334         fitting_max = max_finder + imax*bin_width*n_rebin - max_lowloop*bin_width*n_rebin;
0335         hist->Fit(func, "Q", "", fitting_min, fitting_max);
0336         float chi2ndf = func->GetChisquare()/(float)func->GetNDF();
0337         if (chi2ndf < fitting_best_chi2 && !std::isnan(func->GetParError(0)) && !std::isnan(func->GetParError(1)) && !std::isnan(func->GetParError(2)) && !std::isnan(func->GetParError(3))) {
0338           fitting_best_chi2 = chi2ndf;
0339           fitting_min_recorder = fitting_min;
0340           fitting_max_recorder = fitting_max;
0341         }
0342       }
0343     }
0344     set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0345     hist->Fit(func, "Q", "", fitting_min_recorder, fitting_max_recorder);
0346
0347     // Try to fit again if the error is large. (Usually caused by low fitting range.)

0348     if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0349       float fitting_min_recorder_test1 = fitting_min_recorder - 2*min_lowloop*bin_width;
0350       float fitting_max_recorder_test1 = fitting_max_recorder + 2*min_lowloop*bin_width;
0351       set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0352       hist->Fit(func, "Q", "", fitting_min_recorder_test1, fitting_max_recorder_test1);
0353       if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0354         fitting_min_recorder_test1 = fitting_min_recorder;
0355         fitting_max_recorder_test1 = fitting_max_recorder;
0356         set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0357         hist->Fit(func, "Q", "", fitting_min_recorder_test1, fitting_max_recorder_test1);
0358       }
0359       if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0360         float fitting_min_recorder_test2 = fitting_min_recorder - 4*min_lowloop*bin_width;
0361         float fitting_max_recorder_test2 = fitting_max_recorder + 4*min_lowloop*bin_width;
0362         set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0363         hist->Fit(func, "Q", "", fitting_min_recorder_test2, fitting_max_recorder_test2);
0364         if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0365           fitting_min_recorder_test2 = fitting_min_recorder_test1;
0366           fitting_max_recorder_test2 = fitting_max_recorder_test1;
0367           set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0368           hist->Fit(func, "Q", "", fitting_min_recorder_test2, fitting_max_recorder_test2);
0369         }
0370         if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0371           float fitting_min_recorder_test3 = fitting_min_recorder - 6*min_lowloop*bin_width;
0372           float fitting_max_recorder_test3 = fitting_max_recorder + 6*min_lowloop*bin_width;
0373           set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0374           hist->Fit(func, "Q", "", fitting_min_recorder_test3, fitting_max_recorder_test3);
0375           if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0376             fitting_min_recorder_test3 = fitting_min_recorder_test2;
0377             fitting_max_recorder_test3 = fitting_max_recorder_test2;
0378             set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0379             hist->Fit(func, "Q", "", fitting_min_recorder_test3, fitting_max_recorder_test3);
0380           }
0381         }
0382       }
0383     }
0384   }
0385
0386
0387   // *** Third fitting try (with larger range and certain preset) ***//

0388   if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0389     min_nloop = 15; // loop through min ?bins

0390     min_lowloop = 8; // loop from finder_min - ?bins

0391     max_nloop = 25; // loop through max ?bins

0392     max_lowloop = 5; // loop from finder_max - ?bins

0393
0394     min_lowlimit = 1; // min limit: shift + ?bins

0395     min_highlimit = 3; // min limit: peak - ?bins

0396     max_lowlimit = 4; // max limit: peak + ?bins

0397
0398     fitting_best_chi2 = 1e8;
0399     for (int imin = 0; imin < min_nloop; ++imin) {
0400       bool reach_max = false;
0401       if (min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin <= shift_finder + min_lowlimit*bin_width*n_rebin) continue;
0402       if (reach_max) continue;
0403       if (min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin >= peak_position_finder - min_highlimit*bin_width*n_rebin) {
0404         fitting_min = peak_position_finder - min_highlimit*bin_width*n_rebin;
0405         reach_max = true;
0406       } else {
0407         fitting_min = min_finder + imin*bin_width*n_rebin - min_lowloop*bin_width*n_rebin;
0408       }
0409       for (int imax = 0; imax < max_nloop; ++imax) {
0410       if (max_finder + imax*bin_width*n_rebin - max_lowloop*bin_width*n_rebin <= peak_position_finder + max_lowlimit*bin_width*n_rebin) continue;
0411         set_parameters(func, peak_position_finder, 2, 1.5, N_finder);
0412         fitting_max = max_finder + imax*bin_width*n_rebin - max_lowloop*bin_width*n_rebin;
0413         hist->Fit(func, "Q", "", fitting_min, fitting_max);
0414         float chi2ndf = func->GetChisquare()/(float)func->GetNDF();
0415         if (chi2ndf < fitting_best_chi2 && !std::isnan(func->GetParError(0)) && !std::isnan(func->GetParError(1)) && !std::isnan(func->GetParError(2)) && !std::isnan(func->GetParError(3))) {
0416           fitting_best_chi2 = chi2ndf;
0417           fitting_min_recorder = fitting_min;
0418           fitting_max_recorder = fitting_max;
0419         }
0420       }
0421     }
0422     set_parameters(func, peak_position_finder, 2, 1.5, N_finder);
0423     hist->Fit(func, "Q", "", fitting_min_recorder, fitting_max_recorder);
0424
0425     // Try to fit again if the error is large. (Usually caused by low fitting range.)

0426     if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0427       float fitting_min_recorder_test1 = fitting_min_recorder - 2*min_lowloop*bin_width;
0428       float fitting_max_recorder_test1 = fitting_max_recorder + 2*min_lowloop*bin_width;
0429       set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0430       hist->Fit(func, "Q", "", fitting_min_recorder_test1, fitting_max_recorder_test1);
0431       if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0432         fitting_min_recorder_test1 = fitting_min_recorder;
0433         fitting_max_recorder_test1 = fitting_max_recorder;
0434         set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0435         hist->Fit(func, "Q", "", fitting_min_recorder_test1, fitting_max_recorder_test1);
0436       }
0437       if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0438         float fitting_min_recorder_test2 = fitting_min_recorder - 4*min_lowloop*bin_width;
0439         float fitting_max_recorder_test2 = fitting_max_recorder + 4*min_lowloop*bin_width;
0440         set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0441         hist->Fit(func, "Q", "", fitting_min_recorder_test2, fitting_max_recorder_test2);
0442         if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0443           fitting_min_recorder_test2 = fitting_min_recorder_test1;
0444           fitting_max_recorder_test2 = fitting_max_recorder_test1;
0445           set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0446           hist->Fit(func, "Q", "", fitting_min_recorder_test2, fitting_max_recorder_test2);
0447         }
0448         if (func->GetParError(0) / (float)func->GetParameter(0) > 1.1*(0.0121+9.323/(float)hist->GetEntries()) || func->GetParError(0) / (float)func->GetParameter(0) < 0) {
0449           float fitting_min_recorder_test3 = fitting_min_recorder - 6*min_lowloop*bin_width;
0450           float fitting_max_recorder_test3 = fitting_max_recorder + 6*min_lowloop*bin_width;
0451           set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0452           hist->Fit(func, "Q", "", fitting_min_recorder_test3, fitting_max_recorder_test3);
0453           if (std::isnan(func->GetParError(0)) || std::isnan(func->GetParError(1)) || std::isnan(func->GetParError(2)) || std::isnan(func->GetParError(3))) {
0454             fitting_min_recorder_test3 = fitting_min_recorder_test2;
0455             fitting_max_recorder_test3 = fitting_max_recorder_test2;
0456             set_parameters(func, peak_position_finder, shift_finder, scale_finder, N_finder);
0457             hist->Fit(func, "Q", "", fitting_min_recorder_test3, fitting_max_recorder_test3);
0458           }
0459         }
0460       }
0461     }
0462   }
0463 }