LightFlavorRatios/yield_and_ratios/SidebandFitter.h

0001 #ifndef SIDEBAND_FITTER_H
0002 #define SIDEBAND_FITTER_H
0003
0004 #include <TDatabasePDG.h>
0005
0006 // bin integration methods from ROOT tutorial rf614_binned_fit_problems.C
0007
0008 // Force numeric integration and do this numeric integration with the
0009 // RooBinIntegrator, which sums the function values at the bin centers.
0010 void enableBinIntegrator(RooAbsReal &func, int numBins)
0011 {
0012    RooNumIntConfig customConfig(*func.getIntegratorConfig());
0013    customConfig.method1D().setLabel("RooBinIntegrator");
0014    customConfig.getConfigSection("RooBinIntegrator").setRealValue("numBins", numBins);
0015    func.setIntegratorConfig(customConfig);
0016    func.forceNumInt(true);
0017 }
0018
0019 // Reset the integrator config to disable the RooBinIntegrator.
0020 void disableBinIntegrator(RooAbsReal &func)
0021 {
0022    func.setIntegratorConfig();
0023    func.forceNumInt(false);
0024 }
0025
0026 class SidebandFitter
0027 {
0028   public:
0029
0030   SidebandFitter(const std::string& nametag, const int pdgid,
0031                  const std::pair<double,double>& left_sideband, const std::pair<double,double>& right_sideband, bool binned, int nmassbins);
0032
0033   bool is_valid_fit(RooFitResult* result) const;
0034   void fit_orderN(RooAbsData* data, RooRealVar* mass, int n);
0035   void fit_all_orders_toN(RooAbsData* data, RooRealVar* mass, int n);
0036
0037   RooWorkspace* get_best_fit() const;
0038   double get_best_nsignal() const;
0039   double get_best_nsignal_err() const;
0040   double get_best_nbkg() const;
0041   double get_best_nbkg_err() const;
0042
0043   RooPlot* plot_best_fit(RooAbsData* data, RooRealVar* mass);
0044
0045   double get_yield_syserr() const;
0046
0047   void set_valid_binned_ranges(RooAbsData* data, RooRealVar* mass, double threshold);
0048
0049   private:
0050   std::string m_nametag;
0051   int m_pdgid;
0052   double m_pdg_mass;
0053   std::pair<double,double> m_left_sideband;
0054   std::pair<double,double> m_right_sideband;
0055   bool m_binned;
0056   int m_nmassbins;
0057
0058   // for yield extraction systematic uncertainty
0059   std::vector<double> m_valid_yields;
0060
0061   RooWorkspace* best_fit_ws = nullptr;
0062   double best_AIC = DBL_MAX;
0063 };
0064
0065 SidebandFitter::SidebandFitter(const std::string& nametag, const int pdgid,
0066                                const std::pair<double,double>& left_sideband, const std::pair<double,double>& right_sideband, bool binned, int nmassbins=1)
0067 : m_nametag(nametag), m_pdgid(pdgid), m_left_sideband(left_sideband), m_right_sideband(right_sideband), m_binned(binned), m_nmassbins(nmassbins)
0068 {
0069   m_pdg_mass = TDatabasePDG::Instance()->GetParticle(pdgid)->Mass();
0070 }
0071
0072 void SidebandFitter::set_valid_binned_ranges(RooAbsData* data, RooRealVar* mass, double threshold)
0073 {
0074   TH1* h = data->createHistogram("h_temp",*mass);
0075   int n_bins_rightsideband = h->FindBin(m_right_sideband.second) - h->FindBin(m_right_sideband.first);
0076   std::cout << "n_bins_rightsideband " << n_bins_rightsideband << std::endl;
0077   double avg_background = h->Integral(h->FindBin(m_right_sideband.first),h->FindBin(m_right_sideband.second))/n_bins_rightsideband;
0078   std::cout << "avg_background " << avg_background << std::endl;
0079   int last_left_sideband_bin = h->FindBin(m_left_sideband.second);
0080
0081   int last_bin_below_threshold = 1;
0082   for(int i=1;i<=last_left_sideband_bin;i++)
0083   {
0084     std::cout << "bin content " << h->GetBinContent(i) << std::endl;
0085     if(h->GetBinContent(i)<threshold*avg_background)
0086     {
0087       last_bin_below_threshold = i;
0088     }
0089   }
0090   std::cout << "last bin below threshold: " << last_bin_below_threshold << std::endl;
0091   double left_sideband_newmin = h->GetBinCenter(last_bin_below_threshold);
0092   if(m_left_sideband.second < left_sideband_newmin) m_left_sideband.second = left_sideband_newmin;
0093   mass->setRange("left_sideband",left_sideband_newmin,m_left_sideband.second);
0094   mass->setRange("full_range",left_sideband_newmin,m_right_sideband.second);
0095 }
0096
0097 void SidebandFitter::fit_orderN(RooAbsData* data, RooRealVar* mass, int n)
0098 {
0099   std::vector<RooRealVar> signal_parameters;
0100   std::vector<RooRealVar> sideband_parameters;
0101   std::vector<RooRealVar> background_parameters;
0102
0103   mass->setRange("signal",m_left_sideband.second,m_right_sideband.first);
0104   if(!m_binned)
0105   {
0106     mass->setRange("left_sideband",m_left_sideband.first,m_left_sideband.second);
0107   }
0108   else
0109   {
0110     set_valid_binned_ranges(data,mass,0.02);
0111   }
0112   mass->setRange("right_sideband",m_right_sideband.first,m_right_sideband.second);
0113
0114   std::string suffix = "_order"+std::to_string(n)+"_"+m_nametag;
0115
0116   RooRealVar n_signal("n_signal","nSignal",0.9*data->sumEntries(),0.,data->sumEntries());
0117   RooRealVar n_bkg("n_bkg","nBackground",0.1*data->sumEntries(),0.,data->sumEntries());
0118
0119   for(int order=0; order<n; order++)
0120   {
0121     std::string sideband_name = "q"+std::to_string(order+1);
0122     sideband_parameters.emplace_back(sideband_name.c_str(),sideband_name.c_str(),0.,-10.,10.);
0123     //background_parameters.emplace_back(bkg_name.c_str(),bkg_name.c_str(),0.,-10.,10.);
0124   }
0125
0126   std::string sideband_name = "sideband_order"+std::to_string(n)+"_"+m_nametag;
0127   RooChebychev sideband_pdf(sideband_name.c_str(),"sideband fit",*mass,RooArgList(sideband_parameters.begin(),sideband_parameters.end()));
0128
0129   //enableBinIntegrator(sideband_pdf,m_nmassbins);
0130
0131   if(m_binned)
0132   {
0133     sideband_pdf.fitTo(*data,RooFit::Range("left_sideband,right_sideband"),RooFit::Minimizer("Minuit2","minimize"),RooFit::Offset("bin"),RooFit::Strategy(2),RooFit::PrintLevel(0));
0134   }
0135   else
0136   {
0137     sideband_pdf.fitTo(*data,RooFit::Range("left_sideband,right_sideband"),RooFit::Minimizer("Minuit2","minimize"),RooFit::Offset(true),RooFit::Strategy(2),RooFit::PrintLevel(0));
0138   }
0139
0140   double signal_window_width = m_right_sideband.first-m_left_sideband.second;
0141 /*
0142   RooRealVar mean("mean","mean",m_pdg_mass,(m_pdg_mass-0.5*signal_window_width),(m_pdg_mass+0.5*signal_window_width));
0143   RooRealVar width1("width1","width1",0.2*signal_window_width,0.01*signal_window_width,1.*signal_window_width);
0144   RooRealVar width2("width2","width2",0.2*signal_window_width,0.01*signal_window_width,1.*signal_window_width);
0145
0146   RooGaussian gaus1("gaus1","gaus1",*mass,mean,width1);
0147   RooGaussian gaus2("gaus2","gaus2",*mass,mean,width2);
0148
0149   RooRealVar gaus2_frac("gaus2_frac","gaus2_frac",0.01,0.,1.);
0150
0151   RooAddPdf signal_pdf("signal","signal",RooArgList(gaus1,gaus2),RooArgList(gaus2_frac));
0152 */
0153
0154   RooRealVar mean("mean","mean",m_pdg_mass,(m_pdg_mass-0.5*signal_window_width),(m_pdg_mass+0.5*signal_window_width));
0155   RooRealVar width("width","width",0.5*signal_window_width,0.01*signal_window_width,2.*signal_window_width);
0156   RooRealVar alpha_L("alpha_L","alpha_L",2.,1.,10.);
0157   RooRealVar alpha_R("alpha_R","alpha_R",2.,1.,10.);
0158   RooRealVar n_L("n_L","n_L",2.,1.,100.);
0159   RooRealVar n_R("n_R","n_R",2.,1.,100.);
0160   //alpha_L.setConstant();
0161   //alpha_R.setConstant();
0162
0163   RooCrystalBall symm_signal_pdf("signal","signal",*mass,mean,width,alpha_L,n_L,true);
0164
0165   for(int i=0; i<sideband_parameters.size(); i++)
0166   {
0167     std::string bkg_name = "k"+std::to_string(i+1);
0168     double val = sideband_parameters[i].getVal();
0169     background_parameters.emplace_back(bkg_name.c_str(),bkg_name.c_str(),val,-10.,10.);
0170     //background_parameters[i].setVal(sideband_parameters[i].getVal());
0171     //background_parameters[i].setError(sideband_parameters[i].getError());
0172     background_parameters[i].setConstant();
0173   }
0174
0175   RooChebychev background_pdf("background","background",*mass,RooArgList(background_parameters.begin(),background_parameters.end()));
0176
0177   RooAddPdf symmetric_fit_pdf("symm_fit","symmetric fit",RooArgList(symm_signal_pdf,background_pdf),RooArgList(n_signal,n_bkg));
0178
0179   //enableBinIntegrator(symmetric_fit_pdf,m_nmassbins);
0180
0181   if(m_binned)
0182   {
0183     symmetric_fit_pdf.fitTo(*data,RooFit::Range("full_range"),RooFit::Minimizer("Minuit2","minimize"),RooFit::Strategy(2),RooFit::Offset("bin"),RooFit::PrintLevel(0));
0184   }
0185   else
0186   {
0187     symmetric_fit_pdf.fitTo(*data,RooFit::Minimizer("Minuit2","minimize"),RooFit::Strategy(2),RooFit::Offset(true),RooFit::PrintLevel(0));
0188   }
0189
0190   mean.setConstant();
0191   width.setConstant();
0192   alpha_R.setVal(alpha_L.getVal());
0193   n_R.setVal(n_L.getVal());
0194   //alpha_L.setConstant();
0195   //alpha_R.setConstant();
0196
0197   //for(int i=0;i<sideband_parameters.size();i++)
0198   //{
0199   //  background_parameters[i].setConstant();
0200  // }
0201
0202   RooCrystalBall asymm_signal_pdf("signal","signal",*mass,mean,width,alpha_L,n_L,alpha_R,n_R);
0203
0204   RooAddPdf full_fit_pdf("full_fit","full fit",RooArgList(asymm_signal_pdf,background_pdf),RooArgList(n_signal,n_bkg));
0205
0206   //enableBinIntegrator(full_fit_pdf,m_nmassbins);
0207
0208   RooFitResult* full_result;
0209   if(m_binned)
0210   {
0211     full_result = full_fit_pdf.fitTo(*data,RooFit::Range("full_range"),RooFit::Minimizer("Minuit2","minimize"),RooFit::Strategy(2),RooFit::Offset("bin"),RooFit::Save(true),RooFit::PrintLevel(0));
0212   }
0213   else
0214   {
0215     full_result = full_fit_pdf.fitTo(*data,RooFit::Minimizer("Minuit2","minimize"),RooFit::Strategy(2),RooFit::Offset(true),RooFit::Save(true),RooFit::PrintLevel(0));
0216   }
0217
0218   if(is_valid_fit(full_result))
0219   {
0220     double NLL = full_result->minNll();
0221     std::cout << "minNll = " << NLL << std::endl;
0222
0223     // Akaike information criterion (AIC) counters overfitting:
0224     // AIC = 2NDF - 2NLL
0225     // NDF = {background parameters} + (mean + width + alpha_L + n_L + alpha_R + n_R) + (nsignal + nbackground)
0226     int ndf = background_parameters.size() + 8;
0227     double AIC = 2.*ndf + 2*NLL;
0228
0229     std::cout << "NDF = " << ndf << ", AIC = " << AIC << std::endl;
0230     std::cout << "best AIC = " << best_AIC << std::endl;
0231
0232     m_valid_yields.push_back(n_signal.getVal());
0233
0234     if(AIC < best_AIC)
0235     {
0236       if(best_fit_ws) delete best_fit_ws;
0237
0238       RooWorkspace* ws = new RooWorkspace(("w"+suffix).c_str(),"workspace");
0239       ws->import(full_fit_pdf);
0240       ws->import(*mass->Clone());
0241
0242       best_fit_ws = ws;
0243       best_AIC = AIC;
0244     }
0245   }
0246   else
0247   {
0248     std::cout << "fit for order-" << n << " background failed" << std::endl;
0249   }
0250 }
0251
0252 bool SidebandFitter::is_valid_fit(RooFitResult* result) const
0253 {
0254   int status = result->status();
0255   int cov_status = result->covQual();
0256   std::cout << "status " << status << " covQual " << cov_status << std::endl;
0257   return (status==0 || status==1 || status==5) && (cov_status>=0);
0258 }
0259
0260 void SidebandFitter::fit_all_orders_toN(RooAbsData* data, RooRealVar* mass, int n)
0261 {
0262   for(int order=1; order<=n; order++)
0263   {
0264     fit_orderN(data,mass,order);
0265   }
0266 }
0267
0268 RooWorkspace* SidebandFitter::get_best_fit() const
0269 {
0270   return best_fit_ws;
0271 }
0272
0273 double SidebandFitter::get_best_nsignal() const
0274 {
0275   RooRealVar* nsignal = best_fit_ws->var("n_signal");
0276   return nsignal->getVal();
0277 }
0278
0279 double SidebandFitter::get_best_nsignal_err() const
0280 {
0281   RooRealVar* nsignal = best_fit_ws->var("n_signal");
0282   return nsignal->getError();
0283 }
0284
0285 double SidebandFitter::get_best_nbkg() const
0286 {
0287   RooRealVar* nbkg = best_fit_ws->var("n_bkg");
0288   return nbkg->getVal();
0289 }
0290
0291 double SidebandFitter::get_best_nbkg_err() const
0292 {
0293   RooRealVar* nbkg = best_fit_ws->var("n_bkg");
0294   return nbkg->getError();
0295 }
0296
0297 RooPlot* SidebandFitter::plot_best_fit(RooAbsData* data, RooRealVar* mass)
0298 {
0299   RooAbsPdf* full_fit = best_fit_ws->pdf("full_fit");
0300
0301   RooPlot* plot;
0302   if(m_binned)
0303   {
0304     plot = mass->frame(RooFit::Range("full_range"));
0305   }
0306   else
0307   {
0308     plot = mass->frame();
0309   }
0310   data->plotOn(plot);
0311   if(m_binned)
0312   {
0313     full_fit->plotOn(plot,RooFit::Range("full_range"),RooFit::NormRange("full_range"),RooFit::Components("background"),RooFit::DrawOption("FL"),RooFit::LineStyle(kDashed),RooFit::FillColor(kGray),RooFit::MoveToBack());
0314     full_fit->plotOn(plot,RooFit::Range("full_range"),RooFit::NormRange("full_range"),RooFit::DrawOption("FL"),RooFit::FillColor(kAzure+1),RooFit::MoveToBack());
0315   }
0316   else
0317   {
0318     full_fit->plotOn(plot,RooFit::Components("background"),RooFit::DrawOption("FL"),RooFit::LineStyle(kDashed),RooFit::FillColor(kGray),RooFit::MoveToBack());
0319     full_fit->plotOn(plot,RooFit::DrawOption("FL"),RooFit::FillColor(kAzure+1),RooFit::MoveToBack());
0320   }
0321
0322   return plot;
0323 }
0324
0325 double SidebandFitter::get_yield_syserr() const
0326 {
0327   //auto minmax_yield = std::minmax_element(m_valid_yields.begin(),m_valid_yields.end());
0328   //return (*minmax_yield.second-*minmax_yield.first)/(sqrt(12.)*get_best_nsignal());
0329
0330   const double best_fit_yield = get_best_nsignal();
0331
0332   std::vector<double> relative_diffs;
0333   for(double yield : m_valid_yields)
0334   {
0335     relative_diffs.push_back(fabs(best_fit_yield - yield)/best_fit_yield);
0336   }
0337   std::sort(relative_diffs.begin(),relative_diffs.end());
0338
0339   int median_i = relative_diffs.size() / 2;
0340   if(relative_diffs.size() % 2 == 0)
0341   {
0342     return (relative_diffs[median_i-1] + relative_diffs[median_i])/2.;
0343   }
0344   else
0345   {
0346     return relative_diffs[median_i];
0347   }
0348
0349 }
0350
0351 #endif