OfflineValidation/src/FitWithRooFit.cc

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