HcalRecHitsDQMClient.cc

Go to the documentation of this file.

#include "DQMOffline/Hcal/interface/HcalRecHitsDQMClient.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/Transition.h"
#include "FWCore/Utilities/interface/ESInputTag.h"

#include "DQMServices/Core/interface/DQMStore.h"

HcalRecHitsDQMClient::HcalRecHitsDQMClient(const edm::ParameterSet &iConfig)
    : conf_(iConfig),
      hcalDDDRecConstantsToken_{esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()},
      caloGeometryRunToken_{esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()} {
  outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile", "myfile.root");
  debug_ = false;
  verbose_ = false;
  dirName_ = iConfig.getParameter<std::string>("DQMDirName");
}

HcalRecHitsDQMClient::~HcalRecHitsDQMClient() {}

void HcalRecHitsDQMClient::beginJob() {}

void HcalRecHitsDQMClient::beginRun(const edm::Run &run, const edm::EventSetup &iSetup) {
  HcalDDDRecConstants const &hcons = iSetup.getData(hcalDDDRecConstantsToken_);
  maxDepthHB_ = hcons.getMaxDepth(0);
  maxDepthHE_ = hcons.getMaxDepth(1);
  maxDepthHF_ = std::max(hcons.getMaxDepth(2), 1);
  maxDepthHO_ = hcons.getMaxDepth(3);

  CaloGeometry const &geometry = iSetup.getData(caloGeometryRunToken_);
  const std::vector<DetId> &hbCells = geometry.getValidDetIds(DetId::Hcal, HcalBarrel);
  const std::vector<DetId> &heCells = geometry.getValidDetIds(DetId::Hcal, HcalEndcap);
  const std::vector<DetId> &hoCells = geometry.getValidDetIds(DetId::Hcal, HcalOuter);
  const std::vector<DetId> &hfCells = geometry.getValidDetIds(DetId::Hcal, HcalForward);

  nChannels_[1] = hbCells.size();
  nChannels_[2] = heCells.size();
  nChannels_[3] = hoCells.size();
  nChannels_[4] = hfCells.size();
  nChannels_[0] = nChannels_[1] + nChannels_[2] + nChannels_[3] + nChannels_[4];
  // avoid divide by zero
  for (unsigned i = 0; i < 5; ++i) {
    if (nChannels_[i] == 0)
      nChannels_[i] = 1;
  }

  // std::cout << "Channels HB:" << nChannels_[1] << " HE:" << nChannels_[2] <<
  // " HO:" << nChannels_[3] << " HF:" << nChannels_[4] << std::endl;

  // We hardcode the HF depths because in the dual readout configuration,
  // rechits are not defined for depths 3&4
  maxDepthHF_ = (maxDepthHF_ > 2 ? 2 : maxDepthHF_);  // We reatin the dynamic possibility
                                                      // that HF might have 0 or 1 depths

  maxDepthAll_ = (maxDepthHB_ + maxDepthHO_ > maxDepthHE_ ? maxDepthHB_ + maxDepthHO_ : maxDepthHE_);
  maxDepthAll_ = (maxDepthAll_ > maxDepthHF_ ? maxDepthAll_ : maxDepthHF_);
}

void HcalRecHitsDQMClient::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
  igetter.setCurrentFolder(dirName_);

  if (verbose_)
    std::cout << "\nrunClient" << std::endl;

  std::vector<MonitorElement *> hcalMEs;

  // Since out folders are fixed to three, we can just go over these three
  // folders i.e., CaloTowersD/CaloTowersTask, HcalRecHitsD/HcalRecHitTask,
  // NoiseRatesV/NoiseRatesTask.
  std::vector<std::string> fullPathHLTFolders = igetter.getSubdirs();
  for (unsigned int i = 0; i < fullPathHLTFolders.size(); i++) {
    if (verbose_)
      std::cout << "\nfullPath: " << fullPathHLTFolders[i] << std::endl;
    igetter.setCurrentFolder(fullPathHLTFolders[i]);

    std::vector<std::string> fullSubPathHLTFolders = igetter.getSubdirs();
    for (unsigned int j = 0; j < fullSubPathHLTFolders.size(); j++) {
      if (verbose_)
        std::cout << "fullSub: " << fullSubPathHLTFolders[j] << std::endl;

      if (strcmp(fullSubPathHLTFolders[j].c_str(), "HcalRecHitsD/HcalRecHitTask") == 0) {
        hcalMEs = igetter.getContents(fullSubPathHLTFolders[j]);
        if (verbose_)
          std::cout << "hltMES size : " << hcalMEs.size() << std::endl;
        if (!HcalRecHitsEndjob(hcalMEs))
          std::cout << "\nError in HcalRecHitsEndjob!" << std::endl << std::endl;
      }
    }
  }
}

// called after entering the HcalRecHitsD/HcalRecHitTask directory
// hcalMEs are within that directory
int HcalRecHitsDQMClient::HcalRecHitsEndjob(const std::vector<MonitorElement *> &hcalMEs) {
  MonitorElement *Nhf = nullptr;

  // Search for emap histograms, and collect them into this vector
  // All subdtectors are plotted together in these histograms. We only need to
  // look for different depths
  std::vector<MonitorElement *> emap_depths;

  // This vector is filled occupancy_maps identified by both subdetector and
  // depth
  std::vector<MonitorElement *> occupancy_maps;
  std::vector<std::string> occupancyID;

  // This vector is filled with emean_vs_ieta histograms, they are divided by
  // both subdetector and depth
  std::vector<MonitorElement *> emean_vs_ieta;

  // These are the only histograms filled in this module; however, the
  // histograms are created empty in HcalRecHitsAnalyzer occupancy_vs_ieta,
  // divided by both subdetector and depth
  std::vector<MonitorElement *> occupancy_vs_ieta;
  std::vector<std::string> occupancy_vs_ietaID;

  // RecHit_StatusWord & RecHit_Aux_StatusWord
  // Divided by subdectector
  std::vector<MonitorElement *> RecHit_StatusWord;
  std::vector<float> RecHit_StatusWord_Channels;
  std::vector<MonitorElement *> RecHit_Aux_StatusWord;
  std::vector<float> RecHit_Aux_StatusWord_Channels;

  for (unsigned int ih = 0; ih < hcalMEs.size(); ih++) {
    // N_HF is not special, it is just convient to get the total number of
    // events The number of entries in N_HF is equal to the number of events
    if (hcalMEs[ih]->getName() == "N_HF") {
      Nhf = hcalMEs[ih];
      continue;
    }

    // ***********************
    // * We fill the various MonitorElement vectors by searching for a matching
    // substring
    // * The methods that are used are agnostic to the ordering of vectors
    // ***********************

    if (hcalMEs[ih]->getName().find("emap_depth") != std::string::npos) {
      emap_depths.push_back(hcalMEs[ih]);
      continue;
    }

    if (hcalMEs[ih]->getName().find("occupancy_map_H") != std::string::npos) {
      occupancy_maps.push_back(hcalMEs[ih]);

      // Use occupancyID to save the subdetector and depth information
      // This will help preserve both indifference to vector ordering and
      // specific details of the detector topology The position in occupancyID
      // must correspond to the histogram position in occupancy_maps

      // Save the string after "occupancy_map_"

      std::string prefix = "occupancy_map_";

      occupancyID.push_back(hcalMEs[ih]->getName().substr(prefix.size()));

      continue;
    }

    if (hcalMEs[ih]->getName().find("emean_vs_ieta_H") != std::string::npos) {
      emean_vs_ieta.push_back(hcalMEs[ih]);
      continue;
    }

    if (hcalMEs[ih]->getName().find("occupancy_vs_ieta_H") != std::string::npos) {
      occupancy_vs_ieta.push_back(hcalMEs[ih]);

      // Use occupancy_vs_ietaID to save the subdetector and depth information
      // This will help preserve both indifference to vector ordering and
      // specific details of the detector topology The position in occupancyID
      // must correspond to the histogram position in occupancy_vs_ieta

      // Save the string after "occupancy_vs_ieta_"

      std::string prefix = "occupancy_vs_ieta_";

      occupancy_vs_ietaID.push_back(hcalMEs[ih]->getName().substr(prefix.size()));

      continue;
    }

    if (hcalMEs[ih]->getName().find("HcalRecHitTask_RecHit_StatusWord_H") != std::string::npos) {
      RecHit_StatusWord.push_back(hcalMEs[ih]);

      if (hcalMEs[ih]->getName().find("HB") != std::string::npos) {
        RecHit_StatusWord_Channels.push_back((float)nChannels_[1]);
      } else if (hcalMEs[ih]->getName().find("HE") != std::string::npos) {
        RecHit_StatusWord_Channels.push_back((float)nChannels_[2]);
      } else if (hcalMEs[ih]->getName().find("H0") != std::string::npos) {
        RecHit_StatusWord_Channels.push_back((float)nChannels_[3]);
      } else if (hcalMEs[ih]->getName().find("HF") != std::string::npos) {
        RecHit_StatusWord_Channels.push_back((float)nChannels_[4]);
      } else {
        RecHit_StatusWord_Channels.push_back(1.);
      }

      continue;
    }

    if (hcalMEs[ih]->getName().find("HcalRecHitTask_RecHit_Aux_StatusWord_H") != std::string::npos) {
      RecHit_Aux_StatusWord.push_back(hcalMEs[ih]);

      if (hcalMEs[ih]->getName().find("HB") != std::string::npos) {
        RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[1]);
      } else if (hcalMEs[ih]->getName().find("HE") != std::string::npos) {
        RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[2]);
      } else if (hcalMEs[ih]->getName().find("H0") != std::string::npos) {
        RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[3]);
      } else if (hcalMEs[ih]->getName().find("HF") != std::string::npos) {
        RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[4]);
      } else {
        RecHit_Aux_StatusWord_Channels.push_back(1.);
      }

      continue;
    }
  }

  // mean energies and occupancies evaluation

  double nevtot = Nhf->getEntries();  // Use the number of entries in the Nhf histogram to
                                      // give the total number of events

  if (verbose_)
    std::cout << "nevtot : " << nevtot << std::endl;

  // emap histograms are scaled by the number of events
  float fev = float(nevtot);
  double scaleBynevtot = 1 / fev;

  // In this and the following histogram vectors, recognize that the for-loop
  // index does not have to correspond to any particular depth
  for (unsigned int depthIdx = 0; depthIdx < emap_depths.size(); depthIdx++) {
    int nx = emap_depths[depthIdx]->getNbinsX();
    int ny = emap_depths[depthIdx]->getNbinsY();

    float cnorm;
    float enorm;

    for (int i = 1; i <= nx; i++) {
      for (int j = 1; j <= ny; j++) {
        cnorm = emap_depths[depthIdx]->getBinContent(i, j) * scaleBynevtot;
        enorm = emap_depths[depthIdx]->getBinError(i, j) * scaleBynevtot;
        emap_depths[depthIdx]->setBinContent(i, j, cnorm);
        emap_depths[depthIdx]->setBinError(i, j, enorm);
      }
    }
  }

  // occupancy_maps & matched occupancy_vs_ieta

  bool omatched = false;

  for (unsigned int occupancyIdx = 0; occupancyIdx < occupancy_maps.size(); occupancyIdx++) {
    int nx = occupancy_maps[occupancyIdx]->getNbinsX();
    int ny = occupancy_maps[occupancyIdx]->getNbinsY();

    float cnorm;
    float enorm;

    unsigned int vsIetaIdx = occupancy_vs_ieta.size();
    omatched = false;

    for (vsIetaIdx = 0; vsIetaIdx < occupancy_vs_ieta.size(); vsIetaIdx++) {
      if (occupancyID[occupancyIdx] == occupancy_vs_ietaID[vsIetaIdx]) {
        omatched = true;
        break;
      }
    }  // match occupancy_vs_ieta histogram

    for (int i = 1; i <= nx; i++) {
      for (int j = 1; j <= ny; j++) {
        cnorm = occupancy_maps[occupancyIdx]->getBinContent(i, j) * scaleBynevtot;
        enorm = occupancy_maps[occupancyIdx]->getBinError(i, j) * scaleBynevtot;
        occupancy_maps[occupancyIdx]->setBinContent(i, j, cnorm);
        occupancy_maps[occupancyIdx]->setBinError(i, j, enorm);
      }
    }

    // Fill occupancy_vs_ieta

    if (omatched) {
      // We run over all of the ieta values
      for (int ieta = -41; ieta <= 41; ieta++) {
        float phi_factor = 1.;
        float sumphi = 0.;
        float sumphie = 0.;

        if (ieta == 0)
          continue;  // ieta=0 is not defined

        phi_factor = phifactor(ieta);

        // the rechits occupancy map defines iphi as 0..71
        for (int iphi = 0; iphi <= 71; iphi++) {
          int binIeta = occupancy_maps[occupancyIdx]->getTH2F()->GetXaxis()->FindBin(float(ieta));
          int binIphi = occupancy_maps[occupancyIdx]->getTH2F()->GetYaxis()->FindBin(float(iphi));

          float content = occupancy_maps[occupancyIdx]->getBinContent(binIeta, binIphi);
          float econtent = occupancy_maps[occupancyIdx]->getBinError(binIeta, binIphi);

          sumphi += content;
          sumphie += econtent * econtent;
        }  // for loop over phi

        int ietabin = occupancy_vs_ieta[vsIetaIdx]->getTH1F()->GetXaxis()->FindBin(float(ieta));

        // fill occupancies vs ieta
        cnorm = sumphi / phi_factor;
        enorm = sqrt(sumphie) / phi_factor;
        occupancy_vs_ieta[vsIetaIdx]->setBinContent(ietabin, cnorm);
        occupancy_vs_ieta[vsIetaIdx]->setBinError(ietabin, enorm);

      }  // Fill occupancy_vs_ieta
    }    // if omatched
  }

  // Status Word
  // Normalized by number of events and by number of channels per subdetector as
  // well

  for (unsigned int StatusWordIdx = 0; StatusWordIdx < RecHit_StatusWord.size(); StatusWordIdx++) {
    int nx = RecHit_StatusWord[StatusWordIdx]->getNbinsX();

    float cnorm;
    float enorm;

    for (int i = 1; i <= nx; i++) {
      cnorm = RecHit_StatusWord[StatusWordIdx]->getBinContent(i) * scaleBynevtot /
              RecHit_StatusWord_Channels[StatusWordIdx];
      enorm =
          RecHit_StatusWord[StatusWordIdx]->getBinError(i) * scaleBynevtot / RecHit_StatusWord_Channels[StatusWordIdx];
      RecHit_StatusWord[StatusWordIdx]->setBinContent(i, cnorm);
      RecHit_StatusWord[StatusWordIdx]->setBinError(i, enorm);
    }
  }

  for (unsigned int AuxStatusWordIdx = 0; AuxStatusWordIdx < RecHit_Aux_StatusWord.size(); AuxStatusWordIdx++) {
    int nx = RecHit_Aux_StatusWord[AuxStatusWordIdx]->getNbinsX();

    float cnorm;
    float enorm;

    for (int i = 1; i <= nx; i++) {
      cnorm = RecHit_Aux_StatusWord[AuxStatusWordIdx]->getBinContent(i) * scaleBynevtot /
              RecHit_Aux_StatusWord_Channels[AuxStatusWordIdx];
      enorm = RecHit_Aux_StatusWord[AuxStatusWordIdx]->getBinError(i) * scaleBynevtot /
              RecHit_Aux_StatusWord_Channels[AuxStatusWordIdx];
      RecHit_Aux_StatusWord[AuxStatusWordIdx]->setBinContent(i, cnorm);
      RecHit_Aux_StatusWord[AuxStatusWordIdx]->setBinError(i, enorm);
    }
  }

  return 1;
}

float HcalRecHitsDQMClient::phifactor(int ieta) {
  float phi_factor_;

  if (ieta >= -20 && ieta <= 20) {
    phi_factor_ = 72.;
  } else {
    if (ieta >= 40 || ieta <= -40) {
      phi_factor_ = 18.;
    } else {
      phi_factor_ = 36.;
    }
  }

  return phi_factor_;
}

DEFINE_FWK_MODULE(HcalRecHitsDQMClient);