APVValidationPlots.cc

Go to the documentation of this file.

// system include files
#include <iostream>
#include <memory>
#include <sstream>
#include <vector>

// 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/ServiceRegistry/interface/Service.h"
#include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/Common/interface/DetSetVector.h"
#include "DataFormats/Common/interface/DetSetVectorNew.h"
#include "DQMServices/Core/interface/DQMStore.h"

// ROOT includes
#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TStyle.h"
#include "TTree.h"

//
// class decleration
//

class APVValidationPlots : public edm::one::EDAnalyzer<> {
public:
  typedef dqm::legacy::DQMStore DQMStore;
  typedef dqm::legacy::MonitorElement MonitorElement;

  explicit APVValidationPlots(const edm::ParameterSet&);

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

  const std::string infilename;
  const std::string outfilename;

  // ----------member data ---------------------------
};

//
// constructors and destructor
//
APVValidationPlots::APVValidationPlots(const edm::ParameterSet& iConfig)
    : infilename(iConfig.getUntrackedParameter<std::string>("inputFilename", "in.root")),
      outfilename(iConfig.getUntrackedParameter<std::string>("outputFilename", "out.root"))

{
  //now do what ever initialization is needed
}

//
// member functions
//

// ------------ method called to for each event  ------------
void APVValidationPlots::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {}

// ------------ method called once each job just after ending the event loop  ------------
void APVValidationPlots::endJob() {
  std::ostringstream oss;
  oss << 1;  //runNumber

  DQMStore* dqmStore = edm::Service<DQMStore>().operator->();
  dqmStore->setCurrentFolder("ChannelStatusPlots");

  // Initialize histograms
  std::vector<std::string> subDetName;
  std::vector<unsigned int> nLayers;
  std::vector<std::string> layerName;
  subDetName.push_back("");
  subDetName.push_back("TIB");
  subDetName.push_back("TID");
  subDetName.push_back("TOB");
  subDetName.push_back("TEC");
  nLayers.push_back(0);
  nLayers.push_back(4);
  nLayers.push_back(3);
  nLayers.push_back(6);
  nLayers.push_back(9);
  layerName.push_back("");
  layerName.push_back("Layer");
  layerName.push_back("Disk");
  layerName.push_back("Layer");
  layerName.push_back("Disk");

  std::string histoName;
  std::string histoTitle;
  // Histograms
  // indexes in these arrays are [SubDetId-2][LayerN]
  // histograms for [SubDetId-2][0] are global for the subdetector
  // histogram for [0][0] is global for the tracker
  TH2F* medianVsAbsoluteOccupancy[5][10];
  TH1F* medianOccupancy[5][10];
  TH1F* absoluteOccupancy[5][10];
  for (unsigned int i = 0; i < subDetName.size(); i++) {
    for (unsigned int j = 0; j <= nLayers[i]; j++) {
      histoName = "medianVsAbsoluteOccupancy" + subDetName[i];
      if (j != 0) {
        oss.str("");
        oss << j;
        histoName += layerName[i] + oss.str();
      }
      histoTitle = "Median APV occupancy vs. absolute APV occupancy";
      if (i != 0)
        histoTitle += " in " + subDetName[i];
      if (j != 0) {
        histoTitle += " " + layerName[i] + " " + oss.str();
      }
      MonitorElement* tmp = dqmStore->book2D(histoName.c_str(), histoTitle.c_str(), 1000, 0., 6., 1000, -1., 3.);
      medianVsAbsoluteOccupancy[i][j] = tmp->getTH2F();
      medianVsAbsoluteOccupancy[i][j]->Rebin2D(10, 10);
      medianVsAbsoluteOccupancy[i][j]->GetXaxis()->SetTitle("log_{10}(Abs. Occupancy)");
      medianVsAbsoluteOccupancy[i][j]->GetYaxis()->SetTitle("log_{10}(Median Occupancy)");
      //
      histoName = "medianOccupancy" + subDetName[i];
      if (j != 0) {
        oss.str("");
        oss << j;
        histoName += layerName[i] + oss.str();
      }
      histoTitle = "Median APV occupancy";
      if (i != 0)
        histoTitle += " in " + subDetName[i];
      if (j != 0) {
        histoTitle += " " + layerName[i] + " " + oss.str();
      }
      tmp = dqmStore->book1D(histoName.c_str(), histoTitle.c_str(), 1000, -1., 3.);
      medianOccupancy[i][j] = tmp->getTH1F();
      medianOccupancy[i][j]->GetXaxis()->SetTitle("log_{10}(Occupancy)");
      medianOccupancy[i][j]->GetYaxis()->SetTitle("APVs");
      //
      histoName = "absoluteOccupancy" + subDetName[i];
      if (j != 0) {
        oss.str("");
        oss << j;
        histoName += layerName[i] + oss.str();
      }
      histoTitle = "Absolute APV occupancy";
      if (i != 0)
        histoTitle += " in " + subDetName[i];
      if (j != 0) {
        histoTitle += " " + layerName[i] + " " + oss.str();
      }
      tmp = dqmStore->book1D(histoName.c_str(), histoTitle.c_str(), 1000, 0., 6.);
      absoluteOccupancy[i][j] = tmp->getTH1F();
      absoluteOccupancy[i][j]->GetXaxis()->SetTitle("log_{10}(Occupancy)");
      absoluteOccupancy[i][j]->GetYaxis()->SetTitle("APVs");
    }
  }

  TFile* infile = new TFile(infilename.c_str(), "READ");
  TTree* intree = (TTree*)infile->Get("moduleOccupancy");

  // Declaration of leaf types
  Int_t DetRawId;
  Int_t SubDetId;
  Int_t Layer_Ring;
  Int_t Disc;
  Int_t IsBack;
  Int_t IsExternalString;
  Int_t IsZMinusSide;
  Int_t RodStringPetal;
  Int_t IsStereo;
  Int_t ModulePosition;
  Int_t NumberOfStrips;
  Float_t APVGlobalPositionX;
  Float_t APVGlobalPositionY;
  Float_t APVGlobalPositionZ;
  Int_t APVNumber;
  Int_t APVAbsoluteOccupancy;
  Double_t APVMedianOccupancy;
  intree->SetBranchAddress("DetRawId", &DetRawId);
  intree->SetBranchAddress("SubDetId", &SubDetId);
  intree->SetBranchAddress("Layer_Ring", &Layer_Ring);
  intree->SetBranchAddress("Disc", &Disc);
  intree->SetBranchAddress("IsBack", &IsBack);
  intree->SetBranchAddress("IsExternalString", &IsExternalString);
  intree->SetBranchAddress("IsZMinusSide", &IsZMinusSide);
  intree->SetBranchAddress("RodStringPetal", &RodStringPetal);
  intree->SetBranchAddress("IsStereo", &IsStereo);
  intree->SetBranchAddress("ModuleNumber", &ModulePosition);
  intree->SetBranchAddress("NumberOfStrips", &NumberOfStrips);
  intree->SetBranchAddress("APVGlobalPositionX", &APVGlobalPositionX);
  intree->SetBranchAddress("APVGlobalPositionY", &APVGlobalPositionY);
  intree->SetBranchAddress("APVGlobalPositionZ", &APVGlobalPositionZ);
  intree->SetBranchAddress("APVNumber", &APVNumber);
  intree->SetBranchAddress("APVAbsoluteOccupancy", &APVAbsoluteOccupancy);
  intree->SetBranchAddress("APVMedianOccupancy", &APVMedianOccupancy);

  for (int i = 0; i < intree->GetEntries(); i++) {
    intree->GetEntry(i);

    double logMedianOccupancy = -1;
    double logAbsoluteOccupancy = -1;

    if (APVMedianOccupancy > 0)
      logMedianOccupancy = log10(APVMedianOccupancy);
    if (APVAbsoluteOccupancy > 0)
      logAbsoluteOccupancy = log10(APVAbsoluteOccupancy);

    // The layer/disk information is stored in Layer_Ring for TIB/TOB and in Disc for TID/TEC
    unsigned int layer = 0;
    if (SubDetId == 3 || SubDetId == 5)
      layer = Layer_Ring;
    else
      layer = Disc;

    // Fill histograms for all the tracker
    medianVsAbsoluteOccupancy[0][0]->Fill(logAbsoluteOccupancy, logMedianOccupancy);
    medianOccupancy[0][0]->Fill(logMedianOccupancy);
    absoluteOccupancy[0][0]->Fill(logAbsoluteOccupancy);
    // Fill summary histograms for each subdetector
    medianVsAbsoluteOccupancy[SubDetId - 2][0]->Fill(logAbsoluteOccupancy, logMedianOccupancy);
    medianOccupancy[SubDetId - 2][0]->Fill(logMedianOccupancy);
    absoluteOccupancy[SubDetId - 2][0]->Fill(logAbsoluteOccupancy);
    // Fill histograms for each layer/disk
    medianVsAbsoluteOccupancy[SubDetId - 2][layer]->Fill(logAbsoluteOccupancy, logMedianOccupancy);
    medianOccupancy[SubDetId - 2][layer]->Fill(logMedianOccupancy);
    absoluteOccupancy[SubDetId - 2][layer]->Fill(logAbsoluteOccupancy);
  }

  dqmStore->cd();
  dqmStore->save(outfilename, "ChannelStatusPlots");
}

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