Alignment/src/DiMuonMassBiasClient.cc

0001 // user includes
0002 #include "Alignment/OfflineValidation/interface/FitWithRooFit.h"
0003 #include "DQMOffline/Alignment/interface/DiMuonMassBiasClient.h"
0004 #include "DataFormats/Histograms/interface/DQMToken.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "PhysicsTools/TagAndProbe/interface/RooCMSShape.h"
0008
0009 // RooFit includes
0010 #include "TCanvas.h"
0011 #include "RooAddPdf.h"
0012 #include "RooDataHist.h"
0013 #include "RooExponential.h"
0014 #include "RooGaussian.h"
0015 #include "RooPlot.h"
0016 #include "RooRealVar.h"
0017 #include "RooVoigtian.h"
0018 #include "RooCBShape.h"
0019
0020 //-----------------------------------------------------------------------------------
0021 DiMuonMassBiasClient::DiMuonMassBiasClient(edm::ParameterSet const& iConfig)
0022     : TopFolder_(iConfig.getParameter<std::string>("FolderName")),
0023       useTH1s_(iConfig.getParameter<bool>("useTH1s")),
0024       useBWtimesCB_(iConfig.getParameter<bool>("useBWtimesCB")),
0025       fitBackground_(iConfig.getParameter<bool>("fitBackground")),
0026       useRooCBShape_(iConfig.getParameter<bool>("useRooCBShape")),
0027       useRooCMSShape_(iConfig.getParameter<bool>("useRooCMSShape")),
0028       debugMode_(iConfig.getParameter<bool>("debugMode")),
0029       MEtoHarvest_(iConfig.getParameter<std::vector<std::string>>("MEtoHarvest"))
0030 //-----------------------------------------------------------------------------------
0031 {
0032   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::Constructing DiMuonMassBiasClient ";
0033
0034   // fill the parameters for the fit
0035   edm::ParameterSet fit_par = iConfig.getParameter<edm::ParameterSet>("fit_par");
0036   diMuonMassBias::fillArrayF(meanConfig_, fit_par, "mean_par");
0037   diMuonMassBias::fillArrayF(widthConfig_, fit_par, "width_par");
0038   diMuonMassBias::fillArrayF(sigmaConfig_, fit_par, "sigma_par");
0039
0040   if (debugMode_) {
0041     edm::LogPrint("DiMuonMassBiasClient")
0042         << "mean: " << meanConfig_[0] << " (" << meanConfig_[1] << "," << meanConfig_[2] << ") " << std::endl;
0043     edm::LogPrint("DiMuonMassBiasClient")
0044         << "width: " << widthConfig_[0] << " (" << widthConfig_[1] << "," << widthConfig_[2] << ")" << std::endl;
0045     edm::LogPrint("DiMuonMassBiasClient")
0046         << "sigma: " << sigmaConfig_[0] << " (" << sigmaConfig_[1] << "," << sigmaConfig_[2] << ")" << std::endl;
0047   }
0048 }
0049
0050 //-----------------------------------------------------------------------------------
0051 DiMuonMassBiasClient::~DiMuonMassBiasClient()
0052 //-----------------------------------------------------------------------------------
0053 {
0054   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::Deleting DiMuonMassBiasClient ";
0055 }
0056
0057 //-----------------------------------------------------------------------------------
0058 void DiMuonMassBiasClient::beginJob(void)
0059 //-----------------------------------------------------------------------------------
0060 {
0061   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::beginJob done";
0062 }
0063
0064 //-----------------------------------------------------------------------------------
0065 void DiMuonMassBiasClient::beginRun(edm::Run const& run, edm::EventSetup const& eSetup)
0066 //-----------------------------------------------------------------------------------
0067 {
0068   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient:: Begining of Run";
0069 }
0070
0071 //-----------------------------------------------------------------------------------
0072 void DiMuonMassBiasClient::bookMEs(DQMStore::IBooker& iBooker)
0073 //-----------------------------------------------------------------------------------
0074 {
0075   iBooker.setCurrentFolder(TopFolder_ + "/DiMuonMassBiasMonitor/MassBias/Profiles");
0076   for (const auto& [key, ME] : harvestTargets_) {
0077     if (ME == nullptr) {
0078       edm::LogError("DiMuonMassBiasClient") << "could not find MonitorElement for key: " << key << std::endl;
0079       continue;
0080     }
0081
0082     const auto& title = ME->getTitle();
0083     const auto& xtitle = ME->getAxisTitle(1);
0084     const auto& ytitle = ME->getAxisTitle(2);
0085
0086     const auto& nxbins = ME->getNbinsX();
0087     const auto& xmin = ME->getAxisMin(1);
0088     const auto& xmax = ME->getAxisMax(1);
0089
0090     MonitorElement* meanToBook =
0091         iBooker.book1D(("Mean" + key), (title + ";#LT M_{#mu^{-}#mu^{+}} #GT [GeV];" + ytitle), nxbins, xmin, xmax);
0092     meanHistos_.insert({key, meanToBook});
0093
0094     MonitorElement* sigmaToBook =
0095         iBooker.book1D(("Sigma" + key), (title + ";" + xtitle + ";" + "#sigma of " + ytitle), nxbins, xmin, xmax);
0096     widthHistos_.insert({key, sigmaToBook});
0097   }
0098 }
0099
0100 //-----------------------------------------------------------------------------------
0101 void DiMuonMassBiasClient::getMEsToHarvest(DQMStore::IGetter& iGetter)
0102 //-----------------------------------------------------------------------------------
0103 {
0104   std::string inFolder = TopFolder_ + "/DiMuonMassBiasMonitor/MassBias/";
0105
0106   //loop on the list of histograms to harvest
0107   for (const auto& hname : MEtoHarvest_) {
0108     MonitorElement* toHarvest = iGetter.get(inFolder + hname);
0109
0110     if (toHarvest == nullptr) {
0111       edm::LogError("DiMuonMassBiasClient") << "could not find input MonitorElement: " << inFolder + hname << std::endl;
0112       continue;
0113     }
0114
0115     harvestTargets_.insert({hname, toHarvest});
0116   }
0117 }
0118
0119 //-----------------------------------------------------------------------------------
0120 void DiMuonMassBiasClient::fitAndFillProfile(std::pair<std::string, MonitorElement*> toHarvest,
0121                                              DQMStore::IBooker& iBooker)
0122 //-----------------------------------------------------------------------------------
0123 {
0124   const auto& key = toHarvest.first;
0125   const auto& ME = toHarvest.second;
0126
0127   if (debugMode_)
0128     edm::LogPrint("DiMuonMassBiasClient") << "dealing with key: " << key << std::endl;
0129
0130   if (ME == nullptr) {
0131     edm::LogError("DiMuonMassBiasClient") << "could not find MonitorElement for key: " << key << std::endl;
0132     return;
0133   }
0134
0135   const auto& title = ME->getTitle();
0136   const auto& xtitle = ME->getAxisTitle(1);
0137   const auto& ytitle = ME->getAxisTitle(2);
0138
0139   const auto& nxbins = ME->getNbinsX();
0140   const auto& xmin = ME->getAxisMin(1);
0141   const auto& xmax = ME->getAxisMax(1);
0142
0143   TProfile* p_mean = new TProfile(("Mean" + key).c_str(),
0144                                   (title + ";" + xtitle + ";#LT M_{#mu^{-}#mu^{+}} #GT [GeV]").c_str(),
0145                                   nxbins,
0146                                   xmin,
0147                                   xmax,
0148                                   "g");
0149
0150   TProfile* p_width = new TProfile(
0151       ("Sigma" + key).c_str(), (title + ";" + xtitle + ";#sigma of " + ytitle).c_str(), nxbins, xmin, xmax, "g");
0152
0153   p_mean->Sumw2();
0154   p_width->Sumw2();
0155
0156   TH2F* bareHisto = ME->getTH2F();
0157   for (int bin = 1; bin <= nxbins; bin++) {
0158     const auto& xaxis = bareHisto->GetXaxis();
0159     const auto& low_edge = xaxis->GetBinLowEdge(bin);
0160     const auto& high_edge = xaxis->GetBinUpEdge(bin);
0161
0162     if (debugMode_)
0163       edm::LogPrint("DiMuonMassBiasClient") << "dealing with bin: " << bin << " range: (" << std::setprecision(2)
0164                                             << low_edge << "," << std::setprecision(2) << high_edge << ")";
0165
0166     TH1D* Proj = bareHisto->ProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin);
0167     Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin));
0168
0169     diMuonMassBias::fitOutputs results = useBWtimesCB_ ? fitBWTimesCB(Proj) : fitLineShape(Proj, fitBackground_);
0170
0171     if (results.isInvalid()) {
0172       edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl;
0173       continue;
0174     }
0175
0176     // fill the mean profiles
0177     const Measurement1D& bias = results.getBias();
0178
0179     // ============================================= DISCLAIMER ================================================
0180     // N.B. this is sort of a hack in order to fill arbitrarily both central values and error bars of a TProfile.
0181     // Choosing the option "g" in the constructor the bin error will be 1/sqrt(W(j)), where W(j) is the sum of weights.
0182     // Filling the sum of weights with the 1 / err^2, the bin error automatically becomes "err".
0183     // In order to avoid the central value to be shifted, that's divided by 1 / err^2 as well.
0184     // For more information, please consult the https://root.cern.ch/doc/master/classTProfile.html
0185
0186     p_mean->SetBinContent(bin, bias.value() / (bias.error() * bias.error()));
0187     p_mean->SetBinEntries(bin, 1. / (bias.error() * bias.error()));
0188
0189     if (debugMode_)
0190       LogDebug("DiMuonBassBiasClient") << " Bin: " << bin << " value:  " << bias.value() << " from profile ( "
0191                                        << p_mean->GetBinContent(bin) << ") - error:  " << bias.error()
0192                                        << "  from profile ( " << p_mean->GetBinError(bin) << " )";
0193
0194     // fill the width profiles
0195     const Measurement1D& width = results.getWidth();
0196
0197     // see discussion above
0198     p_width->SetBinContent(bin, width.value() / (width.error() * width.error()));
0199     p_width->SetBinEntries(bin, 1. / (width.error() * width.error()));
0200   }
0201
0202   // now book the profiles
0203   iBooker.setCurrentFolder(TopFolder_ + "/DiMuonMassBiasMonitor/MassBias/Profiles");
0204   MonitorElement* meanToBook = iBooker.bookProfile(p_mean->GetName(), p_mean);
0205   meanProfiles_.insert({key, meanToBook});
0206
0207   MonitorElement* sigmaToBook = iBooker.bookProfile(p_width->GetName(), p_width);
0208   widthProfiles_.insert({key, sigmaToBook});
0209
0210   delete p_mean;
0211   delete p_width;
0212 }
0213
0214 //-----------------------------------------------------------------------------------
0215 void DiMuonMassBiasClient::fitAndFillHisto(std::pair<std::string, MonitorElement*> toHarvest,
0216                                            DQMStore::IBooker& iBooker)
0217 //-----------------------------------------------------------------------------------
0218 {
0219   const auto& key = toHarvest.first;
0220   const auto& ME = toHarvest.second;
0221
0222   if (debugMode_)
0223     edm::LogPrint("DiMuonMassBiasClient") << "dealing with key: " << key << std::endl;
0224
0225   if (ME == nullptr) {
0226     edm::LogError("DiMuonMassBiasClient") << "could not find MonitorElement for key: " << key << std::endl;
0227     return;
0228   }
0229
0230   TH2F* bareHisto = ME->getTH2F();
0231   for (int bin = 1; bin <= ME->getNbinsX(); bin++) {
0232     const auto& xaxis = bareHisto->GetXaxis();
0233     const auto& low_edge = xaxis->GetBinLowEdge(bin);
0234     const auto& high_edge = xaxis->GetBinUpEdge(bin);
0235
0236     if (debugMode_)
0237       edm::LogPrint("DiMuonMassBiasClient") << "dealing with bin: " << bin << " range: (" << std::setprecision(2)
0238                                             << low_edge << "," << std::setprecision(2) << high_edge << ")";
0239     TH1D* Proj = bareHisto->ProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin);
0240     Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin));
0241
0242     diMuonMassBias::fitOutputs results = useBWtimesCB_ ? fitBWTimesCB(Proj) : fitLineShape(Proj, fitBackground_);
0243
0244     if (results.isInvalid()) {
0245       edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl;
0246       continue;
0247     }
0248
0249     // fill the mean profiles
0250     const Measurement1D& bias = results.getBias();
0251     meanHistos_[key]->setBinContent(bin, bias.value());
0252     meanHistos_[key]->setBinError(bin, bias.error());
0253
0254     // fill the width profiles
0255     const Measurement1D& width = results.getWidth();
0256     widthHistos_[key]->setBinContent(bin, width.value());
0257     widthHistos_[key]->setBinError(bin, width.error());
0258   }
0259 }
0260
0261 //-----------------------------------------------------------------------------------
0262 void DiMuonMassBiasClient::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter)
0263 //-----------------------------------------------------------------------------------
0264 {
0265   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::endLuminosityBlock";
0266
0267   getMEsToHarvest(igetter);
0268
0269   // book the histograms upfront
0270   if (useTH1s_) {
0271     bookMEs(ibooker);
0272   }
0273
0274   for (const auto& element : harvestTargets_) {
0275     if (!useTH1s_) {
0276       // if using profiles
0277       this->fitAndFillProfile(element, ibooker);
0278     } else {
0279       // if using histograms
0280       this->fitAndFillHisto(element, ibooker);
0281     }
0282   }
0283 }
0284
0285 //-----------------------------------------------------------------------------------
0286 diMuonMassBias::fitOutputs DiMuonMassBiasClient::fitLineShape(TH1* hist, const bool& fitBackground) const
0287 //-----------------------------------------------------------------------------------
0288 {
0289   if (hist->GetEntries() < diMuonMassBias::minimumHits) {
0290     edm::LogWarning("DiMuonMassBiasClient") << " Input histogram:" << hist->GetName() << " has not enough entries ("
0291                                             << hist->GetEntries() << ") for a meaningful Voigtian fit!\n"
0292                                             << "Skipping!";
0293
0294     return diMuonMassBias::fitOutputs(Measurement1D(0., 0.), Measurement1D(0., 0.));
0295   }
0296
0297   TCanvas* c1 = new TCanvas();
0298   if (debugMode_) {
0299     c1->Clear();
0300     c1->SetLeftMargin(0.15);
0301     c1->SetRightMargin(0.10);
0302   }
0303
0304   // silence messages
0305   RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
0306
0307   Double_t xmin = hist->GetXaxis()->GetXmin();
0308   Double_t xmax = hist->GetXaxis()->GetXmax();
0309
0310   if (debugMode_) {
0311     edm::LogPrint("DiMuonMassBiasClient") << "fitting range: (" << xmin << "-" << xmax << ")" << std::endl;
0312   }
0313
0314   RooRealVar InvMass("InvMass", "di-muon mass M(#mu^{+}#mu^{-}) [GeV]", xmin, xmax);
0315   std::unique_ptr<RooPlot> frame{InvMass.frame()};
0316   RooDataHist datahist("datahist", "datahist", InvMass, RooFit::Import(*hist));
0317   datahist.plotOn(frame.get());
0318
0319   // parameters of the Voigtian
0320   RooRealVar mean("#mu", "mean", meanConfig_[0], meanConfig_[1], meanConfig_[2]);          //90.0, 60.0, 120.0 (for Z)
0321   RooRealVar width("width", "width", widthConfig_[0], widthConfig_[1], widthConfig_[2]);   // 5.0,  0.0, 120.0 (for Z)
0322   RooRealVar sigma("#sigma", "sigma", sigmaConfig_[0], sigmaConfig_[1], sigmaConfig_[2]);  // 5.0,  0.0, 120.0 (for Z)
0323   RooVoigtian voigt("voigt", "voigt", InvMass, mean, width, sigma);
0324
0325   // parameters of the Crystal-ball
0326   RooRealVar peakCB("peakCB", "peakCB", meanConfig_[0], meanConfig_[1], meanConfig_[2]);
0327   RooRealVar sigmaCB("#sigma", "sigma", sigmaConfig_[0], sigmaConfig_[1], sigmaConfig_[2]);
0328   RooRealVar alphaCB("#alpha", "alpha", 1., 0., 10.);
0329   RooRealVar nCB("n", "n", 1., 0., 100.);
0330   RooCBShape crystalball("crystalball", "crystalball", InvMass, peakCB, sigmaCB, alphaCB, nCB);
0331
0332   // for the simple background fit
0333   RooRealVar lambda("#lambda", "slope", 0., -50., 50.);
0334   RooExponential expo("expo", "expo", InvMass, lambda);
0335
0336   // for the more refined background fit
0337   RooRealVar exp_alpha("#alpha", "alpha", 40.0, 20.0, 160.0);
0338   RooRealVar exp_beta("#beta", "beta", 0.05, 0.0, 2.0);
0339   RooRealVar exp_gamma("#gamma", "gamma", 0.02, 0.0, 0.1);
0340   RooRealVar exp_peak("peak", "peak", meanConfig_[0]);
0341   RooCMSShape exp_pdf("exp_pdf", "bkg shape", InvMass, exp_alpha, exp_beta, exp_gamma, exp_peak);
0342
0343   // define the signal and background fractions
0344   RooRealVar b("N_{b}", "Number of background events", 0, hist->GetEntries() / 10.);
0345   RooRealVar s("N_{s}", "Number of signal events", 0, hist->GetEntries());
0346
0347   if (fitBackground) {
0348     RooArgList listPdf;
0349     if (useRooCBShape_) {
0350       if (useRooCMSShape_) {
0351         // crystal-ball + CMS-shape fit
0352         listPdf.add(crystalball);
0353         listPdf.add(exp_pdf);
0354       } else {
0355         // crystal-ball + exponential fit
0356         listPdf.add(crystalball);
0357         listPdf.add(expo);
0358       }
0359     } else {
0360       if (useRooCMSShape_) {
0361         // voigtian + CMS-shape fit
0362         listPdf.add(voigt);
0363         listPdf.add(exp_pdf);
0364       } else {
0365         // voigtian + exponential fit
0366         listPdf.add(voigt);
0367         listPdf.add(expo);
0368       }
0369     }
0370
0371     RooAddPdf fullModel("fullModel", "Signal + Background Model", listPdf, RooArgList(s, b));
0372     fullModel.fitTo(datahist, RooFit::PrintLevel(-1));
0373     fullModel.plotOn(frame.get(), RooFit::LineColor(kRed));
0374     if (useRooCMSShape_) {
0375       fullModel.plotOn(frame.get(), RooFit::Components(exp_pdf), RooFit::LineStyle(kDashed));  //Other option
0376     } else {
0377       fullModel.plotOn(frame.get(), RooFit::Components(expo), RooFit::LineStyle(kDashed));  //Other option
0378     }
0379     fullModel.paramOn(frame.get(), RooFit::Layout(0.65, 0.90, 0.90));
0380   } else {
0381     if (useRooCBShape_) {
0382       // use crystal-ball for a fit-only signal
0383       crystalball.fitTo(datahist, RooFit::PrintLevel(-1));
0384       crystalball.plotOn(frame.get(), RooFit::LineColor(kRed));  //this will show fit overlay on canvas
0385       crystalball.paramOn(frame.get(),
0386                           RooFit::Layout(0.65, 0.90, 0.90));  //this will display the fit parameters on canvas
0387     } else {
0388       // use voigtian for a fit-only signal
0389       voigt.fitTo(datahist, RooFit::PrintLevel(-1));
0390       voigt.plotOn(frame.get(), RooFit::LineColor(kRed));            //this will show fit overlay on canvas
0391       voigt.paramOn(frame.get(), RooFit::Layout(0.65, 0.90, 0.90));  //this will display the fit parameters on canvas
0392     }
0393   }
0394
0395   // Redraw data on top and print / store everything
0396   datahist.plotOn(frame.get());
0397   frame->GetYaxis()->SetTitle("n. of events");
0398   TString histName = hist->GetName();
0399   frame->SetName("frame" + histName);
0400   frame->SetTitle(hist->GetTitle());
0401   frame->Draw();
0402
0403   if (debugMode_) {
0404     c1->Print("fit_debug" + histName + ".pdf");
0405   }
0406   delete c1;
0407
0408   float mass_mean = useRooCBShape_ ? peakCB.getVal() : mean.getVal();
0409   float mass_sigma = useRooCBShape_ ? sigmaCB.getVal() : sigma.getVal();
0410
0411   float mass_mean_err = useRooCBShape_ ? peakCB.getError() : mean.getError();
0412   float mass_sigma_err = useRooCBShape_ ? sigmaCB.getError() : sigma.getError();
0413
0414   Measurement1D resultM(mass_mean, mass_mean_err);
0415   Measurement1D resultW(mass_sigma, mass_sigma_err);
0416
0417   return diMuonMassBias::fitOutputs(resultM, resultW);
0418 }
0419
0420 //-----------------------------------------------------------------------------------
0421 diMuonMassBias::fitOutputs DiMuonMassBiasClient::fitBWTimesCB(TH1* hist) const
0422 //-----------------------------------------------------------------------------------
0423 {
0424   if (hist->GetEntries() < diMuonMassBias::minimumHits) {
0425     edm::LogWarning("DiMuonMassBiasClient") << " Input histogram:" << hist->GetName() << " has not enough entries ("
0426                                             << hist->GetEntries() << ") for a meaningful Voigtian fit!\n"
0427                                             << "Skipping!";
0428
0429     return diMuonMassBias::fitOutputs(Measurement1D(0., 0.), Measurement1D(0., 0.));
0430   }
0431
0432   TCanvas* c1 = new TCanvas();
0433   if (debugMode_) {
0434     c1->Clear();
0435     c1->SetLeftMargin(0.15);
0436     c1->SetRightMargin(0.10);
0437   }
0438
0439   // silence messages
0440   RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
0441
0442   Double_t xMean = 91.1876;
0443   Double_t xMin = hist->GetXaxis()->GetXmin();
0444   Double_t xMax = hist->GetXaxis()->GetXmax();
0445
0446   if (debugMode_) {
0447     edm::LogPrint("DiMuonMassBiasClient") << "fitting range: (" << xMin << "-" << xMax << ")" << std::endl;
0448   }
0449
0450   double sigma(2.);
0451   double sigmaMin(0.1);
0452   double sigmaMax(10.);
0453
0454   double sigma2(0.1);
0455   double sigma2Min(0.);
0456   double sigma2Max(10.);
0457
0458   std::unique_ptr<FitWithRooFit> fitter = std::make_unique<FitWithRooFit>();
0459
0460   bool useChi2(false);
0461
0462   fitter->useChi2_ = useChi2;
0463   fitter->initMean(xMean, xMin, xMax);
0464   fitter->initSigma(sigma, sigmaMin, sigmaMax);
0465   fitter->initSigma2(sigma2, sigma2Min, sigma2Max);
0466   fitter->initAlpha(1.5, 0.05, 10.);
0467   fitter->initN(1, 0.01, 100.);
0468   fitter->initFGCB(0.4, 0., 1.);
0469   fitter->initGamma(2.4952, 0., 10.);
0470   fitter->gamma()->setConstant(kTRUE);
0471   fitter->initMean2(0., -20., 20.);
0472   fitter->mean2()->setConstant(kTRUE);
0473   fitter->initSigma(1.2, 0., 5.);
0474   fitter->initAlpha(1.5, 0.05, 10.);
0475   fitter->initN(1, 0.01, 100.);
0476   fitter->initExpCoeffA0(-1., -10., 10.);
0477   fitter->initExpCoeffA1(0., -10., 10.);
0478   fitter->initExpCoeffA2(0., -2., 2.);
0479   fitter->initFsig(0.9, 0., 1.);
0480   fitter->initA0(0., -10., 10.);
0481   fitter->initA1(0., -10., 10.);
0482   fitter->initA2(0., -10., 10.);
0483   fitter->initA3(0., -10., 10.);
0484   fitter->initA4(0., -10., 10.);
0485   fitter->initA5(0., -10., 10.);
0486   fitter->initA6(0., -10., 10.);
0487
0488   fitter->fit(hist, "breitWignerTimesCB", "exponential", xMin, xMax, false);
0489
0490   TString histName = hist->GetName();
0491
0492   if (debugMode_) {
0493     c1->Print("fit_debug" + histName + ".pdf");
0494   }
0495   delete c1;
0496
0497   Measurement1D resultM(fitter->mean()->getVal(), fitter->mean()->getError());
0498   Measurement1D resultW(fitter->sigma()->getVal(), fitter->sigma()->getError());
0499
0500   return diMuonMassBias::fitOutputs(resultM, resultW);
0501 }
0502
0503 //-----------------------------------------------------------------------------------
0504 void DiMuonMassBiasClient::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
0505 //-----------------------------------------------------------------------------------
0506 {
0507   edm::ParameterSetDescription desc;
0508   desc.add<std::string>("FolderName", "DiMuonMassBiasMonitor");
0509   desc.add<bool>("useTH1s", false);
0510   desc.add<bool>("useBWtimesCB", false);
0511   desc.add<bool>("fitBackground", false);
0512   desc.add<bool>("useRooCMSShape", false);
0513   desc.add<bool>("useRooCBShape", false);
0514   desc.add<bool>("debugMode", false);
0515
0516   edm::ParameterSetDescription fit_par;
0517   fit_par.add<std::vector<double>>("mean_par",
0518                                    {std::numeric_limits<float>::max(),
0519                                     std::numeric_limits<float>::max(),
0520                                     std::numeric_limits<float>::max()});  // par = mean
0521
0522   fit_par.add<std::vector<double>>("width_par",
0523                                    {std::numeric_limits<float>::max(),
0524                                     std::numeric_limits<float>::max(),
0525                                     std::numeric_limits<float>::max()});  // par = width
0526
0527   fit_par.add<std::vector<double>>("sigma_par",
0528                                    {std::numeric_limits<float>::max(),
0529                                     std::numeric_limits<float>::max(),
0530                                     std::numeric_limits<float>::max()});  // par = sigma
0531
0532   desc.add<edm::ParameterSetDescription>("fit_par", fit_par);
0533
0534   desc.add<std::vector<std::string>>("MEtoHarvest",
0535                                      {"DiMuMassVsMuMuPhi",
0536                                       "DiMuMassVsMuMuEta",
0537                                       "DiMuMassVsMuPlusPhi",
0538                                       "DiMuMassVsMuPlusEta",
0539                                       "DiMuMassVsMuMinusPhi",
0540                                       "DiMuMassVsMuMinusEta",
0541                                       "DiMuMassVsMuMuDeltaEta",
0542                                       "DiMuMassVsCosThetaCS"});
0543   descriptions.addWithDefaultLabel(desc);
0544 }
0545
0546 DEFINE_FWK_MODULE(DiMuonMassBiasClient);