// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"

// Geometry
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
#include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"

#include "SimCalorimetry/HGCalSimAlgos/interface/HGCalSiNoiseMap.h"

#include <TH1F.h>
#include <TH2F.h>
#include <TProfile.h>

#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

//STL headers
#include <vector>
#include <sstream>
#include <string>
#include <iomanip>

#include "vdt/vdtMath.h"

using namespace std;

//
// class declaration
//
class HFNoseNoiseMapAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
  explicit HFNoseNoiseMapAnalyzer(const edm::ParameterSet &);
  ~HFNoseNoiseMapAnalyzer() override;

private:
  void beginJob() override {}
  void analyze(const edm::Event &, const edm::EventSetup &) override;
  void endJob() override {}

  // ----------member data ---------------------------
  const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokGeom_;
  edm::Service<TFileService> fs_;
  std::map<DetId::Detector, std::unique_ptr<HGCalSiNoiseMap<HFNoseDetId>>> noiseMaps_;
  std::map<std::pair<DetId::Detector, int>, TH1F *> layerN_, layerCCE_, layerNoise_, layerIleak_, layerSN_, layerF_,
      layerGain_, layerMipPeak_;
  std::map<DetId::Detector, TH2F *> detN_, detCCE_, detNoise_, detIleak_, detSN_, detF_, detGain_, detMipPeak_;
  std::map<std::pair<DetId::Detector, int>, TProfile *> detCCEVsFluence_;

  int aimMIPtoADC_;
  bool ignoreGainSettings_;

  const int plotMargin_ = 20;
};

//
HFNoseNoiseMapAnalyzer::HFNoseNoiseMapAnalyzer(const edm::ParameterSet &iConfig)
    : tokGeom_(esConsumes<CaloGeometry, CaloGeometryRecord>()) {
  usesResource("TFileService");
  fs_->file().cd();

  //configure the dose map
  //  std::string doseMapURL(iConfig.getParameter<std::string>("doseMap"));
  std::string doseMapURL("SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb_fluka_HFNose_3.7.20.12_Eta2.4.txt");
  unsigned int doseMapAlgo(iConfig.getParameter<unsigned int>("doseMapAlgo"));
  double scaleByDoseFactor(iConfig.getParameter<double>("scaleByDoseFactor"));
  std::vector<double> ileakParam(
      iConfig.getParameter<edm::ParameterSet>("ileakParam").template getParameter<std::vector<double>>("ileakParam"));
  std::vector<double> cceParamFine(
      iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamFine"));
  std::vector<double> cceParamThin(
      iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThin"));
  std::vector<double> cceParamThick(
      iConfig.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThick"));

  noiseMaps_[DetId::Forward] = std::unique_ptr<HGCalSiNoiseMap<HFNoseDetId>>(new HGCalSiNoiseMap<HFNoseDetId>);
  noiseMaps_[DetId::Forward]->setDoseMap(doseMapURL, doseMapAlgo);
  noiseMaps_[DetId::Forward]->setFluenceScaleFactor(scaleByDoseFactor);
  noiseMaps_[DetId::Forward]->setIleakParam(ileakParam);
  noiseMaps_[DetId::Forward]->setCceParam(cceParamFine, cceParamThin, cceParamThick);

  aimMIPtoADC_ = iConfig.getParameter<int>("aimMIPtoADC");
  ignoreGainSettings_ = iConfig.getParameter<bool>("ignoreGainSettings");
}

//
HFNoseNoiseMapAnalyzer::~HFNoseNoiseMapAnalyzer() {}

//
void HFNoseNoiseMapAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &es) {
  //get geometry
  const auto &geom = es.getHandle(tokGeom_);

  std::vector<DetId::Detector> dets = {DetId::Forward};
  for (const auto &d : dets) {
    noiseMaps_[d]->setGeometry(geom->getSubdetectorGeometry(d, ForwardSubdetector::HFNose));
    //sub-detector boundaries
    unsigned int nlay = noiseMaps_[d]->ddd()->layers(true);
    std::pair<double, double> ranZ = noiseMaps_[d]->ddd()->rangeZ(true);
    std::pair<double, double> ranRAtZ = noiseMaps_[d]->ddd()->rangeR(ranZ.first, true);
    std::pair<double, double> ranR(ranRAtZ.first - plotMargin_, ranRAtZ.second + plotMargin_);

    const std::vector<DetId> &detIdVec = noiseMaps_[d]->geom()->getValidDetIds();
    /*
    cout << "Subdetector:" << d << " has " << detIdVec.size() << " valid cells" << endl
         << "\t" << ranR.first << "<r<" << ranR.second << "\t" << ranZ.first << "<z<" << ranZ.second << endl;
    */

    //start histos
    TString baseName(Form("d%d_", d));
    TString title("HFNose");
    Int_t nbinsR(100);
    for (unsigned int ilay = 0; ilay < nlay; ilay++) {
      //this layer histos
      int layer(ilay + 1);
      std::pair<DetId::Detector, int> key(d, layer);
      TString layerBaseName(Form("%slayer%d_", baseName.Data(), layer));
      TString layerTitle(Form("%s %d", title.Data(), layer));
      layerTitle += ";Radius [cm];";
      layerN_[key] = fs_->make<TH1F>(layerBaseName + "ncells", layerTitle + "Cells", nbinsR, ranR.first, ranR.second);
      layerCCE_[key] = fs_->make<TH1F>(layerBaseName + "cce", layerTitle + "<CCE>", nbinsR, ranR.first, ranR.second);
      layerNoise_[key] =
          fs_->make<TH1F>(layerBaseName + "noise", layerTitle + "<Noise> [fC]", nbinsR, ranR.first, ranR.second);
      layerIleak_[key] =
          fs_->make<TH1F>(layerBaseName + "ileak", layerTitle + "<I_{leak}> [#muA]", nbinsR, ranR.first, ranR.second);
      layerSN_[key] = fs_->make<TH1F>(layerBaseName + "sn", layerTitle + "<S/N>", nbinsR, ranR.first, ranR.second);
      layerF_[key] = fs_->make<TH1F>(
          layerBaseName + "fluence", layerTitle + "<F> [n_{eq}/cm^{2}]", nbinsR, ranR.first, ranR.second);
      layerGain_[key] = fs_->make<TH1F>(layerBaseName + "gain", layerTitle + "<Gain>", nbinsR, ranR.first, ranR.second);
      layerMipPeak_[key] =
          fs_->make<TH1F>(layerBaseName + "mippeak", layerTitle + "<MIP peak> [ADC]", nbinsR, ranR.first, ranR.second);
    }

    //cce vs fluence
    for (unsigned int wafertype = 0; wafertype < (noiseMaps_[d]->getMipEqfC()).size(); ++wafertype) {
      std::pair<DetId::Detector, int> key2(d, wafertype);
      detCCEVsFluence_[key2] = fs_->make<TProfile>(
          baseName + Form("wt%d_", wafertype) + "cceVsFluence", title + ";<F> [n_{eq}/cm^{2}];<CCE>", 1000, 1e14, 1e16);
    }

    //sub-detector histos
    title += ";Layer;Radius [cm];";
    detN_[d] =
        fs_->make<TH2F>(baseName + "ncells", title + "Cells", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
    detCCE_[d] = fs_->make<TH2F>(baseName + "cce", title + "<CCE>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
    detNoise_[d] =
        fs_->make<TH2F>(baseName + "noise", title + "<Noise> [fC]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
    detIleak_[d] = fs_->make<TH2F>(
        baseName + "ileak", title + "<I_{leak}> [#muA]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
    detSN_[d] = fs_->make<TH2F>(baseName + "sn", title + "<S/N>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
    detF_[d] = fs_->make<TH2F>(
        baseName + "fluence", title + "<F> [n_{eq}/cm^{2}]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
    detGain_[d] =
        fs_->make<TH2F>(baseName + "gain", title + "<Gain>", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);
    detMipPeak_[d] = fs_->make<TH2F>(
        baseName + "mippeak", title + "<MIP peak> [ADC]", nlay, 1, nlay + 1, nbinsR, ranR.first, ranR.second);

    //fill histos
    for (const auto &cellId : detIdVec) {
      HFNoseDetId id(cellId.rawId());
      int layer = std::abs(id.layer());
      GlobalPoint pt = noiseMaps_[d]->geom()->getPosition(id);
      double r(pt.perp());

      HGCalSiNoiseMap<HFNoseDetId>::GainRange_t gainToSet(HGCalSiNoiseMap<HFNoseDetId>::AUTO);
      if (ignoreGainSettings_)
        gainToSet = HGCalSiNoiseMap<HFNoseDetId>::q80fC;
      HGCalSiNoiseMap<HFNoseDetId>::SiCellOpCharacteristics siop =
          noiseMaps_[d]->getSiCellOpCharacteristics(id, gainToSet, aimMIPtoADC_);

      //fill histos (layer,radius)
      detN_[d]->Fill(layer, r, 1);
      detCCE_[d]->Fill(layer, r, siop.core.cce);
      detNoise_[d]->Fill(layer, r, siop.core.noise);
      detSN_[d]->Fill(layer, r, siop.mipfC / siop.core.noise);
      detIleak_[d]->Fill(layer, r, siop.ileak);
      detF_[d]->Fill(layer, r, siop.fluence);
      detGain_[d]->Fill(layer, r, siop.core.gain + 1);
      detMipPeak_[d]->Fill(layer, r, siop.mipADC);

      //per layer histograms
      std::pair<DetId::Detector, int> key(d, layer);
      layerN_[key]->Fill(r, 1);
      layerCCE_[key]->Fill(r, siop.core.cce);
      layerNoise_[key]->Fill(r, siop.core.noise);
      layerSN_[key]->Fill(r, siop.mipfC / siop.core.noise);
      layerIleak_[key]->Fill(r, siop.ileak);
      layerF_[key]->Fill(r, siop.fluence);
      layerGain_[key]->Fill(r, siop.core.gain + 1);
      layerMipPeak_[key]->Fill(r, siop.mipADC);

      std::pair<DetId::Detector, int> key2(d, id.type());
      detCCEVsFluence_[key2]->Fill(siop.fluence, siop.core.cce);
    }

    //normalize histos per cell counts
    detF_[d]->Divide(detN_[d]);
    detCCE_[d]->Divide(detN_[d]);
    detNoise_[d]->Divide(detN_[d]);
    detSN_[d]->Divide(detN_[d]);
    detIleak_[d]->Divide(detN_[d]);
    detGain_[d]->Divide(detN_[d]);
    detMipPeak_[d]->Divide(detN_[d]);
    for (unsigned int ilay = 0; ilay < nlay; ilay++) {
      int layer(ilay + 1);
      std::pair<DetId::Detector, int> key(d, layer);
      layerF_[key]->Divide(layerN_[key]);
      layerCCE_[key]->Divide(layerN_[key]);
      layerNoise_[key]->Divide(layerN_[key]);
      layerSN_[key]->Divide(layerN_[key]);
      layerIleak_[key]->Divide(layerN_[key]);
      layerGain_[key]->Divide(layerN_[key]);
      layerMipPeak_[key]->Divide(layerN_[key]);
    }
  }
}

//define this as a plug-in
DEFINE_FWK_MODULE(HFNoseNoiseMapAnalyzer);
