Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 10:57:30

0001 // user includes
0002 #include "DQMOffline/Alignment/interface/DiMuonMassBiasClient.h"
0003 #include "DataFormats/Histograms/interface/DQMToken.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "PhysicsTools/TagAndProbe/interface/RooCMSShape.h"
0007 
0008 // RooFit includes
0009 #include "TCanvas.h"
0010 #include "RooAddPdf.h"
0011 #include "RooDataHist.h"
0012 #include "RooExponential.h"
0013 #include "RooGaussian.h"
0014 #include "RooPlot.h"
0015 #include "RooRealVar.h"
0016 #include "RooVoigtian.h"
0017 #include "RooCBShape.h"
0018 
0019 //-----------------------------------------------------------------------------------
0020 DiMuonMassBiasClient::DiMuonMassBiasClient(edm::ParameterSet const& iConfig)
0021     : TopFolder_(iConfig.getParameter<std::string>("FolderName")),
0022       fitBackground_(iConfig.getParameter<bool>("fitBackground")),
0023       useRooCBShape_(iConfig.getParameter<bool>("useRooCBShape")),
0024       useRooCMSShape_(iConfig.getParameter<bool>("useRooCMSShape")),
0025       debugMode_(iConfig.getParameter<bool>("debugMode")),
0026       MEtoHarvest_(iConfig.getParameter<std::vector<std::string>>("MEtoHarvest"))
0027 //-----------------------------------------------------------------------------------
0028 {
0029   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::Constructing DiMuonMassBiasClient ";
0030 
0031   // fill the parameters for the fit
0032   edm::ParameterSet fit_par = iConfig.getParameter<edm::ParameterSet>("fit_par");
0033   diMuonMassBias::fillArrayF(meanConfig_, fit_par, "mean_par");
0034   diMuonMassBias::fillArrayF(widthConfig_, fit_par, "width_par");
0035   diMuonMassBias::fillArrayF(sigmaConfig_, fit_par, "sigma_par");
0036 
0037   if (debugMode_) {
0038     edm::LogPrint("DiMuonMassBiasClient")
0039         << "mean: " << meanConfig_[0] << " (" << meanConfig_[1] << "," << meanConfig_[2] << ") " << std::endl;
0040     edm::LogPrint("DiMuonMassBiasClient")
0041         << "width: " << widthConfig_[0] << " (" << widthConfig_[1] << "," << widthConfig_[2] << ")" << std::endl;
0042     edm::LogPrint("DiMuonMassBiasClient")
0043         << "sigma: " << sigmaConfig_[0] << " (" << sigmaConfig_[1] << "," << sigmaConfig_[2] << ")" << std::endl;
0044   }
0045 }
0046 
0047 //-----------------------------------------------------------------------------------
0048 DiMuonMassBiasClient::~DiMuonMassBiasClient()
0049 //-----------------------------------------------------------------------------------
0050 {
0051   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::Deleting DiMuonMassBiasClient ";
0052 }
0053 
0054 //-----------------------------------------------------------------------------------
0055 void DiMuonMassBiasClient::beginJob(void)
0056 //-----------------------------------------------------------------------------------
0057 {
0058   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::beginJob done";
0059 }
0060 
0061 //-----------------------------------------------------------------------------------
0062 void DiMuonMassBiasClient::beginRun(edm::Run const& run, edm::EventSetup const& eSetup)
0063 //-----------------------------------------------------------------------------------
0064 {
0065   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient:: Begining of Run";
0066 }
0067 
0068 //-----------------------------------------------------------------------------------
0069 void DiMuonMassBiasClient::bookMEs(DQMStore::IBooker& iBooker)
0070 //-----------------------------------------------------------------------------------
0071 {
0072   iBooker.setCurrentFolder(TopFolder_ + "/DiMuonMassBiasMonitor/MassBias/Profiles");
0073   for (const auto& [key, ME] : harvestTargets_) {
0074     if (ME == nullptr) {
0075       edm::LogError("DiMuonMassBiasClient") << "could not find MonitorElement for key: " << key << std::endl;
0076       continue;
0077     }
0078 
0079     const auto& title = ME->getTitle();
0080     const auto& xtitle = ME->getAxisTitle(1);
0081     const auto& ytitle = ME->getAxisTitle(2);
0082 
0083     const auto& nxbins = ME->getNbinsX();
0084     const auto& xmin = ME->getAxisMin(1);
0085     const auto& xmax = ME->getAxisMax(1);
0086 
0087     MonitorElement* meanToBook =
0088         iBooker.book1D(("Mean" + key), (title + ";" + xtitle + ";" + ytitle), nxbins, xmin, xmax);
0089     meanProfiles_.insert({key, meanToBook});
0090 
0091     MonitorElement* sigmaToBook =
0092         iBooker.book1D(("Sigma" + key), (title + ";" + xtitle + ";" + "#sigma of " + ytitle), nxbins, xmin, xmax);
0093     widthProfiles_.insert({key, sigmaToBook});
0094   }
0095 }
0096 
0097 //-----------------------------------------------------------------------------------
0098 void DiMuonMassBiasClient::getMEsToHarvest(DQMStore::IGetter& iGetter)
0099 //-----------------------------------------------------------------------------------
0100 {
0101   std::string inFolder = TopFolder_ + "/DiMuonMassBiasMonitor/MassBias/";
0102 
0103   //loop on the list of histograms to harvest
0104   for (const auto& hname : MEtoHarvest_) {
0105     MonitorElement* toHarvest = iGetter.get(inFolder + hname);
0106 
0107     if (toHarvest == nullptr) {
0108       edm::LogError("DiMuonMassBiasClient") << "could not find input MonitorElement: " << inFolder + hname << std::endl;
0109       continue;
0110     }
0111 
0112     harvestTargets_.insert({hname, toHarvest});
0113   }
0114 }
0115 
0116 //-----------------------------------------------------------------------------------
0117 void DiMuonMassBiasClient::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter)
0118 //-----------------------------------------------------------------------------------
0119 {
0120   edm::LogInfo("DiMuonMassBiasClient") << "DiMuonMassBiasClient::endLuminosityBlock";
0121 
0122   getMEsToHarvest(igetter);
0123   bookMEs(ibooker);
0124 
0125   for (const auto& [key, ME] : harvestTargets_) {
0126     if (debugMode_)
0127       edm::LogPrint("DiMuonMassBiasClient") << "dealing with key: " << key << std::endl;
0128     TH2F* bareHisto = ME->getTH2F();
0129     for (int bin = 1; bin <= ME->getNbinsX(); bin++) {
0130       const auto& xaxis = bareHisto->GetXaxis();
0131       const auto& low_edge = xaxis->GetBinLowEdge(bin);
0132       const auto& high_edge = xaxis->GetBinUpEdge(bin);
0133 
0134       if (debugMode_)
0135         edm::LogPrint("DiMuonMassBiasClient") << "dealing with bin: " << bin << " range: (" << std::setprecision(2)
0136                                               << low_edge << "," << std::setprecision(2) << high_edge << ")";
0137       TH1D* Proj = bareHisto->ProjectionY(Form("%s_proj_%i", key.c_str(), bin), bin, bin);
0138       Proj->SetTitle(Form("%s #in (%.2f,%.2f), bin: %i", Proj->GetTitle(), low_edge, high_edge, bin));
0139 
0140       diMuonMassBias::fitOutputs results = fitLineShape(Proj);
0141 
0142       if (results.isInvalid()) {
0143         edm::LogWarning("DiMuonMassBiasClient") << "the current bin has invalid data" << std::endl;
0144         continue;
0145       }
0146 
0147       // fill the mean profiles
0148       const Measurement1D& bias = results.getBias();
0149       meanProfiles_[key]->setBinContent(bin, bias.value());
0150       meanProfiles_[key]->setBinError(bin, bias.error());
0151 
0152       // fill the width profiles
0153       const Measurement1D& width = results.getWidth();
0154       widthProfiles_[key]->setBinContent(bin, width.value());
0155       widthProfiles_[key]->setBinError(bin, width.error());
0156     }
0157   }
0158 }
0159 
0160 //-----------------------------------------------------------------------------------
0161 diMuonMassBias::fitOutputs DiMuonMassBiasClient::fitLineShape(TH1* hist, const bool& fitBackground) const
0162 //-----------------------------------------------------------------------------------
0163 {
0164   if (hist->GetEntries() < diMuonMassBias::minimumHits) {
0165     edm::LogWarning("DiMuonMassBiasClient") << " Input histogram:" << hist->GetName() << " has not enough entries ("
0166                                             << hist->GetEntries() << ") for a meaningful Voigtian fit!\n"
0167                                             << "Skipping!";
0168 
0169     return diMuonMassBias::fitOutputs(Measurement1D(0., 0.), Measurement1D(0., 0.));
0170   }
0171 
0172   TCanvas* c1 = new TCanvas();
0173   if (debugMode_) {
0174     c1->Clear();
0175     c1->SetLeftMargin(0.15);
0176     c1->SetRightMargin(0.10);
0177   }
0178 
0179   // silence messages
0180   RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
0181 
0182   Double_t xmin = hist->GetXaxis()->GetXmin();
0183   Double_t xmax = hist->GetXaxis()->GetXmax();
0184 
0185   if (debugMode_) {
0186     edm::LogPrint("DiMuonMassBiasClient") << "fitting range: (" << xmin << "-" << xmax << ")" << std::endl;
0187   }
0188 
0189   RooRealVar InvMass("InvMass", "di-muon mass M(#mu^{+}#mu^{-}) [GeV]", xmin, xmax);
0190   std::unique_ptr<RooPlot> frame{InvMass.frame()};
0191   RooDataHist datahist("datahist", "datahist", InvMass, RooFit::Import(*hist));
0192   datahist.plotOn(frame.get());
0193 
0194   // parameters of the Voigtian
0195   RooRealVar mean("#mu", "mean", meanConfig_[0], meanConfig_[1], meanConfig_[2]);          //90.0, 60.0, 120.0 (for Z)
0196   RooRealVar width("width", "width", widthConfig_[0], widthConfig_[1], widthConfig_[2]);   // 5.0,  0.0, 120.0 (for Z)
0197   RooRealVar sigma("#sigma", "sigma", sigmaConfig_[0], sigmaConfig_[1], sigmaConfig_[2]);  // 5.0,  0.0, 120.0 (for Z)
0198   RooVoigtian voigt("voigt", "voigt", InvMass, mean, width, sigma);
0199 
0200   // parameters of the Crystal-ball
0201   RooRealVar peakCB("peakCB", "peakCB", meanConfig_[0], meanConfig_[1], meanConfig_[2]);
0202   RooRealVar sigmaCB("#sigma", "sigma", sigmaConfig_[0], sigmaConfig_[1], sigmaConfig_[2]);
0203   RooRealVar alphaCB("#alpha", "alpha", 1., 0., 10.);
0204   RooRealVar nCB("n", "n", 1., 0., 100.);
0205   RooCBShape crystalball("crystalball", "crystalball", InvMass, peakCB, sigmaCB, alphaCB, nCB);
0206 
0207   // for the simple background fit
0208   RooRealVar lambda("#lambda", "slope", 0., -50., 50.);
0209   RooExponential expo("expo", "expo", InvMass, lambda);
0210 
0211   // for the more refined background fit
0212   RooRealVar exp_alpha("#alpha", "alpha", 40.0, 20.0, 160.0);
0213   RooRealVar exp_beta("#beta", "beta", 0.05, 0.0, 2.0);
0214   RooRealVar exp_gamma("#gamma", "gamma", 0.02, 0.0, 0.1);
0215   RooRealVar exp_peak("peak", "peak", meanConfig_[0]);
0216   RooCMSShape exp_pdf("exp_pdf", "bkg shape", InvMass, exp_alpha, exp_beta, exp_gamma, exp_peak);
0217 
0218   // define the signal and background fractions
0219   RooRealVar b("N_{b}", "Number of background events", 0, hist->GetEntries() / 10.);
0220   RooRealVar s("N_{s}", "Number of signal events", 0, hist->GetEntries());
0221 
0222   if (fitBackground_) {
0223     RooArgList listPdf;
0224     if (useRooCBShape_) {
0225       if (useRooCMSShape_) {
0226         // crystal-ball + CMS-shape fit
0227         listPdf.add(crystalball);
0228         listPdf.add(exp_pdf);
0229       } else {
0230         // crystal-ball + exponential fit
0231         listPdf.add(crystalball);
0232         listPdf.add(expo);
0233       }
0234     } else {
0235       if (useRooCMSShape_) {
0236         // voigtian + CMS-shape fit
0237         listPdf.add(voigt);
0238         listPdf.add(exp_pdf);
0239       } else {
0240         // voigtian + exponential fit
0241         listPdf.add(voigt);
0242         listPdf.add(expo);
0243       }
0244     }
0245 
0246     RooAddPdf fullModel("fullModel", "Signal + Background Model", listPdf, RooArgList(s, b));
0247     fullModel.fitTo(datahist, RooFit::PrintLevel(-1));
0248     fullModel.plotOn(frame.get(), RooFit::LineColor(kRed));
0249     if (useRooCMSShape_) {
0250       fullModel.plotOn(frame.get(), RooFit::Components(exp_pdf), RooFit::LineStyle(kDashed));  //Other option
0251     } else {
0252       fullModel.plotOn(frame.get(), RooFit::Components(expo), RooFit::LineStyle(kDashed));  //Other option
0253     }
0254     fullModel.paramOn(frame.get(), RooFit::Layout(0.65, 0.90, 0.90));
0255   } else {
0256     if (useRooCBShape_) {
0257       // use crystal-ball for a fit-only signal
0258       crystalball.fitTo(datahist, RooFit::PrintLevel(-1));
0259       crystalball.plotOn(frame.get(), RooFit::LineColor(kRed));  //this will show fit overlay on canvas
0260       crystalball.paramOn(frame.get(),
0261                           RooFit::Layout(0.65, 0.90, 0.90));  //this will display the fit parameters on canvas
0262     } else {
0263       // use voigtian for a fit-only signal
0264       voigt.fitTo(datahist, RooFit::PrintLevel(-1));
0265       voigt.plotOn(frame.get(), RooFit::LineColor(kRed));            //this will show fit overlay on canvas
0266       voigt.paramOn(frame.get(), RooFit::Layout(0.65, 0.90, 0.90));  //this will display the fit parameters on canvas
0267     }
0268   }
0269 
0270   // Redraw data on top and print / store everything
0271   datahist.plotOn(frame.get());
0272   frame->GetYaxis()->SetTitle("n. of events");
0273   TString histName = hist->GetName();
0274   frame->SetName("frame" + histName);
0275   frame->SetTitle(hist->GetTitle());
0276   frame->Draw();
0277 
0278   if (debugMode_) {
0279     c1->Print("fit_debug" + histName + ".pdf");
0280   }
0281   delete c1;
0282 
0283   float mass_mean = useRooCBShape_ ? peakCB.getVal() : mean.getVal();
0284   float mass_sigma = useRooCBShape_ ? sigmaCB.getVal() : sigma.getVal();
0285 
0286   float mass_mean_err = useRooCBShape_ ? peakCB.getError() : mean.getError();
0287   float mass_sigma_err = useRooCBShape_ ? sigmaCB.getError() : sigma.getError();
0288 
0289   Measurement1D resultM(mass_mean, mass_mean_err);
0290   Measurement1D resultW(mass_sigma, mass_sigma_err);
0291 
0292   return diMuonMassBias::fitOutputs(resultM, resultW);
0293 }
0294 
0295 //-----------------------------------------------------------------------------------
0296 void DiMuonMassBiasClient::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
0297 //-----------------------------------------------------------------------------------
0298 {
0299   edm::ParameterSetDescription desc;
0300   desc.add<std::string>("FolderName", "DiMuonMassBiasMonitor");
0301   desc.add<bool>("fitBackground", false);
0302   desc.add<bool>("useRooCMSShape", false);
0303   desc.add<bool>("useRooCBShape", false);
0304   desc.add<bool>("debugMode", false);
0305 
0306   edm::ParameterSetDescription fit_par;
0307   fit_par.add<std::vector<double>>("mean_par",
0308                                    {std::numeric_limits<float>::max(),
0309                                     std::numeric_limits<float>::max(),
0310                                     std::numeric_limits<float>::max()});  // par = mean
0311 
0312   fit_par.add<std::vector<double>>("width_par",
0313                                    {std::numeric_limits<float>::max(),
0314                                     std::numeric_limits<float>::max(),
0315                                     std::numeric_limits<float>::max()});  // par = width
0316 
0317   fit_par.add<std::vector<double>>("sigma_par",
0318                                    {std::numeric_limits<float>::max(),
0319                                     std::numeric_limits<float>::max(),
0320                                     std::numeric_limits<float>::max()});  // par = sigma
0321 
0322   desc.add<edm::ParameterSetDescription>("fit_par", fit_par);
0323 
0324   desc.add<std::vector<std::string>>("MEtoHarvest",
0325                                      {"DiMuMassVsMuMuPhi",
0326                                       "DiMuMassVsMuMuEta",
0327                                       "DiMuMassVsMuPlusPhi",
0328                                       "DiMuMassVsMuPlusEta",
0329                                       "DiMuMassVsMuMinusPhi",
0330                                       "DiMuMassVsMuMinusEta",
0331                                       "DiMuMassVsMuMuDeltaEta",
0332                                       "DiMuMassVsCosThetaCS"});
0333   descriptions.addWithDefaultLabel(desc);
0334 }
0335 
0336 DEFINE_FWK_MODULE(DiMuonMassBiasClient);