text: // -*- C++ -*-
  //
  // Package:    Validation/RecoParticleFlow
  // Class:      PFJetDQMPostProcessor.cc
  //
  // Main Developer:   "Juska Pekkanen"
  // Original Author:  "Kenichi Hatakeyama"
  //
  
  #include "FWCore/Framework/interface/LuminosityBlock.h"
  #include "FWCore/Framework/interface/Run.h"
  #include "FWCore/Framework/interface/Event.h"
  #include "FWCore/Framework/interface/MakerMacros.h"
  #include "FWCore/ParameterSet/interface/ParameterSet.h"
  #include "FWCore/ServiceRegistry/interface/Service.h"
  
  #include "DQMServices/Core/interface/DQMEDHarvester.h"
  #include "DQMServices/Core/interface/DQMStore.h"
  
  #include "TCanvas.h"
  
  //
  // class declaration
  //
  
  class PFJetDQMPostProcessor : public DQMEDHarvester {
  public:
    explicit PFJetDQMPostProcessor(const edm::ParameterSet&);
    ~PFJetDQMPostProcessor() override;
  
  private:
    void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
    void fitResponse(TH1F* hreso,
                     TH1F* h_genjet_pt,
                     int ptbinlow,
                     int ietahigh,
                     double recoptcut,
                     double& resp,
                     double& resp_err,
                     double& reso,
                     double& reso_err);
    double getRespUnc(double width, double width_err, double mean, double mean_err);
  
    std::vector<std::string> jetResponseDir;
    std::string genjetDir;
    std::vector<double> ptBins;
    std::vector<double> etaBins;
  
    double recoptcut;
  
    bool debug = false;
  
    std::string seta(double eta) {
      std::string seta = std::to_string(int(eta * 10.));
      if (seta.length() < 2)
        seta = "0" + seta;
      return seta;
    }
  
    std::string spt(double ptmin, double ptmax) {
      std::string spt = std::to_string(int(ptmin)) + "_" + std::to_string(int(ptmax));
      return spt;
    }
  };
  
  // Some switches
  //
  // constuctors and destructor
  //
  PFJetDQMPostProcessor::PFJetDQMPostProcessor(const edm::ParameterSet& iConfig) {
    jetResponseDir = iConfig.getParameter<std::vector<std::string>>("jetResponseDir");
    genjetDir = iConfig.getParameter<std::string>("genjetDir");
    ptBins = iConfig.getParameter<std::vector<double>>("ptBins");
    etaBins = iConfig.getParameter<std::vector<double>>("etaBins");
    recoptcut = iConfig.getParameter<double>("recoPtCut");
  }
  
  PFJetDQMPostProcessor::~PFJetDQMPostProcessor() {}
  
  // ------------ method called right after a run ends ------------
  void PFJetDQMPostProcessor::dqmEndJob(DQMStore::IBooker& ibook_, DQMStore::IGetter& iget_) {
    for (unsigned int idir = 0; idir < jetResponseDir.size(); idir++) {
      iget_.setCurrentFolder(genjetDir);
      std::vector<std::string> sME_genjets = iget_.getMEs();
      std::for_each(sME_genjets.begin(), sME_genjets.end(), [&](auto& s) { s.insert(0, genjetDir); });
      //for (unsigned int i=0; i<sME_genjets.size(); i++) std::cout << sME_genjets[i] << std::endl;
  
      iget_.setCurrentFolder(jetResponseDir[idir]);
      std::vector<std::string> sME_response = iget_.getMEs();
      std::for_each(sME_response.begin(), sME_response.end(), [&](auto& s) { s.insert(0, jetResponseDir[idir]); });
      //for (unsigned int i=0; i<sME_response.size(); i++) std::cout << sME_response[i] << std::endl;
  
      iget_.setCurrentFolder(jetResponseDir[idir]);
  
      double ptBinsArray[ptBins.size()];
      unsigned int nPtBins = ptBins.size() - 1;
      std::copy(ptBins.begin(), ptBins.end(), ptBinsArray);
      //for(unsigned int ipt = 0; ipt < ptBins.size(); ++ipt) std::cout << ptBins[ipt] << std::endl;
  
      std::string stitle;
      char ctitle[50];
      std::vector<MonitorElement*> vME_presponse;
      std::vector<MonitorElement*> vME_preso;
      std::vector<MonitorElement*> vME_preso_rms;
  
      MonitorElement* me;
      TH1F* h_resp;
      TH1F* h_genjet_pt;
  
      //
      // Response distributions
      //
      for (unsigned int ieta = 1; ieta < etaBins.size(); ++ieta) {
        stitle = genjetDir + "genjet_pt" + "_eta" + seta(etaBins[ieta]);
        //std::cout << ieta << " " << stitle << std::endl;
  
        std::vector<std::string>::const_iterator it = std::find(sME_genjets.begin(), sME_genjets.end(), stitle);
        if (it == sME_genjets.end())
          continue;
        me = iget_.get(stitle);
        h_genjet_pt = (TH1F*)me->getTH1F();
  
        stitle = "presponse_eta" + seta(etaBins[ieta]);
        // adding "Raw" to the title of raw jet response histograms
        if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
          sprintf(ctitle, "Raw Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
        } else {
          sprintf(ctitle, "Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
        }
        TH1F* h_presponse = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
  
        stitle = "preso_eta" + seta(etaBins[ieta]);
        if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
          sprintf(ctitle, "Raw Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
        } else {
          sprintf(ctitle, "Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
        }
        TH1F* h_preso = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
  
        stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
        if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
          sprintf(ctitle, "Raw Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
        } else {
          sprintf(ctitle, "Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
        }
        TH1F* h_preso_rms = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
  
        for (unsigned int ipt = 0; ipt < ptBins.size() - 1; ++ipt) {
          stitle = jetResponseDir[idir] + "reso_dist_" + spt(ptBins[ipt], ptBins[ipt + 1]) + "_eta" + seta(etaBins[ieta]);
          std::vector<std::string>::const_iterator it = std::find(sME_response.begin(), sME_response.end(), stitle);
          if (it == sME_response.end())
            continue;
          me = iget_.get(stitle);
          h_resp = (TH1F*)me->getTH1F();
  
          // Fit-based
          double resp = 1.0, resp_err = 0.0, reso = 0.0, reso_err = 0.0;
          fitResponse(h_resp, h_genjet_pt, ipt, ieta, recoptcut, resp, resp_err, reso, reso_err);
  
          h_presponse->SetBinContent(ipt + 1, resp);
          h_presponse->SetBinError(ipt + 1, resp_err);
          h_preso->SetBinContent(ipt + 1, reso);
          h_preso->SetBinError(ipt + 1, reso_err);
  
          // RMS-based
          double std = h_resp->GetStdDev();
          double std_error = h_resp->GetStdDevError();
  
          // Scale each bin with mean response
          double mean = 1.0;
          double mean_error = 0.0;
          double err = 0.0;
          if (h_resp->GetMean() > 0) {
            mean = h_resp->GetMean();
            mean_error = h_resp->GetMeanError();
  
            // Scale resolution by response.
            std /= mean;
            std_error /= mean;
  
            err = getRespUnc(std, std_error, mean, mean_error);
          }
  
          h_preso_rms->SetBinContent(ipt + 1, std);
          h_preso_rms->SetBinError(ipt + 1, err);
  
        }  // ipt
  
        stitle = "presponse_eta" + seta(etaBins[ieta]);
        me = ibook_.book1D(stitle.c_str(), h_presponse);
        vME_presponse.push_back(me);
  
        stitle = "preso_eta" + seta(etaBins[ieta]);
        me = ibook_.book1D(stitle.c_str(), h_preso);
        vME_preso.push_back(me);
  
        stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
        me = ibook_.book1D(stitle.c_str(), h_preso_rms);
        vME_preso_rms.push_back(me);
  
      }  // ieta
  
      //
      // Checks
      //
      if (debug) {
        for (std::vector<MonitorElement*>::const_iterator i = vME_presponse.begin(); i != vME_presponse.end(); ++i)
          (*i)->getTH1F()->Print();
        for (std::vector<MonitorElement*>::const_iterator i = vME_preso.begin(); i != vME_preso.end(); ++i)
          (*i)->getTH1F()->Print();
        for (std::vector<MonitorElement*>::const_iterator i = vME_preso_rms.begin(); i != vME_preso_rms.end(); ++i)
          (*i)->getTH1F()->Print();
      }
    }
  }
  
  void PFJetDQMPostProcessor::fitResponse(TH1F* hreso,
                                          TH1F* h_genjet_pt,
                                          int ptbinlow,
                                          int ietahigh,
                                          double recoptcut,
                                          double& resp,
                                          double& resp_err,
                                          double& reso,
                                          double& reso_err) {
    // This 'smartfitter' is converted from the original Python smart_fit() -function
    // implemented in test/helperFunctions.py. See that file for more commentary.
    // Juska 23 May 2019
  
    // Only do plots if needed for debugging
    // NOTE a directory called 'debug' must exist in the working directory
    //
    bool doPlots = false;
  
    double ptlow = ptBins[ptbinlow];
    double pthigh = ptBins[ptbinlow + 1];
  
    // Take range by Mikko's advice: -1.5 and + 1.5 * RMS width
  
    double rmswidth = hreso->GetStdDev();
    double rmsmean = hreso->GetMean();
    double fitlow = rmsmean - 1.5 * rmswidth;
    fitlow = TMath::Max(recoptcut / ptlow, fitlow);
    double fithigh = rmsmean + 1.5 * rmswidth;
  
    TF1* fg = new TF1("mygaus", "gaus", fitlow, fithigh);
    TF1* fg2 = new TF1("fg2", "TMath::Gaus(x,[0],[1],true)*[2]", fitlow, fithigh);
  
    hreso->Fit("mygaus", "RQN");
  
    fg2->SetParameter(0, fg->GetParameter(1));
    fg2->SetParameter(1, fg->GetParameter(2));
  
    // Extract ngenjet in the current pT bin from the genjet histo
    float ngenjet = h_genjet_pt->GetBinContent(ptbinlow + 1);
  
    // Here the fit is forced to take the area of ngenjets.
    // The area is further normalized for the response histogram x-axis length
    // (3) and number of bins (100)
    fg2->FixParameter(2, ngenjet * 3. / 100.);
  
    hreso->Fit("fg2", "RQN");
  
    fitlow = fg2->GetParameter(0) - 1.5 * fg2->GetParameter(1);
    fitlow = TMath::Max(15. / ptlow, fitlow);
    fithigh = fg2->GetParameter(0) + 1.5 * fg2->GetParameter(1);
  
    fg2->SetRange(fitlow, fithigh);
  
    hreso->Fit("fg2", "RQ");
  
    fg->SetRange(0, 3);
    fg2->SetRange(0, 3);
    fg->SetLineWidth(2);
    fg2->SetLineColor(kGreen + 2);
  
    hreso->GetXaxis()->SetRangeUser(0, 2);
  
    // Save plots to a subdirectory if asked (debug-directory must exist!)
    if (doPlots & (hreso->GetEntries() > 0)) {
      TCanvas* cfit = new TCanvas(Form("respofit_%i", int(ptlow)), "respofit", 600, 600);
      hreso->Draw("ehist");
      fg->Draw("same");
      fg2->Draw("same");
      cfit->SaveAs(
          Form("debug/respo_smartfit_%04d_%i_eta%s.pdf", (int)ptlow, (int)pthigh, seta(etaBins[ietahigh]).c_str()));
    }
  
    resp = fg2->GetParameter(0);
    resp_err = fg2->GetParError(0);
  
    // Scale resolution by response. Avoid division by zero.
    if (0 == resp)
      reso = 0;
    else
      reso = fg2->GetParameter(1) / resp;
    reso_err = fg2->GetParError(1) / resp;
  
    reso_err = getRespUnc(reso, reso_err, resp, resp_err);
  }
  
  // Calculate resolution uncertainty
  double PFJetDQMPostProcessor::getRespUnc(double width, double width_err, double mean, double mean_err) {
    if (0 == width || 0 == mean)
      return 0;
    return TMath::Sqrt(pow(width_err, 2) / pow(width, 2) + pow(mean_err, 2) / pow(mean, 2)) * width;
  }
  
  #include "FWCore/Framework/interface/MakerMacros.h"
  DEFINE_FWK_MODULE(PFJetDQMPostProcessor);
  