11 #include "RooAddPdf.h" 12 #include "RooDataHist.h" 13 #include "RooExponential.h" 14 #include "RooGaussian.h" 16 #include "RooRealVar.h" 17 #include "RooVoigtian.h" 18 #include "RooCBShape.h" 22 : TopFolder_(iConfig.getParameter<
std::
string>(
"FolderName")),
23 useTH1s_(iConfig.getParameter<
bool>(
"useTH1s")),
24 useBWtimesCB_(iConfig.getParameter<
bool>(
"useBWtimesCB")),
25 fitBackground_(iConfig.getParameter<
bool>(
"fitBackground")),
26 useRooCBShape_(iConfig.getParameter<
bool>(
"useRooCBShape")),
27 useRooCMSShape_(iConfig.getParameter<
bool>(
"useRooCMSShape")),
28 debugMode_(iConfig.getParameter<
bool>(
"debugMode")),
32 edm::LogInfo(
"DiMuonMassBiasClient") <<
"DiMuonMassBiasClient::Constructing DiMuonMassBiasClient ";
54 edm::LogInfo(
"DiMuonMassBiasClient") <<
"DiMuonMassBiasClient::Deleting DiMuonMassBiasClient ";
61 edm::LogInfo(
"DiMuonMassBiasClient") <<
"DiMuonMassBiasClient::beginJob done";
68 edm::LogInfo(
"DiMuonMassBiasClient") <<
"DiMuonMassBiasClient:: Begining of Run";
78 edm::LogError(
"DiMuonMassBiasClient") <<
"could not find MonitorElement for key: " <<
key << std::endl;
82 const auto&
title =
ME->getTitle();
83 const auto&
xtitle =
ME->getAxisTitle(1);
84 const auto&
ytitle =
ME->getAxisTitle(2);
86 const auto& nxbins =
ME->getNbinsX();
87 const auto&
xmin =
ME->getAxisMin(1);
88 const auto&
xmax =
ME->getAxisMax(1);
110 if (toHarvest ==
nullptr) {
111 edm::LogError(
"DiMuonMassBiasClient") <<
"could not find input MonitorElement: " << inFolder + hname << std::endl;
124 const auto&
key = toHarvest.first;
125 const auto&
ME = toHarvest.second;
128 edm::LogPrint(
"DiMuonMassBiasClient") <<
"dealing with key: " <<
key << std::endl;
131 edm::LogError(
"DiMuonMassBiasClient") <<
"could not find MonitorElement for key: " <<
key << std::endl;
135 const auto&
title =
ME->getTitle();
136 const auto&
xtitle =
ME->getAxisTitle(1);
137 const auto&
ytitle =
ME->getAxisTitle(2);
139 const auto& nxbins =
ME->getNbinsX();
140 const auto&
xmin =
ME->getAxisMin(1);
141 const auto&
xmax =
ME->getAxisMax(1);
143 TProfile* p_mean =
new TProfile((
"Mean" +
key).c_str(),
144 (
title +
";" +
xtitle +
";#LT M_{#mu^{-}#mu^{+}} #GT [GeV]").c_str(),
150 TProfile* p_width =
new TProfile(
156 TH2F* bareHisto =
ME->getTH2F();
157 for (
int bin = 1;
bin <= nxbins;
bin++) {
158 const auto&
xaxis = bareHisto->GetXaxis();
159 const auto& low_edge =
xaxis->GetBinLowEdge(
bin);
160 const auto& high_edge =
xaxis->GetBinUpEdge(
bin);
163 edm::LogPrint(
"DiMuonMassBiasClient") <<
"dealing with bin: " <<
bin <<
" range: (" << std::setprecision(2)
164 << low_edge <<
"," << std::setprecision(2) << high_edge <<
")";
166 TH1D* Proj = bareHisto->ProjectionY(Form(
"%s_proj_%i",
key.c_str(),
bin),
bin,
bin);
167 Proj->SetTitle(Form(
"%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge,
bin));
172 edm::LogWarning(
"DiMuonMassBiasClient") <<
"the current bin has invalid data" << std::endl;
187 p_mean->SetBinEntries(
bin, 1. / (bias.
error() * bias.
error()));
190 LogDebug(
"DiMuonBassBiasClient") <<
" Bin: " <<
bin <<
" value: " << bias.
value() <<
" from profile ( " 191 << p_mean->GetBinContent(
bin) <<
") - error: " << bias.
error()
192 <<
" from profile ( " << p_mean->GetBinError(
bin) <<
" )";
199 p_width->SetBinEntries(
bin, 1. / (
width.error() *
width.error()));
219 const auto&
key = toHarvest.first;
220 const auto&
ME = toHarvest.second;
223 edm::LogPrint(
"DiMuonMassBiasClient") <<
"dealing with key: " <<
key << std::endl;
226 edm::LogError(
"DiMuonMassBiasClient") <<
"could not find MonitorElement for key: " <<
key << std::endl;
230 TH2F* bareHisto =
ME->getTH2F();
232 const auto&
xaxis = bareHisto->GetXaxis();
233 const auto& low_edge =
xaxis->GetBinLowEdge(
bin);
234 const auto& high_edge =
xaxis->GetBinUpEdge(
bin);
237 edm::LogPrint(
"DiMuonMassBiasClient") <<
"dealing with bin: " <<
bin <<
" range: (" << std::setprecision(2)
238 << low_edge <<
"," << std::setprecision(2) << high_edge <<
")";
239 TH1D* Proj = bareHisto->ProjectionY(Form(
"%s_proj_%i",
key.c_str(),
bin),
bin,
bin);
240 Proj->SetTitle(Form(
"%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge,
bin));
245 edm::LogWarning(
"DiMuonMassBiasClient") <<
"the current bin has invalid data" << std::endl;
265 edm::LogInfo(
"DiMuonMassBiasClient") <<
"DiMuonMassBiasClient::endLuminosityBlock";
290 edm::LogWarning(
"DiMuonMassBiasClient") <<
" Input histogram:" <<
hist->GetName() <<
" has not enough entries (" 291 <<
hist->GetEntries() <<
") for a meaningful Voigtian fit!\n" 297 TCanvas*
c1 =
new TCanvas();
300 c1->SetLeftMargin(0.15);
301 c1->SetRightMargin(0.10);
307 Double_t
xmin =
hist->GetXaxis()->GetXmin();
308 Double_t
xmax =
hist->GetXaxis()->GetXmax();
311 edm::LogPrint(
"DiMuonMassBiasClient") <<
"fitting range: (" <<
xmin <<
"-" <<
xmax <<
")" << std::endl;
314 RooRealVar InvMass(
"InvMass",
"di-muon mass M(#mu^{+}#mu^{-}) [GeV]",
xmin,
xmax);
315 std::unique_ptr<RooPlot>
frame{InvMass.frame()};
316 RooDataHist datahist(
"datahist",
"datahist", InvMass, RooFit::Import(*
hist));
317 datahist.plotOn(
frame.get());
320 RooRealVar
mean(
"#mu",
"mean", meanConfig_[0], meanConfig_[1], meanConfig_[2]);
321 RooRealVar
width(
"width",
"width", widthConfig_[0], widthConfig_[1], widthConfig_[2]);
322 RooRealVar sigma(
"#sigma",
"sigma", sigmaConfig_[0], sigmaConfig_[1], sigmaConfig_[2]);
323 RooVoigtian voigt(
"voigt",
"voigt", InvMass,
mean,
width, sigma);
326 RooRealVar peakCB(
"peakCB",
"peakCB", meanConfig_[0], meanConfig_[1], meanConfig_[2]);
327 RooRealVar sigmaCB(
"#sigma",
"sigma", sigmaConfig_[0], sigmaConfig_[1], sigmaConfig_[2]);
328 RooRealVar alphaCB(
"#alpha",
"alpha", 1., 0., 10.);
329 RooRealVar nCB(
"n",
"n", 1., 0., 100.);
330 RooCBShape
crystalball(
"crystalball",
"crystalball", InvMass, peakCB, sigmaCB, alphaCB, nCB);
333 RooRealVar lambda(
"#lambda",
"slope", 0., -50., 50.);
334 RooExponential expo(
"expo",
"expo", InvMass, lambda);
337 RooRealVar exp_alpha(
"#alpha",
"alpha", 40.0, 20.0, 160.0);
338 RooRealVar exp_beta(
"#beta",
"beta", 0.05, 0.0, 2.0);
339 RooRealVar exp_gamma(
"#gamma",
"gamma", 0.02, 0.0, 0.1);
340 RooRealVar exp_peak(
"peak",
"peak", meanConfig_[0]);
341 RooCMSShape exp_pdf(
"exp_pdf",
"bkg shape", InvMass, exp_alpha, exp_beta, exp_gamma, exp_peak);
344 RooRealVar
b(
"N_{b}",
"Number of background events", 0,
hist->GetEntries() / 10.);
345 RooRealVar
s(
"N_{s}",
"Number of signal events", 0,
hist->GetEntries());
349 if (useRooCBShape_) {
350 if (useRooCMSShape_) {
353 listPdf.add(exp_pdf);
360 if (useRooCMSShape_) {
363 listPdf.add(exp_pdf);
371 RooAddPdf fullModel(
"fullModel",
"Signal + Background Model", listPdf, RooArgList(
s,
b));
372 fullModel.fitTo(datahist, RooFit::PrintLevel(-1));
373 fullModel.plotOn(
frame.get(), RooFit::LineColor(kRed));
374 if (useRooCMSShape_) {
375 fullModel.plotOn(
frame.get(), RooFit::Components(exp_pdf), RooFit::LineStyle(kDashed));
377 fullModel.plotOn(
frame.get(), RooFit::Components(expo), RooFit::LineStyle(kDashed));
379 fullModel.paramOn(
frame.get(), RooFit::Layout(0.65, 0.90, 0.90));
381 if (useRooCBShape_) {
383 crystalball.fitTo(datahist, RooFit::PrintLevel(-1));
386 RooFit::Layout(0.65, 0.90, 0.90));
389 voigt.fitTo(datahist, RooFit::PrintLevel(-1));
390 voigt.plotOn(
frame.get(), RooFit::LineColor(kRed));
391 voigt.paramOn(
frame.get(), RooFit::Layout(0.65, 0.90, 0.90));
396 datahist.plotOn(
frame.get());
397 frame->GetYaxis()->SetTitle(
"n. of events");
398 TString histName =
hist->GetName();
399 frame->SetName(
"frame" + histName);
404 c1->Print(
"fit_debug" + histName +
".pdf");
408 float mass_mean = useRooCBShape_ ? peakCB.getVal() :
mean.getVal();
409 float mass_sigma = useRooCBShape_ ? sigmaCB.getVal() : sigma.getVal();
411 float mass_mean_err = useRooCBShape_ ? peakCB.getError() :
mean.getError();
412 float mass_sigma_err = useRooCBShape_ ? sigmaCB.getError() : sigma.getError();
425 edm::LogWarning(
"DiMuonMassBiasClient") <<
" Input histogram:" <<
hist->GetName() <<
" has not enough entries (" 426 <<
hist->GetEntries() <<
") for a meaningful Voigtian fit!\n" 432 TCanvas*
c1 =
new TCanvas();
435 c1->SetLeftMargin(0.15);
436 c1->SetRightMargin(0.10);
442 Double_t xMean = 91.1876;
443 Double_t
xMin =
hist->GetXaxis()->GetXmin();
444 Double_t
xMax =
hist->GetXaxis()->GetXmax();
447 edm::LogPrint(
"DiMuonMassBiasClient") <<
"fitting range: (" <<
xMin <<
"-" <<
xMax <<
")" << std::endl;
451 double sigmaMin(0.1);
452 double sigmaMax(10.);
455 double sigma2Min(0.);
456 double sigma2Max(10.);
458 std::unique_ptr<FitWithRooFit> fitter = std::make_unique<FitWithRooFit>();
462 fitter->useChi2_ = useChi2;
463 fitter->initMean(xMean,
xMin,
xMax);
464 fitter->initSigma(sigma, sigmaMin, sigmaMax);
465 fitter->initSigma2(sigma2, sigma2Min, sigma2Max);
466 fitter->initAlpha(1.5, 0.05, 10.);
467 fitter->initN(1, 0.01, 100.);
468 fitter->initFGCB(0.4, 0., 1.);
469 fitter->initGamma(2.4952, 0., 10.);
470 fitter->gamma()->setConstant(kTRUE);
471 fitter->initMean2(0., -20., 20.);
472 fitter->mean2()->setConstant(kTRUE);
473 fitter->initSigma(1.2, 0., 5.);
474 fitter->initAlpha(1.5, 0.05, 10.);
475 fitter->initN(1, 0.01, 100.);
476 fitter->initExpCoeffA0(-1., -10., 10.);
477 fitter->initExpCoeffA1(0., -10., 10.);
478 fitter->initExpCoeffA2(0., -2., 2.);
479 fitter->initFsig(0.9, 0., 1.);
480 fitter->initA0(0., -10., 10.);
481 fitter->initA1(0., -10., 10.);
482 fitter->initA2(0., -10., 10.);
483 fitter->initA3(0., -10., 10.);
484 fitter->initA4(0., -10., 10.);
485 fitter->initA5(0., -10., 10.);
486 fitter->initA6(0., -10., 10.);
488 fitter->fit(
hist,
"breitWignerTimesCB",
"exponential",
xMin,
xMax,
false);
490 TString histName =
hist->GetName();
493 c1->Print(
"fit_debug" + histName +
".pdf");
497 Measurement1D resultM(fitter->mean()->getValV(), fitter->mean()->getError());
498 Measurement1D resultW(fitter->sigma()->getValV(), fitter->sigma()->getError());
509 desc.add<
bool>(
"useTH1s",
false);
510 desc.add<
bool>(
"useBWtimesCB",
false);
511 desc.add<
bool>(
"fitBackground",
false);
512 desc.add<
bool>(
"useRooCMSShape",
false);
513 desc.add<
bool>(
"useRooCBShape",
false);
514 desc.add<
bool>(
"debugMode",
false);
517 fit_par.add<std::vector<double>>(
"mean_par",
522 fit_par.add<std::vector<double>>(
"width_par",
527 fit_par.add<std::vector<double>>(
"sigma_par",
534 desc.add<std::vector<std::string>>(
"MEtoHarvest",
535 {
"DiMuMassVsMuMuPhi",
537 "DiMuMassVsMuPlusPhi",
538 "DiMuMassVsMuPlusEta",
539 "DiMuMassVsMuMinusPhi",
540 "DiMuMassVsMuMinusEta",
541 "DiMuMassVsMuMuDeltaEta",
542 "DiMuMassVsCosThetaCS"});
void getMEsToHarvest(DQMStore::IGetter &igetter)
void addWithDefaultLabel(ParameterSetDescription const &psetDescription)
T getParameter(std::string const &) const
std::map< std::string, MonitorElement * > widthHistos_
virtual void setCurrentFolder(std::string const &fullpath)
static PFTauRenderPlugin instance
diMuonMassBias::fitOutputs fitLineShape(TH1 *hist, const bool &fitBackground=false) const
~DiMuonMassBiasClient() override
Destructor.
const std::string TopFolder_
void beginRun(edm::Run const &run, edm::EventSetup const &eSetup) override
BeginRun.
Log< level::Error, false > LogError
std::map< std::string, MonitorElement * > harvestTargets_
std::map< std::string, MonitorElement * > meanProfiles_
std::map< std::string, MonitorElement * > widthProfiles_
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
diMuonMassBias::fitOutputs fitBWTimesCB(TH1 *hist) const
#define DEFINE_FWK_MODULE(type)
std::vector< std::string > MEtoHarvest_
Log< level::Warning, true > LogPrint
void dqmEndJob(DQMStore::IBooker &ibooker_, DQMStore::IGetter &igetter_) override
EndJob.
void beginJob(void) override
BeginJob.
Log< level::Info, false > LogInfo
void fitAndFillProfile(std::pair< std::string, MonitorElement *> toHarvest, DQMStore::IBooker &iBooker)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
static constexpr int minimumHits
virtual MonitorElement * get(std::string const &fullpath) const
DiMuonMassBiasClient(const edm::ParameterSet &ps)
Constructor.
std::map< std::string, MonitorElement * > meanHistos_
void bookMEs(DQMStore::IBooker &ibooker)
book MEs
void fitAndFillHisto(std::pair< std::string, MonitorElement *> toHarvest, DQMStore::IBooker &iBooker)
const bool fitBackground_
double crystalball(double *x, double *par)
Log< level::Warning, false > LogWarning
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
void fillArrayF(float *x, const edm::ParameterSet &cfg, const char *name)