From c55d87cf8e48fd98214579c4239d2a3700ecedde Mon Sep 17 00:00:00 2001 From: mmusich Date: Tue, 20 Jun 2023 13:31:27 +0200 Subject: [PATCH 1/2] move FitWitRooFit out of bin and create a class with a public header --- Alignment/OfflineValidation/BuildFile.xml | 1 + .../OfflineValidation/bin/FitWithRooFit.cc | 776 ------------------ Alignment/OfflineValidation/bin/Zmumumerge.cc | 3 +- .../interface/FitWithRooFit.h | 307 +++++++ .../OfflineValidation/src/FitWithRooFit.cc | 578 +++++++++++++ 5 files changed, 888 insertions(+), 777 deletions(-) delete mode 100644 Alignment/OfflineValidation/bin/FitWithRooFit.cc create mode 100644 Alignment/OfflineValidation/interface/FitWithRooFit.h create mode 100644 Alignment/OfflineValidation/src/FitWithRooFit.cc diff --git a/Alignment/OfflineValidation/BuildFile.xml b/Alignment/OfflineValidation/BuildFile.xml index 5f08c1b9948c1..8ed570cb4abba 100644 --- a/Alignment/OfflineValidation/BuildFile.xml +++ b/Alignment/OfflineValidation/BuildFile.xml @@ -26,6 +26,7 @@ + diff --git a/Alignment/OfflineValidation/bin/FitWithRooFit.cc b/Alignment/OfflineValidation/bin/FitWithRooFit.cc deleted file mode 100644 index 4ebb9a47d1e1d..0000000000000 --- a/Alignment/OfflineValidation/bin/FitWithRooFit.cc +++ /dev/null @@ -1,776 +0,0 @@ -#ifndef FitWithRooFit_cc -#define FitWithRooFit_cc - -#ifndef __CINT__ -#include "RooGlobalFunc.h" -#endif -#include "TCanvas.h" -#include "TTree.h" -#include "TH1D.h" -#include "TRandom.h" -#include "RooRealVar.h" -#include "RooDataSet.h" -#include "RooGaussian.h" -#include "RooVoigtian.h" -#include "RooExponential.h" -#include "RooPlot.h" -#include "RooDataHist.h" -#include "RooAddPdf.h" -#include "RooChebychev.h" -#include "RooGenericPdf.h" -#include "RooGaussModel.h" -#include "RooAddModel.h" -#include "RooPolynomial.h" -#include "RooCBShape.h" -#include "RooMinimizer.h" -#include "RooBreitWigner.h" -#include "RooFFTConvPdf.h" - -/** - * This macro allows to use RooFit to perform a fit on a TH1 histogram.
- * The currently implemented functions are:
- * - signal:
- * -- gaussian
- * -- double gaussian
- * -- voigtian
- * - background:
- * -- exponential
- * It is possible to select any combination of signal and background.
- * The fit() method receives the TH1, two strings specifying the signal and background type - * and the min and max x of the histogram over which to perform the fit.
- * The variables of the fit must be initialized separately. For example when doing a gaussian+exponential - * fit the initMean, initSigma, initFSig and initExpCoeff methods must be used to create and initialize - * the corresponding variables.
- * The methods names after the variables return the fit result. - */ - -namespace { - typedef std::pair rooPair; -} - -class FitWithRooFit { -public: - FitWithRooFit() - : useChi2_(false), - mean_(nullptr), - mean2_(nullptr), - mean3_(nullptr), - sigma_(nullptr), - sigma2_(nullptr), - sigma3_(nullptr), - gamma_(nullptr), - gaussFrac_(nullptr), - gaussFrac2_(nullptr), - expCoeffa0_(nullptr), - expCoeffa1_(nullptr), - expCoeffa2_(nullptr), - fsig_(nullptr), - a0_(nullptr), - a1_(nullptr), - a2_(nullptr), - a3_(nullptr), - a4_(nullptr), - a5_(nullptr), - a6_(nullptr), - alpha_(nullptr), - n_(nullptr), - fGCB_(nullptr) {} - - // Import TH1 histogram into a RooDataHist - rooPair importTH1(TH1* histo, const double& inputXmin, const double& inputXmax) { - double xMin = inputXmin; - double xMax = inputXmax; - if ((xMin == xMax) && xMin == 0) { - xMin = histo->GetXaxis()->GetXmin(); - xMax = histo->GetXaxis()->GetXmax(); - } - // Declare observable x - RooRealVar x("x", "x", xMin, xMax); - // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x' - return (std::make_pair(x, new RooDataHist("dh", "dh", x, RooFit::Import(*histo)))); - } - - // Plot and fit a RooDataHist fitting signal and background - void fit(TH1* histo, - const TString signalType, - const TString backgroundType, - const double& xMin = 0., - const double& xMax = 0., - bool sumW2Error = false) { - reinitializeParameters(); - - rooPair imported = importTH1(histo, xMin, xMax); - RooRealVar x(imported.first); - RooDataHist* dh = imported.second; - - // Make plot of binned dataset showing Poisson error bars (RooFit default) - RooPlot* frame = x.frame(RooFit::Title("Imported TH1 with Poisson error bars")); - frame->SetName(TString(histo->GetName()) + "_frame"); - dh->plotOn(frame); - - // Build the composite model - RooAbsPdf* model = buildModel(&x, signalType, backgroundType); - - std::unique_ptr chi2{model->createChi2(*dh, RooFit::DataError(RooAbsData::SumW2))}; - - // Fit the composite model - // ----------------------- - // Fit with likelihood - if (!useChi2_) { - if (sumW2Error) - model->fitTo(*dh, RooFit::Save(), RooFit::SumW2Error(kTRUE)); - else - model->fitTo(*dh); - } - // Fit with chi^2 - else { - std::cout << "FITTING WITH CHI^2" << std::endl; - RooMinimizer m(*chi2); - m.migrad(); - m.hesse(); - // RooFitResult* r_chi2_wgt = m.save(); - } - model->plotOn(frame); - model->plotOn(frame, RooFit::Components(backgroundType), RooFit::LineStyle(kDotted), RooFit::LineColor(kRed)); - model->paramOn(frame, RooFit::Label("fit result"), RooFit::Format("NEU", RooFit::AutoPrecision(2))); - - // TODO: fix next lines to get the prob(chi2) (ndof should be dynamically set according to the choosen pdf) - // double chi2 = xframe.chiSquare("model","data",ndof); - // double ndoff = xframeGetNbinsX(); - // double chi2prob = TMath::Prob(chi2,ndoff); - - // 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 - // --------------------------------------------------------------------------------------------- - - // If histogram has custom error (i.e. its contents is does not originate from a Poisson process - // but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead - // (same error bars as shown by ROOT) - RooPlot* frame2 = x.frame(RooFit::Title("Imported TH1 with internal errors")); - dh->plotOn(frame2, RooFit::DataError(RooAbsData::SumW2)); - model->plotOn(frame2); - model->plotOn(frame2, RooFit::Components(backgroundType), RooFit::LineColor(kRed)); - model->paramOn(frame2, RooFit::Label("fit result"), RooFit::Format("NEU", RooFit::AutoPrecision(2))); - - // Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used - // in a maximum likelihood fit - // - // A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition - // of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly - // fitted with a chi^2 fit (see rf602_chi2fit.C) - - // Draw all frames on a canvas - // if( canvas == 0 ) { - // canvas = new TCanvas("rf102_dataimport","rf102_dataimport",800,800); - // canvas->cd(1); - // } - // canvas->Divide(2,1); - // canvas->cd(1) ; frame->Draw(); - // canvas->cd(2) ; frame2->Draw(); - - frame2->Draw(); - } - - // Initialization methods for all the parameters - void initMean(const double& value, - const double& min, - const double& max, - const TString& name = "mean", - const TString& title = "mean") { - if (mean_ != nullptr) - delete mean_; - mean_ = new RooRealVar(name, title, value, min, max); - initVal_mean = value; - } - void initMean2(const double& value, - const double& min, - const double& max, - const TString& name = "mean2", - const TString& title = "mean2") { - if (mean2_ != nullptr) - delete mean2_; - mean2_ = new RooRealVar(name, title, value, min, max); - initVal_mean2 = value; - } - void initMean3(const double& value, - const double& min, - const double& max, - const TString& name = "mean3", - const TString& title = "mean3") { - if (mean3_ != nullptr) - delete mean3_; - mean3_ = new RooRealVar(name, title, value, min, max); - initVal_mean3 = value; - } - void initSigma(const double& value, - const double& min, - const double& max, - const TString& name = "sigma", - const TString& title = "sigma") { - if (sigma_ != nullptr) - delete sigma_; - sigma_ = new RooRealVar(name, title, value, min, max); - initVal_sigma = value; - } - void initSigma2(const double& value, - const double& min, - const double& max, - const TString& name = "sigma2", - const TString& title = "sigma2") { - if (sigma2_ != nullptr) - delete sigma2_; - sigma2_ = new RooRealVar(name, title, value, min, max); - initVal_sigma2 = value; - } - void initSigma3(const double& value, - const double& min, - const double& max, - const TString& name = "sigma3", - const TString& title = "sigma3") { - if (sigma3_ != nullptr) - delete sigma3_; - sigma3_ = new RooRealVar(name, title, value, min, max); - initVal_sigma3 = value; - } - void initGamma(const double& value, - const double& min, - const double& max, - const TString& name = "gamma", - const TString& title = "gamma") { - if (gamma_ != nullptr) - delete gamma_; - gamma_ = new RooRealVar(name, title, value, min, max); - initVal_gamma = value; - } - void initGaussFrac(const double& value, - const double& min, - const double& max, - const TString& name = "GaussFrac", - const TString& title = "GaussFrac") { - if (gaussFrac_ != nullptr) - delete gaussFrac_; - gaussFrac_ = new RooRealVar(name, title, value, min, max); - initVal_gaussFrac = value; - } - void initGaussFrac2(const double& value, - const double& min, - const double& max, - const TString& name = "GaussFrac2", - const TString& title = "GaussFrac2") { - if (gaussFrac2_ != nullptr) - delete gaussFrac2_; - gaussFrac2_ = new RooRealVar(name, title, value, min, max); - initVal_gaussFrac2 = value; - } - void initExpCoeffA0(const double& value, - const double& min, - const double& max, - const TString& name = "expCoeffa0", - const TString& title = "expCoeffa0") { - if (expCoeffa0_ != nullptr) - delete expCoeffa0_; - expCoeffa0_ = new RooRealVar(name, title, value, min, max); - initVal_expCoeffa0 = value; - } - void initExpCoeffA1(const double& value, - const double& min, - const double& max, - const TString& name = "expCoeffa1", - const TString& title = "expCoeffa1") { - if (expCoeffa1_ != nullptr) - delete expCoeffa1_; - expCoeffa1_ = new RooRealVar(name, title, value, min, max); - initVal_expCoeffa1 = value; - } - void initExpCoeffA2(const double& value, - const double& min, - const double& max, - const TString& name = "expCoeffa2", - const TString& title = "expCoeffa2") { - if (expCoeffa2_ != nullptr) - delete expCoeffa2_; - expCoeffa2_ = new RooRealVar(name, title, value, min, max); - initVal_expCoeffa2 = value; - } - void initFsig(const double& value, - const double& min, - const double& max, - const TString& name = "fsig", - const TString& title = "signal fraction") { - if (fsig_ != nullptr) - delete fsig_; - fsig_ = new RooRealVar(name, title, value, min, max); - initVal_fsig = value; - } - void initA0(const double& value, - const double& min, - const double& max, - const TString& name = "a0", - const TString& title = "a0") { - if (a0_ != nullptr) - delete a0_; - a0_ = new RooRealVar(name, title, value, min, max); - initVal_a0 = value; - } - void initA1(const double& value, - const double& min, - const double& max, - const TString& name = "a1", - const TString& title = "a1") { - if (a1_ != nullptr) - delete a1_; - a1_ = new RooRealVar(name, title, value, min, max); - initVal_a1 = value; - } - void initA2(const double& value, - const double& min, - const double& max, - const TString& name = "a2", - const TString& title = "a2") { - if (a2_ != nullptr) - delete a2_; - a2_ = new RooRealVar(name, title, value, min, max); - initVal_a2 = value; - } - void initA3(const double& value, - const double& min, - const double& max, - const TString& name = "a3", - const TString& title = "a3") { - if (a3_ != nullptr) - delete a3_; - a3_ = new RooRealVar(name, title, value, min, max); - initVal_a3 = value; - } - void initA4(const double& value, - const double& min, - const double& max, - const TString& name = "a4", - const TString& title = "a4") { - if (a4_ != nullptr) - delete a4_; - a4_ = new RooRealVar(name, title, value, min, max); - initVal_a4 = value; - } - void initA5(const double& value, - const double& min, - const double& max, - const TString& name = "a5", - const TString& title = "a5") { - if (a5_ != nullptr) - delete a5_; - a5_ = new RooRealVar(name, title, value, min, max); - initVal_a5 = value; - } - void initA6(const double& value, - const double& min, - const double& max, - const TString& name = "a6", - const TString& title = "a6") { - if (a6_ != nullptr) - delete a6_; - a6_ = new RooRealVar(name, title, value, min, max); - initVal_a6 = value; - } - void initAlpha(const double& value, - const double& min, - const double& max, - const TString& name = "alpha", - const TString& title = "alpha") { - if (alpha_ != nullptr) - delete alpha_; - alpha_ = new RooRealVar(name, title, value, min, max); - initVal_alpha = value; - } - void initN(const double& value, - const double& min, - const double& max, - const TString& name = "n", - const TString& title = "n") { - if (n_ != nullptr) - delete n_; - n_ = new RooRealVar(name, title, value, min, max); - initVal_n = value; - } - void initFGCB(const double& value, - const double& min, - const double& max, - const TString& name = "fGCB", - const TString& title = "fGCB") { - if (fGCB_ != nullptr) - delete fGCB_; - fGCB_ = new RooRealVar(name, title, value, min, max); - initVal_fGCB = value; - } - - void reinitializeParameters() { - if (mean_ != nullptr) - mean_->setVal(initVal_mean); - if (mean2_ != nullptr) - mean2_->setVal(initVal_mean2); - if (mean3_ != nullptr) - mean3_->setVal(initVal_mean3); - if (sigma_ != nullptr) - sigma_->setVal(initVal_sigma); - if (sigma2_ != nullptr) - sigma2_->setVal(initVal_sigma2); - if (sigma3_ != nullptr) - sigma3_->setVal(initVal_sigma3); - if (gamma_ != nullptr) - gamma_->setVal(initVal_gamma); - if (gaussFrac_ != nullptr) - gaussFrac_->setVal(initVal_gaussFrac); - if (gaussFrac2_ != nullptr) - gaussFrac2_->setVal(initVal_gaussFrac2); - if (expCoeffa0_ != nullptr) - expCoeffa0_->setVal(initVal_expCoeffa0); - if (expCoeffa1_ != nullptr) - expCoeffa1_->setVal(initVal_expCoeffa1); - if (expCoeffa2_ != nullptr) - expCoeffa2_->setVal(initVal_expCoeffa2); - if (fsig_ != nullptr) - fsig_->setVal(initVal_fsig); - if (a0_ != nullptr) - a0_->setVal(initVal_a0); - if (a1_ != nullptr) - a1_->setVal(initVal_a1); - if (a2_ != nullptr) - a2_->setVal(initVal_a2); - if (a3_ != nullptr) - a3_->setVal(initVal_a3); - if (a4_ != nullptr) - a4_->setVal(initVal_a4); - if (a5_ != nullptr) - a5_->setVal(initVal_a5); - if (a6_ != nullptr) - a6_->setVal(initVal_a6); - if (alpha_ != nullptr) - alpha_->setVal(initVal_alpha); - if (n_ != nullptr) - n_->setVal(initVal_n); - if (fGCB_ != nullptr) - fGCB_->setVal(initVal_fGCB); - } - - inline RooRealVar* mean() { return mean_; } - inline RooRealVar* mean2() { return mean2_; } - inline RooRealVar* mean3() { return mean3_; } - inline RooRealVar* sigma() { return sigma_; } - inline RooRealVar* sigma2() { return sigma2_; } - inline RooRealVar* sigma3() { return sigma3_; } - inline RooRealVar* gamma() { return gamma_; } - inline RooRealVar* gaussFrac() { return gaussFrac_; } - inline RooRealVar* gaussFrac2() { return gaussFrac2_; } - inline RooRealVar* expCoeffa0() { return expCoeffa0_; } - inline RooRealVar* expCoeffa1() { return expCoeffa1_; } - inline RooRealVar* expCoeffa2() { return expCoeffa2_; } - inline RooRealVar* fsig() { return fsig_; } - inline RooRealVar* a0() { return a0_; } - inline RooRealVar* a1() { return a1_; } - inline RooRealVar* a2() { return a2_; } - inline RooRealVar* a3() { return a3_; } - inline RooRealVar* a4() { return a4_; } - inline RooRealVar* a5() { return a5_; } - inline RooRealVar* a6() { return a6_; } - inline RooRealVar* alpha() { return alpha_; } - inline RooRealVar* n() { return n_; } - inline RooRealVar* fGCB() { return fGCB_; } - - /// Build the model for the specified signal type - RooAbsPdf* buildSignalModel(RooRealVar* x, const TString& signalType) { - RooAbsPdf* signal = nullptr; - if (signalType == "gaussian") { - // Fit a Gaussian p.d.f to the data - if ((mean_ == nullptr) || (sigma_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and sigma" - << std::endl; - exit(1); - } - signal = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma_); - } else if (signalType == "doubleGaussian") { - // Fit with double gaussian - if ((mean_ == nullptr) || (sigma_ == nullptr) || (sigma2_ == nullptr)) { - std::cout - << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and sigma2" - << std::endl; - exit(1); - } - RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_); - RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean_, *sigma2_); - signal = new RooAddModel("doubleGaussian", "double gaussian", RooArgList(*gaussModel, *gaussModel2), *gaussFrac_); - } else if (signalType == "tripleGaussian") { - // Fit with triple gaussian - if ((mean_ == nullptr) || (mean2_ == nullptr) || (mean3_ == nullptr) || (sigma_ == nullptr) || - (sigma2_ == nullptr) || (sigma3_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, " - "mean3, sigma, sigma2, sigma3" - << std::endl; - exit(1); - } - RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_); - RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean2_, *sigma2_); - RooGaussModel* gaussModel3 = new RooGaussModel("gaussModel3", "gaussModel3", *x, *mean3_, *sigma3_); - signal = new RooAddModel("tripleGaussian", - "triple gaussian", - RooArgList(*gaussModel, *gaussModel2, *gaussModel3), - RooArgList(*gaussFrac_, *gaussFrac2_)); - } else if (signalType == "breitWigner") { - // Fit a Breit-Wigner - if ((mean_ == nullptr) || (gamma_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma" - << std::endl; - exit(1); - } - signal = new RooBreitWigner("breiWign", "breitWign", *x, *mean_, *gamma_); - } else if (signalType == "relBreitWigner") { - // Fit a relativistic Breit-Wigner - if ((mean_ == nullptr) || (gamma_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma" - << std::endl; - exit(1); - } - signal = new RooGenericPdf("Relativistic Breit-Wigner", - "RBW", - "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))", - RooArgList(*x, *mean_, *gamma_)); - } else if (signalType == "voigtian") { - // Fit a Voigtian - if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr)) { - std::cout - << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and gamma" - << std::endl; - exit(1); - } - signal = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_); - } else if (signalType == "crystalBall") { - // Fit a CrystalBall - if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, " - "alpha and n" - << std::endl; - exit(1); - } - signal = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_); - } else if (signalType == "breitWignerTimesCB") { - // Fit a Breit Wigner convoluted with a CrystalBall - if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || - (alpha_ == nullptr) || (n_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, " - "sigma, gamma, alpha and n" - << std::endl; - exit(1); - } - RooAbsPdf* bw = new RooBreitWigner("breiWigner", "breitWigner", *x, *mean_, *gamma_); - RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_); - signal = new RooFFTConvPdf("breitWignerTimesCB", "breitWignerTimesCB", *x, *bw, *cb); - } else if (signalType == "relBreitWignerTimesCB") { - // Fit a relativistic Breit Wigner convoluted with a CrystalBall - if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || - (alpha_ == nullptr) || (n_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, " - "sigma, gamma, alpha and n" - << std::endl; - exit(1); - } - RooGenericPdf* bw = new RooGenericPdf("Relativistic Breit-Wigner", - "RBW", - "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))", - RooArgList(*x, *mean_, *gamma_)); - RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_); - signal = new RooFFTConvPdf("relBreitWignerTimesCB", "relBreitWignerTimesCB", *x, *bw, *cb); - } else if (signalType == "gaussianPlusCrystalBall") { - // Fit a Gaussian + CrystalBall with the same mean - if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) || - (fGCB_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, " - "sigma2, alpha, n and fGCB" - << std::endl; - exit(1); - } - RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_); - RooAbsPdf* tempGaussian = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma2_); - - signal = new RooAddPdf( - "gaussianPlusCrystalBall", "gaussianPlusCrystalBall", RooArgList(*tempCB, *tempGaussian), *fGCB_); - } else if (signalType == "voigtianPlusCrystalBall") { - // Fit a Voigtian + CrystalBall with the same mean - if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || - (sigma2_ == nullptr) || (fGCB_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, " - "sigma, sigma2, alpha, n and fGCB" - << std::endl; - exit(1); - } - RooAbsPdf* tempVoigt = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_); - RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_); - - signal = - new RooAddPdf("voigtianPlusCrystalBall", "voigtianPlusCrystalBall", RooArgList(*tempCB, *tempVoigt), *fGCB_); - } else if (signalType == "breitWignerPlusCrystalBall") { - // Fit a Breit-Wigner + CrystalBall with the same mean - if ((mean_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) || - (fGCB_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, " - "sigma, alpha, n and fGCB" - << std::endl; - exit(1); - } - RooAbsPdf* tempBW = new RooBreitWigner("breitWign", "breitWign", *x, *mean_, *gamma_); - RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_); - - signal = new RooAddPdf( - "breitWignerPlusCrystalBall", "breitWignerPlusCrystalBall", RooArgList(*tempCB, *tempBW), *fGCB_); - } - - else if (signalType != "") { - std::cout << "Unknown signal function: " << signalType << ". Signal will not be in the model" << std::endl; - } - return signal; - } - - /// Build the model for the specified background type - RooAbsPdf* buildBackgroundModel(RooRealVar* x, const TString& backgroundType) { - RooAbsPdf* background = nullptr; - if (backgroundType == "exponential") { - // Add an exponential for the background - if ((expCoeffa1_ == nullptr) || (fsig_ == nullptr)) { - std::cout - << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeffa1 and fsig" - << std::endl; - exit(1); - } - background = new RooExponential("exponential", "exponential", *x, *expCoeffa1_); - } - - if (backgroundType == "exponentialpol") { - // Add an exponential for the background - if ((expCoeffa0_ == nullptr) || (expCoeffa1_ == nullptr) || (expCoeffa2_ == nullptr) || (fsig_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeff and fsig" - << std::endl; - exit(1); - } - background = new RooGenericPdf("exponential", - "exponential", - "TMath::Exp(@1+@2*@0+@3*@0*@0)", - RooArgList(*x, *expCoeffa0_, *expCoeffa1_, *expCoeffa2_)); - } - - else if (backgroundType == "chebychev0") { - // Add a linear background - if (a0_ == nullptr) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0" << std::endl; - exit(1); - } - background = new RooChebychev("chebychev0", "chebychev0", *x, *a0_); - } else if (backgroundType == "chebychev1") { - // Add a 2nd order chebychev polynomial background - if ((a0_ == nullptr) || (a1_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0 and a1" - << std::endl; - exit(1); - } - background = new RooChebychev("chebychev1", "chebychev1", *x, RooArgList(*a0_, *a1_)); - } else if (backgroundType == "chebychev3") { - // Add a 3rd order chebychev polynomial background - if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2 and a3" - << std::endl; - exit(1); - } - background = new RooChebychev("3rdOrderPol", "3rdOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_)); - } - - else if (backgroundType == "chebychev6") { - // Add a 6th order chebychev polynomial background - if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr) || (a4_ == nullptr) || - (a5_ == nullptr) || (a6_ == nullptr)) { - std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2, a3, " - "a4, a5 and a6" - << std::endl; - exit(1); - } - background = - new RooChebychev("6thOrderPol", "6thOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_, *a4_, *a5_, *a6_)); - } - - return background; - } - - /// Build the model to fit - RooAbsPdf* buildModel(RooRealVar* x, const TString& signalType, const TString& backgroundType) { - RooAbsPdf* model = nullptr; - - RooAbsPdf* signal = buildSignalModel(x, signalType); - RooAbsPdf* background = buildBackgroundModel(x, backgroundType); - - if ((signal != nullptr) && (background != nullptr)) { - // Combine signal and background pdfs - std::cout << "Building model with signal and backgound" << std::endl; - model = new RooAddPdf("model", "model", RooArgList(*signal, *background), *fsig_); - } else if (signal != nullptr) { - std::cout << "Building model with signal" << std::endl; - model = signal; - } else if (background != nullptr) { - std::cout << "Building model with backgound" << std::endl; - model = background; - } else { - std::cout << "Nothing to fit" << std::endl; - exit(0); - } - return model; - } - - bool useChi2_; - -protected: - // Declare all variables - RooRealVar* mean_; - RooRealVar* mean2_; - RooRealVar* mean3_; - RooRealVar* sigma_; - RooRealVar* sigma2_; - RooRealVar* sigma3_; - RooRealVar* gamma_; - RooRealVar* gaussFrac_; - RooRealVar* gaussFrac2_; - RooRealVar* expCoeffa0_; - RooRealVar* expCoeffa1_; - RooRealVar* expCoeffa2_; - RooRealVar* fsig_; - RooRealVar* a0_; - RooRealVar* a1_; - RooRealVar* a2_; - RooRealVar* a3_; - RooRealVar* a4_; - RooRealVar* a5_; - RooRealVar* a6_; - RooRealVar* alpha_; - RooRealVar* n_; - RooRealVar* fGCB_; - - // Initial values - double initVal_mean; - double initVal_mean2; - double initVal_mean3; - double initVal_sigma; - double initVal_sigma2; - double initVal_sigma3; - double initVal_gamma; - double initVal_gaussFrac; - double initVal_gaussFrac2; - double initVal_expCoeffa0; - double initVal_expCoeffa1; - double initVal_expCoeffa2; - double initVal_fsig; - double initVal_a0; - double initVal_a1; - double initVal_a2; - double initVal_a3; - double initVal_a4; - double initVal_a5; - double initVal_a6; - double initVal_alpha; - double initVal_n; - double initVal_fGCB; -}; - -#endif diff --git a/Alignment/OfflineValidation/bin/Zmumumerge.cc b/Alignment/OfflineValidation/bin/Zmumumerge.cc index cc8479f647a4a..03eb44f9ca9a5 100644 --- a/Alignment/OfflineValidation/bin/Zmumumerge.cc +++ b/Alignment/OfflineValidation/bin/Zmumumerge.cc @@ -5,7 +5,6 @@ #include #include "toolbox.h" -#include "FitWithRooFit.cc" #include "Options.h" #include "TAxis.h" @@ -39,6 +38,8 @@ #include #include +#include "Alignment/OfflineValidation/interface/FitWithRooFit.h" + using namespace RooFit; using namespace std; using namespace AllInOneConfig; diff --git a/Alignment/OfflineValidation/interface/FitWithRooFit.h b/Alignment/OfflineValidation/interface/FitWithRooFit.h new file mode 100644 index 0000000000000..015b2284879ca --- /dev/null +++ b/Alignment/OfflineValidation/interface/FitWithRooFit.h @@ -0,0 +1,307 @@ +#ifndef Alignment_OfflineValidation_FitWithRooFit_h +#define Alignment_OfflineValidation_FitWithRooFit_h + +#ifndef __CINT__ +#include "RooGlobalFunc.h" +#endif +#include "TCanvas.h" +#include "TTree.h" +#include "TH1D.h" +#include "TRandom.h" +#include "RooRealVar.h" +#include "RooDataSet.h" +#include "RooGaussian.h" +#include "RooVoigtian.h" +#include "RooExponential.h" +#include "RooPlot.h" +#include "RooDataHist.h" +#include "RooAddPdf.h" +#include "RooChebychev.h" +#include "RooGenericPdf.h" +#include "RooGaussModel.h" +#include "RooAddModel.h" +#include "RooPolynomial.h" +#include "RooCBShape.h" +#include "RooMinimizer.h" +#include "RooBreitWigner.h" +#include "RooFFTConvPdf.h" + +/** + * This class allows to use RooFit to perform a fit on a TH1 histogram.
+ * The currently implemented functions are:
+ * - signal:
+ * -- gaussian
+ * -- double gaussian
+ * -- voigtian
+ * - background:
+ * -- exponential
+ * It is possible to select any combination of signal and background.
+ * The fit() method receives the TH1, two strings specifying the signal and background type + * and the min and max x of the histogram over which to perform the fit.
+ * The variables of the fit must be initialized separately. For example when doing a gaussian+exponential + * fit the initMean, initSigma, initFSig and initExpCoeff methods must be used to create and initialize + * the corresponding variables.
+ * The methods names after the variables return the fit result. + */ + +namespace { + typedef std::pair rooPair; +} + +class FitWithRooFit { +public: + FitWithRooFit() + : useChi2_(false), + mean_(nullptr), + mean2_(nullptr), + mean3_(nullptr), + sigma_(nullptr), + sigma2_(nullptr), + sigma3_(nullptr), + gamma_(nullptr), + gaussFrac_(nullptr), + gaussFrac2_(nullptr), + expCoeffa0_(nullptr), + expCoeffa1_(nullptr), + expCoeffa2_(nullptr), + fsig_(nullptr), + a0_(nullptr), + a1_(nullptr), + a2_(nullptr), + a3_(nullptr), + a4_(nullptr), + a5_(nullptr), + a6_(nullptr), + alpha_(nullptr), + n_(nullptr), + fGCB_(nullptr) {} + + rooPair importTH1(TH1* histo, const double& inputXmin, const double& inputXmax); + + void fit(TH1* histo, + const TString signalType, + const TString backgroundType, + const double& xMin = 0., + const double& xMax = 0., + bool sumW2Error = false); + + void initMean(const double& value, + const double& min, + const double& max, + const TString& name = "mean", + const TString& title = "mean"); + + void initMean2(const double& value, + const double& min, + const double& max, + const TString& name = "mean2", + const TString& title = "mean2"); + + void initMean3(const double& value, + const double& min, + const double& max, + const TString& name = "mean3", + const TString& title = "mean3"); + + void initSigma(const double& value, + const double& min, + const double& max, + const TString& name = "sigma", + const TString& title = "sigma"); + + void initSigma2(const double& value, + const double& min, + const double& max, + const TString& name = "sigma2", + const TString& title = "sigma2"); + + void initSigma3(const double& value, + const double& min, + const double& max, + const TString& name = "sigma3", + const TString& title = "sigma3"); + + void initGamma(const double& value, + const double& min, + const double& max, + const TString& name = "gamma", + const TString& title = "gamma"); + + void initGaussFrac(const double& value, + const double& min, + const double& max, + const TString& name = "GaussFrac", + const TString& title = "GaussFrac"); + + void initGaussFrac2(const double& value, + const double& min, + const double& max, + const TString& name = "GaussFrac2", + const TString& title = "GaussFrac2"); + + void initExpCoeffA0(const double& value, + const double& min, + const double& max, + const TString& name = "expCoeffa0", + const TString& title = "expCoeffa0"); + + void initExpCoeffA1(const double& value, + const double& min, + const double& max, + const TString& name = "expCoeffa1", + const TString& title = "expCoeffa1"); + + void initExpCoeffA2(const double& value, + const double& min, + const double& max, + const TString& name = "expCoeffa2", + const TString& title = "expCoeffa2"); + + void initFsig(const double& value, + const double& min, + const double& max, + const TString& name = "fsig", + const TString& title = "signal fraction"); + + void initA0(const double& value, + const double& min, + const double& max, + const TString& name = "a0", + const TString& title = "a0"); + + void initA1(const double& value, + const double& min, + const double& max, + const TString& name = "a1", + const TString& title = "a1"); + + void initA2(const double& value, + const double& min, + const double& max, + const TString& name = "a2", + const TString& title = "a2"); + + void initA3(const double& value, + const double& min, + const double& max, + const TString& name = "a3", + const TString& title = "a3"); + + void initA4(const double& value, + const double& min, + const double& max, + const TString& name = "a4", + const TString& title = "a4"); + + void initA5(const double& value, + const double& min, + const double& max, + const TString& name = "a5", + const TString& title = "a5"); + + void initA6(const double& value, + const double& min, + const double& max, + const TString& name = "a6", + const TString& title = "a6"); + + void initAlpha(const double& value, + const double& min, + const double& max, + const TString& name = "alpha", + const TString& title = "alpha"); + + void initN( + const double& value, const double& min, const double& max, const TString& name = "n", const TString& title = "n"); + + void initFGCB(const double& value, + const double& min, + const double& max, + const TString& name = "fGCB", + const TString& title = "fGCB"); + + inline RooRealVar* mean() { return mean_; } + inline RooRealVar* mean2() { return mean2_; } + inline RooRealVar* mean3() { return mean3_; } + inline RooRealVar* sigma() { return sigma_; } + inline RooRealVar* sigma2() { return sigma2_; } + inline RooRealVar* sigma3() { return sigma3_; } + inline RooRealVar* gamma() { return gamma_; } + inline RooRealVar* gaussFrac() { return gaussFrac_; } + inline RooRealVar* gaussFrac2() { return gaussFrac2_; } + inline RooRealVar* expCoeffa0() { return expCoeffa0_; } + inline RooRealVar* expCoeffa1() { return expCoeffa1_; } + inline RooRealVar* expCoeffa2() { return expCoeffa2_; } + inline RooRealVar* fsig() { return fsig_; } + inline RooRealVar* a0() { return a0_; } + inline RooRealVar* a1() { return a1_; } + inline RooRealVar* a2() { return a2_; } + inline RooRealVar* a3() { return a3_; } + inline RooRealVar* a4() { return a4_; } + inline RooRealVar* a5() { return a5_; } + inline RooRealVar* a6() { return a6_; } + inline RooRealVar* alpha() { return alpha_; } + inline RooRealVar* n() { return n_; } + inline RooRealVar* fGCB() { return fGCB_; } + + void reinitializeParameters(); + + RooAbsPdf* buildSignalModel(RooRealVar* x, const TString& signalType); + RooAbsPdf* buildBackgroundModel(RooRealVar* x, const TString& backgroundType); + RooAbsPdf* buildModel(RooRealVar* x, const TString& signalType, const TString& backgroundType); + + bool useChi2_; + +protected: + // Declare all variables + RooRealVar* mean_; + RooRealVar* mean2_; + RooRealVar* mean3_; + RooRealVar* sigma_; + RooRealVar* sigma2_; + RooRealVar* sigma3_; + RooRealVar* gamma_; + RooRealVar* gaussFrac_; + RooRealVar* gaussFrac2_; + RooRealVar* expCoeffa0_; + RooRealVar* expCoeffa1_; + RooRealVar* expCoeffa2_; + RooRealVar* fsig_; + RooRealVar* a0_; + RooRealVar* a1_; + RooRealVar* a2_; + RooRealVar* a3_; + RooRealVar* a4_; + RooRealVar* a5_; + RooRealVar* a6_; + RooRealVar* alpha_; + RooRealVar* n_; + RooRealVar* fGCB_; + + // Initial values + double initVal_mean; + double initVal_mean2; + double initVal_mean3; + double initVal_sigma; + double initVal_sigma2; + double initVal_sigma3; + double initVal_gamma; + double initVal_gaussFrac; + double initVal_gaussFrac2; + double initVal_expCoeffa0; + double initVal_expCoeffa1; + double initVal_expCoeffa2; + double initVal_fsig; + double initVal_a0; + double initVal_a1; + double initVal_a2; + double initVal_a3; + double initVal_a4; + double initVal_a5; + double initVal_a6; + double initVal_alpha; + double initVal_n; + double initVal_fGCB; +}; + +#endif diff --git a/Alignment/OfflineValidation/src/FitWithRooFit.cc b/Alignment/OfflineValidation/src/FitWithRooFit.cc new file mode 100644 index 0000000000000..a277c98d7fba1 --- /dev/null +++ b/Alignment/OfflineValidation/src/FitWithRooFit.cc @@ -0,0 +1,578 @@ +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "Alignment/OfflineValidation/interface/FitWithRooFit.h" + +// Import TH1 histogram into a RooDataHist +rooPair FitWithRooFit::importTH1(TH1* histo, const double& inputXmin, const double& inputXmax) { + double xMin = inputXmin; + double xMax = inputXmax; + if ((xMin == xMax) && xMin == 0) { + xMin = histo->GetXaxis()->GetXmin(); + xMax = histo->GetXaxis()->GetXmax(); + } + // Declare observable x + RooRealVar x("InvMass", "di-muon mass M(#mu^{+}#mu^{-}) [GeV]", xMin, xMax); + // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x' + return (std::make_pair(x, new RooDataHist("dh", "dh", x, RooFit::Import(*histo)))); +} + +// Plot and fit a RooDataHist fitting signal and background +void FitWithRooFit::fit(TH1* histo, + const TString signalType, + const TString backgroundType, + const double& xMin, + const double& xMax, + bool sumW2Error) { + reinitializeParameters(); + + rooPair imported = importTH1(histo, xMin, xMax); + RooRealVar x(imported.first); + RooDataHist* dh = imported.second; + + // Make plot of binned dataset showing Poisson error bars (RooFit default) + RooPlot* frame = x.frame(RooFit::Title("di-muon mass M(#mu^{+}#mu^{-}) [GeV]")); + frame->SetName(TString(histo->GetName()) + "_frame"); + dh->plotOn(frame); + + // Build the composite model + RooAbsPdf* model = buildModel(&x, signalType, backgroundType); + + std::unique_ptr chi2{model->createChi2(*dh, RooFit::DataError(RooAbsData::SumW2))}; + + // Fit the composite model + // ----------------------- + // Fit with likelihood + if (!useChi2_) { + if (sumW2Error) + model->fitTo(*dh, RooFit::Save(), RooFit::SumW2Error(kTRUE)); + else + model->fitTo(*dh); + } + // Fit with chi^2 + else { + edm::LogWarning("FitWithRooFit") << "FITTING WITH CHI^2"; + RooMinimizer m(*chi2); + m.migrad(); + m.hesse(); + // RooFitResult* r_chi2_wgt = m.save(); + } + model->plotOn(frame, RooFit::LineColor(kRed)); + model->plotOn(frame, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed)); + model->paramOn(frame, + RooFit::Label("fit result"), + RooFit::Layout(0.65, 0.90, 0.90), + RooFit::Format("NEU", RooFit::AutoPrecision(2))); + + // TODO: fix next lines to get the prob(chi2) (ndof should be dynamically set according to the choosen pdf) + // double chi2 = xframe.chiSquare("model","data",ndof); + // double ndoff = xframeGetNbinsX(); + // double chi2prob = TMath::Prob(chi2,ndoff); + + // 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 + // --------------------------------------------------------------------------------------------- + + // If histogram has custom error (i.e. its contents is does not originate from a Poisson process + // but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead + // (same error bars as shown by ROOT) + RooPlot* frame2 = x.frame(RooFit::Title("di-muon mass M(#mu^{+}#mu^{-}) [GeV]")); + dh->plotOn(frame2, RooFit::DataError(RooAbsData::SumW2)); + model->plotOn(frame2, RooFit::LineColor(kRed)); + model->plotOn(frame2, RooFit::Components(backgroundType), RooFit::LineStyle(kDashed)); + model->paramOn(frame2, + RooFit::Label("fit result"), + RooFit::Layout(0.65, 0.90, 0.90), + RooFit::Format("NEU", RooFit::AutoPrecision(2))); + + // Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used + // in a maximum likelihood fit + // + // A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition + // of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly + // fitted with a chi^2 fit (see rf602_chi2fit.C) + + // Draw all frames on a canvas + // if( canvas == 0 ) { + // canvas = new TCanvas("rf102_dataimport","rf102_dataimport",800,800); + // canvas->cd(1); + // } + // canvas->Divide(2,1); + // canvas->cd(1) ; frame->Draw(); + // canvas->cd(2) ; frame2->Draw(); + + frame2->Draw(); +} + +// Initialization methods for all the parameters +void FitWithRooFit::initMean( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (mean_ != nullptr) + delete mean_; + mean_ = new RooRealVar(name, title, value, min, max); + initVal_mean = value; +} + +void FitWithRooFit::initMean2( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (mean2_ != nullptr) + delete mean2_; + mean2_ = new RooRealVar(name, title, value, min, max); + initVal_mean2 = value; +} + +void FitWithRooFit::initMean3( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (mean3_ != nullptr) + delete mean3_; + mean3_ = new RooRealVar(name, title, value, min, max); + initVal_mean3 = value; +} + +void FitWithRooFit::initSigma( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (sigma_ != nullptr) + delete sigma_; + sigma_ = new RooRealVar(name, title, value, min, max); + initVal_sigma = value; +} + +void FitWithRooFit::initSigma2( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (sigma2_ != nullptr) + delete sigma2_; + sigma2_ = new RooRealVar(name, title, value, min, max); + initVal_sigma2 = value; +} + +void FitWithRooFit::initSigma3( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (sigma3_ != nullptr) + delete sigma3_; + sigma3_ = new RooRealVar(name, title, value, min, max); + initVal_sigma3 = value; +} + +void FitWithRooFit::initGamma( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (gamma_ != nullptr) + delete gamma_; + gamma_ = new RooRealVar(name, title, value, min, max); + initVal_gamma = value; +} + +void FitWithRooFit::initGaussFrac( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (gaussFrac_ != nullptr) + delete gaussFrac_; + gaussFrac_ = new RooRealVar(name, title, value, min, max); + initVal_gaussFrac = value; +} + +void FitWithRooFit::initGaussFrac2( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (gaussFrac2_ != nullptr) + delete gaussFrac2_; + gaussFrac2_ = new RooRealVar(name, title, value, min, max); + initVal_gaussFrac2 = value; +} + +void FitWithRooFit::initExpCoeffA0( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (expCoeffa0_ != nullptr) + delete expCoeffa0_; + expCoeffa0_ = new RooRealVar(name, title, value, min, max); + initVal_expCoeffa0 = value; +} + +void FitWithRooFit::initExpCoeffA1( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (expCoeffa1_ != nullptr) + delete expCoeffa1_; + expCoeffa1_ = new RooRealVar(name, title, value, min, max); + initVal_expCoeffa1 = value; +} + +void FitWithRooFit::initExpCoeffA2( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (expCoeffa2_ != nullptr) + delete expCoeffa2_; + expCoeffa2_ = new RooRealVar(name, title, value, min, max); + initVal_expCoeffa2 = value; +} + +void FitWithRooFit::initFsig( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (fsig_ != nullptr) + delete fsig_; + fsig_ = new RooRealVar(name, title, value, min, max); + initVal_fsig = value; +} + +void FitWithRooFit::initA0( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (a0_ != nullptr) + delete a0_; + a0_ = new RooRealVar(name, title, value, min, max); + initVal_a0 = value; +} + +void FitWithRooFit::initA1( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (a1_ != nullptr) + delete a1_; + a1_ = new RooRealVar(name, title, value, min, max); + initVal_a1 = value; +} + +void FitWithRooFit::initA2( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (a2_ != nullptr) + delete a2_; + a2_ = new RooRealVar(name, title, value, min, max); + initVal_a2 = value; +} + +void FitWithRooFit::initA3( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (a3_ != nullptr) + delete a3_; + a3_ = new RooRealVar(name, title, value, min, max); + initVal_a3 = value; +} + +void FitWithRooFit::initA4( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (a4_ != nullptr) + delete a4_; + a4_ = new RooRealVar(name, title, value, min, max); + initVal_a4 = value; +} + +void FitWithRooFit::initA5( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (a5_ != nullptr) + delete a5_; + a5_ = new RooRealVar(name, title, value, min, max); + initVal_a5 = value; +} + +void FitWithRooFit::initA6( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (a6_ != nullptr) + delete a6_; + a6_ = new RooRealVar(name, title, value, min, max); + initVal_a6 = value; +} + +void FitWithRooFit::initAlpha( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (alpha_ != nullptr) + delete alpha_; + alpha_ = new RooRealVar(name, title, value, min, max); + initVal_alpha = value; +} + +void FitWithRooFit::initN( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (n_ != nullptr) + delete n_; + n_ = new RooRealVar(name, title, value, min, max); + initVal_n = value; +} + +void FitWithRooFit::initFGCB( + const double& value, const double& min, const double& max, const TString& name, const TString& title) { + if (fGCB_ != nullptr) + delete fGCB_; + fGCB_ = new RooRealVar(name, title, value, min, max); + initVal_fGCB = value; +} + +void FitWithRooFit::reinitializeParameters() { + if (mean_ != nullptr) + mean_->setVal(initVal_mean); + if (mean2_ != nullptr) + mean2_->setVal(initVal_mean2); + if (mean3_ != nullptr) + mean3_->setVal(initVal_mean3); + if (sigma_ != nullptr) + sigma_->setVal(initVal_sigma); + if (sigma2_ != nullptr) + sigma2_->setVal(initVal_sigma2); + if (sigma3_ != nullptr) + sigma3_->setVal(initVal_sigma3); + if (gamma_ != nullptr) + gamma_->setVal(initVal_gamma); + if (gaussFrac_ != nullptr) + gaussFrac_->setVal(initVal_gaussFrac); + if (gaussFrac2_ != nullptr) + gaussFrac2_->setVal(initVal_gaussFrac2); + if (expCoeffa0_ != nullptr) + expCoeffa0_->setVal(initVal_expCoeffa0); + if (expCoeffa1_ != nullptr) + expCoeffa1_->setVal(initVal_expCoeffa1); + if (expCoeffa2_ != nullptr) + expCoeffa2_->setVal(initVal_expCoeffa2); + if (fsig_ != nullptr) + fsig_->setVal(initVal_fsig); + if (a0_ != nullptr) + a0_->setVal(initVal_a0); + if (a1_ != nullptr) + a1_->setVal(initVal_a1); + if (a2_ != nullptr) + a2_->setVal(initVal_a2); + if (a3_ != nullptr) + a3_->setVal(initVal_a3); + if (a4_ != nullptr) + a4_->setVal(initVal_a4); + if (a5_ != nullptr) + a5_->setVal(initVal_a5); + if (a6_ != nullptr) + a6_->setVal(initVal_a6); + if (alpha_ != nullptr) + alpha_->setVal(initVal_alpha); + if (n_ != nullptr) + n_->setVal(initVal_n); + if (fGCB_ != nullptr) + fGCB_->setVal(initVal_fGCB); +} + +/// Build the model for the specified signal type +RooAbsPdf* FitWithRooFit::buildSignalModel(RooRealVar* x, const TString& signalType) { + RooAbsPdf* signal = nullptr; + if (signalType == "gaussian") { + // Fit a Gaussian p.d.f to the data + if ((mean_ == nullptr) || (sigma_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean and sigma"; + exit(1); + } + signal = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma_); + } else if (signalType == "doubleGaussian") { + // Fit with double gaussian + if ((mean_ == nullptr) || (sigma_ == nullptr) || (sigma2_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and sigma2"; + exit(1); + } + RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_); + RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean_, *sigma2_); + signal = new RooAddModel("doubleGaussian", "double gaussian", RooArgList(*gaussModel, *gaussModel2), *gaussFrac_); + } else if (signalType == "tripleGaussian") { + // Fit with triple gaussian + if ((mean_ == nullptr) || (mean2_ == nullptr) || (mean3_ == nullptr) || (sigma_ == nullptr) || + (sigma2_ == nullptr) || (sigma3_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, " + "mean3, sigma, sigma2, sigma3"; + exit(1); + } + RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_); + RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean2_, *sigma2_); + RooGaussModel* gaussModel3 = new RooGaussModel("gaussModel3", "gaussModel3", *x, *mean3_, *sigma3_); + signal = new RooAddModel("tripleGaussian", + "triple gaussian", + RooArgList(*gaussModel, *gaussModel2, *gaussModel3), + RooArgList(*gaussFrac_, *gaussFrac2_)); + } else if (signalType == "breitWigner") { + // Fit a Breit-Wigner + if ((mean_ == nullptr) || (gamma_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma"; + exit(1); + } + signal = new RooBreitWigner("breiWign", "breitWign", *x, *mean_, *gamma_); + } else if (signalType == "relBreitWigner") { + // Fit a relativistic Breit-Wigner + if ((mean_ == nullptr) || (gamma_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma"; + exit(1); + } + signal = new RooGenericPdf("Relativistic Breit-Wigner", + "RBW", + "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))", + RooArgList(*x, *mean_, *gamma_)); + } else if (signalType == "voigtian") { + // Fit a Voigtian + if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and gamma"; + exit(1); + } + signal = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_); + } else if (signalType == "crystalBall") { + // Fit a CrystalBall + if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, " + "alpha and n"; + exit(1); + } + signal = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_); + } else if (signalType == "breitWignerTimesCB") { + // Fit a Breit Wigner convoluted with a CrystalBall + if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || + (alpha_ == nullptr) || (n_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, " + "sigma, gamma, alpha and n"; + exit(1); + } + RooAbsPdf* bw = new RooBreitWigner("breiWigner", "breitWigner", *x, *mean_, *gamma_); + RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_); + signal = new RooFFTConvPdf("breitWignerTimesCB", "breitWignerTimesCB", *x, *bw, *cb); + } else if (signalType == "relBreitWignerTimesCB") { + // Fit a relativistic Breit Wigner convoluted with a CrystalBall + if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || + (alpha_ == nullptr) || (n_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, " + "sigma, gamma, alpha and n"; + exit(1); + } + RooGenericPdf* bw = new RooGenericPdf("Relativistic Breit-Wigner", + "RBW", + "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))", + RooArgList(*x, *mean_, *gamma_)); + RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_); + signal = new RooFFTConvPdf("relBreitWignerTimesCB", "relBreitWignerTimesCB", *x, *bw, *cb); + } else if (signalType == "gaussianPlusCrystalBall") { + // Fit a Gaussian + CrystalBall with the same mean + if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) || + (fGCB_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, " + "sigma2, alpha, n and fGCB"; + exit(1); + } + RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_); + RooAbsPdf* tempGaussian = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma2_); + + signal = + new RooAddPdf("gaussianPlusCrystalBall", "gaussianPlusCrystalBall", RooArgList(*tempCB, *tempGaussian), *fGCB_); + } else if (signalType == "voigtianPlusCrystalBall") { + // Fit a Voigtian + CrystalBall with the same mean + if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || + (sigma2_ == nullptr) || (fGCB_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, " + "sigma, sigma2, alpha, n and fGCB"; + exit(1); + } + RooAbsPdf* tempVoigt = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_); + RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_); + + signal = + new RooAddPdf("voigtianPlusCrystalBall", "voigtianPlusCrystalBall", RooArgList(*tempCB, *tempVoigt), *fGCB_); + } else if (signalType == "breitWignerPlusCrystalBall") { + // Fit a Breit-Wigner + CrystalBall with the same mean + if ((mean_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) || + (fGCB_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, " + "sigma, alpha, n and fGCB"; + exit(1); + } + RooAbsPdf* tempBW = new RooBreitWigner("breitWign", "breitWign", *x, *mean_, *gamma_); + RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_); + + signal = + new RooAddPdf("breitWignerPlusCrystalBall", "breitWignerPlusCrystalBall", RooArgList(*tempCB, *tempBW), *fGCB_); + } + + else if (signalType != "") { + edm::LogError("FitWithRooFit") << "Unknown signal function: " << signalType << ". Signal will not be in the model"; + } + return signal; +} + +/// Build the model for the specified background type +RooAbsPdf* FitWithRooFit::buildBackgroundModel(RooRealVar* x, const TString& backgroundType) { + RooAbsPdf* background = nullptr; + if (backgroundType == "exponential") { + // Add an exponential for the background + if ((expCoeffa1_ == nullptr) || (fsig_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeffa1 and fsig"; + exit(1); + } + background = new RooExponential("exponential", "exponential", *x, *expCoeffa1_); + } + + if (backgroundType == "exponentialpol") { + // Add an exponential for the background + if ((expCoeffa0_ == nullptr) || (expCoeffa1_ == nullptr) || (expCoeffa2_ == nullptr) || (fsig_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeff and fsig"; + exit(1); + } + background = new RooGenericPdf("exponential", + "exponential", + "TMath::Exp(@1+@2*@0+@3*@0*@0)", + RooArgList(*x, *expCoeffa0_, *expCoeffa1_, *expCoeffa2_)); + } + + else if (backgroundType == "chebychev0") { + // Add a linear background + if (a0_ == nullptr) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize a0"; + exit(1); + } + background = new RooChebychev("chebychev0", "chebychev0", *x, *a0_); + } else if (backgroundType == "chebychev1") { + // Add a 2nd order chebychev polynomial background + if ((a0_ == nullptr) || (a1_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize a0 and a1"; + exit(1); + } + background = new RooChebychev("chebychev1", "chebychev1", *x, RooArgList(*a0_, *a1_)); + } else if (backgroundType == "chebychev3") { + // Add a 3rd order chebychev polynomial background + if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2 and a3"; + exit(1); + } + background = new RooChebychev("3rdOrderPol", "3rdOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_)); + } + + else if (backgroundType == "chebychev6") { + // Add a 6th order chebychev polynomial background + if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr) || (a4_ == nullptr) || + (a5_ == nullptr) || (a6_ == nullptr)) { + edm::LogError("FitWithRooFit") + << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2, a3, " + "a4, a5 and a6"; + exit(1); + } + background = + new RooChebychev("6thOrderPol", "6thOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_, *a4_, *a5_, *a6_)); + } + + return background; +} + +/// Build the model to fit +RooAbsPdf* FitWithRooFit::buildModel(RooRealVar* x, const TString& signalType, const TString& backgroundType) { + RooAbsPdf* model = nullptr; + + RooAbsPdf* signal = buildSignalModel(x, signalType); + RooAbsPdf* background = buildBackgroundModel(x, backgroundType); + + if ((signal != nullptr) && (background != nullptr)) { + // Combine signal and background pdfs + edm::LogPrint("FitWithRooFit") << "Building model with signal and backgound"; + model = new RooAddPdf("model", "model", RooArgList(*signal, *background), *fsig_); + } else if (signal != nullptr) { + edm::LogPrint("FitWithRooFit") << "Building model with signal"; + model = signal; + } else if (background != nullptr) { + edm::LogPrint("FitWithRooFit") << "Building model with backgound"; + model = background; + } else { + edm::LogWarning("FitWithRooFit") << "Nothing to fit"; + exit(0); + } + return model; +} From 267ef5a0869551ed6ebcb620e95aeb522d24c573 Mon Sep 17 00:00:00 2001 From: mmusich Date: Tue, 20 Jun 2023 16:36:42 +0200 Subject: [PATCH 2/2] use FitWithRooFit in DiMuonMassBiasClient --- DQMOffline/Alignment/BuildFile.xml | 1 + .../interface/DiMuonMassBiasClient.h | 2 + .../Alignment/src/DiMuonMassBiasClient.cc | 92 ++++++++++++++++++- 3 files changed, 92 insertions(+), 3 deletions(-) diff --git a/DQMOffline/Alignment/BuildFile.xml b/DQMOffline/Alignment/BuildFile.xml index 97b183bef88e3..60ff2194046e9 100644 --- a/DQMOffline/Alignment/BuildFile.xml +++ b/DQMOffline/Alignment/BuildFile.xml @@ -1,4 +1,5 @@ + diff --git a/DQMOffline/Alignment/interface/DiMuonMassBiasClient.h b/DQMOffline/Alignment/interface/DiMuonMassBiasClient.h index 6edef3b8b6c8b..8781fc70a7136 100644 --- a/DQMOffline/Alignment/interface/DiMuonMassBiasClient.h +++ b/DQMOffline/Alignment/interface/DiMuonMassBiasClient.h @@ -81,12 +81,14 @@ class DiMuonMassBiasClient : public DQMEDHarvester { void bookMEs(DQMStore::IBooker& ibooker); void getMEsToHarvest(DQMStore::IGetter& igetter); diMuonMassBias::fitOutputs fitLineShape(TH1* hist, const bool& fitBackground = false) const; + diMuonMassBias::fitOutputs fitBWTimesCB(TH1* hist) const; void fitAndFillProfile(std::pair toHarvest, DQMStore::IBooker& iBooker); void fitAndFillHisto(std::pair toHarvest, DQMStore::IBooker& iBooker); // data members const std::string TopFolder_; const bool useTH1s_; + const bool useBWtimesCB_; const bool fitBackground_; const bool useRooCBShape_; const bool useRooCMSShape_; diff --git a/DQMOffline/Alignment/src/DiMuonMassBiasClient.cc b/DQMOffline/Alignment/src/DiMuonMassBiasClient.cc index 108c73a46ef81..26c1dbb825ba7 100644 --- a/DQMOffline/Alignment/src/DiMuonMassBiasClient.cc +++ b/DQMOffline/Alignment/src/DiMuonMassBiasClient.cc @@ -1,4 +1,5 @@ // user includes +#include "Alignment/OfflineValidation/interface/FitWithRooFit.h" #include "DQMOffline/Alignment/interface/DiMuonMassBiasClient.h" #include "DataFormats/Histograms/interface/DQMToken.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -20,6 +21,7 @@ DiMuonMassBiasClient::DiMuonMassBiasClient(edm::ParameterSet const& iConfig) : TopFolder_(iConfig.getParameter("FolderName")), useTH1s_(iConfig.getParameter("useTH1s")), + useBWtimesCB_(iConfig.getParameter("useBWtimesCB")), fitBackground_(iConfig.getParameter("fitBackground")), useRooCBShape_(iConfig.getParameter("useRooCBShape")), useRooCMSShape_(iConfig.getParameter("useRooCMSShape")), @@ -164,7 +166,7 @@ void DiMuonMassBiasClient::fitAndFillProfile(std::pairProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin); Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin)); - diMuonMassBias::fitOutputs results = fitLineShape(Proj); + diMuonMassBias::fitOutputs results = useBWtimesCB_ ? fitBWTimesCB(Proj) : fitLineShape(Proj, fitBackground_); if (results.isInvalid()) { edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl; @@ -237,7 +239,7 @@ void DiMuonMassBiasClient::fitAndFillHisto(std::pairProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin); Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin)); - diMuonMassBias::fitOutputs results = fitLineShape(Proj); + diMuonMassBias::fitOutputs results = useBWtimesCB_ ? fitBWTimesCB(Proj) : fitLineShape(Proj, fitBackground_); if (results.isInvalid()) { edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl; @@ -342,7 +344,7 @@ diMuonMassBias::fitOutputs DiMuonMassBiasClient::fitLineShape(TH1* hist, const b RooRealVar b("N_{b}", "Number of background events", 0, hist->GetEntries() / 10.); RooRealVar s("N_{s}", "Number of signal events", 0, hist->GetEntries()); - if (fitBackground_) { + if (fitBackground) { RooArgList listPdf; if (useRooCBShape_) { if (useRooCMSShape_) { @@ -415,6 +417,89 @@ diMuonMassBias::fitOutputs DiMuonMassBiasClient::fitLineShape(TH1* hist, const b return diMuonMassBias::fitOutputs(resultM, resultW); } +//----------------------------------------------------------------------------------- +diMuonMassBias::fitOutputs DiMuonMassBiasClient::fitBWTimesCB(TH1* hist) const +//----------------------------------------------------------------------------------- +{ + if (hist->GetEntries() < diMuonMassBias::minimumHits) { + edm::LogWarning("DiMuonMassBiasClient") << " Input histogram:" << hist->GetName() << " has not enough entries (" + << hist->GetEntries() << ") for a meaningful Voigtian fit!\n" + << "Skipping!"; + + return diMuonMassBias::fitOutputs(Measurement1D(0., 0.), Measurement1D(0., 0.)); + } + + TCanvas* c1 = new TCanvas(); + if (debugMode_) { + c1->Clear(); + c1->SetLeftMargin(0.15); + c1->SetRightMargin(0.10); + } + + // silence messages + RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL); + + Double_t xMean = 91.1876; + Double_t xMin = hist->GetXaxis()->GetXmin(); + Double_t xMax = hist->GetXaxis()->GetXmax(); + + if (debugMode_) { + edm::LogPrint("DiMuonMassBiasClient") << "fitting range: (" << xMin << "-" << xMax << ")" << std::endl; + } + + double sigma(2.); + double sigmaMin(0.1); + double sigmaMax(10.); + + double sigma2(0.1); + double sigma2Min(0.); + double sigma2Max(10.); + + std::unique_ptr fitter = std::make_unique(); + + bool useChi2(false); + + fitter->useChi2_ = useChi2; + fitter->initMean(xMean, xMin, xMax); + fitter->initSigma(sigma, sigmaMin, sigmaMax); + fitter->initSigma2(sigma2, sigma2Min, sigma2Max); + fitter->initAlpha(1.5, 0.05, 10.); + fitter->initN(1, 0.01, 100.); + fitter->initFGCB(0.4, 0., 1.); + fitter->initGamma(2.4952, 0., 10.); + fitter->gamma()->setConstant(kTRUE); + fitter->initMean2(0., -20., 20.); + fitter->mean2()->setConstant(kTRUE); + fitter->initSigma(1.2, 0., 5.); + fitter->initAlpha(1.5, 0.05, 10.); + fitter->initN(1, 0.01, 100.); + fitter->initExpCoeffA0(-1., -10., 10.); + fitter->initExpCoeffA1(0., -10., 10.); + fitter->initExpCoeffA2(0., -2., 2.); + fitter->initFsig(0.9, 0., 1.); + fitter->initA0(0., -10., 10.); + fitter->initA1(0., -10., 10.); + fitter->initA2(0., -10., 10.); + fitter->initA3(0., -10., 10.); + fitter->initA4(0., -10., 10.); + fitter->initA5(0., -10., 10.); + fitter->initA6(0., -10., 10.); + + fitter->fit(hist, "breitWignerTimesCB", "exponential", xMin, xMax, false); + + TString histName = hist->GetName(); + + if (debugMode_) { + c1->Print("fit_debug" + histName + ".pdf"); + } + delete c1; + + Measurement1D resultM(fitter->mean()->getValV(), fitter->mean()->getError()); + Measurement1D resultW(fitter->sigma()->getValV(), fitter->sigma()->getError()); + + return diMuonMassBias::fitOutputs(resultM, resultW); +} + //----------------------------------------------------------------------------------- void DiMuonMassBiasClient::fillDescriptions(edm::ConfigurationDescriptions& descriptions) //----------------------------------------------------------------------------------- @@ -422,6 +507,7 @@ void DiMuonMassBiasClient::fillDescriptions(edm::ConfigurationDescriptions& desc edm::ParameterSetDescription desc; desc.add("FolderName", "DiMuonMassBiasMonitor"); desc.add("useTH1s", false); + desc.add("useBWtimesCB", false); desc.add("fitBackground", false); desc.add("useRooCMSShape", false); desc.add("useRooCBShape", false);