CTPPS/plugins/TotemRPDQMHarvester.cc

0001 /****************************************************************************
0002 *
0003 * This is a part of TotemDQM and TOTEM offline software.
0004 * Authors:
0005 *   Jan Kašpar (jan.kaspar@gmail.com)
0006 *
0007 ****************************************************************************/
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0014 #include "DQMServices/Core/interface/DQMStore.h"
0015
0016 #include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h"
0017
0018 //----------------------------------------------------------------------------------------------------
0019
0020 class TotemRPDQMHarvester : public DQMEDHarvester {
0021 public:
0022   TotemRPDQMHarvester(const edm::ParameterSet &ps);
0023   ~TotemRPDQMHarvester() override;
0024
0025 protected:
0026   void dqmEndJob(DQMStore::IBooker &, DQMStore::IGetter &) override {}
0027
0028   void dqmEndLuminosityBlock(DQMStore::IBooker &,
0029                              DQMStore::IGetter &,
0030                              const edm::LuminosityBlock &,
0031                              const edm::EventSetup &) override;
0032
0033 private:
0034   void MakeHitNumberRatios(unsigned int id, DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter);
0035
0036   void MakePlaneEfficiencyHistograms(unsigned int id,
0037                                      DQMStore::IBooker &ibooker,
0038                                      DQMStore::IGetter &igetter,
0039                                      bool &rpPlotInitialized);
0040 };
0041
0042 //----------------------------------------------------------------------------------------------------
0043 //----------------------------------------------------------------------------------------------------
0044
0045 using namespace std;
0046 using namespace edm;
0047
0048 //----------------------------------------------------------------------------------------------------
0049
0050 TotemRPDQMHarvester::TotemRPDQMHarvester(const edm::ParameterSet &ps) {}
0051
0052 //----------------------------------------------------------------------------------------------------
0053
0054 TotemRPDQMHarvester::~TotemRPDQMHarvester() {}
0055
0056 //----------------------------------------------------------------------------------------------------
0057
0058 void TotemRPDQMHarvester::MakeHitNumberRatios(unsigned int id, DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0059   // get source histogram
0060   string path;
0061   TotemRPDetId(id).rpName(path, TotemRPDetId::nPath);
0062
0063   MonitorElement *activity = igetter.get(path + "/activity in planes (2D)");
0064
0065   if (!activity)
0066     return;
0067
0068   // book new histogram, if not yet done
0069   const string hit_ratio_name = "hit ratio in hot spot";
0070   MonitorElement *hit_ratio = igetter.get(path + "/" + hit_ratio_name);
0071
0072   if (hit_ratio == nullptr) {
0073     ibooker.setCurrentFolder(path);
0074     string title;
0075     TotemRPDetId(id).rpName(title, TotemRPDetId::nFull);
0076     hit_ratio = ibooker.book1D(hit_ratio_name, title + ";plane;N_hits(320<strip<440) / N_hits(all)", 10, -0.5, 9.5);
0077   } else {
0078     hit_ratio->getTH1F()->Reset();
0079   }
0080
0081   // calculate ratios
0082   TAxis *y_axis = activity->getTH2F()->GetYaxis();
0083   for (int bix = 1; bix <= activity->getNbinsX(); ++bix) {
0084     double S_full = 0., S_sel = 0.;
0085     for (int biy = 1; biy <= activity->getNbinsY(); ++biy) {
0086       double c = activity->getBinContent(bix, biy);
0087       double s = y_axis->GetBinCenter(biy);
0088
0089       S_full += c;
0090
0091       if (s > 320. && s < 440.)
0092         S_sel += c;
0093     }
0094
0095     double r = (S_full > 0.) ? S_sel / S_full : 0.;
0096
0097     hit_ratio->setBinContent(bix, r);
0098   }
0099 }
0100
0101 //----------------------------------------------------------------------------------------------------
0102
0103 void TotemRPDQMHarvester::MakePlaneEfficiencyHistograms(unsigned int id,
0104                                                         DQMStore::IBooker &ibooker,
0105                                                         DQMStore::IGetter &igetter,
0106                                                         bool &rpPlotInitialized) {
0107   TotemRPDetId detId(id);
0108
0109   // get source histograms
0110   string path;
0111   detId.planeName(path, TotemRPDetId::nPath);
0112
0113   MonitorElement *efficiency_num = igetter.get(path + "/efficiency num");
0114   MonitorElement *efficiency_den = igetter.get(path + "/efficiency den");
0115
0116   if (!efficiency_num || !efficiency_den)
0117     return;
0118
0119   // book new plane histogram, if not yet done
0120   const string efficiency_name = "efficiency";
0121   MonitorElement *efficiency = igetter.get(path + "/" + efficiency_name);
0122
0123   if (efficiency == nullptr) {
0124     string title;
0125     detId.planeName(title, TotemRPDetId::nFull);
0126     TAxis *axis = efficiency_den->getTH1()->GetXaxis();
0127     ibooker.setCurrentFolder(path);
0128     efficiency = ibooker.book1D(
0129         efficiency_name, title + ";track position   (mm)", axis->GetNbins(), axis->GetXmin(), axis->GetXmax());
0130   } else {
0131     efficiency->getTH1F()->Reset();
0132   }
0133
0134   // book new RP histogram, if not yet done
0135   CTPPSDetId rpId = detId.rpId();
0136   rpId.rpName(path, TotemRPDetId::nPath);
0137   const string rp_efficiency_name = "plane efficiency";
0138   MonitorElement *rp_efficiency = igetter.get(path + "/" + rp_efficiency_name);
0139
0140   if (rp_efficiency == nullptr) {
0141     string title;
0142     rpId.rpName(title, TotemRPDetId::nFull);
0143     TAxis *axis = efficiency_den->getTH1()->GetXaxis();
0144     ibooker.setCurrentFolder(path);
0145     rp_efficiency = ibooker.book2D(rp_efficiency_name,
0146                                    title + ";plane;track position   (mm)",
0147                                    10,
0148                                    -0.5,
0149                                    9.5,
0150                                    axis->GetNbins(),
0151                                    axis->GetXmin(),
0152                                    axis->GetXmax());
0153     rpPlotInitialized = true;
0154   } else {
0155     if (!rpPlotInitialized)
0156       rp_efficiency->getTH2F()->Reset();
0157     rpPlotInitialized = true;
0158   }
0159
0160   // calculate and fill efficiencies
0161   for (signed int bi = 1; bi <= efficiency->getNbinsX(); bi++) {
0162     double num = efficiency_num->getBinContent(bi);
0163     double den = efficiency_den->getBinContent(bi);
0164
0165     if (den > 0) {
0166       double p = num / den;
0167       double p_unc = sqrt(p * (1. - p) / den);
0168       efficiency->setBinContent(bi, p);
0169       efficiency->setBinError(bi, p_unc);
0170
0171       int pl_bi = detId.plane() + 1;
0172       rp_efficiency->setBinContent(pl_bi, bi, p);
0173     } else {
0174       efficiency->setBinContent(bi, 0.);
0175     }
0176   }
0177 }
0178
0179 //----------------------------------------------------------------------------------------------------
0180
0181 void TotemRPDQMHarvester::dqmEndLuminosityBlock(DQMStore::IBooker &ibooker,
0182                                                 DQMStore::IGetter &igetter,
0183                                                 const edm::LuminosityBlock &,
0184                                                 const edm::EventSetup &) {
0185   // loop over arms
0186   for (unsigned int arm = 0; arm < 2; arm++) {
0187     // loop over stations
0188     for (unsigned int st = 0; st < 3; st += 2) {
0189       // loop over RPs
0190       for (unsigned int rp = 0; rp < 6; ++rp) {
0191         if (st == 2) {
0192           // unit 220-nr is not equipped
0193           if (rp <= 2)
0194             continue;
0195
0196           // RP 220-fr-hr contains pixels
0197           if (rp == 3)
0198             continue;
0199         }
0200
0201         TotemRPDetId rpId(arm, st, rp);
0202
0203         MakeHitNumberRatios(rpId, ibooker, igetter);
0204
0205         bool rpPlotInitialized = false;
0206
0207         // loop over planes
0208         for (unsigned int pl = 0; pl < 10; ++pl) {
0209           TotemRPDetId plId(arm, st, rp, pl);
0210
0211           MakePlaneEfficiencyHistograms(plId, ibooker, igetter, rpPlotInitialized);
0212         }
0213       }
0214     }
0215   }
0216 }
0217
0218 //----------------------------------------------------------------------------------------------------
0219
0220 DEFINE_FWK_MODULE(TotemRPDQMHarvester);