HcalObjRepresent.h

Go to the documentation of this file.

#ifndef HcalObjRepresent_h
#define HcalObjRepresent_h

#include "CondFormats/HcalObjects/interface/HcalGains.h"

//#include "CondCore/Utilities/interface/PayLoadInspector.h"
//#include "CondCore/Utilities/interface/InspectorPythonWrapper.h"
#include "CondFormats/HcalObjects/interface/HcalDetIdRelationship.h"
#include "CondFormats/HcalObjects/interface/HcalCondObjectContainer.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/HcalDetId/interface/HcalOtherDetId.h"
#include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
#include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
#include "FWCore/Utilities/interface/Exception.h"

#include <string>
#include <fstream>
#include <sstream>
#include <vector>

#include "TH1F.h"
#include "TH2F.h"
#include "DataFormats/DetId/interface/DetId.h"

#include "TROOT.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TColor.h"
#include "TLine.h"
#include "TLatex.h"
#include "TProfile.h"
#include "TPaveLabel.h"

//functions for correct representation of data in summary and plot
namespace HcalObjRepresent {

  //used to produce all display objects for payload inspector
  template <class Items, class Item>
  class HcalDataContainer {
  public:
    HcalDataContainer(std::shared_ptr<Items> payload, unsigned int run) : payload_(payload), run_(run) {
      PlotMode_ = "Map";
    }

    virtual ~HcalDataContainer(){};
    // For easier channel mapping
    typedef std::tuple<int, int, int> Coord;
    typedef std::map<Coord, Item> tHcalValCont;
    // mapping of pair of subdetector name (e.g, "HE") and depth number (e.g. 3) to Histogram of data for that subdetector/depth pair
    typedef std::map<std::pair<std::string, int>, TH2F*> DepthMap;

    unsigned int GetRun() { return run_; }

    std::string GetTopoMode() { return TopoMode_; }

    std::map<std::string, int> GetSubDetDepths() { return subDetDepths_; }

    DepthMap GetDepths() {
      fillValConts();
      //std::cout << "Got depths with run number = " << std::to_string(GetRun()) << std::endl;
      return depths_;
    }


    // Fills a tHcalValCont for each subdetector, setting Topology Mode along the way
    void fillValConts() {
      if (!depths_.empty())
        return;
      int iphi, ieta, depth;
      HcalDetId detId;

      std::string subDetName;
      std::vector<Item> itemsVec;
      std::pair<std::string, int> depthKey;
      const char* histLabel;
      for (std::pair<std::string, std::vector<Item> > cont : (*payload_).getAllContainers()) {
        subDetName = std::get<0>(cont);
        itemsVec = std::get<1>(cont);

        auto valContainer = getContFromString(subDetName);

        for (Item item : itemsVec) {
          detId = HcalDetId(item.rawId());
          iphi = detId.iphi();
          ieta = detId.ieta();
          depth = detId.depth();
          Coord coord = std::make_tuple(depth, ieta, iphi);
          //Add hist if it's not there, AND if it is one we care about: HO,HB,HE,HF, AND it's not an empty entry (depth = 0 in HO)
          if (subDetName[0] == 'H' && depth != 0) {
            valContainer->insert(std::pair<std::tuple<int, int, int>, Item>(coord, item));

            depthKey = std::make_pair(subDetName, depth);

            auto depthVal = depths_.find(depthKey);
            if (depthVal == depths_.end()) {
              histLabel = ("run" + std::to_string(run_) + "_" + subDetName + "_d" + std::to_string(depth)).c_str();
              depths_.insert(std::make_pair(std::make_pair(subDetName, depth),
                                            new TH2F(histLabel, histLabel, 83, -42.5, 41.5, 71, 0.5, 71.5)));
            }
            depths_[depthKey]->Fill(ieta, iphi, getValue(&item));
          }
        }
      }

      //Still need to know which hists to take when done; decide now
      setTopoModeFromValConts();
    }

    const std::vector<std::pair<std::string, std::vector<Item> > > getAllItems() {
      return (*payload_).getAllContainers();
    }

    virtual float getValue(const Item* item) {
      throw cms::Exception("Value definition not found") << "getValue definition not found for " << payload_->myname();
    }

    //Gets Hcal Object at given coordinate
    //Currently unused but remains as a potentially useful function
    Item* getItemFromValCont(std::string subDetName, int depth, int ieta, int iphi, bool throwOnFail) {
      Item* cell = nullptr;
      Coord coord = std::make_tuple(depth, ieta, iphi);
      tHcalValCont* valContainer = getContFromString(subDetName);

      auto it = valContainer->find(coord);
      if (it != valContainer->end())
        cell = &it->second;

      if ((!cell)) {
        if (throwOnFail) {
          throw cms::Exception("Conditions not found")
              << "Unavailable Conditions of type " << payload_->myname() << " for cell " << subDetName << " (" << depth
              << ", " << ieta << ", " << iphi << ")";
        }
      }
      return cell;
    }


    //TODO: remove zero entries from doing divide and subtract

    // To generate Ratios of two IOVs
    void Divide(HcalDataContainer* dataCont2) {
      //TODO: Do something like looping over this and that depths setting every empty bin (which I think means content=0) to -999. Then replacing the divide call with another manual loop that sets all bins with content -999 to 0

      PlotMode_ = "Ratio";
      DepthMap::iterator depth1;
      std::pair<std::string, int> key;
      DepthMap myDepths = this->GetDepths();
      DepthMap theirDepths = dataCont2->GetDepths();
      for (depth1 = myDepths.begin(); depth1 != myDepths.end(); depth1++) {
        key = depth1->first;
        if (theirDepths.count(key) != 0) {
          myDepths.at(key)->Divide((const TH1*)theirDepths.at(key));
        } else {
          throw cms::Exception("Unaligned Conditions")
              << "trying to plot conditions for " << payload_->myname() << "; found value for " << std::get<0>(key)
              << " depth " << std::to_string(std::get<1>(key)) << " in run " << GetRun() << " but not in run "
              << dataCont2->GetRun();
        }
      }
    }

    // To generate Diffs of two IOVs
    void Subtract(HcalDataContainer* dataCont2) {
      PlotMode_ = "Diff";
      DepthMap::iterator depth1;
      std::pair<std::string, int> key;
      DepthMap myDepths = this->GetDepths();
      DepthMap theirDepths = dataCont2->GetDepths();
      for (depth1 = myDepths.begin(); depth1 != myDepths.end(); depth1++) {
        key = depth1->first;
        if (theirDepths.count(key) != 0) {
          myDepths.at(key)->Add(myDepths.at(key), theirDepths.at(key), 1, -1);
        } else {
          throw cms::Exception("Unaligned Conditions")
              << "trying to plot conditions for " << payload_->myname() << "; found value for " << std::get<0>(key)
              << " depth " << std::to_string(std::get<1>(key)) << " in run " << GetRun() << " but not in run "
              << dataCont2->GetRun();
        }
      }
    }

    // wisely determines what range to set histogram axis to
    std::pair<float, float> GetRange(TH1* hist) {
      if (PlotMode_ == "Ratio") {
        float amp;
        Double_t adjustMin = 1;
        Double_t tempMin;
        int nBinsX = hist->GetXaxis()->GetNbins();
        int nBinsY = hist->GetYaxis()->GetNbins();
        for (int i = 0; i < nBinsX; i++) {
          for (int j = 0; j < nBinsY; j++) {
            tempMin = hist->GetBinContent(i, j);
            if ((tempMin != 0) && (tempMin < adjustMin))
              adjustMin = tempMin;
          }
        }
        amp = std::max((1 - adjustMin), (hist->GetMaximum() - 1));
        //amp = std::max((1 - hist->GetMinimum()),(hist->GetMaximum() - 1) );
        return std::make_pair(1 - amp, 1 + amp);
      } else if (PlotMode_ == "Diff") {
        float amp;
        amp = std::max((0 - hist->GetMinimum()), hist->GetMaximum());
        return std::make_pair((-1 * amp), amp);
      } else {
        Double_t adjustMin = 10000;
        Double_t tempMin;
        int nBinsX = hist->GetXaxis()->GetNbins();
        int nBinsY = hist->GetYaxis()->GetNbins();
        for (int i = 0; i < nBinsX; i++) {
          for (int j = 0; j < nBinsY; j++) {
            tempMin = hist->GetBinContent(i, j);
            if ((tempMin != 0) && (tempMin < adjustMin))
              adjustMin = tempMin;
          }
        }
        return std::make_pair(((adjustMin == 10000) ? hist->GetMinimum() : adjustMin), hist->GetMaximum());
      }
    }

    // set style
    void initGraphics() {
      gStyle->SetOptStat(0);
      gStyle->SetPalette(1);
      gStyle->SetOptFit(0);
      gStyle->SetLabelFont(42);
      gStyle->SetLabelFont(42);
      gStyle->SetTitleFont(42);
      gStyle->SetTitleFont(42);
      gStyle->SetMarkerSize(0);
      gStyle->SetTitleOffset(1.3, "Y");
      gStyle->SetTitleOffset(1.0, "X");
      gStyle->SetNdivisions(510);
      gStyle->SetStatH(0.11);
      gStyle->SetStatW(0.33);
      gStyle->SetTitleW(0.4);
      gStyle->SetTitleX(0.13);
      gStyle->SetPadTickX(1);
      gStyle->SetPadTickY(1);
    }

    TH1D* GetProjection(TH2F* hist, std::string plotType, const char* newName, std::string subDetName, int depth) {
      //TODO: Also want average for standard projection of 2DHist (not ratio or diff)?
      //if (PlotMode_ != "Ratio") return (plotType=="EtaProfile") ? ((TH2F*)(hist->Clone("temp")))->ProjectionX(newName) : ((TH2F*)(hist->Clone("temp")))->ProjectionY(newName);

      //TH1D* projection = ((TH2F*)(depths_[std::make_pair(subDetName,depth)]->Clone("temp")))->ProjectionX(newName);

      int xBins = (plotType == "EtaProfile") ? 83 : 71;
      int etaMin = -42, etaMax = 42, phiMin = 1, phiMax = 72;
      int xMin = (plotType == "EtaProfile") ? etaMin : phiMin;
      int xMax = (plotType == "EtaProfile") ? etaMax : phiMax;
      int otherMin = (plotType == "EtaProfile") ? phiMin : etaMin;
      int otherMax = (plotType == "EtaProfile") ? phiMax : etaMax;
      TH1D* retHist = new TH1D(newName, newName, xBins, xMin, xMax);
      int numChannels;
      Double_t sumVal;
      Double_t channelVal;
      int ieta, iphi;
      int bin = 0;
      for (int i = xMin; i <= xMax; i++ && bin++) {
        numChannels = 0;
        sumVal = 0;
        for (int j = otherMin; j <= otherMax; j++) {
          ieta = (plotType == "EtaProfile") ? i : j;
          ieta += 42;
          iphi = (plotType == "EtaProfile") ? j : i;
          iphi += -1;
          channelVal = hist->GetBinContent(ieta, iphi);
          //std::cout << "(ieta, iphi)    :       (" << std::to_string(ieta) << ", " << std::to_string(iphi) << ")" << std::endl;
          if (channelVal != 0) {
            sumVal += channelVal;
            numChannels++;
          }
        }
        //if(sumVal !=0) projection->SetBinContent(i, sumVal/((Double_t)numChannels));//retHist->Fill(i,sumVal/((Double_t)numChannels));
        if (sumVal != 0)
          retHist->Fill(i, sumVal / ((Double_t)numChannels));
      }
      return retHist;
      //return projection;
    }

    // fills a canvas with given subdetector information, plotting all depths
    void FillCanv(TCanvas* canvas,
                  std::string subDetName,
                  int startDepth = 1,
                  int startCanv = 1,
                  std::string plotForm = "2DHist") {
      const char* newName;
      std::pair<float, float> range;
      int padNum;
      int maxDepth = (subDetName == "HO") ? 4 : subDetDepths_[subDetName];
      TH1D* projection;
      TLatex label;
      for (int i = startDepth; i <= maxDepth; i++) {
        //skip if data not obtained; TODO: Maybe add text on plot saying data not found?
        if (depths_.count(std::make_pair(subDetName, i)) == 0) {
          return;
        }

        padNum = i + startCanv - 1;
        if (subDetName == "HO")
          padNum = padNum - 3;
        canvas->cd(padNum);
        canvas->GetPad(padNum)->SetGridx(1);
        canvas->GetPad(padNum)->SetGridy(1);

        if (plotForm == "2DHist") {
          canvas->GetPad(padNum)->SetRightMargin(0.13);
          range = GetRange(depths_[std::make_pair(subDetName, i)]);
          depths_[std::make_pair(subDetName, i)]->Draw("colz");
          depths_[std::make_pair(subDetName, i)]->SetContour(100);
          depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetTitle("ieta");
          depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetTitle("iphi");
          depths_[std::make_pair(subDetName, i)]->GetXaxis()->CenterTitle();
          depths_[std::make_pair(subDetName, i)]->GetYaxis()->CenterTitle();
          depths_[std::make_pair(subDetName, i)]->GetZaxis()->SetRangeUser(std::get<0>(range), std::get<1>(range));
          depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetTitleSize(0.06);
          depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetTitleOffset(0.80);
          depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetTitleSize(0.06);
          depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetTitleOffset(0.80);
          depths_[std::make_pair(subDetName, i)]->GetYaxis()->SetLabelSize(0.055);
          depths_[std::make_pair(subDetName, i)]->GetXaxis()->SetLabelSize(0.055);
        } else {
          canvas->GetPad(padNum)->SetLeftMargin(0.152);
          canvas->GetPad(padNum)->SetRightMargin(0.02);
          //gStyle->SetTitleOffset(1.6,"Y");
          newName = ("run_" + std::to_string(run_) + "_" + subDetName + "_d" + std::to_string(i) + "_" +
                     (plotForm == "EtaProfile" ? "ieta" : "iphi"))
                        .c_str();
          //projection = ((TH2F*)(depths_[std::make_pair(subDetName,i)]->Clone("temp")))->ProjectionX(newName);
          projection = GetProjection(depths_[std::make_pair(subDetName, i)], plotForm, newName, subDetName, i);
          range = GetRange(projection);
          projection->Draw("hist");
          projection->GetXaxis()->SetTitle((plotForm == "EtaProfile" ? "ieta" : "iphi"));
          projection->GetXaxis()->CenterTitle();
          projection->GetYaxis()->SetTitle(
              (payload_->myname() + " " + (PlotMode_ == "Map" ? "" : PlotMode_) + " " + GetUnit(payload_->myname()))
                  .c_str());
          label.SetNDC();
          label.SetTextAlign(26);
          label.SetTextSize(0.05);
          label.DrawLatex(
              0.3, 0.95, ("run_" + std::to_string(run_) + "_" + subDetName + "_d" + std::to_string(i)).c_str());
          projection->GetYaxis()->CenterTitle();
          projection->GetXaxis()->SetTitleSize(0.06);
          projection->GetYaxis()->SetTitleSize(0.06);
          projection->GetXaxis()->SetTitleOffset(0.80);
          projection->GetXaxis()->SetLabelSize(0.055);
          projection->GetYaxis()->SetTitleOffset(1.34);
          projection->GetYaxis()->SetLabelSize(0.055);
        }
      }
    }


    // profile = "EtaProfile" || "PhiProfile"
    TCanvas* getCanvasAll(std::string profile = "2DHist") {
      fillValConts();
      initGraphics();
      TCanvas* HAll = new TCanvas("HAll", "HAll", 1680, (GetTopoMode() == "2015/2016") ? 1680 : 2500);
      HAll->Divide(3, (GetTopoMode() == "2015/2016") ? 3 : 6, 0.02, 0.01);
      FillCanv(HAll, "HB", 1, 1, profile);
      FillCanv(HAll, "HO", 4, 3, profile);
      FillCanv(HAll, "HF", 1, 4, profile);
      FillCanv(HAll, "HE", 1, (GetTopoMode() == "2015/2016") ? 7 : 10, profile);
      return HAll;
    }

    TCanvas* getCanvasHF() {
      fillValConts();
      initGraphics();
      TCanvas* HF = new TCanvas("HF", "HF", 1600, 1000);
      HF->Divide(3, 2, 0.02, 0.01);
      FillCanv(HF, "HF");
      return HF;
    }
    TCanvas* getCanvasHE() {
      fillValConts();
      initGraphics();
      TCanvas* HE = new TCanvas("HE", "HE", 1680, 1680);
      HE->Divide(3, 3, 0.02, 0.01);
      FillCanv(HE, "HE");
      return HE;
    }
    TCanvas* getCanvasHBHO() {
      fillValConts();
      initGraphics();
      TCanvas* HBHO = new TCanvas("HBHO", "HBHO", 1680, 1680);
      HBHO->Divide(3, 3, 0.02, 0.01);
      FillCanv(HBHO, "HB");
      FillCanv(HBHO, "HO", 4, 3);
      FillCanv(HBHO, "HB", 1, 4, "EtaProfile");
      FillCanv(HBHO, "HO", 4, 6, "EtaProfile");
      FillCanv(HBHO, "HB", 1, 7, "PhiProfile");
      FillCanv(HBHO, "HO", 4, 9, "PhiProfile");
      return HBHO;
    }

    std::string GetUnit(std::string type) {
      std::string unit = units_[type];
      if (unit.empty())
        return "";
      else
        return "(" + unit + ")";
    }

  private:
    DepthMap depths_;
    std::shared_ptr<Items> payload_;
    unsigned int run_;
    std::string TopoMode_;
    // "Map", "Ratio", or "Diff"
    std::string PlotMode_;

    tHcalValCont HBvalContainer;
    tHcalValCont HEvalContainer;
    tHcalValCont HOvalContainer;
    tHcalValCont HFvalContainer;
    tHcalValCont HTvalContainer;
    tHcalValCont ZDCvalContainer;
    tHcalValCont CALIBvalContainer;
    tHcalValCont CASTORvalContainer;
    std::map<std::string, int> subDetDepths_;
    std::map<std::string, std::string> units_ = {{"HcalPedestals", "ADC"},
                                                 {"HcalGains", ""},             //dimensionless TODO: verify
                                                 {"HcalL1TriggerObjects", ""},  //dimensionless TODO: Verify
                                                 {"HcalPedestalWidths", "ADC"},
                                                 {"HcalRespCorrs", ""},  //dimensionless TODO: verify
                                                 {"Dark Current", ""},
                                                 {"fcByPE", ""},
                                                 {"crossTalk", ""},
                                                 {"parLin", ""}};

    tHcalValCont* getContFromString(std::string subDetString) {
      if (subDetString == "HB")
        return &HBvalContainer;
      else if (subDetString == "HE")
        return &HEvalContainer;
      else if (subDetString == "HF")
        return &HFvalContainer;
      else if (subDetString == "HO")
        return &HOvalContainer;
      else if (subDetString == "HT")
        return &HTvalContainer;
      else if (subDetString == "CALIB")
        return &CALIBvalContainer;
      else if (subDetString == "CASTOR")
        return &CASTORvalContainer;
      //else return &ZDCvalContainer;
      else if (subDetString == "ZDC_EM" || subDetString == "ZDC" || subDetString == "ZDC_HAD" ||
               subDetString == "ZDC_LUM")
        return &ZDCvalContainer;
      else
        throw cms::Exception("subDetString " + subDetString + " not found in Item");
    }

    void setTopoModeFromValConts(bool throwOnFail = false) {
      // Check HEP17 alternate channel for 2017, just by checking if the 7th depth is there, or if HF has 4 depths
      if (depths_.count(std::make_pair("HF", 4)) != 0 || depths_.count(std::make_pair("HE", 7)) != 0)
        TopoMode_ = "2017";
      // Check endcap depth unique to 2018
      else if (HEvalContainer.count(std::make_tuple(7, -26, 63)) != 0)
        TopoMode_ = "2018";
      // if not 2017 or 2018, 2015 and 2016 are the same
      else
        TopoMode_ = "2015/2016";
      //NOTE: HO's one depth is labeled depth 4
      if (TopoMode_ == "2018" || TopoMode_ == "2017") {
        subDetDepths_.insert(std::pair<std::string, int>("HB", 2));
        subDetDepths_.insert(std::pair<std::string, int>("HE", 7));
        subDetDepths_.insert(std::pair<std::string, int>("HF", 4));
        subDetDepths_.insert(std::pair<std::string, int>("HO", 1));
      } else if (TopoMode_ == "2015/2016") {
        subDetDepths_.insert(std::pair<std::string, int>("HB", 2));
        subDetDepths_.insert(std::pair<std::string, int>("HE", 3));
        subDetDepths_.insert(std::pair<std::string, int>("HF", 2));
        subDetDepths_.insert(std::pair<std::string, int>("HO", 1));
      }
    }
  };

  inline void drawTable(int nbRows, int nbColumns) {
    TLine* l = new TLine;
    l->SetLineWidth(1);
    for (int i = 1; i < nbRows; i++) {
      double y = (double)i;
      l = new TLine(0., y, nbColumns, y);
      l->Draw();
    }

    for (int i = 1; i < nbColumns; i++) {
      double x = (double)i;
      double y = (double)nbRows;
      l = new TLine(x, 0., x, y);
      l->Draw();
    }
  }

  inline std::string SciNotatStr(float num) {
    // Create an output string stream
    std::ostringstream streamObj2;

    if (num == -1)
      return "NOT FOUND";

    //Add double to stream
    streamObj2 << num;
    // Get string from output string stream
    std::string strObj2 = streamObj2.str();

    return strObj2;
  }

  inline std::string IntToBinary(unsigned int number) {
    std::stringstream ss;
    unsigned int mask = 1 << 31;
    for (unsigned short int i = 0; i < 32; ++i) {
      //if (!(i % 4))
      //    ss << "_";
      if (mask & number)
        ss << "1";
      else
        ss << "0";
      mask = mask >> 1;
    }
    return ss.str();
  }

  inline const bool isBitSet(unsigned int bitnumber, unsigned int status) {
    unsigned int statadd = 0x1 << (bitnumber);
    return (status & statadd) ? (true) : (false);
  }

  inline std::string getBitsSummary(uint32_t bits, std::string statusBitArray[], short unsigned int bitMap[]) {
    std::stringstream ss;
    for (unsigned int i = 0; i < 9; ++i) {
      if (isBitSet(bitMap[i], bits)) {
        ss << "[" << bitMap[i] << "]" << statusBitArray[bitMap[i]] << "; ";
      }
    }
    ss << std::endl;
    return ss.str();
  }

  //functions for making plot:
  inline void setBinLabels(std::vector<TH2F>& depth) {
    // Set labels for all depth histograms
    for (unsigned int i = 0; i < depth.size(); ++i) {
      depth[i].SetXTitle("i#eta");
      depth[i].SetYTitle("i#phi");
    }

    std::stringstream label;

    // set label on every other bin
    for (int i = -41; i <= -29; i = i + 2) {
      label << i;
      depth[0].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
      depth[1].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
      label.str("");
    }
    depth[0].GetXaxis()->SetBinLabel(14, "-29HE");
    depth[1].GetXaxis()->SetBinLabel(14, "-29HE");

    // offset by one for HE
    for (int i = -27; i <= 27; i = i + 2) {
      label << i;
      depth[0].GetXaxis()->SetBinLabel(i + 43, label.str().c_str());
      label.str("");
    }
    depth[0].GetXaxis()->SetBinLabel(72, "29HE");
    for (int i = 29; i <= 41; i = i + 2) {
      label << i;
      depth[0].GetXaxis()->SetBinLabel(i + 44, label.str().c_str());
      label.str("");
    }
    for (int i = 16; i <= 28; i = i + 2) {
      label << i - 43;
      depth[1].GetXaxis()->SetBinLabel(i, label.str().c_str());
      label.str("");
    }
    depth[1].GetXaxis()->SetBinLabel(29, "NULL");
    for (int i = 15; i <= 27; i = i + 2) {
      label << i;
      depth[1].GetXaxis()->SetBinLabel(i + 15, label.str().c_str());
      label.str("");
    }

    depth[1].GetXaxis()->SetBinLabel(44, "29HE");
    for (int i = 29; i <= 41; i = i + 2) {
      label << i;
      depth[1].GetXaxis()->SetBinLabel(i + 16, label.str().c_str());
      label.str("");
    }

    // HE depth 3 labels;
    depth[2].GetXaxis()->SetBinLabel(1, "-28");
    depth[2].GetXaxis()->SetBinLabel(2, "-27");
    depth[2].GetXaxis()->SetBinLabel(3, "Null");
    depth[2].GetXaxis()->SetBinLabel(4, "-16");
    depth[2].GetXaxis()->SetBinLabel(5, "Null");
    depth[2].GetXaxis()->SetBinLabel(6, "16");
    depth[2].GetXaxis()->SetBinLabel(7, "Null");
    depth[2].GetXaxis()->SetBinLabel(8, "27");
    depth[2].GetXaxis()->SetBinLabel(9, "28");
  }
  // Now define functions that can be used in conjunction with EtaPhi histograms

  // This arrays the eta binning for depth 2 histograms (with a gap between -15 -> +15)
  const int binmapd2[] = {-42, -41, -40, -39, -38, -37, -36, -35, -34, -33, -32, -31, -30, -29,   -28,
                          -27, -26, -25, -24, -23, -22, -21, -20, -19, -18, -17, -16, -15, -9999, 15,
                          16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,  28,  29,    30,
                          31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,  42};

  // This stores eta binning in depth 3 (where HE is only present at a few ieta values)

  const int binmapd3[] = {-28, -27, -9999, -16, -9999, 16, -9999, 27, 28};

  inline int CalcEtaBin(int subdet, int ieta, int depth) {
    // This takes the eta value from a subdetector and return an eta counter value as used by eta-phi array
    // (ieta=-41 corresponds to bin 0, +41 to bin 85 -- there are two offsets to deal with the overlap at |ieta|=29).
    // For HO, ieta = -15 corresponds to bin 0, and ieta=15 is bin 30
    // For HE depth 3, things are more complicated, but feeding the ieta value will give back the corresponding counter eta value

    // The CalcEtaBin value is the value as used within our array counters, and thus starts at 0.
    // If you are using it with getBinContent or setBinContent, you will need to add +1 to the result of this function

    int etabin = -9999;  // default invalid value

    if (depth == 1) {
      // Depth 1 is fairly straightforward -- just shift HF-, HF+ by -/+1
      etabin = ieta + 42;
      if (subdet == HcalForward) {
        ieta < 0 ? etabin-- : etabin++;
      }
    }

    else if (depth == 2) {
      // Depth 2 is more complicated, given that there are no cells in the range |ieta|<15
      if (ieta < -14) {
        etabin = ieta + 42;
        if (subdet == HcalForward)
          etabin--;
      } else if (ieta > 14) {
        etabin = ieta + 14;
        if (subdet == HcalForward)
          etabin++;
      }

    }
    // HO is also straightforward; a simple offset to the ieta value is applied
    else if (subdet == HcalOuter && abs(ieta) < 16)
      etabin = ieta + 15;
    else if (subdet == HcalEndcap) {
      // HE depth 3 has spotty coverage; hard-code the bin response
      if (depth == 3) {
        if (ieta == -28)
          etabin = 0;
        else if (ieta == -27)
          etabin = 1;
        else if (ieta == -16)
          etabin = 3;
        else if (ieta == 16)
          etabin = 5;
        else if (ieta == 27)
          etabin = 7;
        else if (ieta == 28)
          etabin = 8;
      }
    }
    return etabin;
  }

  inline int CalcIeta(int subdet, int eta, int depth) {
    // This function returns the 'true' ieta value given subdet, eta, and depth
    // Here 'eta' is the index from our arrays (it starts at 0);
    // remember that histogram bins start with bin 1, so there's an offset of 1
    // to consider if using getBinContent(eta,phi)

    // eta runs from 0...X  (X depends on depth)
    int ieta = -9999;  // default value is nonsensical
    if (subdet == HcalBarrel) {
      if (depth == 1) {
        ieta = eta - 42;
        if (ieta == 0)
          return -9999;
        return ieta;
      } else if (depth == 2) {
        ieta = binmapd2[eta];
        if (ieta == 0)
          return -9999;
        if (ieta == 17 || ieta == -17)
          return -9999;  // no depth 2 cells at |ieta| = 17
        return ieta;
      } else
        return -9999;  // non-physical value
    } else if (subdet == HcalForward) {
      if (depth == 1) {
        ieta = eta - 42;
        if (eta < 13)
          ieta++;
        else if (eta > 71)
          ieta--;
        else
          return -9999;  // if outside forward range, return dummy
        return ieta;
      } else if (depth == 2) {
        ieta = binmapd2[eta];  // special map for depth 2
        if (ieta <= -30)
          ieta++;
        else if (ieta >= 30)
          ieta--;
        else
          return -9999;
        return ieta;
      } else
        return -9999;
    }

    else if (subdet == HcalEndcap) {
      if (depth == 1)
        ieta = eta - 42;
      else if (depth == 2) {
        ieta = binmapd2[eta];
        if (abs(ieta) > 29 || abs(ieta) < 18)
          return -9999;  // outside HE
        if (ieta == 0)
          return -9999;
        return ieta;
      } else if (depth == 3) {
        if (eta < 0 || eta > 8)
          return -9999;
        else
          ieta = binmapd3[eta];  // special map for depth 3
        if (ieta == 0)
          return -9999;
        return ieta;
      } else
        return -9999;
    }  // HcalEndcap
    else if (subdet == HcalOuter) {
      if (depth != 4)
        return -9999;
      else {
        ieta = eta - 15;  // bin 0 is ieta=-15, all bins increment normally from there
        if (abs(ieta) > 15)
          return -9999;
        if (ieta == 0)
          return -9999;
        return ieta;
      }
    }  // HcalOuter
    if (ieta == 0)
      return -9999;
    return ieta;
  }

  inline int CalcIeta(int eta, int depth) {
    // This version of CalcIeta does the same as the function above,
    // but does not require that 'subdet' be specified.

    // returns ieta value give an eta counter.
    // eta runs from 0...X  (X depends on depth)
    int ieta = -9999;
    if (eta < 0)
      return ieta;
    if (depth == 1) {
      ieta =
          eta -
          42;  // default shift: bin 0 corresponds to a histogram ieta of -42 (which is offset by 1 from true HF value of -41)
      if (eta < 13)
        ieta++;
      else if (eta > 71)
        ieta--;
      if (ieta == 0)
        ieta = -9999;
      return ieta;
    } else if (depth == 2) {
      if (eta > 57)
        return -9999;
      else {
        ieta = binmapd2[eta];
        if (ieta == -9999)
          return ieta;
        if (ieta == 0)
          return -9999;
        if (ieta == 17 || ieta == -17)
          return -9999;  // no depth 2 cells at |ieta| = 17
        else if (ieta <= -30)
          ieta++;
        else if (ieta >= 30)
          ieta--;
        return ieta;
      }
    } else if (depth == 3) {
      if (eta > 8)
        return -9999;
      else
        ieta = binmapd3[eta];
      if (ieta == 0)
        return -9999;
      return ieta;
    } else if (depth == 4) {
      ieta = eta - 15;  // bin 0 is ieta=-15, all bins increment normally from there
      if (abs(ieta) > 15)
        return -9999;
      if (ieta == 0)
        return -9999;
      return ieta;
    }
    return ieta;  // avoids compilation warning
  }

  // Functions to check whether a given (eta,depth) value is valid for a given subdetector

  inline std::vector<std::string> HcalEtaPhiHistNames() {
    std::vector<std::string> name;
    name.push_back("HB HE HF Depth 1 ");
    name.push_back("HB HE HF Depth 2 ");
    name.push_back("HE Depth 3 ");
    name.push_back("HO Depth 4 ");
    return name;
  }

  inline bool isHB(int etabin, int depth) {
    if (depth > 2)
      return false;
    else if (depth < 1)
      return false;
    else {
      int ieta = CalcIeta(etabin, depth);
      if (ieta == -9999)
        return false;
      if (depth == 1) {
        if (abs(ieta) <= 16)
          return true;
        else
          return false;
      } else if (depth == 2) {
        if (abs(ieta) == 15 || abs(ieta) == 16)
          return true;
        else
          return false;
      }
    }
    return false;
  }

  inline bool isHE(int etabin, int depth) {
    if (depth > 3)
      return false;
    else if (depth < 1)
      return false;
    else {
      int ieta = CalcIeta(etabin, depth);
      if (ieta == -9999)
        return false;
      if (depth == 1) {
        if (abs(ieta) >= 17 && abs(ieta) <= 28)
          return true;
        if (ieta == -29 && etabin == 13)
          return true;  // HE -29
        if (ieta == 29 && etabin == 71)
          return true;  // HE +29
      } else if (depth == 2) {
        if (abs(ieta) >= 17 && abs(ieta) <= 28)
          return true;
        if (ieta == -29 && etabin == 13)
          return true;  // HE -29
        if (ieta == 29 && etabin == 43)
          return true;  // HE +29
      } else if (depth == 3)
        return true;
    }
    return false;
  }

  inline bool isHF(int etabin, int depth) {
    if (depth > 2)
      return false;
    else if (depth < 1)
      return false;
    else {
      int ieta = CalcIeta(etabin, depth);
      if (ieta == -9999)
        return false;
      if (depth == 1) {
        if (ieta == -29 && etabin == 13)
          return false;  // HE -29
        else if (ieta == 29 && etabin == 71)
          return false;  // HE +29
        else if (abs(ieta) >= 29)
          return true;
      } else if (depth == 2) {
        if (ieta == -29 && etabin == 13)
          return false;  // HE -29
        else if (ieta == 29 && etabin == 43)
          return false;  // HE +29
        else if (abs(ieta) >= 29)
          return true;
      }
    }
    return false;
  }

  inline bool isHO(int etabin, int depth) {
    if (depth != 4)
      return false;
    int ieta = CalcIeta(etabin, depth);
    if (ieta != -9999)
      return true;
    return false;
  }

  // Checks whether HO region contains SiPM

  inline bool isSiPM(int ieta, int iphi, int depth) {
    if (depth != 4)
      return false;
    // HOP1
    if (ieta >= 5 && ieta <= 10 && iphi >= 47 && iphi <= 58)
      return true;
    // HOP2
    if (ieta >= 11 && ieta <= 15 && iphi >= 59 && iphi <= 70)
      return true;
    return false;
  }  // bool isSiPM

  // Checks whether (subdet, ieta, iphi, depth) value is a valid Hcal cell

  inline bool validDetId(HcalSubdetector sd, int ies, int ip, int dp) {
    // inputs are (subdetector, ieta, iphi, depth)
    // stolen from latest version of DataFormats/HcalDetId/src/HcalDetId.cc (not yet available in CMSSW_2_1_9)

    const int ie(abs(ies));

    return (
        (ip >= 1) && (ip <= 72) && (dp >= 1) && (ie >= 1) &&
        (((sd == HcalBarrel) && (((ie <= 14) && (dp == 1)) || (((ie == 15) || (ie == 16)) && (dp <= 2)))) ||
         ((sd == HcalEndcap) &&
          (((ie == 16) && (dp == 3)) || ((ie == 17) && (dp == 1)) || ((ie >= 18) && (ie <= 20) && (dp <= 2)) ||
           ((ie >= 21) && (ie <= 26) && (dp <= 2) && (ip % 2 == 1)) ||
           ((ie >= 27) && (ie <= 28) && (dp <= 3) && (ip % 2 == 1)) || ((ie == 29) && (dp <= 2) && (ip % 2 == 1)))) ||
         ((sd == HcalOuter) && (ie <= 15) && (dp == 4)) ||
         ((sd == HcalForward) && (dp <= 2) &&
          (((ie >= 29) && (ie <= 39) && (ip % 2 == 1)) || ((ie >= 40) && (ie <= 41) && (ip % 4 == 3))))));

  }  // bool validDetId(HcalSubdetector sd, int ies, int ip, int dp)

  // Sets eta, phi labels for 'summary' eta-phi plots (identical to Depth 1 Eta-Phi labelling)

  inline void SetEtaPhiLabels(TH2F& h) {
    std::stringstream label;
    for (int i = -41; i <= -29; i = i + 2) {
      label << i;
      h.GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
      label.str("");
    }
    h.GetXaxis()->SetBinLabel(14, "-29HE");

    // offset by one for HE
    for (int i = -27; i <= 27; i = i + 2) {
      label << i;
      h.GetXaxis()->SetBinLabel(i + 43, label.str().c_str());
      label.str("");
    }
    h.GetXaxis()->SetBinLabel(72, "29HE");
    for (int i = 29; i <= 41; i = i + 2) {
      label << i;
      h.GetXaxis()->SetBinLabel(i + 44, label.str().c_str());
      label.str("");
    }
    return;
  }

  // Fill Unphysical bins in histograms
  inline void FillUnphysicalHEHFBins(std::vector<TH2F>& hh) {
    int ieta = 0;
    int iphi = 0;
    // First 2 depths have 5-10-20 degree corrections
    for (unsigned int d = 0; d < 3; ++d) {
      //BUG CAN BE HERE:
      //if (hh[d] != 0) continue;

      for (int eta = 0; eta < hh[d].GetNbinsX(); ++eta) {
        ieta = CalcIeta(eta, d + 1);
        if (ieta == -9999 || abs(ieta) < 21)
          continue;
        for (int phi = 0; phi < hh[d].GetNbinsY(); ++phi) {
          iphi = phi + 1;
          if (iphi % 2 == 1 && abs(ieta) < 40 && iphi < 73) {
            hh[d].SetBinContent(eta + 1, iphi + 1, hh[d].GetBinContent(eta + 1, iphi));
          }
          // last two eta strips span 20 degrees in phi
          // Fill the phi cell above iphi, and the 2 below it
          else if (abs(ieta) > 39 && iphi % 4 == 3 && iphi < 73) {
            //ieta=40, iphi=3 covers iphi 3,4,5,6
            hh[d].SetBinContent(eta + 1, (iphi) % 72 + 1, hh[d].GetBinContent(eta + 1, iphi));
            hh[d].SetBinContent(eta + 1, (iphi + 1) % 72 + 1, hh[d].GetBinContent(eta + 1, iphi));
            hh[d].SetBinContent(eta + 1, (iphi + 2) % 72 + 1, hh[d].GetBinContent(eta + 1, iphi));
          }
        }  // for (int phi...)
      }    // for (int eta...)
    }      // for (int d=0;...)
    // no corrections needed for HO (depth 4)
    return;
  }  // FillUnphysicalHEHFBins(MonitorElement* hh)

  //Fill unphysical bins for single ME
  inline void FillUnphysicalHEHFBins(TH2F& hh) {
    // Fills unphysical HE/HF bins for Summary Histogram
    // Summary Histogram is binned with the same binning as the Depth 1 EtaPhiHists

    //CAN BE BUG HERE:
    //if (hh==0) return;

    int ieta = 0;
    int iphi = 0;
    int etabins = hh.GetNbinsX();
    int phibins = hh.GetNbinsY();
    float binval = 0;
    for (int eta = 0; eta < etabins; ++eta)  // loop over eta bins
    {
      ieta = CalcIeta(eta, 1);
      if (ieta == -9999 || abs(ieta) < 21)
        continue;  // ignore etas that don't exist, or that have 5 degree phi binning

      for (int phi = 0; phi < phibins; ++phi) {
        iphi = phi + 1;
        if (iphi % 2 == 1 && abs(ieta) < 40 && iphi < 73)  // 10 degree phi binning condition
        {
          binval = hh.GetBinContent(eta + 1, iphi);
          hh.SetBinContent(eta + 1, iphi + 1, binval);
        }                                                       // if (iphi%2==1...)
        else if (abs(ieta) > 39 && iphi % 4 == 3 && iphi < 73)  // 20 degree phi binning condition
        {
          // Set last two eta strips where each cell spans 20 degrees in phi
          // Set next phi cell above iphi, and 2 cells below the actual cell
          hh.SetBinContent(eta + 1, (iphi) % 72 + 1, hh.GetBinContent(eta + 1, iphi));
          hh.SetBinContent(eta + 1, (iphi + 1) % 72 + 1, hh.GetBinContent(eta + 1, iphi));
          hh.SetBinContent(eta + 1, (iphi + 2) % 72 + 1, hh.GetBinContent(eta + 1, iphi));
        }  // else if (abs(ieta)>39 ...)
      }    // for (int phi=0;phi<72;++phi)

    }  // for (int eta=0; eta< (etaBins_-2);++eta)

    return;
  }  // FillUnphysicalHEHFBins(std::vector<MonitorElement*> &hh)

  // special fill call based on detid -- eventually will need special treatment
  inline void Fill(HcalDetId& id, double val /*=1*/, std::vector<TH2F>& depth) {
    // If in HF, need to shift by 1 bin (-1 bin lower in -HF, +1 bin higher in +HF)
    if (id.subdet() == HcalForward)
      depth[id.depth() - 1].Fill(id.ieta() < 0 ? id.ieta() - 1 : id.ieta() + 1, id.iphi(), val);
    else
      depth[id.depth() - 1].Fill(id.ieta(), id.iphi(), val);
  }

  inline void Reset(std::vector<TH2F>& depth) {
    for (unsigned int d = 0; d < depth.size(); d++)
      //BUG CAN BE HERE:
      //if(depth[d])
      depth[d].Reset();
  }  // void Reset(void)

  inline void setup(std::vector<TH2F>& depth, std::string name, std::string units = "") {
    std::string unittitle, unitname;
    if (units.empty()) {
      unitname = units;
      unittitle = "No Units";
    } else {
      unitname = " " + units;
      unittitle = units;
    }

    // Push back depth plots
    //1. create first plot
    depth.push_back(TH2F(("HB HE HF Depth 1 " + name + unitname).c_str(),
                         (name + " Depth 1 -- HB HE HF (" + unittitle + ")").c_str(),
                         85,
                         -42.5,
                         42.5,
                         72,
                         0.5,
                         72.5));

    //2.1 prepare second plot
    float ybins[73];
    for (int i = 0; i <= 72; i++)
      ybins[i] = (float)(i + 0.5);
    float xbinsd2[] = {-42.5, -41.5, -40.5, -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5,
                       -30.5, -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5, -19.5,
                       -18.5, -17.5, -16.5, -15.5, -14.5, 14.5,  15.5,  16.5,  17.5,  18.5,  19.5,  20.5,
                       21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,  30.5,  31.5,  32.5,
                       33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,  40.5,  41.5,  42.5};

    //2.2 create second plot
    depth.push_back(TH2F(("HB HE HF Depth 2 " + name + unitname).c_str(),
                         (name + " Depth 2 -- HB HE HF (" + unittitle + ")").c_str(),
                         57,
                         xbinsd2,
                         72,
                         ybins));

    //3.1 Set up variable-sized bins for HE depth 3 (MonitorElement also requires phi bins to be entered in array format)
    float xbins[] = {-28.5, -27.5, -26.5, -16.5, -15.5, 15.5, 16.5, 26.5, 27.5, 28.5};
    //3.2
    depth.push_back(TH2F(("HE Depth 3 " + name + unitname).c_str(),
                         (name + " Depth 3 -- HE (" + unittitle + ")").c_str(),
                         // Use variable-sized eta bins
                         9,
                         xbins,
                         72,
                         ybins));

    //4.1 HO bins are fixed width, but cover a smaller eta range (-15 -> 15)
    depth.push_back(TH2F(("HO Depth 4 " + name + unitname).c_str(),
                         (name + " Depth 4 -- HO (" + unittitle + ")").c_str(),
                         31,
                         -15.5,
                         15.5,
                         72,
                         0.5,
                         72.5));

    for (unsigned int i = 0; i < depth.size(); ++i)
      depth[i].Draw("colz");

    setBinLabels(depth);  // set axis titles, special bins
  }

  inline void fillOneGain(std::vector<TH2F>& graphData,
                          HcalGains::tAllContWithNames& allContainers,
                          std::string name,
                          int id,
                          std::string units = "") {
    setup(graphData, name);

    std::stringstream x;
    // Change the titles of each individual histogram
    for (unsigned int d = 0; d < graphData.size(); ++d) {
      graphData[d].Reset();
      x << "Gain " << id << " for HCAL depth " << d + 1;

      //BUG CAN BE HERE:
      //if (ChannelStatus->depth[d])
      graphData[d].SetTitle(
          x.str()
              .c_str());  // replace "setTitle" with "SetTitle", since you are using TH2F objects instead of MonitorElements
      x.str("");
    }

    HcalDetId hcal_id;
    int ieta, depth, iphi;

    //main loop
    // get all containers with names
    //HcalGains::tAllContWithNames allContainers = object().getAllContainers();

    //ITERATORS AND VALUES:
    HcalGains::tAllContWithNames::const_iterator iter;
    std::vector<HcalGain>::const_iterator contIter;
    float gain = 0.0;

    //Run trough given id gain:
    int i = id;

    //run trough all pair containers
    for (iter = allContainers.begin(); iter != allContainers.end(); ++iter) {
      //Run trough all values:
      for (contIter = (*iter).second.begin(); contIter != (*iter).second.end(); ++contIter) {
        hcal_id = HcalDetId((uint32_t)(*contIter).rawId());

        depth = hcal_id.depth();
        if (depth < 1 || depth > 4)
          continue;

        ieta = hcal_id.ieta();
        iphi = hcal_id.iphi();

        if (hcal_id.subdet() == HcalForward)
          ieta > 0 ? ++ieta : --ieta;

        //GET VALUE:
        gain = (*contIter).getValue(i);
        //logstatus = log2(1.*channelBits)+1;

        //FILLING GOES HERE:
        graphData[depth - 1].Fill(ieta, iphi, gain);

        //FOR DEBUGGING:
        //std::cout << "ieta: " << ieta << "; iphi: " << iphi << "; logstatus: " << logstatus << "; channelBits: " << channelBits<< std::endl;
      }
    }
  }

  //FillUnphysicalHEHFBins(graphData);
  //return ("kasdasd");

  class ADataRepr  //Sample base class for c++ inheritance tutorial
  {
  public:
    ADataRepr(unsigned int d) : m_total(d){};
    virtual ~ADataRepr(){};
    unsigned int nr, id;
    std::stringstream filename, rootname, plotname;

    void fillOneGain(std::vector<TH2F>& graphData, std::string units = "") {
      std::stringstream ss("");

      if (m_total == 1)
        ss << rootname.str() << " for HCAL depth ";
      else
        ss << rootname.str() << nr << " for HCAL depth ";

      setup(graphData, ss.str());

      // Change the titles of each individual histogram
      for (unsigned int d = 0; d < graphData.size(); ++d) {
        graphData[d].Reset();

        ss.str("");
        if (m_total == 1)
          ss << plotname.str() << " for HCAL depth " << d + 1;
        else
          ss << plotname.str() << nr << " for HCAL depth " << d + 1;

        //BUG CAN BE HERE:
        //if (ChannelStatus->depth[d])
        graphData[d].SetTitle(
            ss.str()
                .c_str());  // replace "setTitle" with "SetTitle", since you are using TH2F objects instead of MonitorElements
        ss.str("");
      }
      //overload this function:
      doFillIn(graphData);

      FillUnphysicalHEHFBins(graphData);

      ss.str("");
      if (m_total == 1)
        ss << filename.str() << ".png";
      else
        ss << filename.str() << nr << ".png";
      draw(graphData, ss.str());
      //FOR DEBUGGING:
      //std::cout << "ieta: " << ieta << "; iphi: " << iphi << "; logstatus: " << logstatus << "; channelBits: " << channelBits<< std::endl;
    }

  protected:
    unsigned int m_total;
    HcalDetId hcal_id;
    int ieta, depth, iphi;

    virtual void doFillIn(std::vector<TH2F>& graphData) = 0;

  private:
    void draw(std::vector<TH2F>& graphData, std::string filename) {
      //Drawing...
      // use David's palette
      gStyle->SetPalette(1);
      const Int_t NCont = 999;
      gStyle->SetNumberContours(NCont);
      TCanvas canvas("CC map", "CC map", 840, 369 * 4);

      TPad pad1("pad1", "pad1", 0.0, 0.75, 1.0, 1.0);
      pad1.Draw();
      TPad pad2("pad2", "pad2", 0.0, 0.5, 1.0, 0.75);
      pad2.Draw();
      TPad pad3("pad3", "pad3", 0.0, 0.25, 1.0, 0.5);
      pad3.Draw();
      TPad pad4("pad4", "pad4", 0.0, 0.0, 1.0, 0.25);
      pad4.Draw();

      pad1.cd();
      graphData[0].SetStats(false);
      graphData[0].Draw("colz");

      pad2.cd();
      graphData[1].SetStats(false);
      graphData[1].Draw("colz");

      pad3.cd();
      graphData[2].SetStats(false);
      graphData[2].Draw("colz");

      pad4.cd();
      graphData[3].SetStats(false);
      graphData[3].Draw("colz");

      canvas.SaveAs(filename.c_str());
    }

    void setup(std::vector<TH2F>& depth, std::string name, std::string units = "") {
      std::string unittitle, unitname;
      if (units.empty()) {
        unitname = units;
        unittitle = "No Units";
      } else {
        unitname = " " + units;
        unittitle = units;
      }

      // Push back depth plots
      //1. create first plot
      depth.push_back(TH2F(("HB HE HF Depth 1 " + name + unitname).c_str(),
                           (name + " Depth 1 -- HB HE HF (" + unittitle + ")").c_str(),
                           85,
                           -42.5,
                           42.5,
                           72,
                           0.5,
                           72.5));

      //2.1 prepare second plot
      float ybins[73];
      for (int i = 0; i <= 72; i++)
        ybins[i] = (float)(i + 0.5);
      float xbinsd2[] = {-42.5, -41.5, -40.5, -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5,
                         -30.5, -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5, -19.5,
                         -18.5, -17.5, -16.5, -15.5, -14.5, 14.5,  15.5,  16.5,  17.5,  18.5,  19.5,  20.5,
                         21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,  30.5,  31.5,  32.5,
                         33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,  40.5,  41.5,  42.5};

      //2.2 create second plot
      depth.push_back(TH2F(("HB HE HF Depth 2 " + name + unitname).c_str(),
                           (name + " Depth 2 -- HB HE HF (" + unittitle + ")").c_str(),
                           57,
                           xbinsd2,
                           72,
                           ybins));

      //3.1 Set up variable-sized bins for HE depth 3 (MonitorElement also requires phi bins to be entered in array format)
      float xbins[] = {-28.5, -27.5, -26.5, -16.5, -15.5, 15.5, 16.5, 26.5, 27.5, 28.5};
      //3.2
      depth.push_back(TH2F(("HE Depth 3 " + name + unitname).c_str(),
                           (name + " Depth 3 -- HE (" + unittitle + ")").c_str(),
                           // Use variable-sized eta bins
                           9,
                           xbins,
                           72,
                           ybins));

      //4.1 HO bins are fixed width, but cover a smaller eta range (-15 -> 15)
      depth.push_back(TH2F(("HO Depth 4 " + name + unitname).c_str(),
                           (name + " Depth 4 -- HO (" + unittitle + ")").c_str(),
                           31,
                           -15.5,
                           15.5,
                           72,
                           0.5,
                           72.5));

      for (unsigned int i = 0; i < depth.size(); ++i)
        depth[i].Draw("colz");

      setBinLabels(depth);  // set axis titles, special bins
    }

    //functions for making plot:
    void setBinLabels(std::vector<TH2F>& depth) {
      // Set labels for all depth histograms
      for (unsigned int i = 0; i < depth.size(); ++i) {
        depth[i].SetXTitle("i#eta");
        depth[i].SetYTitle("i#phi");
      }

      std::stringstream label;

      // set label on every other bin
      for (int i = -41; i <= -29; i = i + 2) {
        label << i;
        depth[0].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
        depth[1].GetXaxis()->SetBinLabel(i + 42, label.str().c_str());
        label.str("");
      }
      depth[0].GetXaxis()->SetBinLabel(14, "-29HE");
      depth[1].GetXaxis()->SetBinLabel(14, "-29HE");

      // offset by one for HE
      for (int i = -27; i <= 27; i = i + 2) {
        label << i;
        depth[0].GetXaxis()->SetBinLabel(i + 43, label.str().c_str());
        label.str("");
      }
      depth[0].GetXaxis()->SetBinLabel(72, "29HE");
      for (int i = 29; i <= 41; i = i + 2) {
        label << i;
        depth[0].GetXaxis()->SetBinLabel(i + 44, label.str().c_str());
        label.str("");
      }
      for (int i = 16; i <= 28; i = i + 2) {
        label << i - 43;
        depth[1].GetXaxis()->SetBinLabel(i, label.str().c_str());
        label.str("");
      }
      depth[1].GetXaxis()->SetBinLabel(29, "NULL");
      for (int i = 15; i <= 27; i = i + 2) {
        label << i;
        depth[1].GetXaxis()->SetBinLabel(i + 15, label.str().c_str());
        label.str("");
      }

      depth[1].GetXaxis()->SetBinLabel(44, "29HE");
      for (int i = 29; i <= 41; i = i + 2) {
        label << i;
        depth[1].GetXaxis()->SetBinLabel(i + 16, label.str().c_str());
        label.str("");
      }

      // HE depth 3 labels;
      depth[2].GetXaxis()->SetBinLabel(1, "-28");
      depth[2].GetXaxis()->SetBinLabel(2, "-27");
      depth[2].GetXaxis()->SetBinLabel(3, "Null");
      depth[2].GetXaxis()->SetBinLabel(4, "-16");
      depth[2].GetXaxis()->SetBinLabel(5, "Null");
      depth[2].GetXaxis()->SetBinLabel(6, "16");
      depth[2].GetXaxis()->SetBinLabel(7, "Null");
      depth[2].GetXaxis()->SetBinLabel(8, "27");
      depth[2].GetXaxis()->SetBinLabel(9, "28");
    }
  };
}  // namespace HcalObjRepresent
#endif