Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:42 

0001 
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author G. Mila - INFN Torino
0006  */
0007 
0008 #include <DQMOffline/Muon/interface/MuonTestSummary.h>
0009 
0010 // Framework
0011 #include <FWCore/Framework/interface/Event.h>
0012 #include <FWCore/Framework/interface/EventSetup.h>
0013 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0014 
0015 #include "DQMServices/Core/interface/DQMStore.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/Framework/interface/Run.h"
0018 
0019 #include "DQMOffline/Muon/test/langauFit.C"
0020 #include <string>
0021 
0022 using namespace edm;
0023 using namespace std;
0024 
0025 MuonTestSummary::MuonTestSummary(const edm::ParameterSet &ps) {
0026   // parameter initialization for kinematics test
0027   etaExpected = ps.getParameter<double>("etaExpected");
0028   phiExpected = ps.getParameter<double>("phiExpected");
0029   chi2Fraction = ps.getParameter<double>("chi2Fraction");
0030   chi2Spread = ps.getParameter<double>("chi2Spread");
0031   resEtaSpread_tkGlb = ps.getParameter<double>("resEtaSpread_tkGlb");
0032   resEtaSpread_glbSta = ps.getParameter<double>("resEtaSpread_glbSta");
0033   resPhiSpread_tkGlb = ps.getParameter<double>("resPhiSpread_tkGlb");
0034   resPhiSpread_glbSta = ps.getParameter<double>("resPhiSpread_glbSta");
0035   resOneOvPSpread_tkGlb = ps.getParameter<double>("resOneOvPSpread_tkGlb");
0036   resOneOvPSpread_glbSta = ps.getParameter<double>("resOneOvPSpread_glbSta");
0037   pullEtaSpread = ps.getParameter<double>("pullEtaSpread");
0038   pullPhiSpread = ps.getParameter<double>("pullPhiSpread");
0039   pullOneOvPSpread = ps.getParameter<double>("pullOneOvPSpread");
0040   resChargeLimit_tkGlb = ps.getParameter<double>("resChargeLimit_tkGlb");
0041   resChargeLimit_glbSta = ps.getParameter<double>("resChargeLimit_glbSta");
0042   resChargeLimit_tkSta = ps.getParameter<double>("resChargeLimit_tkSta");
0043   numMatchedExpected_min = ps.getParameter<double>("numMatchedExpected_min");
0044   numMatchedExpected_max = ps.getParameter<double>("numMatchedExpected_max");
0045   matchesFractionDt_min = ps.getParameter<double>("matchesFractionDt_min");
0046   matchesFractionDt_max = ps.getParameter<double>("matchesFractionDt_max");
0047   matchesFractionCsc_min = ps.getParameter<double>("matchesFractionCsc_min");
0048   matchesFractionCsc_max = ps.getParameter<double>("matchesFractionCsc_max");
0049   resSegmTrack_rms_min = ps.getParameter<double>("resSegmTrack_rms_min");
0050   resSegmTrack_rms_max = ps.getParameter<double>("resSegmTrack_rms_max");
0051   resSegmTrack_mean_min = ps.getParameter<double>("resSegmTrack_mean_min");
0052   resSegmTrack_mean_max = ps.getParameter<double>("resSegmTrack_mean_max");
0053   expPeakEcalS9_min = ps.getParameter<double>("expPeakEcalS9_min");
0054   expPeakEcalS9_max = ps.getParameter<double>("expPeakEcalS9_max");
0055   expPeakHadS9_min = ps.getParameter<double>("expPeakHadS9_min");
0056   expPeakHadS9_max = ps.getParameter<double>("expPeakHadS9_max");
0057   expMultiplicityGlb_max = ps.getParameter<double>("expMultiplicityGlb_max");
0058   expMultiplicityTk_max = ps.getParameter<double>("expMultiplicityTk_max");
0059   expMultiplicitySta_max = ps.getParameter<double>("expMultiplicitySta_max");
0060   expMultiplicityGlb_min = ps.getParameter<double>("expMultiplicityGlb_min");
0061   expMultiplicityTk_min = ps.getParameter<double>("expMultiplicityTk_min");
0062   expMultiplicitySta_min = ps.getParameter<double>("expMultiplicitySta_min");
0063 }
0064 
0065 MuonTestSummary::~MuonTestSummary() {}
0066 
0067 void MuonTestSummary::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0068   // BOOKING NEW HISTOGRAMS
0069   ibooker.cd();
0070   ibooker.setCurrentFolder("Muons/TestSummary");
0071 
0072   // kinematics test report
0073   kinematicsSummaryMap = ibooker.book2D("kinematicsSummaryMap", "Kinematics test summary", 5, 1, 6, 3, 1, 4);
0074   kinematicsSummaryMap->setAxisTitle("track monitored", 1);
0075   kinematicsSummaryMap->setBinLabel(1, "GLB", 1);
0076   kinematicsSummaryMap->setBinLabel(2, "TKfromGLB", 1);
0077   kinematicsSummaryMap->setBinLabel(3, "STAfromGLB", 1);
0078   kinematicsSummaryMap->setBinLabel(4, "TK", 1);
0079   kinematicsSummaryMap->setBinLabel(5, "STA", 1);
0080   kinematicsSummaryMap->setAxisTitle("parameter tested", 2);
0081   kinematicsSummaryMap->setBinLabel(1, "#chi^{2}", 2);
0082   kinematicsSummaryMap->setBinLabel(2, "#eta", 2);
0083   kinematicsSummaryMap->setBinLabel(3, "#phi", 2);
0084 
0085   //chi2 kinematics quality test report
0086   chi2TestSummaryMap = ibooker.book2D("chi2TestSummaryMap", "#chi2 quality test summary", 5, 1, 6, 5, 1, 6);
0087   chi2TestSummaryMap->setAxisTitle("track monitored", 1);
0088   chi2TestSummaryMap->setBinLabel(1, "GLB", 1);
0089   chi2TestSummaryMap->setBinLabel(2, "TKfromGLB", 1);
0090   chi2TestSummaryMap->setBinLabel(3, "STAfromGLB", 1);
0091   chi2TestSummaryMap->setBinLabel(4, "TK", 1);
0092   chi2TestSummaryMap->setBinLabel(5, "STA", 1);
0093   chi2TestSummaryMap->setAxisTitle("parameter tested", 2);
0094   chi2TestSummaryMap->setBinLabel(1, "#chi^{2}", 2);
0095   chi2TestSummaryMap->setBinLabel(2, "#eta", 2);
0096   chi2TestSummaryMap->setBinLabel(3, "#phi", 2);
0097   chi2TestSummaryMap->setBinLabel(4, "#pt", 2);
0098   chi2TestSummaryMap->setBinLabel(5, "#q", 2);
0099 
0100   // residuals test report
0101   residualsSummaryMap = ibooker.book2D("residualsSummaryMap", "Residuals test summary", 4, 1, 5, 4, 1, 5);
0102   residualsSummaryMap->setAxisTitle("residuals", 1);
0103   residualsSummaryMap->setBinLabel(1, "TK-GLB", 1);
0104   residualsSummaryMap->setBinLabel(2, "GLB-STA", 1);
0105   residualsSummaryMap->setBinLabel(3, "TK-STA", 1);
0106   residualsSummaryMap->setBinLabel(4, "TK-STA Pull", 1);
0107   residualsSummaryMap->setAxisTitle("parameter tested", 2);
0108   residualsSummaryMap->setBinLabel(1, "#eta", 2);
0109   residualsSummaryMap->setBinLabel(2, "#phi", 2);
0110   residualsSummaryMap->setBinLabel(3, "1/p", 2);
0111   residualsSummaryMap->setBinLabel(4, "q", 2);
0112 
0113   // muonId test report
0114   muonIdSummaryMap = ibooker.book2D("muonIdSummaryMap", "muonId test summary", 4, 1, 5, 5, 1, 6);
0115   muonIdSummaryMap->setAxisTitle("muons", 1);
0116   muonIdSummaryMap->setBinLabel(1, "GLB DT", 1);
0117   muonIdSummaryMap->setBinLabel(2, "GLB CSC", 1);
0118   muonIdSummaryMap->setBinLabel(3, "TK DT", 1);
0119   muonIdSummaryMap->setBinLabel(4, "TK CSC", 1);
0120   muonIdSummaryMap->setAxisTitle("tests", 2);
0121   muonIdSummaryMap->setBinLabel(1, "#assSeg", 2);
0122   muonIdSummaryMap->setBinLabel(2, "x mean", 2);
0123   muonIdSummaryMap->setBinLabel(3, "x rms", 2);
0124   muonIdSummaryMap->setBinLabel(4, "y mean", 2);
0125   muonIdSummaryMap->setBinLabel(5, "y rms", 2);
0126 
0127   // energy test report
0128   energySummaryMap = ibooker.book2D("energySummaryMap", "Energy deposits test summary", 3, 1, 4, 3, 1, 4);
0129   energySummaryMap->setAxisTitle("muons", 1);
0130   energySummaryMap->setBinLabel(1, "GLB", 1);
0131   energySummaryMap->setBinLabel(2, "TK", 1);
0132   energySummaryMap->setBinLabel(3, "STA", 1);
0133   energySummaryMap->setAxisTitle("calorimeter tested", 2);
0134   energySummaryMap->setBinLabel(1, "ECAL", 2);
0135   energySummaryMap->setBinLabel(2, "HAD", 2);
0136   energySummaryMap->setBinLabel(3, "H0", 2);
0137 
0138   // multiplicity tests report
0139   multiplicitySummaryMap = ibooker.book1D("multiplicitySummaryMap", "muon multiplicity test summary", 3, 1, 4);
0140   multiplicitySummaryMap->setAxisTitle("muon");
0141   multiplicitySummaryMap->setBinLabel(1, "GLB");
0142   multiplicitySummaryMap->setBinLabel(2, "TK");
0143   multiplicitySummaryMap->setBinLabel(3, "STA");
0144 
0145   // summary test report
0146   ibooker.setCurrentFolder("Muons/EventInfo");
0147   summaryReport = ibooker.bookFloat("reportSummary");
0148 
0149   summaryReportMap = ibooker.book2D("reportSummaryMap", "Muon Report Summary Map", 3, 1, 4, 7, 1, 8);
0150   summaryReportMap->setAxisTitle("muons", 1);
0151   summaryReportMap->setBinLabel(1, "GLB", 1);
0152   summaryReportMap->setBinLabel(2, "TK", 1);
0153   summaryReportMap->setBinLabel(3, "STA", 1);
0154   summaryReportMap->setAxisTitle("test", 2);
0155   summaryReportMap->setBinLabel(1, "#chi^{2}/Df", 2);
0156   summaryReportMap->setBinLabel(2, "#eta", 2);
0157   summaryReportMap->setBinLabel(3, "#phi", 2);
0158   summaryReportMap->setBinLabel(4, "residuals", 2);
0159   summaryReportMap->setBinLabel(5, "muonId", 2);
0160   summaryReportMap->setBinLabel(6, "energyDeposits", 2);
0161   summaryReportMap->setBinLabel(7, "multiplicity", 2);
0162 
0163   ibooker.setCurrentFolder("Muons/EventInfo/reportSummaryContents");
0164   theSummaryContents.push_back(ibooker.bookFloat("kinematics_GLB"));
0165   theSummaryContents.push_back(ibooker.bookFloat("muonId_GLB"));
0166   theSummaryContents.push_back(ibooker.bookFloat("residuals_GLB"));
0167   theSummaryContents.push_back(ibooker.bookFloat("GLB"));
0168   theSummaryContents.push_back(ibooker.bookFloat("kinematics_TK"));
0169   theSummaryContents.push_back(ibooker.bookFloat("muonId_TK"));
0170   theSummaryContents.push_back(ibooker.bookFloat("residuals_TK"));
0171   theSummaryContents.push_back(ibooker.bookFloat("TK"));
0172   theSummaryContents.push_back(ibooker.bookFloat("kinematics_STA"));
0173   theSummaryContents.push_back(ibooker.bookFloat("residuals_STA"));
0174   theSummaryContents.push_back(ibooker.bookFloat("STA"));
0175   theSummaryContents.push_back(ibooker.bookFloat("energyDeposits"));
0176   theSummaryContents.push_back(ibooker.bookFloat("multiplicity"));
0177 
0178   // certification report
0179   ibooker.setCurrentFolder("Muons/EventInfo");
0180   summaryCertification = ibooker.bookFloat("CertificationSummary");
0181   summaryCertification->Fill(-1);
0182 
0183   summaryCertificationMap =
0184       ibooker.book2D("CertificationSummaryMap", "Muon Certification Summary Map", 9, 1, 10, 7, 1, 8);
0185   summaryCertificationMap->setAxisTitle("muons", 1);
0186   summaryCertificationMap->setBinLabel(1, "GLB_Tot", 1);
0187   summaryCertificationMap->setBinLabel(2, "TK_Tot", 1);
0188   summaryCertificationMap->setBinLabel(3, "STA_tot", 1);
0189   summaryCertificationMap->setBinLabel(4, "GLB_B", 1);
0190   summaryCertificationMap->setBinLabel(5, "TK_B", 1);
0191   summaryCertificationMap->setBinLabel(6, "STA_B", 1);
0192   summaryCertificationMap->setBinLabel(7, "GLB_EC", 1);
0193   summaryCertificationMap->setBinLabel(8, "TK_EC", 1);
0194   summaryCertificationMap->setBinLabel(9, "STA_EC", 1);
0195   summaryCertificationMap->setAxisTitle("test", 2);
0196   summaryCertificationMap->setBinLabel(1, "#chi^{2}/Df", 2);
0197   summaryCertificationMap->setBinLabel(2, "#eta", 2);
0198   summaryCertificationMap->setBinLabel(3, "#phi", 2);
0199   summaryCertificationMap->setBinLabel(4, "residuals", 2);
0200   summaryCertificationMap->setBinLabel(5, "muonId", 2);
0201   summaryCertificationMap->setBinLabel(6, "energyDeposits", 2);
0202   summaryCertificationMap->setBinLabel(7, "multiplicity", 2);
0203 
0204   ibooker.setCurrentFolder("Muons/EventInfo/CertificationContents");
0205   theCertificationContents.push_back(ibooker.bookFloat("GLB_Tot"));
0206   theCertificationContents.push_back(ibooker.bookFloat("STA_Tot"));
0207   theCertificationContents.push_back(ibooker.bookFloat("TK_Tot"));
0208   theCertificationContents.push_back(ibooker.bookFloat("GLB_B"));
0209   theCertificationContents.push_back(ibooker.bookFloat("STA_B"));
0210   theCertificationContents.push_back(ibooker.bookFloat("TK_B"));
0211   theCertificationContents.push_back(ibooker.bookFloat("GLB_EC"));
0212   theCertificationContents.push_back(ibooker.bookFloat("STA_EC"));
0213   theCertificationContents.push_back(ibooker.bookFloat("TK_EC"));
0214 
0215   for (unsigned int icert = 0; icert < theCertificationContents.size(); icert++) {
0216     theCertificationContents[icert]->Fill(-1);
0217   }
0218 
0219   // initialisation of histo bins
0220   for (int xBin = 1; xBin <= 5; xBin++) {
0221     for (int yBin = 1; yBin <= 3; yBin++) {
0222       kinematicsSummaryMap->Fill(xBin, yBin, 1);
0223     }
0224     for (int yBin = 1; yBin <= 5; yBin++) {
0225       chi2TestSummaryMap->Fill(xBin, yBin, 1);
0226     }
0227   }
0228   for (int xBin = 1; xBin <= residualsSummaryMap->getNbinsX(); xBin++) {
0229     for (int yBin = 1; yBin <= residualsSummaryMap->getNbinsY(); yBin++) {
0230       residualsSummaryMap->Fill(xBin, yBin, 1);
0231     }
0232   }
0233   residualsSummaryMap->setBinContent(4, 4, 1);  //not used for now
0234 
0235   for (int xBin = 1; xBin <= muonIdSummaryMap->getNbinsX(); xBin++) {
0236     for (int yBin = 1; yBin <= muonIdSummaryMap->getNbinsY(); yBin++) {
0237       muonIdSummaryMap->Fill(xBin, yBin, 1);
0238     }
0239   }
0240   for (int xBin = 1; xBin <= 3; xBin++) {
0241     for (int yBin = 1; yBin <= 3; yBin++) {
0242       energySummaryMap->Fill(xBin, yBin, 1);
0243     }
0244   }
0245   for (int xBin = 1; xBin <= 3; xBin++) {
0246     multiplicitySummaryMap->Fill(xBin, 1);
0247   }
0248   for (int xBin = 1; xBin <= 3; xBin++) {
0249     for (int yBin = 1; yBin <= 7; yBin++) {
0250       summaryReportMap->Fill(xBin, yBin, 1);
0251     }
0252   }
0253   for (int xBin = 1; xBin <= 9; xBin++) {
0254     for (int yBin = 1; yBin <= 7; yBin++) {
0255       summaryCertificationMap->Fill(xBin, yBin, 1);
0256     }
0257   }
0258   ////////////////////////////////////////////////////////////
0259   /// DO OTHER OPERATIONS WITH HISTOGRAMS
0260 
0261   // fill the kinematics report summary
0262   doKinematicsTests(igetter, "GlbMuon_Glb_", 1);
0263   doKinematicsTests(igetter, "GlbMuon_Tk_", 2);
0264   doKinematicsTests(igetter, "GlbMuon_Sta_", 3);
0265   doKinematicsTests(igetter, "TkMuon_", 4);
0266   doKinematicsTests(igetter, "StaMuon_", 5);
0267 
0268   // fill the residuals report summary
0269   doResidualsTests(igetter, "TkGlb", "eta", 1);
0270   doResidualsTests(igetter, "GlbSta", "eta", 2);
0271   doResidualsTests(igetter, "TkSta", "eta", 3);
0272   doResidualsTests(igetter, "TkGlb", "phi", 1);
0273   doResidualsTests(igetter, "GlbSta", "phi", 2);
0274   doResidualsTests(igetter, "TkSta", "phi", 3);
0275   doResidualsTests(igetter, "TkGlb", "oneOverp", 1);
0276   doResidualsTests(igetter, "GlbSta", "oneOverp", 2);
0277   doResidualsTests(igetter, "TkSta", "oneOverp", 3);
0278   doResidualsTests(igetter, "GlbMuon", "qComparison", -1);
0279 
0280   // fill the muonID report summary
0281   doMuonIDTests(igetter);
0282 
0283   // fill the energy report summary
0284   doEnergyTests(igetter, "ecalS9PointingMuDepositedEnergy_", "Glb_muons", 1);
0285   doEnergyTests(igetter, "hadS9PointingMuDepositedEnergy_", "Glb_muons", 1);
0286   doEnergyTests(igetter, "hoS9PointingMuDepositedEnergy_", "Glb_muons", 1);
0287   doEnergyTests(igetter, "ecalS9PointingMuDepositedEnergy_", "Tk_muons", 2);
0288   doEnergyTests(igetter, "hadS9PointingMuDepositedEnergy_", "Tk_muons", 2);
0289   doEnergyTests(igetter, "hoS9PointingMuDepositedEnergy_", "Tk_muons", 2);
0290   doEnergyTests(igetter, "ecalS9PointingMuDepositedEnergy_", "Sta_muons", 3);
0291   doEnergyTests(igetter, "hadS9PointingMuDepositedEnergy_", "Sta_muons", 3);
0292   doEnergyTests(igetter, "hoS9PointingMuDepositedEnergy_", "Sta_muons", 3);
0293 
0294   // fill the multiplicity test summary
0295   doMultiplicityTests(igetter);
0296 
0297   // fill the final report summary
0298 
0299   //-- modified GH
0300   double residualsSummary = 0;
0301   //put the TRK-STA resid & pulls in the first bin ("GLB")
0302   //then the GLB-TRK and GLB-STA residuals in the 2nd and 3rd
0303   for (int i = 3; i <= residualsSummaryMap->getNbinsX(); i++)
0304     for (int j = 1; j <= residualsSummaryMap->getNbinsY(); j++)
0305       residualsSummary += residualsSummaryMap->getBinContent(i, j);
0306   residualsSummary /= 2 * residualsSummaryMap->getNbinsY();
0307   summaryReportMap->setBinContent(1, 4, residualsSummary);
0308 
0309   residualsSummary = 0;
0310   for (int i = 1; i <= 1; i++)
0311     for (int j = 1; j <= residualsSummaryMap->getNbinsY(); j++)
0312       residualsSummary += residualsSummaryMap->getBinContent(i, j);
0313   residualsSummary /= 1 * residualsSummaryMap->getNbinsY();
0314   summaryReportMap->setBinContent(2, 4, residualsSummary);
0315 
0316   residualsSummary = 0;
0317   for (int i = 2; i <= 2; i++)
0318     for (int j = 1; j <= residualsSummaryMap->getNbinsY(); j++)
0319       residualsSummary += residualsSummaryMap->getBinContent(i, j);
0320   residualsSummary /= 1 * residualsSummaryMap->getNbinsY();
0321   summaryReportMap->setBinContent(3, 4, residualsSummary);
0322 
0323   //--
0324 
0325   //-- modified GH
0326   float idtest = 0;
0327   for (int i = 1; i <= 2; i++)
0328     for (int j = 1; j <= 5; j++) {
0329       if (j == 3 || j == 5)
0330         continue;  //ignore pull widths for now
0331       idtest += muonIdSummaryMap->getBinContent(i, j);
0332     }
0333   //  idtest/=10.;
0334   idtest /= 6.;
0335   summaryReportMap->setBinContent(1, 5, idtest);
0336   idtest = 0;
0337   for (int i = 3; i <= 4; i++)
0338     for (int j = 1; j <= 5; j++) {
0339       if (j == 3 || j == 5)
0340         continue;  //ignore pull widths for now
0341       idtest += muonIdSummaryMap->getBinContent(i, j);
0342     }
0343   //  idtest/=10.;
0344   idtest /= 6.;
0345   summaryReportMap->setBinContent(2, 5, idtest);
0346   summaryReportMap->setBinContent(3, 5, -1.0 / 6.0);
0347   //--
0348 
0349   summaryReportMap->setBinContent(
0350       1, 6, double(energySummaryMap->getBinContent(1, 1) + energySummaryMap->getBinContent(1, 2)) / 2.0);
0351   summaryReportMap->setBinContent(
0352       2, 6, double(energySummaryMap->getBinContent(2, 1) + energySummaryMap->getBinContent(2, 2)) / 2.0);
0353   summaryReportMap->setBinContent(
0354       3, 6, double(energySummaryMap->getBinContent(3, 1) + energySummaryMap->getBinContent(3, 2)) / 2.0);
0355   summaryReportMap->setBinContent(1, 7, multiplicitySummaryMap->getBinContent(1));
0356   summaryReportMap->setBinContent(2, 7, multiplicitySummaryMap->getBinContent(2));
0357   summaryReportMap->setBinContent(3, 7, multiplicitySummaryMap->getBinContent(3));
0358 
0359   double kinematics_GLB = double(summaryReportMap->getBinContent(1, 1) + summaryReportMap->getBinContent(1, 2) +
0360                                  summaryReportMap->getBinContent(1, 3)) /
0361                           3.0;
0362   theSummaryContents[0]->Fill(kinematics_GLB);
0363   double muonId_GLB = double(summaryReportMap->getBinContent(1, 5));
0364   theSummaryContents[1]->Fill(muonId_GLB);
0365   double residuals_GLB = double(summaryReportMap->getBinContent(1, 4));
0366   theSummaryContents[2]->Fill(residuals_GLB);
0367   double GLB = (kinematics_GLB + muonId_GLB + residuals_GLB) / 3.0;
0368   theSummaryContents[3]->Fill(GLB);
0369 
0370   double kinematics_TK = double(summaryReportMap->getBinContent(2, 1) + summaryReportMap->getBinContent(2, 2) +
0371                                 summaryReportMap->getBinContent(2, 3)) /
0372                          3.0;
0373   theSummaryContents[4]->Fill(kinematics_TK);
0374   double muonId_TK = double(summaryReportMap->getBinContent(2, 5));
0375   theSummaryContents[5]->Fill(muonId_TK);
0376   double residuals_TK = double(summaryReportMap->getBinContent(2, 4));
0377   theSummaryContents[6]->Fill(residuals_TK);
0378   double TK = double(kinematics_TK + muonId_TK + residuals_TK) / 3.0;
0379   theSummaryContents[7]->Fill(TK);
0380 
0381   double kinematics_STA = double(summaryReportMap->getBinContent(3, 1) + summaryReportMap->getBinContent(3, 2) +
0382                                  summaryReportMap->getBinContent(3, 3)) /
0383                           3.0;
0384   theSummaryContents[8]->Fill(kinematics_STA);
0385   double residuals_STA = double(summaryReportMap->getBinContent(3, 4));
0386   theSummaryContents[9]->Fill(residuals_STA);
0387   double STA = double(kinematics_STA + residuals_STA) / 2.0;
0388   theSummaryContents[10]->Fill(STA);
0389   double energyDeposits = double(summaryReportMap->getBinContent(1, 6) + summaryReportMap->getBinContent(2, 6) +
0390                                  summaryReportMap->getBinContent(3, 6)) /
0391                           3.0;
0392   theSummaryContents[11]->Fill(energyDeposits);
0393   double multiplicity = double(summaryReportMap->getBinContent(1, 7) + summaryReportMap->getBinContent(2, 7) +
0394                                summaryReportMap->getBinContent(3, 7)) /
0395                         3.0;
0396   theSummaryContents[12]->Fill(multiplicity);
0397 
0398   summaryReport->Fill((GLB + TK + STA + energyDeposits + multiplicity) / 5.0);
0399 
0400   //global barrel:
0401   float muonIDsummary = 0;
0402   //  for(int i=2; i<=5; i++)
0403   //     muonIDsummary += muonIdSummaryMap->getBinContent(2, i);
0404   //  summaryCertificationMap->setBinContent(4, 5, muonIDsummary/4.);
0405   //for now, just report the mean:
0406   muonIDsummary += muonIdSummaryMap->getBinContent(1, 2);
0407   muonIDsummary += muonIdSummaryMap->getBinContent(1, 4);
0408   summaryCertificationMap->setBinContent(4, 5, muonIDsummary / 2.);
0409 
0410   //global EC:
0411   muonIDsummary = 0;
0412   //  for(int i=2; i<=5; i++)
0413   //  muonIDsummary += muonIdSummaryMap->getBinContent(2, i);
0414   // summaryCertificationMap->setBinContent(7, 5, muonIDsummary/4.);
0415   muonIDsummary += muonIdSummaryMap->getBinContent(2, 2);
0416   muonIDsummary += muonIdSummaryMap->getBinContent(2, 4);
0417   summaryCertificationMap->setBinContent(7, 5, muonIDsummary / 2.);
0418 
0419   //tracker barrel:
0420   muonIDsummary = 0;
0421   //  for(int i=2; i<=5; i++)
0422   //    muonIDsummary += muonIdSummaryMap->getBinContent(3, i);
0423   //  summaryCertificationMap->setBinContent(5, 5, muonIDsummary/4.);
0424   muonIDsummary += muonIdSummaryMap->getBinContent(3, 2);
0425   muonIDsummary += muonIdSummaryMap->getBinContent(3, 4);
0426   summaryCertificationMap->setBinContent(5, 5, muonIDsummary / 2.);
0427 
0428   //tracker EC:
0429   muonIDsummary = 0;
0430   //  for(int i=2; i<=5; i++)
0431   //    muonIDsummary += muonIdSummaryMap->getBinContent(4, i);
0432   //  summaryCertificationMap->setBinContent(8, 5, muonIDsummary/4.);
0433   muonIDsummary += muonIdSummaryMap->getBinContent(4, 2);
0434   muonIDsummary += muonIdSummaryMap->getBinContent(4, 4);
0435   summaryCertificationMap->setBinContent(8, 5, muonIDsummary / 2.);
0436 
0437   double muonId_GLB_B = double(summaryCertificationMap->getBinContent(4, 5));
0438   theCertificationContents[3]->Fill(muonId_GLB_B);
0439   double muonId_GLB_EC = double(summaryCertificationMap->getBinContent(7, 5));
0440   theCertificationContents[6]->Fill(muonId_GLB_EC);
0441 
0442   double muonId_TK_B = double(summaryCertificationMap->getBinContent(5, 5));
0443   theCertificationContents[5]->Fill(muonId_TK_B);
0444   double muonId_TK_EC = double(summaryCertificationMap->getBinContent(8, 5));
0445   theCertificationContents[8]->Fill(muonId_TK_EC);
0446 }
0447 
0448 void MuonTestSummary::doKinematicsTests(DQMStore::IGetter &igetter, string muonType, int bin) {
0449   // chi2 test
0450   string path = "Muons/MuonRecoAnalyzer/" + muonType + "chi2OverDf";
0451   MonitorElement *chi2Histo = igetter.get(path);
0452 
0453   if (chi2Histo) {
0454     TH1F *chi2Histo_root = chi2Histo->getTH1F();
0455     if (chi2Histo_root->GetEntries() > 20) {
0456       //Standard QT based on fraction of events below and above a cut
0457       LogTrace(metname) << "chi2 kin test based on fraction for " << muonType << endl;
0458       int maxBin = chi2Histo_root->GetMaximumBin();
0459       if (chi2Histo_root->Integral(maxBin + 1, chi2Histo_root->GetNbinsX()) != 0) {
0460         double fraction = double(chi2Histo_root->Integral(1, maxBin)) /
0461                           double(chi2Histo_root->Integral(maxBin + 1, chi2Histo_root->GetNbinsX()));
0462         LogTrace(metname) << "chi2 fraction for " << muonType << " : " << fraction << " must be within "
0463                           << chi2Fraction - chi2Spread << "," << chi2Fraction + chi2Spread << endl;
0464         if (fraction > (chi2Fraction - chi2Spread) && fraction < (chi2Fraction + chi2Spread))
0465           kinematicsSummaryMap->setBinContent(bin, 1, 1);
0466         else
0467           kinematicsSummaryMap->setBinContent(bin, 1, 0);
0468       }
0469     } else {
0470       LogTrace(metname) << "[MuonTestSummary]: Test of Chi2 Kin not performed for " << muonType
0471                         << " because # entries < 20 ";
0472     }
0473   }
0474 
0475   // pseudorapidity test
0476   path = "Muons/MuonRecoAnalyzer/" + muonType + "eta";
0477   MonitorElement *etaHisto = igetter.get(path);
0478 
0479   if (etaHisto) {
0480     TH1F *etaHisto_root = etaHisto->getTH1F();
0481     if (etaHisto_root->GetEntries() > 20) {
0482       //Standard QT based on fraction of events below and above a cut
0483       LogTrace(metname) << "eta kin test based on fraction for " << muonType << endl;
0484       double binSize =
0485           (etaHisto_root->GetXaxis()->GetXmax() - etaHisto_root->GetXaxis()->GetXmin()) / etaHisto_root->GetNbinsX();
0486       int binZero = int((0 - etaHisto_root->GetXaxis()->GetXmin()) / binSize);
0487       if (etaHisto_root->Integral(1, binZero - 1) != 0 &&
0488           etaHisto_root->Integral(binZero, etaHisto_root->GetNbinsX()) != 0) {
0489         double symmetryFactor = double(etaHisto_root->Integral(1, binZero - 1)) /
0490                                 double(etaHisto_root->Integral(binZero, etaHisto_root->GetNbinsX()));
0491         double errSymmetryFactor =
0492             symmetryFactor * sqrt(1.0 / double(etaHisto_root->Integral(1, binZero - 1)) +
0493                                   1.0 / double(etaHisto_root->Integral(binZero, etaHisto_root->GetNbinsX())));
0494         LogTrace(metname) << "eta symmetryFactor for " << muonType << " : " << symmetryFactor
0495                           << " (expected :" << etaExpected << ")" << endl;
0496         LogTrace(metname) << "eta errSymmetryFactor for " << muonType << " : " << errSymmetryFactor << endl;
0497         double tParameter;
0498         if ((symmetryFactor - etaExpected) > 0)
0499           tParameter = double(symmetryFactor - etaExpected) / errSymmetryFactor;
0500         else
0501           tParameter = double(-symmetryFactor + etaExpected) / errSymmetryFactor;
0502         LogTrace(metname) << "eta tParameter for " << muonType << " : " << tParameter << " (expected < 1.95)" << endl;
0503         if (tParameter < 1.95)  //2sigma rejection
0504           kinematicsSummaryMap->setBinContent(bin, 2, 1);
0505         else
0506           kinematicsSummaryMap->setBinContent(bin, 2, 0);
0507       }
0508 
0509     } else {
0510       LogTrace(metname) << "[MuonTestSummary]: Test of Eta Kin not performed for " << muonType
0511                         << " because # entries < 20 ";
0512     }
0513   }
0514 
0515   // phi test
0516   path = "Muons/MuonRecoAnalyzer/" + muonType + "phi";
0517   MonitorElement *phiHisto = igetter.get(path);
0518 
0519   if (phiHisto) {
0520     TH1F *phiHisto_root = phiHisto->getTH1F();
0521     if (phiHisto_root->GetEntries() > 20) {
0522       //Standard QT based on fraction of events below and above a cut
0523       LogTrace(metname) << "phi kin test based on fraction for " << muonType << endl;
0524       double binSize =
0525           (phiHisto_root->GetXaxis()->GetXmax() - phiHisto_root->GetXaxis()->GetXmin()) / phiHisto_root->GetNbinsX();
0526       int binZero = int((0 - phiHisto_root->GetXaxis()->GetXmin()) / binSize);
0527       if (phiHisto_root->Integral(binZero + 1, phiHisto_root->GetNbinsX()) != 0 &&
0528           phiHisto_root->Integral(1, binZero) != 0) {
0529         double symmetryFactor = double(phiHisto_root->Integral(binZero + 1, phiHisto_root->GetNbinsX())) /
0530                                 double(phiHisto_root->Integral(1, binZero));
0531         double errSymmetryFactor =
0532             symmetryFactor * sqrt(1.0 / double(phiHisto_root->Integral(binZero + 1, phiHisto_root->GetNbinsX())) +
0533                                   1.0 / double(phiHisto_root->Integral(1, binZero)));
0534         LogTrace(metname) << "phi symmetryFactor for " << muonType << " : " << symmetryFactor
0535                           << "(phi expected :" << phiExpected << ")" << endl;
0536         LogTrace(metname) << "phi errSymmetryFactor for " << muonType << " : " << errSymmetryFactor << endl;
0537         double tParameter;
0538         if ((symmetryFactor - phiExpected) > 0)
0539           tParameter = double(symmetryFactor - phiExpected) / errSymmetryFactor;
0540         else
0541           tParameter = double(-symmetryFactor + phiExpected) / errSymmetryFactor;
0542         LogTrace(metname) << "phi tParameter for " << muonType << " : " << tParameter << " (expected < 1.95)" << endl;
0543         if (tParameter < 1.95)  //2sigma rejection
0544           kinematicsSummaryMap->setBinContent(bin, 3, 1);
0545         else
0546           kinematicsSummaryMap->setBinContent(bin, 3, 0);
0547       }
0548 
0549     } else {
0550       LogTrace(metname) << "[MuonTestSummary]: Test of Phi Kin not performed for " << muonType
0551                         << " because # entries < 20 ";
0552     }
0553   }
0554 }
0555 //--GH new
0556 void MuonTestSummary::GaussFit(
0557     string type, string parameter, MonitorElement *Histo, float &mean, float &mean_err, float &sigma, float &sigma_err) {
0558   // Gaussian Fit
0559   float statMean = Histo->getMean(1);
0560   float statSigma = Histo->getRMS(1);
0561   TH1F *histo_root = Histo->getTH1F();
0562   if (histo_root->GetEntries() > 20) {
0563     TF1 *gfit = new TF1("Gaussian", "gaus", (statMean - (2 * statSigma)), (statMean + (2 * statSigma)));
0564     try {
0565       histo_root->Fit(gfit, "Q0");
0566     } catch (cms::Exception &iException) {
0567       edm::LogError(metname) << "[MuonTestSummary]: Exception when fitting Res_" << type << "_" << parameter;
0568       mean = 1;
0569       mean_err = 1;
0570       sigma = 1;
0571       sigma_err = 1;
0572       return;
0573     }
0574     if (gfit) {
0575       mean = gfit->GetParameter(1);
0576       mean_err = gfit->GetParErrors()[2];
0577       sigma = gfit->GetParameter(2);
0578       sigma_err = gfit->GetParErrors()[2];
0579       LogTrace(metname) << "Gaussian fit mean: " << mean << " +- " << mean_err << " for " << type << "_" << parameter
0580                         << endl;
0581       LogTrace(metname) << "Gaussina fit sigma: " << sigma << " +- " << sigma_err << " for " << type << "_" << parameter
0582                         << endl;
0583     }
0584   } else {
0585     LogTrace(metname) << "[MuonTestSummary]: Test of  Res_" << type << "_" << parameter
0586                       << " not performed because # entries < 20 ";
0587     //auto-pass if not enough events.
0588     mean = 1;
0589     mean_err = 1;
0590     sigma = 1;
0591     sigma_err = 1;
0592   }
0593 }
0594 void MuonTestSummary::doResidualsTests(DQMStore::IGetter &igetter, string type, string parameter, int bin) {
0595   // residuals test
0596   if (type != "GlbMuon") {
0597     string path = "Muons/MuonRecoAnalyzer/Res_" + type + "_" + parameter;
0598     MonitorElement *residualsHisto = igetter.get(path);
0599 
0600     float mean = -1;
0601     float mean_err = -1;
0602     float sigma = -1;
0603     float sigma_err = -1;
0604 
0605     if (residualsHisto) {
0606       LogTrace(metname) << "[MuonTestSummary]: Starting Gaussian fit for Test of  Res_" << type << "_" << parameter
0607                         << endl;
0608       GaussFit(type, parameter, residualsHisto, mean, mean_err, sigma, sigma_err);
0609 
0610       if (sigma != -1 && parameter == "eta" && type == "TkGlb") {
0611         if (sigma - sigma_err < resEtaSpread_tkGlb)
0612           residualsSummaryMap->setBinContent(bin, 1, 1);
0613         else
0614           residualsSummaryMap->setBinContent(bin, 1, 0);
0615       }
0616       if (sigma != -1 && parameter == "eta" && (type == "GlbSta" || type == "TkSta")) {
0617         if (sigma - sigma_err < resEtaSpread_glbSta)
0618           residualsSummaryMap->setBinContent(bin, 1, 1);
0619         else
0620           residualsSummaryMap->setBinContent(bin, 1, 0);
0621       }
0622       if (sigma != -1 && parameter == "phi" && type == "TkGlb") {
0623         if (sigma - sigma_err < resPhiSpread_tkGlb)
0624           residualsSummaryMap->setBinContent(bin, 2, 1);
0625         else
0626           residualsSummaryMap->setBinContent(bin, 2, 0);
0627       }
0628       if (sigma != -1 && parameter == "phi" && (type == "GlbSta" || type == "TkSta")) {
0629         if (sigma - sigma_err < resPhiSpread_glbSta)
0630           residualsSummaryMap->setBinContent(bin, 2, 1);
0631         else
0632           residualsSummaryMap->setBinContent(bin, 2, 0);
0633       }
0634       if (sigma != -1 && parameter == "oneOverp" && type == "TkGlb") {
0635         if (sigma - sigma_err < resOneOvPSpread_tkGlb)
0636           residualsSummaryMap->setBinContent(bin, 3, 1);
0637         else
0638           residualsSummaryMap->setBinContent(bin, 3, 0);
0639       }
0640       if (sigma != -1 && parameter == "oneOverp" && (type == "GlbSta" || type == "TkSta")) {
0641         if (sigma - sigma_err < resOneOvPSpread_glbSta)
0642           residualsSummaryMap->setBinContent(bin, 3, 1);
0643         else
0644           residualsSummaryMap->setBinContent(bin, 3, 0);
0645       }
0646     }
0647 
0648     //--GH modified
0649     if (type == "TkSta") {
0650       //look at the pull:
0651       string path = "Muons/MuonRecoAnalyzer/Pull_" + type + "_" + parameter;
0652       MonitorElement *pullHisto = igetter.get(path);
0653 
0654       if (pullHisto) {
0655         LogTrace(metname) << "[MuonTestSummary]: Starting Gaussian fit for Test of  Pull_" << type << "_" << parameter
0656                           << endl;
0657         GaussFit(type, parameter, pullHisto, mean, mean_err, sigma, sigma_err);
0658 
0659         if (sigma != -1 && parameter == "eta") {
0660           if (sigma - sigma_err < pullEtaSpread)
0661             residualsSummaryMap->setBinContent(4, 1, 1);
0662           else
0663             residualsSummaryMap->setBinContent(4, 1, 0);
0664         }
0665         if (sigma != -1 && parameter == "phi") {
0666           if (sigma - sigma_err < pullPhiSpread)
0667             residualsSummaryMap->setBinContent(4, 2, 1);
0668           else
0669             residualsSummaryMap->setBinContent(4, 2, 0);
0670         }
0671         if (sigma != -1 && parameter == "oneOverp") {
0672           if (sigma - sigma_err < pullOneOvPSpread)
0673             residualsSummaryMap->setBinContent(4, 3, 1);
0674           else
0675             residualsSummaryMap->setBinContent(4, 3, 0);
0676         }
0677 
0678       }  //have pull histo
0679     }    //TkSta muons
0680   }
0681 
0682   //this part for Global Muons:
0683   else {
0684     string path = "Muons/MuonRecoAnalyzer/" + type + "_" + parameter;
0685     MonitorElement *residualsHisto = igetter.get(path);
0686 
0687     if (residualsHisto) {
0688       LogTrace(metname) << "[MuonTestSummary]: Test of  Charge Comparison " << type << "_" << parameter << endl;
0689       if ((residualsHisto->getBinContent(3) + residualsHisto->getBinContent(4)) != 0) {
0690         LogTrace(metname) << "Charge comparison TkGlb: "
0691                           << residualsHisto->getBinContent(4) /
0692                                  double(residualsHisto->getBinContent(3) + residualsHisto->getBinContent(4))
0693                           << endl;
0694         if (residualsHisto->getBinContent(4) /
0695                 double(residualsHisto->getBinContent(3) + residualsHisto->getBinContent(4)) <
0696             resChargeLimit_tkGlb)
0697           residualsSummaryMap->setBinContent(1, 4, 1);
0698         else
0699           residualsSummaryMap->setBinContent(1, 4, 0);
0700       }
0701       if ((residualsHisto->getBinContent(1) + residualsHisto->getBinContent(2)) != 0) {
0702         LogTrace(metname) << "charge comparison GlbSta: "
0703                           << residualsHisto->getBinContent(2) /
0704                                  double(residualsHisto->getBinContent(1) + residualsHisto->getBinContent(2))
0705                           << endl;
0706         if (residualsHisto->getBinContent(2) /
0707                 double(residualsHisto->getBinContent(1) + residualsHisto->getBinContent(2)) <
0708             resChargeLimit_glbSta)
0709           residualsSummaryMap->setBinContent(2, 4, 1);
0710         else
0711           residualsSummaryMap->setBinContent(2, 4, 0);
0712       }
0713       if (residualsHisto->getBinContent(5) + residualsHisto->getBinContent(6) != 0) {
0714         LogTrace(metname) << "charge comparison TkSta: "
0715                           << residualsHisto->getBinContent(6) /
0716                                  double(residualsHisto->getBinContent(5) + residualsHisto->getBinContent(6))
0717                           << endl;
0718         if (residualsHisto->getBinContent(6) /
0719                 double(residualsHisto->getBinContent(5) + residualsHisto->getBinContent(6)) <
0720             resChargeLimit_tkSta)
0721           residualsSummaryMap->setBinContent(3, 4, 1);
0722         else
0723           residualsSummaryMap->setBinContent(3, 4, 0);
0724       }
0725     }
0726   }
0727 }
0728 
0729 void MuonTestSummary::doMuonIDTests(DQMStore::IGetter &igetter) {
0730   vector<string> muType;
0731   muType.push_back("GlobalMuons");
0732   muType.push_back("TrackerMuons");
0733 
0734   for (int i = 0; i <= 1; i++) {
0735     // num matches test
0736     string path = "Muons/MuonIdDQM/" + muType[i] + "/hNumMatches";
0737     MonitorElement *matchesHisto = igetter.get(path);
0738 
0739     if (matchesHisto) {
0740       TH1F *matchesHisto_root = matchesHisto->getTH1F();
0741       if (matchesHisto_root->GetMaximumBin() >= numMatchedExpected_min &&
0742           matchesHisto_root->GetMaximumBin() <= numMatchedExpected_max)
0743         muonIdSummaryMap->setBinContent(i + 1, 1, 1);
0744       else
0745         muonIdSummaryMap->setBinContent(i + 1, 1, 0);
0746     }
0747 
0748     // num of 0 associated segments
0749     double numOneSegm_dt = 0;
0750     int numHistos_dt = 0;
0751     int numHistos_csc = 0;
0752     MonitorElement *DT1Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hDT1NumSegments");
0753     if (DT1Histo) {
0754       numHistos_dt++;
0755       if (DT1Histo->getEntries() != 0)
0756         numOneSegm_dt += double(DT1Histo->getBinContent(2)) / double(DT1Histo->getEntries());
0757     }
0758     MonitorElement *DT2Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hDT2NumSegments");
0759     if (DT2Histo) {
0760       numHistos_dt++;
0761       if (DT2Histo->getEntries() != 0)
0762         numOneSegm_dt += double(DT2Histo->getBinContent(2)) / double(DT2Histo->getEntries());
0763     }
0764     MonitorElement *DT3Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hDT3NumSegments");
0765     if (DT3Histo) {
0766       numHistos_dt++;
0767       if (DT3Histo->getEntries() != 0)
0768         numOneSegm_dt += double(DT3Histo->getBinContent(2)) / double(DT3Histo->getEntries());
0769     }
0770     MonitorElement *DT4Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hDT4NumSegments");
0771     if (DT4Histo) {
0772       numHistos_dt++;
0773       if (DT4Histo->getEntries() != 0)
0774         numOneSegm_dt += double(DT4Histo->getBinContent(2)) / double(DT4Histo->getEntries());
0775     }
0776     double fraction_dt = 0;
0777     if (numOneSegm_dt != 0) {
0778       fraction_dt = numOneSegm_dt / double(numHistos_dt);
0779       LogTrace(metname) << "fraction_dt: " << fraction_dt << " for " << muType[i] << endl;
0780     }
0781 
0782     double numOneSegm_csc = 0;
0783     MonitorElement *CSC1Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hCSC1NumSegments");
0784     if (CSC1Histo) {
0785       numHistos_csc++;
0786       if (CSC1Histo->getEntries() != 0)
0787         numOneSegm_csc += double(CSC1Histo->getBinContent(2)) / double(CSC1Histo->getEntries());
0788     }
0789     MonitorElement *CSC2Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hCSC2NumSegments");
0790     if (CSC2Histo) {
0791       numHistos_csc++;
0792       if (CSC2Histo->getEntries() != 0)
0793         numOneSegm_csc += double(CSC2Histo->getBinContent(2)) / double(CSC2Histo->getEntries());
0794     }
0795     MonitorElement *CSC3Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hCSC3NumSegments");
0796     if (CSC3Histo) {
0797       numHistos_csc++;
0798       if (CSC3Histo->getEntries() != 0)
0799         numOneSegm_csc += double(CSC3Histo->getBinContent(2)) / double(CSC3Histo->getEntries());
0800     }
0801     MonitorElement *CSC4Histo = igetter.get("Muons/MuonIdDQM/" + muType[i] + "/hCSC4NumSegments");
0802     if (CSC4Histo) {
0803       numHistos_csc++;
0804       if (CSC4Histo->getEntries() != 0)
0805         numOneSegm_csc += double(CSC4Histo->getBinContent(2)) / double(CSC4Histo->getEntries());
0806     }
0807     double fraction_csc = 0;
0808     if (numOneSegm_csc != 0) {
0809       fraction_csc = numOneSegm_csc / double(numHistos_csc);
0810       LogTrace(metname) << "fraction_csc: " << fraction_csc << " for " << muType[i] << endl;
0811     }
0812 
0813     //--GH modified
0814 
0815     if (fraction_dt > matchesFractionDt_min && fraction_dt < matchesFractionDt_max)
0816       muonIdSummaryMap->setBinContent(2 * i + 1, 1, 1);
0817     else
0818       muonIdSummaryMap->setBinContent(2 * i + 1, 1, 0);
0819 
0820     if (fraction_csc > matchesFractionCsc_min && fraction_csc < matchesFractionCsc_max)
0821       muonIdSummaryMap->setBinContent(2 * i + 2, 1, 1);
0822     else
0823       muonIdSummaryMap->setBinContent(2 * i + 2, 1, 0);
0824 
0825     //--GH modified
0826 
0827     // residuals test
0828     vector<string> DTXresHistos, DTYresHistos, CSCXresHistos, CSCYresHistos;
0829     DTXresHistos.push_back("hDT1Pullx");
0830     DTXresHistos.push_back("hDT2Pullx");
0831     DTXresHistos.push_back("hDT3Pullx");
0832     DTXresHistos.push_back("hDT4Pullx");
0833 
0834     DTYresHistos.push_back("hDT1Pully");
0835     DTYresHistos.push_back("hDT2Pully");
0836     DTYresHistos.push_back("hDT3Pully");
0837 
0838     CSCXresHistos.push_back("hCSC1Pullx");
0839     CSCXresHistos.push_back("hCSC2Pullx");
0840     CSCXresHistos.push_back("hCSC3Pullx");
0841     CSCXresHistos.push_back("hCSC4Pullx");
0842 
0843     CSCYresHistos.push_back("hCSC1Pully");
0844     CSCYresHistos.push_back("hCSC2Pully");
0845     CSCYresHistos.push_back("hCSC3Pully");
0846     CSCYresHistos.push_back("hCSC4Pully");
0847 
0848     int numPlot_dtX, numPlot_dtY, numPlot_cscX, numPlot_cscY;
0849     double dtSigmaX, dtSigmaY, cscSigmaX, cscSigmaY;
0850     double dtSigmaX_err, dtSigmaY_err, cscSigmaX_err, cscSigmaY_err;
0851     double dtMeanX, dtMeanY, cscMeanX, cscMeanY;
0852     double dtMeanX_err, dtMeanY_err, cscMeanX_err, cscMeanY_err;
0853     MuonTestSummary::ResidualCheck(
0854         igetter, muType[i], DTXresHistos, numPlot_dtX, dtMeanX, dtMeanX_err, dtSigmaX, dtSigmaX_err);
0855     MuonTestSummary::ResidualCheck(
0856         igetter, muType[i], DTYresHistos, numPlot_dtY, dtMeanY, dtMeanY_err, dtSigmaY, dtSigmaY_err);
0857     MuonTestSummary::ResidualCheck(
0858         igetter, muType[i], CSCXresHistos, numPlot_cscX, cscMeanX, cscMeanX_err, cscSigmaX, cscSigmaX_err);
0859     MuonTestSummary::ResidualCheck(
0860         igetter, muType[i], CSCYresHistos, numPlot_cscY, cscMeanY, cscMeanY_err, cscSigmaY, cscSigmaY_err);
0861 
0862     LogTrace(metname) << "DT mean must be between: " << resSegmTrack_mean_min << " and " << resSegmTrack_mean_max
0863                       << endl;
0864     LogTrace(metname) << "DT rms must be between: " << resSegmTrack_rms_min << " and " << resSegmTrack_rms_max << endl;
0865     LogTrace(metname) << "DT X residual " << muType[i] << " mean: " << dtMeanX << " +- " << dtMeanX_err
0866                       << ", sigma: " << dtSigmaX << " +- " << dtSigmaX_err << endl;
0867     LogTrace(metname) << "DT Y residual " << muType[i] << " mean: " << dtMeanY << " +- " << dtMeanY_err
0868                       << ", sigma: " << dtSigmaY << " +- " << dtSigmaY_err << endl;
0869     LogTrace(metname) << "CSC X residual " << muType[i] << " mean: " << cscMeanX << " +- " << cscMeanX_err
0870                       << ", sigma: " << cscSigmaX << " +- " << cscSigmaX_err << endl;
0871     LogTrace(metname) << "CSC Y residual " << muType[i] << " mean: " << cscMeanY << " +- " << cscMeanY_err
0872                       << ", sigma: " << cscSigmaY << " +- " << cscSigmaY_err << endl;
0873 
0874     //require the mean and rms to be within nsig sigma of preferred range;
0875     const int nsig = 2;
0876     if (numPlot_dtX > 0) {
0877       if (dtMeanX + nsig * dtMeanX_err > resSegmTrack_mean_min && dtMeanX - nsig * dtMeanX_err < resSegmTrack_mean_max)
0878         muonIdSummaryMap->setBinContent(2 * i + 1, 2, 1);
0879       else
0880         muonIdSummaryMap->setBinContent(2 * i + 1, 2, 0);
0881 
0882       if (dtSigmaX + nsig * dtSigmaX_err > resSegmTrack_rms_min &&
0883           dtSigmaX - nsig * dtSigmaX_err < resSegmTrack_rms_max)
0884         muonIdSummaryMap->setBinContent(2 * i + 1, 3, 1);
0885       else
0886         muonIdSummaryMap->setBinContent(2 * i + 1, 3, 0);
0887     }
0888     if (numPlot_dtY > 0) {
0889       if (dtMeanY + nsig * dtMeanY_err > resSegmTrack_mean_min && dtMeanY - nsig * dtMeanY_err < resSegmTrack_mean_max)
0890         muonIdSummaryMap->setBinContent(2 * i + 1, 4, 1);
0891       else
0892         muonIdSummaryMap->setBinContent(2 * i + 1, 4, 0);
0893 
0894       if (dtSigmaY + nsig * dtSigmaY_err > resSegmTrack_rms_min &&
0895           dtSigmaY - nsig * dtSigmaX_err < resSegmTrack_rms_max)
0896         muonIdSummaryMap->setBinContent(2 * i + 1, 5, 1);
0897       else
0898         muonIdSummaryMap->setBinContent(2 * i + 1, 5, 0);
0899     }
0900 
0901     if (numPlot_cscX > 0) {
0902       if (cscMeanX + nsig * cscMeanX_err > resSegmTrack_mean_min &&
0903           cscMeanX - nsig * cscMeanX_err < resSegmTrack_mean_max)
0904         muonIdSummaryMap->setBinContent(2 * i + 2, 2, 1);
0905       else
0906         muonIdSummaryMap->setBinContent(2 * i + 2, 2, 0);
0907 
0908       if (cscSigmaX + nsig * cscSigmaX_err > resSegmTrack_rms_min &&
0909           cscSigmaX - nsig * cscSigmaX_err < resSegmTrack_rms_max)
0910         muonIdSummaryMap->setBinContent(2 * i + 2, 3, 1);
0911       else
0912         muonIdSummaryMap->setBinContent(2 * i + 2, 3, 0);
0913     }
0914     if (numPlot_cscY > 0) {
0915       if (cscMeanY + nsig * cscMeanY_err > resSegmTrack_mean_min &&
0916           cscMeanY - nsig * cscMeanY_err < resSegmTrack_mean_max)
0917         muonIdSummaryMap->setBinContent(2 * i + 2, 4, 1);
0918       else
0919         muonIdSummaryMap->setBinContent(2 * i + 2, 4, 0);
0920 
0921       if (cscSigmaY + nsig * cscSigmaY_err > resSegmTrack_rms_min &&
0922           cscSigmaY - nsig * cscSigmaY_err < resSegmTrack_rms_max)
0923         muonIdSummaryMap->setBinContent(2 * i + 2, 5, 1);
0924       else
0925         muonIdSummaryMap->setBinContent(2 * i + 2, 5, 0);
0926     }
0927 
0928     //---- end of modification
0929   }
0930 }
0931 
0932 void MuonTestSummary::ResidualCheck(DQMStore::IGetter &igetter,
0933                                     std::string muType,
0934                                     const std::vector<std::string> &resHistos,
0935                                     int &numPlot,
0936                                     double &Mean,
0937                                     double &Mean_err,
0938                                     double &Sigma,
0939                                     double &Sigma_err) {
0940   numPlot = 0;
0941   Mean = 0;
0942   Mean_err = 0;
0943   Sigma = 0;
0944   Sigma_err = 0;
0945   for (uint name = 0; name < resHistos.size(); name++) {
0946     MonitorElement *resHisto = igetter.get("Muons/MuonIdDQM/" + muType + "/" + resHistos[name]);
0947 
0948     if (resHisto) {
0949       TH1F *resHisto_root = resHisto->getTH1F();
0950       if (resHisto_root->GetEntries() < 20) {
0951         LogTrace(metname) << "[MuonTestSummary]: Test of " << muType << " for " << resHistos[name]
0952                           << " not performed because # entries < 20 ";
0953         continue;
0954       }
0955 
0956       //we also want to check if the peak is away from zero.
0957       //so, try fitting in 3 sigma around the histogram mean.
0958       //alternatively, could use the maximum bin (less sensitive to 1-sided tails).
0959       //  float mean = resHisto_root->GetMean();
0960       float mean = resHisto_root->GetBinLowEdge(resHisto_root->GetMaximumBin());
0961       TF1 *gfit = new TF1("Gaussian", "gaus", mean - 3, mean + 3);
0962 
0963       try {
0964         resHisto_root->Fit(gfit, "Q0");
0965       } catch (cms::Exception &iException) {
0966         edm::LogError(metname) << "[MuonTestSummary]: Exception when fitting " << resHistos[name];
0967         continue;
0968       }
0969       if (gfit) {
0970         double mean = gfit->GetParameter(1);
0971         double mean_err = gfit->GetParError(1);
0972         double sigma = gfit->GetParameter(2);
0973         double sigma_err = gfit->GetParError(2);
0974         LogTrace(metname) << "meanRes: " << mean << " +- " << mean_err << " for " << resHistos[name] << endl;
0975         LogTrace(metname) << "sigmaRes: " << sigma << " +- " << sigma_err << " for " << resHistos[name] << endl;
0976 
0977         Mean += mean;
0978         Mean_err += mean_err * mean_err;
0979         Sigma += sigma;
0980         Sigma_err += sigma_err * sigma_err;
0981         numPlot++;
0982       }  //if gfit? why would we not have gfit?
0983 
0984     }  //histogram exists...
0985   }    // loop over residuals histos
0986 
0987   if (numPlot == 0) {  //eg no stats
0988     Mean_err = 1;
0989     Mean = 1;
0990     Sigma_err = 1;
0991     Sigma = 1;
0992   } else {
0993     Mean_err = sqrt(Mean_err);
0994     Mean_err /= numPlot;
0995     Mean /= numPlot;
0996 
0997     Sigma_err = sqrt(Sigma_err);
0998     Sigma_err /= numPlot;
0999     Sigma /= numPlot;
1000   }
1001   return;
1002 }
1003 void MuonTestSummary::doEnergyTests(DQMStore::IGetter &igetter, string histname, string muonType, int binNumber) {
1004   // num matches test
1005   string path = "Muons/MuonEnergyDepositAnalyzer/" + histname + muonType;
1006   MonitorElement *energyHisto = igetter.get(path);
1007   Double_t hPeak = -1, hFWHM = -1;
1008   if (energyHisto) {
1009     TH1F *energyHisto_root = energyHisto->getTH1F();
1010 
1011     // Setting fit range and start values
1012     Double_t fitRange[2];
1013     Double_t startValues[4], parlimitslo[4], parlimitshi[4], fitPar[4], fitParErr[4];
1014 
1015     if (energyHisto->getEntries() > 20) {
1016       if (histname == "ecalS9PointingMuDepositedEnergy_") {
1017         fitRange[0] = 0.04;
1018         fitRange[1] = 3.0;
1019 
1020         startValues[0] = 0.036;
1021         startValues[1] = 0.193;
1022         startValues[2] = 110.0;
1023         startValues[3] = 0.06;
1024         parlimitslo[0] = 0.0;
1025         parlimitslo[1] = 0.;
1026         parlimitslo[2] = 1.0;
1027         parlimitslo[3] = 0.;
1028         parlimitshi[0] = 0.05;
1029         parlimitshi[1] = 0.5;
1030         parlimitshi[2] = 80000.0;
1031         parlimitshi[3] = 0.1;
1032 
1033         Double_t chisqr;
1034         Int_t ndf;
1035         TF1 *fit = langaufit(
1036             energyHisto_root, fitRange, startValues, parlimitslo, parlimitshi, fitPar, fitParErr, &chisqr, &ndf);
1037         if (fit) {
1038           langaupro(fitPar, hPeak, hFWHM);
1039           LogTrace(metname) << "hPeak from langau fit: " << hPeak << " for: " << histname + muonType << endl;
1040           LogTrace(metname) << "hFWHM from langau fit: " << hFWHM << " for: " << histname + muonType << endl;
1041         }
1042       }
1043 
1044       if (histname == "hadS9PointingMuDepositedEnergy_") {
1045         fitRange[0] = 0.0;
1046         fitRange[1] = 7.0;
1047 
1048         startValues[0] = 2.0;
1049         startValues[1] = 2.4;
1050         startValues[2] = 110.0;
1051         startValues[3] = 4.0;
1052         parlimitslo[0] = 0.0;
1053         parlimitslo[1] = 0.;
1054         parlimitslo[2] = 1.0;
1055         parlimitslo[3] = 0.;
1056         parlimitshi[0] = 4.0;
1057         parlimitshi[1] = 4.0;
1058         parlimitshi[2] = 80000.0;
1059         parlimitshi[3] = 8.0;
1060 
1061         Double_t chisqr;
1062         Int_t ndf;
1063         TF1 *fit = langaufit(
1064             energyHisto_root, fitRange, startValues, parlimitslo, parlimitshi, fitPar, fitParErr, &chisqr, &ndf);
1065         if (fit) {
1066           langaupro(fitPar, hPeak, hFWHM);
1067           LogTrace(metname) << "hPeak from langau fit: " << hPeak << " for: " << histname + muonType << endl;
1068           LogTrace(metname) << "hFWHM from langau fit: " << hFWHM << " for: " << histname + muonType << endl;
1069         }
1070       }
1071     } else {
1072       LogTrace(metname) << "[MuonTestSummary]: Test of  Energy for " << histname + muonType
1073                         << " not performed because # entries < 20 ";
1074     }
1075   }
1076 
1077   if (histname == "ecalS9PointingMuDepositedEnergy_" && hPeak > expPeakEcalS9_min && hPeak < expPeakEcalS9_max)
1078     energySummaryMap->setBinContent(binNumber, 1, 1);
1079   if (histname == "ecalS9PointingMuDepositedEnergy_" && (hPeak != -1) &&
1080       !(hPeak > expPeakEcalS9_min && hPeak < expPeakEcalS9_max))
1081     energySummaryMap->setBinContent(binNumber, 1, 0);
1082 
1083   if (histname == "hadS9PointingMuDepositedEnergy_" && hPeak > expPeakHadS9_min && hPeak < expPeakHadS9_max)
1084     energySummaryMap->setBinContent(binNumber, 2, 1);
1085   if (histname == "hadS9PointingMuDepositedEnergy_" && (hPeak != -1) &&
1086       !(hPeak > expPeakHadS9_min && hPeak < expPeakHadS9_max))
1087     energySummaryMap->setBinContent(binNumber, 2, 0);
1088 
1089   //missing test on ho distributions
1090 }
1091 void MuonTestSummary::doMultiplicityTests(DQMStore::IGetter &igetter) {
1092   MonitorElement *multiplicityHisto = igetter.get("Muons/MuonRecoAnalyzer/muReco");
1093 
1094   if (multiplicityHisto) {
1095     if (multiplicityHisto->getEntries() > 20) {
1096       double multiplicity_GLB = double(multiplicityHisto->getBinContent(1) + multiplicityHisto->getBinContent(2)) /
1097                                 double(multiplicityHisto->getEntries());
1098       LogTrace(metname) << "multiplicity_GLB: " << multiplicity_GLB << " ExpMultiplicityGlb_min "
1099                         << expMultiplicityGlb_min << " ExpMultiplicityGlb_max " << expMultiplicityGlb_max << endl;
1100       double multiplicity_TK = double(multiplicityHisto->getBinContent(3) + multiplicityHisto->getBinContent(4)) /
1101                                double(multiplicityHisto->getEntries());
1102       LogTrace(metname) << "multiplicity_TK: " << multiplicity_TK << " ExpMultiplicityTk_min " << expMultiplicityTk_min
1103                         << " ExpMultiplicityTk_max " << expMultiplicityTk_max << endl;
1104       double multiplicity_STA = double(multiplicityHisto->getBinContent(3) + multiplicityHisto->getBinContent(5)) /
1105                                 double(multiplicityHisto->getEntries());
1106       LogTrace(metname) << "multiplicity_STA: " << multiplicity_STA << " ExpMultiplicitySta_min "
1107                         << expMultiplicitySta_min << " ExpMultiplicitySta_max " << expMultiplicitySta_max << endl;
1108 
1109       if (multiplicity_GLB > expMultiplicityGlb_min && multiplicity_GLB < expMultiplicityGlb_max)
1110         multiplicitySummaryMap->setBinContent(1, 1);
1111       else
1112         multiplicitySummaryMap->setBinContent(1, 0);
1113 
1114       if (multiplicity_TK > expMultiplicityTk_min && multiplicity_TK < expMultiplicityTk_max)
1115         multiplicitySummaryMap->setBinContent(2, 1);
1116       else
1117         multiplicitySummaryMap->setBinContent(2, 0);
1118 
1119       if (multiplicity_STA > expMultiplicitySta_min && multiplicity_STA < expMultiplicitySta_max)
1120         multiplicitySummaryMap->setBinContent(3, 1);
1121       else
1122         multiplicitySummaryMap->setBinContent(3, 0);
1123     } else {
1124       LogTrace(metname) << "[MuonTestSummary]: Test of  Multiplicity not performed because # entries < 20 ";
1125     }
1126   }
1127 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!