LightFlavorRatios/yield_and_ratios/ResonanceRatio.h

0001 #ifndef RESONANCE_RATIO_H
0002 #define RESONANCE_RATIO_H
0003
0004 #include "TFile.h"
0005 #include "TH1F.h"
0006
0007 #include "RooDataSet.h"
0008 #include "RooDataHist.h"
0009 #include "RooPlot.h"
0010
0011 #include <phool/PHRandomSeed.h>
0012 #include <gsl/gsl_rng.h>
0013
0014 #include "../util/binning.h"
0015 #include "../util/DifferentialContainer.h"
0016 #include "ParticleModel.h"
0017
0018 class ResonanceRatio
0019 {
0020   public:
0021   ResonanceRatio(ParticleModel& numerator_model, ParticleModel& denominator_model,
0022                  TFile* outf, std::string rationame, std::string ratiotitle, float scalefactor, bool blind,
0023                  std::vector<HistogramInfo> variables, std::vector<std::vector<std::shared_ptr<CorrectionHistogram1D>>> corrections)
0024   : _numerator_model(numerator_model), _denominator_model(denominator_model),
0025     outfile(outf), _rationame(rationame), _ratiotitle(ratiotitle), _scalefactor(scalefactor), _blind(blind),
0026     _variables(variables), _corrections(corrections)
0027   {
0028     setup_yield_histograms();
0029   }
0030
0031   // main analysis functions
0032   void calculate_ratios_unbinned(RooAbsData& numerator_data, RooAbsData& denominator_data);
0033   void calculate_ratios_binned(TH1F* integrated_numerator_data, std::vector<DifferentialContainer>& diff_numerator_data,
0034                                TH1F* integrated_denominator_data, std::vector<DifferentialContainer>& diff_denominator_data);
0035
0036   // yield histograms
0037   TH1F* numerator_integrated_yield;
0038   TH1F* denominator_integrated_yield;
0039   std::vector<TH1F*> numerator_diff_yields;
0040   std::vector<TH1F*> denominator_diff_yields;
0041   // ratio histograms
0042   std::vector<TH1F*> integrated_ratio_w_corrections;
0043   std::vector<std::vector<TH1F*>> diff_ratios_w_corrections;
0044   // output file
0045   TFile* outfile;
0046
0047   protected:
0048   // workflow functions
0049   void setup_yield_histograms();
0050   void get_yield(TH1F* h_yield, int i, RooAbsData& ds, ParticleModel& model);
0051   void get_diff_yield_unbinned(TH1F* h_yield, HistogramInfo& hinfo, RooAbsData& data, ParticleModel& model);
0052   void get_diff_yield_binned(TH1F* h_yield, HistogramInfo& hinfo, DifferentialContainer& data, ParticleModel& model);
0053   std::string get_corrected_title(std::string current_title, std::string correction_title);
0054   void calculate_ratios();
0055   void save_results();
0056
0057   // fit models
0058   ParticleModel _numerator_model;
0059   ParticleModel _denominator_model;
0060   // differential variables and corresponding slates of corrections to apply
0061   std::vector<HistogramInfo> _variables;
0062   std::vector<std::vector<std::shared_ptr<CorrectionHistogram1D>>> _corrections;
0063   // naming, scale factor, and blind settings
0064   std::string _rationame;
0065   std::string _ratiotitle;
0066   float _scalefactor;
0067   bool _blind;
0068 };
0069
0070 void ResonanceRatio::setup_yield_histograms()
0071 {
0072   HistogramInfo numerator_massbins = BinInfo::mass_bins.at(_numerator_model.name);
0073   HistogramInfo denominator_massbins = BinInfo::mass_bins.at(_denominator_model.name);
0074
0075   for(HistogramInfo& hinfo : _variables)
0076   {
0077     numerator_diff_yields.push_back(makeHistogram(_numerator_model.name+"_yield",_numerator_model.name+" yield",hinfo));
0078     denominator_diff_yields.push_back(makeHistogram(_denominator_model.name+"_yield",_denominator_model.name+" yield",hinfo));
0079   }
0080
0081   numerator_integrated_yield = new TH1F(("all_"+_numerator_model.name+"_yield").c_str(),("All "+_numerator_model.name+" yield").c_str(),1,0.,1.);
0082   denominator_integrated_yield = new TH1F(("all_"+_denominator_model.name+"_yield").c_str(),("All "+_denominator_model.name+" yield").c_str(),1,0.,1.);
0083
0084   numerator_integrated_yield->SetTitle((numerator_massbins.title+";"+numerator_massbins.axis_label+";Candidates").c_str());
0085   denominator_integrated_yield->SetTitle((denominator_massbins.title+";"+denominator_massbins.axis_label+";Candidates").c_str());
0086 }
0087
0088 void ResonanceRatio::get_yield(TH1F* h_yield, int i, RooAbsData& ds, ParticleModel& model)
0089 {
0090   model.fitTo(ds);
0091
0092   double nsignal = model.n_signal->getVal();
0093   double nsignal_err = model.n_signal->getError();
0094
0095   double nbkg = model.n_background->getVal();
0096   double nbkg_err = model.n_background->getError();
0097
0098   double yield = ds.sumEntries() - nbkg;
0099   double yield_err = nbkg_err;
0100
0101   std::cout << "nsignal val " << nsignal << std::endl;
0102   std::cout << "nsignal err " << nsignal_err << std::endl;
0103   std::cout << "nbkg " << nbkg << std::endl;
0104   std::cout << "nbkg err " << nbkg_err << std::endl;
0105   std::cout << "sum entries " << ds.sumEntries() << std::endl;
0106   std::cout << "yield " << yield << std::endl;
0107
0108   h_yield->SetBinContent(i,yield);
0109   h_yield->SetBinError(i,yield_err);
0110
0111   std::string name;
0112   if(i>=0) name = std::string(h_yield->GetName())+"_"+std::to_string(i);
0113   else name = std::string(h_yield->GetName())+"_";
0114   std::string title;
0115   if(i>=0) title = std::string(h_yield->GetTitle())+" bin "+std::to_string(i);
0116   else title = std::string(h_yield->GetTitle());
0117
0118   RooPlot* plot = model.mass->frame(RooFit::Title(title.c_str()));
0119   plot->SetName(name.c_str());
0120   ds.plotOn(plot);
0121   model.fit_function->plotOn(plot,RooFit::Components(model.background_function->GetName()),RooFit::DrawOption("FL"),RooFit::LineStyle(kDashed),RooFit::FillColor(kGray),RooFit::MoveToBack());
0122   model.fit_function->plotOn(plot,RooFit::DrawOption("FL"),RooFit::FillColor(kAzure+1),RooFit::MoveToBack());
0123   plot->Write();
0124 }
0125
0126 void ResonanceRatio::get_diff_yield_unbinned(TH1F* h_yield, HistogramInfo& hinfo, RooAbsData& data, ParticleModel& model)
0127 {
0128   for(int i=1; i<=h_yield->GetNbinsX(); i++)
0129   {
0130     std::cout << "bin " << i << " of " << h_yield->GetNbinsX() << std::endl;
0131
0132     std::string selection = hinfo.get_bin_selection(std::string(data.GetName())+"_"+hinfo.name,i);
0133     std::cout << "selection: " << selection << std::endl;
0134
0135     RooDataSet ds_selected = static_cast<RooDataSet&>(*(data.reduce({*(model.mass)},selection.c_str())));
0136
0137     get_yield(h_yield,i,ds_selected,model);
0138   }
0139 }
0140
0141 void ResonanceRatio::get_diff_yield_binned(TH1F* h_yield, HistogramInfo& hinfo, DifferentialContainer& data, ParticleModel& model)
0142 {
0143   for(int i=1; i<=h_yield->GetNbinsX(); i++)
0144   {
0145     std::cout << "bin " << i << " of " << h_yield->GetNbinsX() << std::endl;
0146     RooDataHist dh("binned_massfit","binned_massfit",*(model.mass),RooFit::Import(*(data.hists[i])));
0147     get_yield(h_yield,i,dh,model);
0148   }
0149 }
0150
0151 std::string ResonanceRatio::get_corrected_title(std::string current_title,std::string correction_title)
0152 {
0153   std::string new_title;
0154   int semicolon_pos = current_title.find(";");
0155   if(semicolon_pos<current_title.length())
0156   {
0157     std::string main_title = current_title.substr(0,current_title.find(";"));
0158     std::string axis_labels = current_title.substr(current_title.find(";"),current_title.length());
0159     std::string new_main_title = main_title + ", " + correction_title + " corrected";
0160     new_title = new_main_title + axis_labels;
0161   }
0162   else
0163   {
0164     new_title = current_title + ", " + correction_title + " corrected";
0165   }
0166   return new_title;
0167 }
0168
0169 void ResonanceRatio::calculate_ratios()
0170 {
0171   // calculate integrated ratio and apply corrections
0172
0173   TH1F* integrated_ratio = new TH1F("integrated_ratio",("Integrated "+_rationame).c_str(),1,0.,1.);
0174   float integrated_ratio_val = numerator_integrated_yield->GetBinContent(1)/denominator_integrated_yield->GetBinContent(1)*_scalefactor;
0175   float integrated_ratio_err = integrated_ratio_val * sqrt(pow(numerator_integrated_yield->GetBinError(1)/numerator_integrated_yield->GetBinContent(1),2.)+
0176                                                            pow(denominator_integrated_yield->GetBinError(1)/denominator_integrated_yield->GetBinContent(1),2.));
0177   integrated_ratio->SetBinContent(1,integrated_ratio_val);
0178   integrated_ratio->SetBinError(1,integrated_ratio_err);
0179
0180   integrated_ratio_w_corrections.push_back(integrated_ratio);
0181
0182   for(std::shared_ptr<CorrectionHistogram1D> correction : _corrections.at(0)) // retrieve one set of corrections to integrate over
0183   {
0184     TH1F* current_h = integrated_ratio_w_corrections.back();
0185     std::string name = std::string(current_h->GetName()) + "_" + correction->name + "corrected";
0186
0187     std::string current_title = std::string(current_h->GetTitle());
0188     std::string new_title = get_corrected_title(current_title,correction->title);
0189
0190     TH1F* new_h = (TH1F*)current_h->Clone(name.c_str());
0191     new_h->SetTitle(new_title.c_str());
0192
0193     correction->apply_correction(correction->get_xmin(),correction->get_xmax(),new_h,1);
0194
0195     integrated_ratio_w_corrections.push_back(new_h);
0196   }
0197
0198   // calculate differential ratios
0199
0200   for(int ivar=0; ivar<_variables.size(); ivar++)
0201   {
0202     TH1F* diffratio = makeHistogram(_rationame,_ratiotitle,_variables[ivar]);
0203     for(int ibin=1; ibin<=diffratio->GetNbinsX(); ibin++)
0204     {
0205       float ratio_val = numerator_diff_yields[ivar]->GetBinContent(ibin) * _scalefactor / denominator_diff_yields[ivar]->GetBinContent(ibin);
0206       float ratio_err = ratio_val * sqrt(pow(numerator_diff_yields[ivar]->GetBinError(ibin)/numerator_diff_yields[ivar]->GetBinContent(ibin),2.)+pow(denominator_diff_yields[ivar]->GetBinError(ibin)/denominator_diff_yields[ivar]->GetBinContent(ibin),2.));
0207       diffratio->SetBinContent(ibin,ratio_val);
0208       diffratio->SetBinError(ibin,ratio_err);
0209     }
0210     diff_ratios_w_corrections.push_back({diffratio});
0211   }
0212
0213   // apply corrections to differential ratios
0214
0215   for(int ivar=0; ivar<_variables.size(); ivar++)
0216   {
0217     for(int icorr=0; icorr<_corrections[ivar].size(); icorr++)
0218     {
0219       TH1F* current_h = diff_ratios_w_corrections[ivar][icorr];
0220       std::shared_ptr<CorrectionHistogram1D> correction = _corrections[ivar][icorr];
0221       std::string name = std::string(current_h->GetName()) + "_" + correction->name + "corrected";
0222
0223       // transforming title is a bit complicated due to axis labels being in the way
0224       std::string current_title = std::string(current_h->GetTitle());
0225       std::string new_title = get_corrected_title(current_title,correction->title);
0226
0227       TH1F* new_h = (TH1F*)current_h->Clone(name.c_str());
0228       new_h->SetTitle(new_title.c_str());
0229
0230       for(int ibin=1; ibin<=new_h->GetNbinsX(); ibin++)
0231       {
0232         float xlow = new_h->GetBinLowEdge(ibin);
0233         float xhigh = xlow + new_h->GetBinWidth(ibin);
0234         std::cout << "ivar " << ivar << " icorr " << icorr << " ibin " << ibin << std::endl;
0235         std::cout << "before correction: " << new_h->GetBinContent(ibin) << " +- " << new_h->GetBinError(ibin) << std::endl;
0236         correction->apply_correction(xlow,xhigh,new_h,ibin);
0237         std::cout << "after " << correction->title << " correction: " << new_h->GetBinContent(ibin) << " +- " << new_h->GetBinError(ibin) << std::endl;
0238       }
0239       diff_ratios_w_corrections[ivar].push_back(new_h);
0240     }
0241   }
0242
0243   // apply blinding
0244
0245   if (_blind)
0246   {
0247     const uint seed = PHRandomSeed();
0248     std::unique_ptr<gsl_rng> m_rng;
0249     m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
0250     gsl_rng_set(m_rng.get(), seed);
0251     float blind_par = 9.9*gsl_rng_uniform_pos(m_rng.get())+0.1; //Take a value anywhere between 0.1 and 10.
0252
0253     for(TH1F* h : integrated_ratio_w_corrections)
0254     {
0255       float content = h->GetBinContent(1);
0256       h->SetBinContent(1,content * blind_par);
0257       h->SetBinError(1,h->GetBinError(1) * blind_par);
0258     }
0259
0260     for(std::vector<TH1F*>& v_corr : diff_ratios_w_corrections)
0261     {
0262       for(TH1F* h : v_corr)
0263       {
0264         for(int ibin=0; ibin<=h->GetNbinsX(); ibin++)
0265         {
0266           float content = h->GetBinContent(ibin);
0267           h->SetBinContent(ibin,content * blind_par);
0268           h->SetBinError(ibin,h->GetBinError(ibin) * blind_par);
0269         }
0270       }
0271     }
0272   }
0273 }
0274
0275 void ResonanceRatio::save_results()
0276 {
0277   // Save everything to file
0278
0279   outfile->cd();
0280
0281   for(TH1F* h : integrated_ratio_w_corrections)
0282   {
0283     h->Write();
0284   }
0285
0286   for(std::vector<TH1F*>& v_corr : diff_ratios_w_corrections)
0287   {
0288     for(TH1F* h : v_corr)
0289     {
0290       h->Write();
0291     }
0292   }
0293 }
0294
0295 void ResonanceRatio::calculate_ratios_unbinned(RooAbsData& numerator_data, RooAbsData& denominator_data)
0296 {
0297   get_yield(numerator_integrated_yield,-1,numerator_data,_numerator_model);
0298   get_yield(denominator_integrated_yield,-1,denominator_data,_denominator_model);
0299
0300   // extract differential yields
0301
0302   for(size_t i=0; i<_variables.size(); i++)
0303   {
0304     std::cout << "======= Differential " << numerator_diff_yields[i]->GetName() << " =======" << std::endl;
0305     get_diff_yield_unbinned(numerator_diff_yields[i],_variables[i],numerator_data,_numerator_model);
0306     std::cout << "======= Differential " << denominator_diff_yields[i]->GetName() << " =======" << std::endl;
0307     get_diff_yield_unbinned(denominator_diff_yields[i],_variables[i],denominator_data,_denominator_model);
0308   }
0309
0310   calculate_ratios();
0311
0312   save_results();
0313 }
0314
0315 void ResonanceRatio::calculate_ratios_binned(TH1F* integrated_numerator_data, std::vector<DifferentialContainer>& diff_numerator_data,
0316                                              TH1F* integrated_denominator_data, std::vector<DifferentialContainer>& diff_denominator_data)
0317 {
0318   RooDataHist integrated_numerator_dh("integrated_numerator_dh","integrated_numerator_dh",RooArgList(*(_numerator_model.mass)),RooFit::Import(*integrated_numerator_data));
0319   RooDataHist integrated_denominator_dh("integrated_denominator_dh","integrated_denominator_dh",RooArgList(*(_denominator_model.mass)),RooFit::Import(*integrated_denominator_data));
0320
0321   get_yield(numerator_integrated_yield,-1,integrated_numerator_dh,_numerator_model);
0322   get_yield(denominator_integrated_yield,-1,integrated_denominator_dh,_denominator_model);
0323
0324   for(int i=0; i<_variables.size(); i++)
0325   {
0326     get_diff_yield_binned(numerator_diff_yields[i],_variables[i],diff_numerator_data[i],_numerator_model);
0327     get_diff_yield_binned(denominator_diff_yields[i],_variables[i],diff_denominator_data[i],_denominator_model);
0328   }
0329 }
0330
0331 #endif