Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:32:57

0001 #include "Alignment/OfflineValidation/interface/FitWithRooFit.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 
0005 // Import TH1 histogram into a RooDataHist
0006 std::unique_ptr<RooDataHist> FitWithRooFit::importTH1(TH1* histo, double xMin, double xMax) {
0007   if ((xMin == xMax) && xMin == 0) {
0008     xMin = histo->GetXaxis()->GetXmin();
0009     xMax = histo->GetXaxis()->GetXmax();
0010   }
0011   // Declare observable x
0012   RooRealVar x("InvMass", "di-muon mass M(#mu^{+}#mu^{-}) [GeV]", xMin, xMax);
0013   // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
0014   return std::make_unique<RooDataHist>("dh", "dh", x, RooFit::Import(*histo));
0015 }
0016 
0017 // Plot and fit a RooDataHist fitting signal and background
0018 void FitWithRooFit::fit(
0019     TH1* histo, const TString signalType, const TString backgroundType, double xMin, double xMax, bool sumW2Error) {
0020   reinitializeParameters();
0021 
0022   std::unique_ptr<RooDataHist> dh = importTH1(histo, xMin, xMax);
0023   RooRealVar x(*static_cast<RooRealVar*>(dh->get()->find("x")));
0024 
0025   // Make plot of binned dataset showing Poisson error bars (RooFit default)
0026   RooPlot* frame = x.frame(RooFit::Title("di-muon mass M(#mu^{+}#mu^{-}) [GeV]"));
0027   frame->SetName(TString(histo->GetName()) + "_frame");
0028   dh->plotOn(frame);
0029 
0030   // Build the composite model
0031   std::unique_ptr<RooAbsPdf> model = buildModel(&x, signalType, backgroundType);
0032 
0033   std::unique_ptr<RooAbsReal> chi2{model->createChi2(*dh, RooFit::DataError(RooAbsData::SumW2))};
0034 
0035   // Fit the composite model
0036   // -----------------------
0037   // Fit with likelihood
0038   if (!useChi2_) {
0039     model->fitTo(*dh, RooFit::SumW2Error(sumW2Error));
0040   }
0041   // Fit with chi^2
0042   else {
0043     edm::LogWarning("FitWithRooFit") << "FITTING WITH CHI^2";
0044     RooMinimizer m(*chi2);
0045     m.migrad();
0046     m.hesse();
0047     // RooFitResult* r_chi2_wgt = m.save();
0048   }
0049   model->plotOn(frame, RooFit::LineColor(kRed));
0050   model->plotOn(frame, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed));
0051   model->paramOn(frame,
0052                  RooFit::Label("fit result"),
0053                  RooFit::Layout(0.65, 0.90, 0.90),
0054                  RooFit::Format("NEU", RooFit::AutoPrecision(2)));
0055 
0056   // TODO: fix next lines to get the prob(chi2) (ndof should be dynamically set according to the choosen pdf)
0057   // double chi2 = xframe.chiSquare("model","data",ndof);
0058   // double ndoff = xframeGetNbinsX();
0059   // double chi2prob = TMath::Prob(chi2,ndoff);
0060 
0061   // P l o t   a n d   f i t   a   R o o D a t a H i s t   w i t h   i n t e r n a l   e r r o r s
0062   // ---------------------------------------------------------------------------------------------
0063 
0064   // If histogram has custom error (i.e. its contents is does not originate from a Poisson process
0065   // but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead
0066   // (same error bars as shown by ROOT)
0067   RooPlot* frame2 = x.frame(RooFit::Title("di-muon mass M(#mu^{+}#mu^{-}) [GeV]"));
0068   dh->plotOn(frame2, RooFit::DataError(RooAbsData::SumW2));
0069   model->plotOn(frame2, RooFit::LineColor(kRed));
0070   model->plotOn(frame2, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed));
0071   model->paramOn(frame2,
0072                  RooFit::Label("fit result"),
0073                  RooFit::Layout(0.65, 0.90, 0.90),
0074                  RooFit::Format("NEU", RooFit::AutoPrecision(2)));
0075 
0076   // Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used
0077   // in a maximum likelihood fit
0078   //
0079   // A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition
0080   // of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly
0081   // fitted with a chi^2 fit (see rf602_chi2fit.C)
0082 
0083   // Draw all frames on a canvas
0084   // if( canvas == 0 ) {
0085   //   canvas = new TCanvas("rf102_dataimport","rf102_dataimport",800,800);
0086   //   canvas->cd(1);
0087   // }
0088   // canvas->Divide(2,1);
0089   // canvas->cd(1) ; frame->Draw();
0090   // canvas->cd(2) ; frame2->Draw();
0091 
0092   frame2->Draw();
0093 }
0094 
0095 // Initialization methods for all the parameters
0096 void FitWithRooFit::initMean(double value, double min, double max, const TString& name, const TString& title) {
0097   mean_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0098   initVal_mean = value;
0099 }
0100 
0101 void FitWithRooFit::initMean2(double value, double min, double max, const TString& name, const TString& title) {
0102   mean2_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0103   initVal_mean2 = value;
0104 }
0105 
0106 void FitWithRooFit::initMean3(double value, double min, double max, const TString& name, const TString& title) {
0107   mean3_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0108   initVal_mean3 = value;
0109 }
0110 
0111 void FitWithRooFit::initSigma(double value, double min, double max, const TString& name, const TString& title) {
0112   sigma_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0113   initVal_sigma = value;
0114 }
0115 
0116 void FitWithRooFit::initSigma2(double value, double min, double max, const TString& name, const TString& title) {
0117   sigma2_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0118   initVal_sigma2 = value;
0119 }
0120 
0121 void FitWithRooFit::initSigma3(double value, double min, double max, const TString& name, const TString& title) {
0122   sigma3_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0123   initVal_sigma3 = value;
0124 }
0125 
0126 void FitWithRooFit::initGamma(double value, double min, double max, const TString& name, const TString& title) {
0127   gamma_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0128   initVal_gamma = value;
0129 }
0130 
0131 void FitWithRooFit::initGaussFrac(double value, double min, double max, const TString& name, const TString& title) {
0132   gaussFrac_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0133   initVal_gaussFrac = value;
0134 }
0135 
0136 void FitWithRooFit::initGaussFrac2(double value, double min, double max, const TString& name, const TString& title) {
0137   gaussFrac2_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0138   initVal_gaussFrac2 = value;
0139 }
0140 
0141 void FitWithRooFit::initExpCoeffA0(double value, double min, double max, const TString& name, const TString& title) {
0142   expCoeffa0_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0143   initVal_expCoeffa0 = value;
0144 }
0145 
0146 void FitWithRooFit::initExpCoeffA1(double value, double min, double max, const TString& name, const TString& title) {
0147   expCoeffa1_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0148   initVal_expCoeffa1 = value;
0149 }
0150 
0151 void FitWithRooFit::initExpCoeffA2(double value, double min, double max, const TString& name, const TString& title) {
0152   expCoeffa2_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0153   initVal_expCoeffa2 = value;
0154 }
0155 
0156 void FitWithRooFit::initFsig(double value, double min, double max, const TString& name, const TString& title) {
0157   fsig_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0158   initVal_fsig = value;
0159 }
0160 
0161 void FitWithRooFit::initA0(double value, double min, double max, const TString& name, const TString& title) {
0162   a0_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0163   initVal_a0 = value;
0164 }
0165 
0166 void FitWithRooFit::initA1(double value, double min, double max, const TString& name, const TString& title) {
0167   a1_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0168   initVal_a1 = value;
0169 }
0170 
0171 void FitWithRooFit::initA2(double value, double min, double max, const TString& name, const TString& title) {
0172   a2_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0173   initVal_a2 = value;
0174 }
0175 
0176 void FitWithRooFit::initA3(double value, double min, double max, const TString& name, const TString& title) {
0177   a3_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0178   initVal_a3 = value;
0179 }
0180 
0181 void FitWithRooFit::initA4(double value, double min, double max, const TString& name, const TString& title) {
0182   a4_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0183   initVal_a4 = value;
0184 }
0185 
0186 void FitWithRooFit::initA5(double value, double min, double max, const TString& name, const TString& title) {
0187   a5_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0188   initVal_a5 = value;
0189 }
0190 
0191 void FitWithRooFit::initA6(double value, double min, double max, const TString& name, const TString& title) {
0192   a6_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0193   initVal_a6 = value;
0194 }
0195 
0196 void FitWithRooFit::initAlpha(double value, double min, double max, const TString& name, const TString& title) {
0197   alpha_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0198   initVal_alpha = value;
0199 }
0200 
0201 void FitWithRooFit::initN(double value, double min, double max, const TString& name, const TString& title) {
0202   n_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0203   initVal_n = value;
0204 }
0205 
0206 void FitWithRooFit::initFGCB(double value, double min, double max, const TString& name, const TString& title) {
0207   fGCB_ = std::make_unique<RooRealVar>(name, title, value, min, max);
0208   initVal_fGCB = value;
0209 }
0210 
0211 void FitWithRooFit::reinitializeParameters() {
0212   auto initParam = [](std::unique_ptr<RooRealVar>& var, double val) {
0213     if (var)
0214       var->setVal(val);
0215   };
0216 
0217   initParam(mean_, initVal_mean);
0218   initParam(mean2_, initVal_mean2);
0219   initParam(mean3_, initVal_mean3);
0220   initParam(sigma_, initVal_sigma);
0221   initParam(sigma2_, initVal_sigma2);
0222   initParam(sigma3_, initVal_sigma3);
0223   initParam(gamma_, initVal_gamma);
0224   initParam(gaussFrac_, initVal_gaussFrac);
0225   initParam(gaussFrac2_, initVal_gaussFrac2);
0226   initParam(expCoeffa0_, initVal_expCoeffa0);
0227   initParam(expCoeffa1_, initVal_expCoeffa1);
0228   initParam(expCoeffa2_, initVal_expCoeffa2);
0229   initParam(fsig_, initVal_fsig);
0230   initParam(a0_, initVal_a0);
0231   initParam(a1_, initVal_a1);
0232   initParam(a2_, initVal_a2);
0233   initParam(a3_, initVal_a3);
0234   initParam(a4_, initVal_a4);
0235   initParam(a5_, initVal_a5);
0236   initParam(a6_, initVal_a6);
0237   initParam(alpha_, initVal_alpha);
0238   initParam(n_, initVal_n);
0239   initParam(fGCB_, initVal_fGCB);
0240 }
0241 
0242 /// Build the model for the specified signal type
0243 std::unique_ptr<RooAbsPdf> FitWithRooFit::buildSignalModel(RooRealVar* x, const TString& signalType) {
0244   if (signalType == "gaussian") {
0245     // Fit a Gaussian p.d.f to the data
0246     if ((mean_ == nullptr) || (sigma_ == nullptr)) {
0247       edm::LogError("FitWithRooFit")
0248           << "Error: one or more parameters are not initialized. Please be sure to initialize mean and sigma";
0249       exit(1);
0250     }
0251     return std::make_unique<RooGaussian>("gauss", "gauss", *x, *mean_, *sigma_);
0252   } else if (signalType == "doubleGaussian") {
0253     // Fit with double gaussian
0254     if ((mean_ == nullptr) || (sigma_ == nullptr) || (sigma2_ == nullptr)) {
0255       edm::LogError("FitWithRooFit")
0256           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and sigma2";
0257       exit(1);
0258     }
0259     RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_);
0260     RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean_, *sigma2_);
0261     RooArgList components{*gaussModel, *gaussModel2};
0262     auto out = std::make_unique<RooAddModel>("doubleGaussian", "double gaussian", components, *gaussFrac_);
0263     out->addOwnedComponents(components);
0264     return out;
0265   } else if (signalType == "tripleGaussian") {
0266     // Fit with triple gaussian
0267     if ((mean_ == nullptr) || (mean2_ == nullptr) || (mean3_ == nullptr) || (sigma_ == nullptr) ||
0268         (sigma2_ == nullptr) || (sigma3_ == nullptr)) {
0269       edm::LogError("FitWithRooFit")
0270           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0271              "mean3, sigma, sigma2, sigma3";
0272       exit(1);
0273     }
0274     RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_);
0275     RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean2_, *sigma2_);
0276     RooGaussModel* gaussModel3 = new RooGaussModel("gaussModel3", "gaussModel3", *x, *mean3_, *sigma3_);
0277     RooArgList components{*gaussModel, *gaussModel2, *gaussModel3};
0278     auto out = std::make_unique<RooAddModel>(
0279         "tripleGaussian", "triple gaussian", components, RooArgList{*gaussFrac_, *gaussFrac2_});
0280     out->addOwnedComponents(components);
0281     return out;
0282   } else if (signalType == "breitWigner") {
0283     // Fit a Breit-Wigner
0284     if ((mean_ == nullptr) || (gamma_ == nullptr)) {
0285       edm::LogError("FitWithRooFit")
0286           << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma";
0287       exit(1);
0288     }
0289     return std::make_unique<RooBreitWigner>("breiWign", "breitWign", *x, *mean_, *gamma_);
0290   } else if (signalType == "relBreitWigner") {
0291     // Fit a relativistic Breit-Wigner
0292     if ((mean_ == nullptr) || (gamma_ == nullptr)) {
0293       edm::LogError("FitWithRooFit")
0294           << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma";
0295       exit(1);
0296     }
0297     return std::make_unique<RooGenericPdf>("Relativistic Breit-Wigner",
0298                                            "RBW",
0299                                            "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",
0300                                            RooArgList{*x, *mean_, *gamma_});
0301   } else if (signalType == "voigtian") {
0302     // Fit a Voigtian
0303     if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr)) {
0304       edm::LogError("FitWithRooFit")
0305           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and gamma";
0306       exit(1);
0307     }
0308     return std::make_unique<RooVoigtian>("voigt", "voigt", *x, *mean_, *gamma_, *sigma_);
0309   } else if (signalType == "crystalBall") {
0310     // Fit a CrystalBall
0311     if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr)) {
0312       edm::LogError("FitWithRooFit")
0313           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, "
0314              "alpha and n";
0315       exit(1);
0316     }
0317     return std::make_unique<RooCBShape>("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_);
0318   } else if (signalType == "breitWignerTimesCB") {
0319     // Fit a Breit Wigner convoluted with a CrystalBall
0320     if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) ||
0321         (alpha_ == nullptr) || (n_ == nullptr)) {
0322       edm::LogError("FitWithRooFit")
0323           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0324              "sigma, gamma, alpha and n";
0325       exit(1);
0326     }
0327     RooAbsPdf* bw = new RooBreitWigner("breiWigner", "breitWigner", *x, *mean_, *gamma_);
0328     RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0329     auto out = std::make_unique<RooFFTConvPdf>("breitWignerTimesCB", "breitWignerTimesCB", *x, *bw, *cb);
0330     out->addOwnedComponents({*bw, *cb});
0331     return out;
0332   } else if (signalType == "relBreitWignerTimesCB") {
0333     // Fit a relativistic Breit Wigner convoluted with a CrystalBall
0334     if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) ||
0335         (alpha_ == nullptr) || (n_ == nullptr)) {
0336       edm::LogError("FitWithRooFit")
0337           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0338              "sigma, gamma, alpha and n";
0339       exit(1);
0340     }
0341     RooGenericPdf* bw = new RooGenericPdf("Relativistic Breit-Wigner",
0342                                           "RBW",
0343                                           "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",
0344                                           {*x, *mean_, *gamma_});
0345     RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0346     auto out = std::make_unique<RooFFTConvPdf>("relBreitWignerTimesCB", "relBreitWignerTimesCB", *x, *bw, *cb);
0347     out->addOwnedComponents({*bw, *cb});
0348     return out;
0349   } else if (signalType == "gaussianPlusCrystalBall") {
0350     // Fit a Gaussian + CrystalBall with the same mean
0351     if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) ||
0352         (fGCB_ == nullptr)) {
0353       edm::LogError("FitWithRooFit")
0354           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, "
0355              "sigma2, alpha, n and fGCB";
0356       exit(1);
0357     }
0358     RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_);
0359     RooAbsPdf* tempGaussian = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma2_);
0360     RooArgList components{*tempCB, *tempGaussian};
0361 
0362     auto out = std::make_unique<RooAddPdf>("gaussianPlusCrystalBall", "gaussianPlusCrystalBall", components, *fGCB_);
0363     out->addOwnedComponents(components);
0364     return out;
0365   } else if (signalType == "voigtianPlusCrystalBall") {
0366     // Fit a Voigtian + CrystalBall with the same mean
0367     if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) ||
0368         (sigma2_ == nullptr) || (fGCB_ == nullptr)) {
0369       edm::LogError("FitWithRooFit")
0370           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0371              "sigma, sigma2, alpha, n and fGCB";
0372       exit(1);
0373     }
0374     RooAbsPdf* tempVoigt = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_);
0375     RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0376     RooArgList components{*tempVoigt, *tempCB};
0377 
0378     auto out = std::make_unique<RooAddPdf>("voigtianPlusCrystalBall", "voigtianPlusCrystalBall", components, *fGCB_);
0379     out->addOwnedComponents(components);
0380     return out;
0381   } else if (signalType == "breitWignerPlusCrystalBall") {
0382     // Fit a Breit-Wigner + CrystalBall with the same mean
0383     if ((mean_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) ||
0384         (fGCB_ == nullptr)) {
0385       edm::LogError("FitWithRooFit")
0386           << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0387              "sigma, alpha, n and fGCB";
0388       exit(1);
0389     }
0390     RooAbsPdf* tempBW = new RooBreitWigner("breitWign", "breitWign", *x, *mean_, *gamma_);
0391     RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0392     RooArgList components{*tempCB, *tempBW};
0393 
0394     auto out =
0395         std::make_unique<RooAddPdf>("breitWignerPlusCrystalBall", "breitWignerPlusCrystalBall", components, *fGCB_);
0396     out->addOwnedComponents(components);
0397     return out;
0398   }
0399 
0400   else if (signalType != "") {
0401     edm::LogError("FitWithRooFit") << "Unknown signal function: " << signalType << ". Signal will not be in the model";
0402   }
0403   return nullptr;
0404 }
0405 
0406 /// Build the model for the specified background type
0407 std::unique_ptr<RooAbsPdf> FitWithRooFit::buildBackgroundModel(RooRealVar* x, const TString& backgroundType) {
0408   if (backgroundType == "exponential") {
0409     // Add an exponential for the background
0410     if ((expCoeffa1_ == nullptr) || (fsig_ == nullptr)) {
0411       edm::LogError("FitWithRooFit")
0412           << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeffa1 and fsig";
0413       exit(1);
0414     }
0415     return std::make_unique<RooExponential>("exponential", "exponential", *x, *expCoeffa1_);
0416   }
0417 
0418   if (backgroundType == "exponentialpol") {
0419     // Add an exponential for the background
0420     if ((expCoeffa0_ == nullptr) || (expCoeffa1_ == nullptr) || (expCoeffa2_ == nullptr) || (fsig_ == nullptr)) {
0421       edm::LogError("FitWithRooFit")
0422           << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeff and fsig";
0423       exit(1);
0424     }
0425     return std::make_unique<RooGenericPdf>("exponential",
0426                                            "exponential",
0427                                            "TMath::Exp(@1+@2*@0+@3*@0*@0)",
0428                                            RooArgList{*x, *expCoeffa0_, *expCoeffa1_, *expCoeffa2_});
0429   }
0430 
0431   else if (backgroundType == "chebychev0") {
0432     // Add a linear background
0433     if (a0_ == nullptr) {
0434       edm::LogError("FitWithRooFit")
0435           << "Error: one or more parameters are not initialized. Please be sure to initialize a0";
0436       exit(1);
0437     }
0438     return std::make_unique<RooChebychev>("chebychev0", "chebychev0", *x, *a0_);
0439   } else if (backgroundType == "chebychev1") {
0440     // Add a 2nd order chebychev polynomial background
0441     if ((a0_ == nullptr) || (a1_ == nullptr)) {
0442       edm::LogError("FitWithRooFit")
0443           << "Error: one or more parameters are not initialized. Please be sure to initialize a0 and a1";
0444       exit(1);
0445     }
0446     return std::make_unique<RooChebychev>("chebychev1", "chebychev1", *x, RooArgList{*a0_, *a1_});
0447   } else if (backgroundType == "chebychev3") {
0448     // Add a 3rd order chebychev polynomial background
0449     if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr)) {
0450       edm::LogError("FitWithRooFit")
0451           << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2 and a3";
0452       exit(1);
0453     }
0454     return std::make_unique<RooChebychev>("3rdOrderPol", "3rdOrderPol", *x, RooArgList{*a0_, *a1_, *a2_, *a3_});
0455   }
0456 
0457   else if (backgroundType == "chebychev6") {
0458     // Add a 6th order chebychev polynomial background
0459     if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr) || (a4_ == nullptr) ||
0460         (a5_ == nullptr) || (a6_ == nullptr)) {
0461       edm::LogError("FitWithRooFit")
0462           << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2, a3, "
0463              "a4, a5 and a6";
0464       exit(1);
0465     }
0466     return std::make_unique<RooChebychev>(
0467         "6thOrderPol", "6thOrderPol", *x, RooArgList{*a0_, *a1_, *a2_, *a3_, *a4_, *a5_, *a6_});
0468   }
0469   return nullptr;
0470 }
0471 
0472 /// Build the model to fit
0473 std::unique_ptr<RooAbsPdf> FitWithRooFit::buildModel(RooRealVar* x,
0474                                                      const TString& signalType,
0475                                                      const TString& backgroundType) {
0476   std::unique_ptr<RooAbsPdf> model;
0477 
0478   std::unique_ptr<RooAbsPdf> signal = buildSignalModel(x, signalType);
0479   std::unique_ptr<RooAbsPdf> background = buildBackgroundModel(x, backgroundType);
0480 
0481   if ((signal != nullptr) && (background != nullptr)) {
0482     // Combine signal and background pdfs
0483     edm::LogPrint("FitWithRooFit") << "Building model with signal and backgound";
0484     RooArgList components{*signal.release(), *background.release()};
0485     model = std::make_unique<RooAddPdf>("model", "model", components, *fsig_);
0486     model->addOwnedComponents(components);
0487   } else if (signal != nullptr) {
0488     edm::LogPrint("FitWithRooFit") << "Building model with signal";
0489     model = std::move(signal);
0490   } else if (background != nullptr) {
0491     edm::LogPrint("FitWithRooFit") << "Building model with backgound";
0492     model = std::move(background);
0493   } else {
0494     edm::LogWarning("FitWithRooFit") << "Nothing to fit";
0495     exit(0);
0496   }
0497   return model;
0498 }