Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-22 07:34:00 

0001 /** \class PFAnalyzer
0002  *
0003  *  DQM ParticleFlow analysis monitoring
0004  *
0005  *  \author J. Roloff - Brown University
0006  *
0007  */
0008 
0009 #include "DQMOffline/ParticleFlow/plugins/PFAnalyzer.h"
0010 
0011 // ***********************************************************
0012 PFAnalyzer::PFAnalyzer(const edm::ParameterSet& pSet) {
0013   m_directory = "ParticleFlow";
0014   parameters_ = pSet.getParameter<edm::ParameterSet>("pfAnalysis");
0015 
0016   thePfCandidateCollection_ = consumes<reco::PFCandidateCollection>(pSet.getParameter<edm::InputTag>("pfCandidates"));
0017   pfJetsToken_ = consumes<reco::PFJetCollection>(pSet.getParameter<edm::InputTag>("pfJetCollection"));
0018 
0019   theTriggerResultsLabel_ = pSet.getParameter<edm::InputTag>("TriggerResultsLabel");
0020   triggerResultsToken_ = consumes<edm::TriggerResults>(edm::InputTag(theTriggerResultsLabel_));
0021   highPtJetExpr_ = pSet.getParameter<edm::InputTag>("TriggerName");
0022 
0023   srcWeights = pSet.getParameter<edm::InputTag>("srcWeights");
0024   weightsToken_ = consumes<edm::ValueMap<float>>(srcWeights);
0025 
0026   m_pfNames = {"allPFC", "neutralHadPFC", "chargedHadPFC", "electronPFC", "muonPFC", "gammaPFC", "hadHFPFC", "emHFPFC"};
0027   vertexTag_ = pSet.getParameter<edm::InputTag>("PVCollection");
0028   vertexToken_ = consumes<std::vector<reco::Vertex>>(edm::InputTag(vertexTag_));
0029 
0030   tok_ew_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
0031 
0032   m_observables = parameters_.getParameter<vstring>("observables");
0033   m_eventObservables = parameters_.getParameter<vstring>("eventObservables");
0034   m_pfInJetObservables = parameters_.getParameter<vstring>("pfInJetObservables");
0035   m_npvBins = parameters_.getParameter<vDouble>("NPVBins");
0036 
0037   // List of cuts applied to PFCs that we want to plot
0038   m_cutList = parameters_.getParameter<vstring>("cutList");
0039   // List of jet cuts that we apply for the case of plotting PFCs in jets
0040   m_jetCutList = parameters_.getParameter<vstring>("jetCutList");
0041 
0042   // Link observable strings to the static functions defined in the header file
0043   // Many of these are quite trivial, but this enables a simple way to include a
0044   // variety of observables on-the-fly.
0045   m_funcMap["pt"] = &getPt;
0046   m_funcMap["energy"] = getEnergy;
0047   m_funcMap["eta"] = getEta;
0048   m_funcMap["phi"] = getPhi;
0049 
0050   m_funcMap["HCalE_depth1"] = getHcalEnergy_depth1;
0051   m_funcMap["HCalE_depth2"] = getHcalEnergy_depth2;
0052   m_funcMap["HCalE_depth3"] = getHcalEnergy_depth3;
0053   m_funcMap["HCalE_depth4"] = getHcalEnergy_depth4;
0054   m_funcMap["HCalE_depth5"] = getHcalEnergy_depth5;
0055   m_funcMap["HCalE_depth6"] = getHcalEnergy_depth6;
0056   m_funcMap["HCalE_depth7"] = getHcalEnergy_depth7;
0057 
0058   m_funcMap["ECal_E"] = getEcalEnergy;
0059   m_funcMap["RawECal_E"] = getRawEcalEnergy;
0060   m_funcMap["HCal_E"] = getHcalEnergy;
0061   m_funcMap["RawHCal_E"] = getRawHcalEnergy;
0062   m_funcMap["HO_E"] = getHOEnergy;
0063   m_funcMap["RawHO_E"] = getRawHOEnergy;
0064   m_funcMap["PFHad_calibration"] = getHadCalibration;
0065 
0066   m_funcMap["MVAIsolated"] = getMVAIsolated;
0067   m_funcMap["MVAEPi"] = getMVAEPi;
0068   m_funcMap["MVAEMu"] = getMVAEMu;
0069   m_funcMap["MVAPiMu"] = getMVAPiMu;
0070   m_funcMap["MVANothingGamma"] = getMVANothingGamma;
0071   m_funcMap["MVANothingNH"] = getMVANothingNH;
0072   m_funcMap["MVAGammaNH"] = getMVAGammaNH;
0073 
0074   m_funcMap["DNNESigIsolated"] = getDNNESigIsolated;
0075   m_funcMap["DNNESigNonIsolated"] = getDNNESigNonIsolated;
0076   m_funcMap["DNNEBkgNonIsolated"] = getDNNEBkgNonIsolated;
0077   m_funcMap["DNNEBkgTauIsolated"] = getDNNEBkgTauIsolated;
0078   m_funcMap["DNNEBkgPhotonIsolated"] = getDNNEBkgPhotonIsolated;
0079 
0080   m_funcMap["hcalE"] = getHCalEnergy;
0081   m_funcMap["eOverP"] = getEoverP;
0082   m_funcMap["nTrkInBlock"] = getNTracksInBlock;
0083 
0084   m_eventFuncMap["NPFC"] = getNPFC;
0085   m_jetWideFuncMap["NPFC"] = getNPFCinJet;
0086 
0087   m_pfInJetFuncMap["PFSpectrum"] = getEnergySpectrum;
0088 
0089   // Link jet observables to static functions in the header file.
0090   // This is very similar to m_funcMap, but for jets instead.
0091   m_jetFuncMap["pt"] = getJetPt;
0092 
0093   // Convert the cutList strings into real cuts that can be applied
0094   // The format should be three comma separated values
0095   // with the first number being the name of the observable
0096   // (corresponding to a key in m_funcMap),
0097   // the second being the minimum value, and the last being the max.
0098   //
0099   for (unsigned int i = 0; i < m_cutList.size(); i++) {
0100     m_fullCutList.push_back(std::vector<std::string>());
0101     while (m_cutList[i].find(']') != std::string::npos) {
0102       size_t pos = m_cutList[i].find(']');
0103       m_fullCutList[i].push_back(m_cutList[i].substr(1, pos));
0104       m_cutList[i].erase(0, pos + 1);
0105     }
0106   }
0107 
0108   for (unsigned int i = 0; i < m_fullCutList.size(); i++) {
0109     m_binList.push_back(std::vector<std::vector<double>>());
0110     for (unsigned int j = 0; j < m_fullCutList[i].size(); j++) {
0111       size_t pos = m_fullCutList[i][j].find(';');
0112       std::string observableName = m_fullCutList[i][j].substr(0, pos);
0113       m_fullCutList[i][j].erase(0, pos + 1);
0114 
0115       m_binList[i].push_back(getBinList(m_fullCutList[i][j]));
0116       m_fullCutList[i][j] = observableName;
0117     }
0118   }
0119 
0120   // Convert the jetCutList strings into real cuts that can be applied
0121   // The format should be three comma separated values,
0122   // with the first number being the name of the observable
0123   // (corresponding to a key in m_jetFuncMap),
0124   // the second being the minimum value, and the last being the max value.
0125   //
0126   for (unsigned int i = 0; i < m_jetCutList.size(); i++) {
0127     m_fullJetCutList.push_back(std::vector<std::string>());
0128     while (m_jetCutList[i].find(']') != std::string::npos) {
0129       size_t pos = m_jetCutList[i].find(']');
0130       m_fullJetCutList[i].push_back(m_jetCutList[i].substr(1, pos));
0131       m_jetCutList[i].erase(0, pos + 1);
0132     }
0133   }
0134 
0135   for (unsigned int i = 0; i < m_fullJetCutList.size(); i++) {
0136     m_jetBinList.push_back(std::vector<std::vector<double>>());
0137     for (unsigned int j = 0; j < m_fullJetCutList[i].size(); j++) {
0138       size_t pos = m_fullJetCutList[i][j].find(';');
0139       std::string observableName = m_fullJetCutList[i][j].substr(0, pos);
0140       m_fullJetCutList[i][j].erase(0, pos + 1);
0141 
0142       m_jetBinList[i].push_back(getBinList(m_fullJetCutList[i][j]));
0143       m_fullJetCutList[i][j] = observableName;
0144     }
0145   }
0146 }
0147 
0148 // ***********************************************************
0149 PFAnalyzer::~PFAnalyzer() { LogTrace("PFAnalyzer") << "[PFAnalyzer] Saving the histos"; }
0150 
0151 // ***********************************************************
0152 void PFAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const&) {
0153   ibooker.setCurrentFolder(m_directory);
0154 
0155   for (unsigned int i = 0; i < m_fullCutList.size(); i++) {
0156     m_allSuffixes.push_back(getAllSuffixes(m_fullCutList[i], m_binList[i]));
0157   }
0158 
0159   for (unsigned int i = 0; i < m_fullJetCutList.size(); i++) {
0160     m_allJetSuffixes.push_back(getAllSuffixes(m_fullJetCutList[i], m_jetBinList[i]));
0161   }
0162 
0163   for (unsigned int npv = 0; npv < m_npvBins.size() - 1; npv++) {
0164     std::string npvString = Form("npv_%.0f_%.0f", m_npvBins[npv], m_npvBins[npv + 1]);
0165     // TODO: Make it possible to use an arbitrary list of bins instead of evenly space bins?
0166     // It is not clear if this is straightforward to do with these classes and CMSSW.
0167     // If it is, it should be an easy change to the code.
0168     //
0169     //
0170     // Books a histogram for each histogram in the config file.
0171     // The format for the observables should be four comma separated values,
0172     // with the first being the observable name (corresponding to one of
0173     // the keys in m_funcMap), the second being the number of bins,
0174     // and the last two being the min and max value for the histogram respectively.
0175     for (unsigned int i = 0; i < m_observables.size(); i++) {
0176       std::string cObservable = m_observables[i];
0177       PFAnalyzer::binInfo obsInfo = getBinInfo(cObservable);
0178 
0179       if (npv == 0)
0180         m_observableNames.push_back(obsInfo.observable);
0181 
0182       for (unsigned int j = 0; j < m_allSuffixes.size(); j++) {
0183         for (unsigned int n = 0; n < m_allSuffixes[j].size(); n++) {
0184           // Loop over all of the different types of PF candidates
0185           for (unsigned int m = 0; m < m_pfNames.size(); m++) {
0186             // For each observable, we make a couple histograms based on a few generic categorizations.
0187             // In all cases, the PFCs that go into these histograms must pass the PFC selection from m_cutList.
0188             std::string histName = Form("%s_%s%s_%s",
0189                                         m_pfNames[m].c_str(),
0190                                         obsInfo.observable.c_str(),
0191                                         m_allSuffixes[j][n].c_str(),
0192                                         npvString.c_str());
0193             MonitorElement* mHist = ibooker.book1D(
0194                 histName, Form(";%s;", obsInfo.axisName.c_str()), obsInfo.nBins, obsInfo.binMin, obsInfo.binMax);
0195             map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHist));
0196           }
0197 
0198           for (unsigned int k = 0; k < m_allJetSuffixes.size(); k++) {
0199             for (unsigned int p = 0; p < m_allJetSuffixes[k].size(); p++) {
0200               for (unsigned int m = 0; m < m_pfNames.size(); m++) {
0201                 // These histograms are for PFCs passing the basic selection, and which are matched to jets
0202                 // that pass the jet selection
0203                 std::string histName = Form("%s_jetMatched_%s%s_jetCuts%s_%s",
0204                                             m_pfNames[m].c_str(),
0205                                             obsInfo.observable.c_str(),
0206                                             m_allSuffixes[j][n].c_str(),
0207                                             m_allJetSuffixes[k][p].c_str(),
0208                                             npvString.c_str());
0209                 MonitorElement* mHistInJet = ibooker.book1D(
0210                     histName, Form(";%s;", obsInfo.axisName.c_str()), obsInfo.nBins, obsInfo.binMin, obsInfo.binMax);
0211                 map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHistInJet));
0212               }
0213             }
0214           }
0215         }
0216       }
0217     }
0218 
0219     // Do the same for global observables (things like the number of PFCs in an event, or in a jet, etc)
0220     for (unsigned int i = 0; i < m_eventObservables.size(); i++) {
0221       std::string cEventObservable = m_eventObservables[i];
0222       size_t pos = cEventObservable.find(';');
0223       std::string observableName = cEventObservable.substr(0, pos);
0224       cEventObservable.erase(0, pos + 1);
0225 
0226       pos = cEventObservable.find(';');
0227       std::string axisString = cEventObservable.substr(0, pos);
0228       cEventObservable.erase(0, pos + 1);
0229 
0230       pos = cEventObservable.find(';');
0231       int nBins = atoi(cEventObservable.substr(0, pos).c_str());
0232       cEventObservable.erase(0, pos + 1);
0233 
0234       pos = cEventObservable.find(';');
0235       float binMin = atof(cEventObservable.substr(0, pos).c_str());
0236       cEventObservable.erase(0, pos + 1);
0237 
0238       pos = cEventObservable.find(';');
0239       float binMax = atof(cEventObservable.substr(0, pos).c_str());
0240       cEventObservable.erase(0, pos + 1);
0241 
0242       pos = cEventObservable.find(';');
0243       int nBinsJet = atoi(cEventObservable.substr(0, pos).c_str());
0244       cEventObservable.erase(0, pos + 1);
0245 
0246       pos = cEventObservable.find(';');
0247       float binMinJet = atof(cEventObservable.substr(0, pos).c_str());
0248       cEventObservable.erase(0, pos + 1);
0249 
0250       float binMaxJet = atof(cEventObservable.c_str());
0251       if (npv == 0)
0252         m_eventObservableNames.push_back(observableName);
0253 
0254       for (unsigned int m = 0; m < m_pfNames.size(); m++) {
0255         std::string histName = Form("%s_%s_%s", m_pfNames[m].c_str(), observableName.c_str(), npvString.c_str());
0256         MonitorElement* mHist = ibooker.book1D(histName, Form(";%s;", axisString.c_str()), nBins, binMin, binMax);
0257         map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHist));
0258       }
0259 
0260       for (unsigned int k = 0; k < m_allJetSuffixes.size(); k++) {
0261         for (unsigned int p = 0; p < m_allJetSuffixes[k].size(); p++) {
0262           for (unsigned int m = 0; m < m_pfNames.size(); m++) {
0263             // These histograms are for PFCs passing the basic selection, and which are matched to jets
0264             // that pass the jet selection
0265             std::string histName = Form("%s_jetMatched_%s_jetCuts%s_%s",
0266                                         m_pfNames[m].c_str(),
0267                                         observableName.c_str(),
0268                                         m_allJetSuffixes[k][p].c_str(),
0269                                         npvString.c_str());
0270             MonitorElement* mHistInJet =
0271                 ibooker.book1D(histName, Form(";%s;", axisString.c_str()), nBinsJet, binMinJet, binMaxJet);
0272             map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHistInJet));
0273           }
0274         }
0275       }
0276     }
0277 
0278     for (unsigned int i = 0; i < m_pfInJetObservables.size(); i++) {
0279       std::string cPfInJetObservable = m_pfInJetObservables[i];
0280       PFAnalyzer::binInfo pfInJetInfo = getBinInfo(cPfInJetObservable);
0281       if (npv == 0)
0282         m_pfInJetObservableNames.push_back(pfInJetInfo.observable);
0283 
0284       for (unsigned int j = 0; j < m_allSuffixes.size(); j++) {
0285         for (unsigned int n = 0; n < m_allSuffixes[j].size(); n++) {
0286           for (unsigned int k = 0; k < m_allJetSuffixes.size(); k++) {
0287             for (unsigned int p = 0; p < m_allJetSuffixes[k].size(); p++) {
0288               for (unsigned int m = 0; m < m_pfNames.size(); m++) {
0289                 // These histograms are for PFCs passing the basic selection, and which are matched to jets
0290                 // that pass the jet selection
0291                 std::string histName = Form("%s_jetMatched_%s%s_jetCuts%s_%s",
0292                                             m_pfNames[m].c_str(),
0293                                             pfInJetInfo.observable.c_str(),
0294                                             m_allSuffixes[j][n].c_str(),
0295                                             m_allJetSuffixes[k][p].c_str(),
0296                                             npvString.c_str());
0297                 MonitorElement* mHistInJet = ibooker.book1D(histName,
0298                                                             Form(";%s;", pfInJetInfo.axisName.c_str()),
0299                                                             pfInJetInfo.nBins,
0300                                                             pfInJetInfo.binMin,
0301                                                             pfInJetInfo.binMax);
0302                 map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHistInJet));
0303               }
0304             }
0305           }
0306         }
0307       }
0308     }
0309 
0310     // Extra histograms for basic validation of the selection etc.
0311     std::string histName = Form("jetPt_%s", npvString.c_str());
0312     MonitorElement* mHist = ibooker.book1D(histName, Form(";%s;", "p_{T,jet}"), 2000, 0, 2000);
0313     map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHist));
0314 
0315     histName = Form("jetPtLead_%s", npvString.c_str());
0316     mHist = ibooker.book1D(histName, Form(";%s;", "p_{T, leading jet}"), 2000, 0, 2000);
0317     map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHist));
0318 
0319     histName = Form("jetEta_%s", npvString.c_str());
0320     mHist = ibooker.book1D(histName, Form(";%s;", "#eta_{jet}"), 200, -5, 5);
0321     map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHist));
0322 
0323     histName = Form("jetEtaLead_%s", npvString.c_str());
0324     mHist = ibooker.book1D(histName, Form(";%s;", "#eta_{leading jet}"), 200, -5, 5);
0325     map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHist));
0326   }
0327 
0328   std::string histName = Form("NPV");
0329   MonitorElement* mHist = ibooker.book1D(histName, Form(";%s;", "N_PV"), 100, 0, 100);
0330   map_of_MEs.insert(std::pair<std::string, MonitorElement*>(m_directory + "/" + histName, mHist));
0331 }
0332 
0333 PFAnalyzer::binInfo PFAnalyzer::getBinInfo(std::string observableString) {
0334   PFAnalyzer::binInfo binningDetails;
0335 
0336   size_t pos = observableString.find(';');
0337   binningDetails.observable = observableString.substr(0, pos);
0338   observableString.erase(0, pos + 1);
0339 
0340   std::vector<double> binList = getBinList(observableString);
0341   pos = observableString.find(';');
0342   binningDetails.axisName = observableString.substr(0, pos);
0343   observableString.erase(0, pos + 1);
0344 
0345   pos = observableString.find(';');
0346   binningDetails.nBins = atoi(observableString.substr(0, pos).c_str());
0347   observableString.erase(0, pos + 1);
0348 
0349   pos = observableString.find(';');
0350   binningDetails.binMin = atof(observableString.substr(0, pos).c_str());
0351   observableString.erase(0, pos + 1);
0352 
0353   binningDetails.binMax = atof(observableString.c_str());
0354 
0355   return binningDetails;
0356 }
0357 
0358 void PFAnalyzer::bookMESetSelection(std::string DirName, DQMStore::IBooker& ibooker) {
0359   ibooker.setCurrentFolder(DirName);
0360 }
0361 
0362 // ***********************************************************
0363 void PFAnalyzer::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0364 
0365 bool PFAnalyzer::passesEventSelection(const edm::Event& iEvent) { return true; }
0366 
0367 // How many significant digits do we need to save for the values to be distinct?
0368 std::string PFAnalyzer::stringWithDecimals(int bin, std::vector<double> bins) {
0369   double diff = bins[bin + 1] - bins[bin];
0370   double sigFigs = log10(diff);
0371 
0372   // We only want to save as many significant digits as we need to.
0373   // Currently, we might lose some information, so we should think about
0374   // if we want to specify more digits
0375   if (sigFigs >= 1) {
0376     return Form("%.0f_%.0f", bins[bin], bins[bin + 1]);
0377   }
0378 
0379   int nDecimals = int(-1 * sigFigs) + 1;
0380   // We do not want to use decimals since these can mess up histogram retrieval in some cases.
0381   // Instead, we use a 'p' to indicate the decimal.
0382   double newDigit = abs((bins[bin] - int(bins[bin])) * pow(10, nDecimals));
0383   double newDigit2 = (bins[bin + 1] - int(bins[bin + 1])) * pow(10, nDecimals);
0384   std::string signStringLow = "";
0385   std::string signStringHigh = "";
0386   if (bins[bin] < 0)
0387     signStringLow = "m";
0388   if (bins[bin + 1] < 0)
0389     signStringHigh = "m";
0390   return Form("%s%.0fp%.0f_%s%.0fp%.0f",
0391               signStringLow.c_str(),
0392               abs(bins[bin]),
0393               newDigit,
0394               signStringHigh.c_str(),
0395               abs(bins[bin + 1]),
0396               newDigit2);
0397 }
0398 
0399 std::vector<double> PFAnalyzer::getBinList(std::string binString) {
0400   std::vector<double> binList;
0401 
0402   while (binString.find(';') != std::string::npos) {
0403     size_t pos = binString.find(';');
0404     binList.push_back(atof(binString.substr(0, pos).c_str()));
0405     binString.erase(0, pos + 1);
0406   }
0407   binList.push_back(atof(binString.c_str()));
0408 
0409   if (binList.size() == 3) {
0410     int nBins = int(binList[0]);
0411     double minVal = binList[1];
0412     double maxVal = binList[2];
0413     binList.clear();
0414 
0415     for (int i = 0; i <= nBins; i++) {
0416       binList.push_back(minVal + i * (maxVal - minVal) / nBins);
0417     }
0418   }
0419 
0420   return binList;
0421 }
0422 
0423 std::vector<std::string> PFAnalyzer::getAllSuffixes(std::vector<std::string> observables,
0424                                                     std::vector<std::vector<double>> binnings) {
0425   int nTotalBins = 1;
0426   std::vector<int> nBins;
0427   for (unsigned int i = 0; i < binnings.size(); i++) {
0428     nTotalBins = (binnings[i].size() - 1) * nTotalBins;
0429     nBins.push_back(binnings[i].size() - 1);
0430   }
0431 
0432   std::vector<std::vector<int>> binList;
0433 
0434   binList.reserve(nTotalBins);
0435   for (int i = 0; i < nTotalBins; i++) {
0436     binList.push_back(std::vector<int>());
0437   }
0438 
0439   int factor = nTotalBins;
0440   int otherFactor = 1;
0441   for (unsigned int i = 0; i < binnings.size(); i++) {
0442     factor = factor / nBins[i];
0443 
0444     for (int j = 0; j < nBins[i]; j++) {
0445       for (int k = 0; k < factor; k++) {
0446         for (int m = 0; m < otherFactor; m++) {
0447           binList[m * otherFactor + j * factor + k].push_back(j);
0448         }
0449       }
0450     }
0451     otherFactor = otherFactor * nBins[i];
0452   }
0453 
0454   std::vector<std::string> allSuffixes;
0455   allSuffixes.reserve(nTotalBins);
0456   for (int i = 0; i < nTotalBins; i++) {
0457     allSuffixes.push_back(getSuffix(binList[i], observables, binnings));
0458   }
0459 
0460   return allSuffixes;
0461 }
0462 
0463 // Get a unique string corresponding to the selection cuts
0464 std::string PFAnalyzer::getSuffix(std::vector<int> binList,
0465                                   std::vector<std::string> observables,
0466                                   std::vector<std::vector<double>> binnings) {
0467   std::string suffix = "";
0468   for (unsigned int i = 0; i < binList.size(); i++) {
0469     if (binList[i] < 0)
0470       return "";
0471     std::string digitString = stringWithDecimals(binList[i], binnings[i]);
0472 
0473     suffix = Form("%s_%s_%s", suffix.c_str(), observables[i].c_str(), digitString.c_str());
0474   }
0475 
0476   return suffix;
0477 }
0478 
0479 int PFAnalyzer::getBinNumber(double binVal, std::vector<double> bins) {
0480   if (binVal < bins[0])
0481     return -1;
0482   for (unsigned int i = 0; i < bins.size(); i++) {
0483     if (binVal < bins[i])
0484       return i - 1;
0485   }
0486 
0487   return -1;
0488 }
0489 
0490 int PFAnalyzer::getBinNumbers(std::vector<double> binVal, std::vector<std::vector<double>> bins) {
0491   std::vector<int> cbins;
0492   std::vector<int> nBins;
0493   for (unsigned int i = 0; i < binVal.size(); i++) {
0494     int cbin = getBinNumber(binVal[i], bins[i]);
0495     if (cbin < 0)
0496       return -1;
0497     nBins.push_back(bins[i].size() - 1);
0498     cbins.push_back(cbin);
0499   }
0500 
0501   int bin = 0;
0502   int factor = 1;
0503   for (unsigned int i = 0; i < binVal.size(); i++) {
0504     bin += cbins[i] * factor;
0505     factor = factor * nBins[i];
0506   }
0507 
0508   return bin;
0509 }
0510 
0511 int PFAnalyzer::getPFBin(const reco::PFCandidate pfCand, int i) {
0512   std::vector<double> binVals;
0513   for (unsigned int j = 0; j < m_fullCutList[i].size(); j++) {
0514     binVals.push_back(m_funcMap[m_fullCutList[i][j]](pfCand));
0515   }
0516 
0517   return getBinNumbers(binVals, m_binList[i]);
0518 }
0519 
0520 int PFAnalyzer::getJetBin(const reco::PFJet jetCand, int i) {
0521   std::vector<double> binVals;
0522   for (unsigned int j = 0; j < m_fullJetCutList[i].size(); j++) {
0523     binVals.push_back(m_jetFuncMap[m_fullJetCutList[i][j]](jetCand));
0524   }
0525 
0526   return getBinNumbers(binVals, m_jetBinList[i]);
0527 }
0528 
0529 // ***********************************************************
0530 void PFAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0531   const edm::Handle<GenEventInfoProduct> genEventInfo = iEvent.getHandle(tok_ew_);
0532   double eventWeight = 1;
0533   if (genEventInfo.isValid()) {
0534     eventWeight = genEventInfo->weight();
0535   }
0536 
0537   weights_ = &iEvent.get(weightsToken_);
0538 
0539   // **** Get the TriggerResults container
0540   edm::Handle<edm::TriggerResults> triggerResults;
0541   iEvent.getByToken(triggerResultsToken_, triggerResults);
0542 
0543   // Hack to make it pass the lowest unprescaled HLT?
0544   Int_t JetHiPass = 0;
0545 
0546   if (triggerResults.isValid()) {
0547     const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0548 
0549     const unsigned int nTrig(triggerNames.size());
0550     for (unsigned int i = 0; i < nTrig; ++i) {
0551       if (triggerNames.triggerName(i).find(highPtJetExpr_.label()) != std::string::npos && triggerResults->accept(i)) {
0552         JetHiPass = 1;
0553       }
0554     }
0555   }
0556 
0557   //Vertex information
0558   edm::Handle<reco::VertexCollection> vertexHandle;
0559   iEvent.getByToken(vertexToken_, vertexHandle);
0560 
0561   if (!vertexHandle.isValid()) {
0562     LogDebug("") << "PFAnalyzer: Could not find vertex collection" << std::endl;
0563   }
0564   int numPV = 0;
0565 
0566   if (vertexHandle.isValid()) {
0567     reco::VertexCollection vertex = *(vertexHandle.product());
0568     for (reco::VertexCollection::const_iterator v = vertex.begin(); v != vertex.end(); ++v) {
0569       if (v->isFake())
0570         continue;
0571       if (v->ndof() < 4)
0572         continue;
0573       if (fabs(v->z()) > 24.0)
0574         continue;
0575       ++numPV;
0576     }
0577   }
0578 
0579   int npvBin = getBinNumber(numPV, m_npvBins);
0580   if (npvBin < 0)
0581     return;
0582   std::string npvString = Form("npv_%.0f_%.0f", m_npvBins[npvBin], m_npvBins[npvBin + 1]);
0583 
0584   if (!JetHiPass)
0585     return;
0586 
0587   // Retrieve the PFCs
0588   edm::Handle<reco::PFCandidateCollection> pfCollection;
0589   iEvent.getByToken(thePfCandidateCollection_, pfCollection);
0590   if (!pfCollection.isValid()) {
0591     edm::LogError("PFAnalyzer") << "invalid collection: PF candidate \n";
0592     return;
0593   }
0594 
0595   edm::Handle<reco::PFJetCollection> pfJets;
0596   iEvent.getByToken(pfJetsToken_, pfJets);
0597   if (!pfJets.isValid()) {
0598     edm::LogError("PFAnalyzer") << "invalid collection: PF jets \n";
0599     return;
0600   }
0601 
0602   // Probably we want to define a few different options for how the selection will work
0603   // Currently it is just a dummy function, and we hardcode the other cuts.
0604   if (pfJets->size() < 2)
0605     return;
0606   if (pfJets->at(0).pt() < 450)
0607     return;
0608   if (pfJets->at(0).pt() / pfJets->at(1).pt() > 2)
0609     return;
0610 
0611   if (!passesEventSelection(iEvent))
0612     return;
0613 
0614   for (reco::PFCandidateCollection::const_iterator recoPF = pfCollection->begin(); recoPF != pfCollection->end();
0615        ++recoPF) {
0616     for (unsigned int j = 0; j < m_fullCutList.size(); j++) {
0617       int binNumber = getPFBin(*recoPF, j);
0618       if (binNumber < 0)
0619         continue;
0620       if (binNumber >= int(m_allSuffixes[j].size())) {
0621         continue;
0622       }
0623       std::string binString = m_allSuffixes[j][binNumber];
0624 
0625       // Eventually, we might want the hist name to include the cuts that we are applying,
0626       // so I am keepking it as a separate string for now, even though it is redundant.
0627       // Make plots of all observables
0628       for (unsigned int i = 0; i < m_observables.size(); i++) {
0629         std::string histName = Form("%s%s_%s", m_observableNames[i].c_str(), binString.c_str(), npvString.c_str());
0630         map_of_MEs[m_directory + "/allPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0631 
0632         switch (recoPF->particleId()) {
0633           case reco::PFCandidate::ParticleType::h:
0634             map_of_MEs[m_directory + "/chargedHadPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0635                                                                          eventWeight);
0636             break;
0637           case reco::PFCandidate::ParticleType::h0:
0638             map_of_MEs[m_directory + "/neutralHadPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0639                                                                          eventWeight);
0640             break;
0641           case reco::PFCandidate::ParticleType::e:
0642             map_of_MEs[m_directory + "/electronPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0643                                                                        eventWeight);
0644             break;
0645           case reco::PFCandidate::ParticleType::mu:
0646             map_of_MEs[m_directory + "/muonPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0647                                                                    eventWeight);
0648             break;
0649           case reco::PFCandidate::ParticleType::gamma:
0650             map_of_MEs[m_directory + "/gammaPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0651                                                                     eventWeight);
0652             break;
0653           case reco::PFCandidate::ParticleType::h_HF:
0654             map_of_MEs[m_directory + "/hadHFPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0655                                                                     eventWeight);
0656             break;
0657           case reco::PFCandidate::ParticleType::egamma_HF:
0658             map_of_MEs[m_directory + "/emHFPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0659                                                                    eventWeight);
0660             break;
0661           default:
0662             break;
0663         }
0664       }
0665     }
0666   }
0667 
0668   for (unsigned int i = 0; i < m_eventObservableNames.size(); i++) {
0669     std::string histName = Form("%s_%s", m_eventObservableNames[i].c_str(), npvString.c_str());
0670     map_of_MEs[m_directory + "/allPFC_" + histName]->Fill(
0671         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::X), eventWeight);
0672     map_of_MEs[m_directory + "/chargedHadPFC_" + histName]->Fill(
0673         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::h), eventWeight);
0674     map_of_MEs[m_directory + "/neutralHadPFC_" + histName]->Fill(
0675         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::h0), eventWeight);
0676     map_of_MEs[m_directory + "/electronPFC_" + histName]->Fill(
0677         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::e), eventWeight);
0678     map_of_MEs[m_directory + "/muonPFC_" + histName]->Fill(
0679         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::mu), eventWeight);
0680     map_of_MEs[m_directory + "/gammaPFC_" + histName]->Fill(
0681         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::gamma), eventWeight);
0682     map_of_MEs[m_directory + "/hadHFPFC_" + histName]->Fill(
0683         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::h_HF), eventWeight);
0684     map_of_MEs[m_directory + "/emHFPFC_" + histName]->Fill(
0685         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::egamma_HF),
0686         eventWeight);
0687   }
0688 
0689   // Plots for generic debugging
0690   map_of_MEs[m_directory + "/NPV"]->Fill(numPV, eventWeight);
0691   map_of_MEs[m_directory + Form("/jetPtLead_%s", npvString.c_str())]->Fill(pfJets->begin()->pt(), eventWeight);
0692   map_of_MEs[m_directory + Form("/jetEtaLead_%s", npvString.c_str())]->Fill(pfJets->begin()->eta(), eventWeight);
0693 
0694   // Make plots of all observables, this time for PF candidates within jets
0695   for (reco::PFJetCollection::const_iterator cjet = pfJets->begin(); cjet != pfJets->end(); ++cjet) {
0696     map_of_MEs[m_directory + Form("/jetPt_%s", npvString.c_str())]->Fill(cjet->pt(), eventWeight);
0697     map_of_MEs[m_directory + Form("/jetEta_%s", npvString.c_str())]->Fill(cjet->eta(), eventWeight);
0698 
0699     for (unsigned int k = 0; k < m_fullJetCutList.size(); k++) {
0700       int jetBinNumber = getJetBin(*cjet, k);
0701       if (jetBinNumber < 0)
0702         continue;
0703       std::string jetBinString = m_allJetSuffixes[k][jetBinNumber];
0704 
0705       std::vector<reco::PFCandidatePtr> pfConstits = cjet->getPFConstituents();
0706 
0707       for (const auto& recoPF : pfConstits) {
0708         for (unsigned int j = 0; j < m_fullCutList.size(); j++) {
0709           int binNumber = getPFBin(*recoPF, j);
0710           if (binNumber < 0)
0711             continue;
0712           if (binNumber >= int(m_allSuffixes[j].size())) {
0713             continue;
0714           }
0715           std::string binString = m_allSuffixes[j][binNumber];
0716 
0717           for (unsigned int i = 0; i < m_observableNames.size(); i++) {
0718             std::string histName = Form("%s%s_jetCuts%s_%s",
0719                                         m_observableNames[i].c_str(),
0720                                         binString.c_str(),
0721                                         jetBinString.c_str(),
0722                                         npvString.c_str());
0723             map_of_MEs[m_directory + "/allPFC_jetMatched_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0724                                                                              eventWeight);
0725 
0726             switch (recoPF->particleId()) {
0727               case reco::PFCandidate::ParticleType::h:
0728                 map_of_MEs[m_directory + "/chargedHadPFC_jetMatched_" + histName]->Fill(
0729                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0730                 break;
0731               case reco::PFCandidate::ParticleType::h0:
0732                 map_of_MEs[m_directory + "/neutralHadPFC_jetMatched_" + histName]->Fill(
0733                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0734                 break;
0735               case reco::PFCandidate::ParticleType::e:
0736                 map_of_MEs[m_directory + "/electronPFC_jetMatched_" + histName]->Fill(
0737                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0738                 break;
0739               case reco::PFCandidate::ParticleType::mu:
0740                 map_of_MEs[m_directory + "/muonPFC_jetMatched_" + histName]->Fill(
0741                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0742                 break;
0743               case reco::PFCandidate::ParticleType::gamma:
0744                 map_of_MEs[m_directory + "/gammaPFC_jetMatched_" + histName]->Fill(
0745                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0746                 break;
0747               case reco::PFCandidate::ParticleType::h_HF:
0748                 map_of_MEs[m_directory + "/hadHFPFC_jetMatched_" + histName]->Fill(
0749                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0750                 break;
0751               case reco::PFCandidate::ParticleType::egamma_HF:
0752                 map_of_MEs[m_directory + "/emHFPFC_jetMatched_" + histName]->Fill(
0753                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0754                 break;
0755               default:
0756                 break;
0757             }
0758           }
0759 
0760           for (unsigned int i = 0; i < m_pfInJetObservableNames.size(); i++) {
0761             std::string histName = Form("%s%s_jetCuts%s_%s",
0762                                         m_pfInJetObservableNames[i].c_str(),
0763                                         binString.c_str(),
0764                                         jetBinString.c_str(),
0765                                         npvString.c_str());
0766             map_of_MEs[m_directory + "/allPFC_jetMatched_" + histName]->Fill(
0767                 m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0768 
0769             switch (recoPF->particleId()) {
0770               case reco::PFCandidate::ParticleType::h:
0771                 map_of_MEs[m_directory + "/chargedHadPFC_jetMatched_" + histName]->Fill(
0772                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0773                 break;
0774               case reco::PFCandidate::ParticleType::h0:
0775                 map_of_MEs[m_directory + "/neutralHadPFC_jetMatched_" + histName]->Fill(
0776                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0777                 break;
0778               case reco::PFCandidate::ParticleType::e:
0779                 map_of_MEs[m_directory + "/electronPFC_jetMatched_" + histName]->Fill(
0780                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0781                 break;
0782               case reco::PFCandidate::ParticleType::mu:
0783                 map_of_MEs[m_directory + "/muonPFC_jetMatched_" + histName]->Fill(
0784                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0785                 break;
0786               case reco::PFCandidate::ParticleType::gamma:
0787                 map_of_MEs[m_directory + "/gammaPFC_jetMatched_" + histName]->Fill(
0788                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0789                 break;
0790               case reco::PFCandidate::ParticleType::h_HF:
0791                 map_of_MEs[m_directory + "/hadHFPFC_jetMatched_" + histName]->Fill(
0792                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0793                 break;
0794               case reco::PFCandidate::ParticleType::egamma_HF:
0795                 map_of_MEs[m_directory + "/emHFPFC_jetMatched_" + histName]->Fill(
0796                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0797                 break;
0798               default:
0799                 break;
0800             }
0801           }
0802         }
0803 
0804         for (unsigned int i = 0; i < m_eventObservableNames.size(); i++) {
0805           std::string histName =
0806               Form("%s_jetCuts%s_%s", m_eventObservableNames[i].c_str(), jetBinString.c_str(), npvString.c_str());
0807           map_of_MEs[m_directory + "/allPFC_jetMatched_" + histName]->Fill(
0808               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::X), eventWeight);
0809           map_of_MEs[m_directory + "/chargedHadPFC_jetMatched_" + histName]->Fill(
0810               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::h), eventWeight);
0811           map_of_MEs[m_directory + "/neutralHadPFC_jetMatched_" + histName]->Fill(
0812               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::h0),
0813               eventWeight);
0814           map_of_MEs[m_directory + "/electronPFC_jetMatched_" + histName]->Fill(
0815               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::e), eventWeight);
0816           map_of_MEs[m_directory + "/muonPFC_jetMatched_" + histName]->Fill(
0817               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::mu),
0818               eventWeight);
0819           map_of_MEs[m_directory + "/gammaPFC_jetMatched_" + histName]->Fill(
0820               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::gamma),
0821               eventWeight);
0822           map_of_MEs[m_directory + "/hadHFPFC_jetMatched_" + histName]->Fill(
0823               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::h_HF),
0824               eventWeight);
0825           map_of_MEs[m_directory + "/emHFPFC_jetMatched_" + histName]->Fill(
0826               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::egamma_HF),
0827               eventWeight);
0828         }
0829       }
0830     }
0831   }
0832 }
[ 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!