MuonME0Validation/plugins/MuonME0DigisHarvestor.cc

0001 // system include files
0002 #include <memory>
0003
0004 // user include files
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ServiceRegistry/interface/Service.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "TGraphAsymmErrors.h"
0012
0013 /// Geometry
0014 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0015 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0016
0017 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0018 #include "Geometry/GEMGeometry/interface/ME0EtaPartition.h"
0019 #include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
0020 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0021
0022 #include "DQMServices/Core/interface/DQMStore.h"
0023
0024 /// Log messages
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027
0028 #include "Validation/MuonME0Validation/plugins/MuonME0DigisHarvestor.h"
0029
0030 MuonME0DigisHarvestor::MuonME0DigisHarvestor(const edm::ParameterSet &ps) {
0031   dbe_path_ = std::string("MuonME0DigisV/ME0DigisTask/");
0032 }
0033
0034 MuonME0DigisHarvestor::~MuonME0DigisHarvestor() {}
0035
0036 TProfile *MuonME0DigisHarvestor::ComputeEff(TH1F *num, TH1F *denum, std::string nameHist) {
0037   std::string name = "eff_" + nameHist;
0038   std::string title = "Digi Efficiency" + std::string(num->GetTitle());
0039   TProfile *efficHist = new TProfile(name.c_str(),
0040                                      title.c_str(),
0041                                      denum->GetXaxis()->GetNbins(),
0042                                      denum->GetXaxis()->GetXmin(),
0043                                      denum->GetXaxis()->GetXmax());
0044
0045   for (int i = 1; i <= denum->GetNbinsX(); i++) {
0046     double nNum = num->GetBinContent(i);
0047     double nDenum = denum->GetBinContent(i);
0048     if (nDenum == 0 || nNum == 0) {
0049       continue;
0050     }
0051     if (nNum > nDenum) {
0052       double temp = nDenum;
0053       nDenum = nNum;
0054       nNum = temp;
0055       edm::LogWarning("MuonME0DigisHarvestor")
0056           << "Alert! specific bin's num is bigger than denum " << i << " " << nNum << " " << nDenum;
0057     }
0058     const double effVal = nNum / nDenum;
0059     efficHist->SetBinContent(i, effVal);
0060     efficHist->SetBinEntries(i, 1);
0061     efficHist->SetBinError(i, 0);
0062     const double errLo = TEfficiency::ClopperPearson((int)nDenum, (int)nNum, 0.683, false);
0063     const double errUp = TEfficiency::ClopperPearson((int)nDenum, (int)nNum, 0.683, true);
0064     const double errVal = (effVal - errLo > errUp - effVal) ? effVal - errLo : errUp - effVal;
0065     efficHist->SetBinError(i, sqrt(effVal * effVal + errVal * errVal));
0066   }
0067   return efficHist;
0068 }
0069
0070 void MuonME0DigisHarvestor::ProcessBooking(
0071     DQMStore::IBooker &ibooker, DQMStore::IGetter &ig, std::string nameHist, TH1F *num, TH1F *den) {
0072   if (num != nullptr && den != nullptr) {
0073     TProfile *profile = ComputeEff(num, den, nameHist);
0074
0075     TString x_axis_title = TString(num->GetXaxis()->GetTitle());
0076     TString title = TString::Format("Digi Efficiency;%s;Eff.", x_axis_title.Data());
0077
0078     profile->SetTitle(title.Data());
0079     ibooker.bookProfile(profile->GetName(), profile);
0080
0081     delete profile;
0082
0083   } else {
0084     edm::LogWarning("MuonME0DigisHarvestor") << "Can not find histograms";
0085     if (num == nullptr)
0086       edm::LogWarning("MuonME0DigisHarvestor") << "num not found";
0087     if (den == nullptr)
0088       edm::LogWarning("MuonME0DigisHarvestor") << "den not found";
0089   }
0090   return;
0091 }
0092
0093 TH1F *MuonME0DigisHarvestor::ComputeBKG(TH1F *hist1, TH1F *hist2, std::string nameHist) {
0094   std::string name = "rate_" + nameHist;
0095   hist1->SetName(name.c_str());
0096   for (int bin = 1; bin <= hist1->GetNbinsX(); ++bin) {
0097     double R_min = hist1->GetBinCenter(bin) - 0.5 * hist1->GetBinWidth(bin);
0098     double R_max = hist1->GetBinCenter(bin) + 0.5 * hist1->GetBinWidth(bin);
0099
0100     double Area = TMath::Pi() * (R_max * R_max - R_min * R_min);
0101     hist1->SetBinContent(bin, (hist1->GetBinContent(bin)) / Area);
0102     hist1->SetBinError(bin, (hist1->GetBinError(bin)) / Area);
0103   }
0104
0105   int nEvts = hist2->GetEntries();
0106   float scale = 6 * 2 * nEvts * 3 * 25e-9;  // New redigitizer saves hits only in the BX range: [-1,+1], so the
0107                                             // number of background hits has to be divided by 3
0108   hist1->Scale(1.0 / scale);
0109   return hist1;
0110 }
0111
0112 void MuonME0DigisHarvestor::ProcessBookingBKG(
0113     DQMStore::IBooker &ibooker, DQMStore::IGetter &ig, std::string nameHist, TH1F *hist1, TH1F *hist2) {
0114   if (hist1 != nullptr && hist2 != nullptr) {
0115     TH1F *rate = ComputeBKG(hist1, hist2, nameHist);
0116
0117     TString x_axis_title = TString(hist1->GetXaxis()->GetTitle());
0118     TString origTitle = TString(hist1->GetTitle());
0119     TString title = TString::Format((origTitle + ";%s;Rate [Hz/cm^{2}]").Data(), x_axis_title.Data());
0120
0121     rate->SetTitle(title.Data());
0122     ibooker.book1D(rate->GetName(), rate);
0123
0124   } else {
0125     edm::LogWarning("MuonME0DigisHarvestor") << "Can not find histograms";
0126     if (hist1 == nullptr)
0127       edm::LogWarning("MuonME0DigisHarvestor") << "num not found";
0128     if (hist2 == nullptr)
0129       edm::LogWarning("MuonME0DigisHarvestor") << "den not found";
0130   }
0131   return;
0132 }
0133
0134 void MuonME0DigisHarvestor::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &ig) {
0135   ig.setCurrentFolder(dbe_path_);
0136
0137   const char *l_suffix[6] = {"_l1", "_l2", "_l3", "_l4", "_l5", "_l6"};
0138   const char *r_suffix[2] = {"-1", "1"};
0139
0140   TString eta_label_den_tot = TString(dbe_path_) + "me0_strip_dg_den_eta_tot";
0141   TString eta_label_num_tot = TString(dbe_path_) + "me0_strip_dg_num_eta_tot";
0142   if (ig.get(eta_label_num_tot.Data()) != nullptr && ig.get(eta_label_den_tot.Data()) != nullptr) {
0143     TH1F *num_vs_eta_tot = (TH1F *)ig.get(eta_label_num_tot.Data())->getTH1F()->Clone();
0144     num_vs_eta_tot->Sumw2();
0145     TH1F *den_vs_eta_tot = (TH1F *)ig.get(eta_label_den_tot.Data())->getTH1F()->Clone();
0146     den_vs_eta_tot->Sumw2();
0147
0148     ProcessBooking(ibooker, ig, "me0_strip_dg_eta_tot", num_vs_eta_tot, den_vs_eta_tot);
0149
0150     delete num_vs_eta_tot;
0151     delete den_vs_eta_tot;
0152
0153   } else
0154     edm::LogWarning("MuonME0DigisHarvestor")
0155         << "Can not find histograms: " << eta_label_num_tot << " or " << eta_label_den_tot;
0156
0157   for (int i = 0; i < 2; i++) {
0158     for (int j = 0; j < 6; j++) {
0159       TString eta_label_den = TString(dbe_path_) + "me0_strip_dg_den_eta" + r_suffix[i] + l_suffix[j];
0160       TString eta_label_num = TString(dbe_path_) + "me0_strip_dg_num_eta" + r_suffix[i] + l_suffix[j];
0161
0162       if (ig.get(eta_label_num.Data()) != nullptr && ig.get(eta_label_den.Data()) != nullptr) {
0163         TH1F *num_vs_eta = (TH1F *)ig.get(eta_label_num.Data())->getTH1F()->Clone();
0164         num_vs_eta->Sumw2();
0165         TH1F *den_vs_eta = (TH1F *)ig.get(eta_label_den.Data())->getTH1F()->Clone();
0166         den_vs_eta->Sumw2();
0167
0168         std::string r_s = r_suffix[i];
0169         std::string l_s = l_suffix[j];
0170         std::string name = "me0_strip_dg_eta" + r_s + l_s;
0171         ProcessBooking(ibooker, ig, name, num_vs_eta, den_vs_eta);
0172
0173         delete num_vs_eta;
0174         delete den_vs_eta;
0175
0176       } else
0177         edm::LogWarning("MuonME0DigisHarvestor")
0178             << "Can not find histograms: " << eta_label_num << " " << eta_label_den;
0179     }
0180   }
0181
0182   TString label_eleBkg = TString(dbe_path_) + "me0_strip_dg_bkgElePos_radius";
0183   TString label_neuBkg = TString(dbe_path_) + "me0_strip_dg_bkgNeutral_radius";
0184   TString label_totBkg = TString(dbe_path_) + "me0_strip_dg_bkg_radius_tot";
0185   TString label_evts = TString(dbe_path_) + "num_evts";
0186
0187   if (ig.get(label_evts.Data()) != nullptr) {
0188     TH1F *numEvts = (TH1F *)ig.get(label_evts.Data())->getTH1F()->Clone();
0189
0190     if (ig.get(label_eleBkg.Data()) != nullptr) {
0191       TH1F *eleBkg = (TH1F *)ig.get(label_eleBkg.Data())->getTH1F()->Clone();
0192       eleBkg->Sumw2();
0193       ProcessBookingBKG(ibooker, ig, "me0_strip_dg_elePosBkg_rad", eleBkg, numEvts);
0194
0195       delete eleBkg;
0196     }
0197     if (ig.get(label_neuBkg.Data()) != nullptr) {
0198       TH1F *neuBkg = (TH1F *)ig.get(label_neuBkg.Data())->getTH1F()->Clone();
0199       neuBkg->Sumw2();
0200       ProcessBookingBKG(ibooker, ig, "me0_strip_dg_neuBkg_rad", neuBkg, numEvts);
0201
0202       delete neuBkg;
0203     }
0204     if (ig.get(label_totBkg.Data()) != nullptr) {
0205       TH1F *totBkg = (TH1F *)ig.get(label_totBkg.Data())->getTH1F()->Clone();
0206       totBkg->Sumw2();
0207       ProcessBookingBKG(ibooker, ig, "me0_strip_dg_totBkg_rad", totBkg, numEvts);
0208
0209       delete totBkg;
0210     }
0211
0212     delete numEvts;
0213   }
0214 }
0215
0216 // define this as a plug-in
0217 DEFINE_FWK_MODULE(MuonME0DigisHarvestor);