TotemRPDQMHarvester.cc

Go to the documentation of this file.

/****************************************************************************
*
* This is a part of TotemDQM and TOTEM offline software.
* Authors:
*   Jan Kašpar (jan.kaspar@gmail.com)
*
****************************************************************************/

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

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

#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h"

//----------------------------------------------------------------------------------------------------
 
class TotemRPDQMHarvester: public DQMEDHarvester
{
  public:
    TotemRPDQMHarvester(const edm::ParameterSet& ps);
    ~TotemRPDQMHarvester() override;
  
  protected:
    void dqmEndJob(DQMStore::IBooker &, DQMStore::IGetter &) override {}

    void dqmEndLuminosityBlock(DQMStore::IBooker &, DQMStore::IGetter &, const edm::LuminosityBlock &, const edm::EventSetup &) override;

  private:
    void MakeHitNumberRatios(unsigned int id, DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter);

    void MakePlaneEfficiencyHistograms(unsigned int id, DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter, bool &rpPlotInitialized);
};

//----------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------

using namespace std;
using namespace edm;

//----------------------------------------------------------------------------------------------------

TotemRPDQMHarvester::TotemRPDQMHarvester(const edm::ParameterSet& ps)
{
}

//----------------------------------------------------------------------------------------------------

TotemRPDQMHarvester::~TotemRPDQMHarvester()
{
}
    
//----------------------------------------------------------------------------------------------------

void TotemRPDQMHarvester::MakeHitNumberRatios(unsigned int id, DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter)
{
  // get source histogram
  string path;
  TotemRPDetId(id).rpName(path, TotemRPDetId::nPath);

  MonitorElement *activity = igetter.get(path + "/activity in planes (2D)");

  if (!activity)
    return;

  // book new histogram, if not yet done
  const string hit_ratio_name = "hit ratio in hot spot";
  MonitorElement *hit_ratio = igetter.get(path + "/" + hit_ratio_name);

  if (hit_ratio == nullptr)
  {
    ibooker.setCurrentFolder(path);
    string title;
    TotemRPDetId(id).rpName(title, TotemRPDetId::nFull);
    hit_ratio = ibooker.book1D(hit_ratio_name, title+";plane;N_hits(320<strip<440) / N_hits(all)", 10, -0.5, 9.5);
  } else {
    hit_ratio->getTH1F()->Reset();
  }

  // calculate ratios
  TAxis *y_axis = activity->getTH2F()->GetYaxis();
  for (int bix = 1; bix <= activity->getNbinsX(); ++bix)
  {
    double S_full = 0., S_sel = 0.;
    for (int biy = 1; biy <= activity->getNbinsY(); ++biy)
    {
      double c = activity->getBinContent(bix, biy);
      double s = y_axis->GetBinCenter(biy);
      
      S_full += c;

      if (s > 320. && s < 440.)
        S_sel += c;
    }

    double r = (S_full > 0.) ? S_sel / S_full : 0.;

    hit_ratio->setBinContent(bix, r);
  }
}

//----------------------------------------------------------------------------------------------------

void TotemRPDQMHarvester::MakePlaneEfficiencyHistograms(unsigned int id, DQMStore::IBooker& ibooker,
  DQMStore::IGetter& igetter, bool &rpPlotInitialized)
{
  TotemRPDetId detId(id);

  // get source histograms
  string path;
  detId.planeName(path, TotemRPDetId::nPath);

  MonitorElement *efficiency_num = igetter.get(path + "/efficiency num");
  MonitorElement *efficiency_den = igetter.get(path + "/efficiency den");

  if (!efficiency_num || !efficiency_den)
    return;

  // book new plane histogram, if not yet done
  const string efficiency_name = "efficiency";
  MonitorElement *efficiency = igetter.get(path + "/" + efficiency_name);

  if (efficiency == nullptr)
  {
    string title;
    detId.planeName(title, TotemRPDetId::nFull);
    TAxis *axis = efficiency_den->getTH1()->GetXaxis();
    ibooker.setCurrentFolder(path);
    efficiency = ibooker.book1D(efficiency_name, title+";track position   (mm)", axis->GetNbins(), axis->GetXmin(), axis->GetXmax());
  } else {
    efficiency->getTH1F()->Reset();
  }

  // book new RP histogram, if not yet done
  CTPPSDetId rpId = detId.getRPId();
  rpId.rpName(path, TotemRPDetId::nPath);
  const string rp_efficiency_name = "plane efficiency";
  MonitorElement *rp_efficiency = igetter.get(path + "/" + rp_efficiency_name);
  
  if (rp_efficiency == nullptr)
  {
    string title;
    rpId.rpName(title, TotemRPDetId::nFull);
    TAxis *axis = efficiency_den->getTH1()->GetXaxis();
    ibooker.setCurrentFolder(path);
    rp_efficiency = ibooker.book2D(rp_efficiency_name, title+";plane;track position   (mm)",
      10, -0.5, 9.5, axis->GetNbins(), axis->GetXmin(), axis->GetXmax());
    rpPlotInitialized = true;
  } else {
    if (!rpPlotInitialized)
      rp_efficiency->getTH2F()->Reset();
    rpPlotInitialized = true;
  }

  // calculate and fill efficiencies
  for (signed int bi = 1; bi <= efficiency->getNbinsX(); bi++)
  {
    double num = efficiency_num->getBinContent(bi);
    double den = efficiency_den->getBinContent(bi);

    if (den > 0)
    {
      double p = num / den;
      double p_unc = sqrt(p * (1. - p) / den);
      efficiency->setBinContent(bi, p);
      efficiency->setBinError(bi, p_unc);

      int pl_bi = detId.plane() + 1;
      rp_efficiency->setBinContent(pl_bi, bi, p);
    } else {
      efficiency->setBinContent(bi, 0.);
    }
  }
}

//----------------------------------------------------------------------------------------------------

void TotemRPDQMHarvester::dqmEndLuminosityBlock(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter,
  const edm::LuminosityBlock &, const edm::EventSetup &)
{
  // loop over arms
  for (unsigned int arm = 0; arm < 2; arm++)
  {
    // loop over stations
    for (unsigned int st = 0; st < 3; st += 2)
    {
      // loop over RPs
      for (unsigned int rp = 0; rp < 6; ++rp)
      {
        if (st == 2)
        {
          // unit 220-nr is not equipped
          if (rp <= 2)
            continue;

          // RP 220-fr-hr contains pixels
          if (rp == 3)
            continue;
        }

        TotemRPDetId rpId(arm, st, rp);

        MakeHitNumberRatios(rpId, ibooker, igetter);

        bool rpPlotInitialized = false;
        
        // loop over planes
        for (unsigned int pl = 0; pl < 10; ++pl)
        {
          TotemRPDetId plId(arm, st, rp, pl);

          MakePlaneEfficiencyHistograms(plId, ibooker, igetter, rpPlotInitialized);
        }
      }
    }
  }
}

//----------------------------------------------------------------------------------------------------

DEFINE_FWK_MODULE(TotemRPDQMHarvester);