Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-22 22:37:05 

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   for (int i = 0; i < nTotalBins; i++) {
0435     binList.push_back(std::vector<int>());
0436   }
0437 
0438   int factor = nTotalBins;
0439   int otherFactor = 1;
0440   for (unsigned int i = 0; i < binnings.size(); i++) {
0441     factor = factor / nBins[i];
0442 
0443     for (int j = 0; j < nBins[i]; j++) {
0444       for (int k = 0; k < factor; k++) {
0445         for (int m = 0; m < otherFactor; m++) {
0446           binList[m * otherFactor + j * factor + k].push_back(j);
0447         }
0448       }
0449     }
0450     otherFactor = otherFactor * nBins[i];
0451   }
0452 
0453   std::vector<std::string> allSuffixes;
0454   for (int i = 0; i < nTotalBins; i++) {
0455     allSuffixes.push_back(getSuffix(binList[i], observables, binnings));
0456   }
0457 
0458   return allSuffixes;
0459 }
0460 
0461 // Get a unique string corresponding to the selection cuts
0462 std::string PFAnalyzer::getSuffix(std::vector<int> binList,
0463                                   std::vector<std::string> observables,
0464                                   std::vector<std::vector<double>> binnings) {
0465   std::string suffix = "";
0466   for (unsigned int i = 0; i < binList.size(); i++) {
0467     if (binList[i] < 0)
0468       return "";
0469     std::string digitString = stringWithDecimals(binList[i], binnings[i]);
0470 
0471     suffix = Form("%s_%s_%s", suffix.c_str(), observables[i].c_str(), digitString.c_str());
0472   }
0473 
0474   return suffix;
0475 }
0476 
0477 int PFAnalyzer::getBinNumber(double binVal, std::vector<double> bins) {
0478   if (binVal < bins[0])
0479     return -1;
0480   for (unsigned int i = 0; i < bins.size(); i++) {
0481     if (binVal < bins[i])
0482       return i - 1;
0483   }
0484 
0485   return -1;
0486 }
0487 
0488 int PFAnalyzer::getBinNumbers(std::vector<double> binVal, std::vector<std::vector<double>> bins) {
0489   std::vector<int> cbins;
0490   std::vector<int> nBins;
0491   for (unsigned int i = 0; i < binVal.size(); i++) {
0492     int cbin = getBinNumber(binVal[i], bins[i]);
0493     if (cbin < 0)
0494       return -1;
0495     nBins.push_back(bins[i].size() - 1);
0496     cbins.push_back(cbin);
0497   }
0498 
0499   int bin = 0;
0500   int factor = 1;
0501   for (unsigned int i = 0; i < binVal.size(); i++) {
0502     bin += cbins[i] * factor;
0503     factor = factor * nBins[i];
0504   }
0505 
0506   return bin;
0507 }
0508 
0509 int PFAnalyzer::getPFBin(const reco::PFCandidate pfCand, int i) {
0510   std::vector<double> binVals;
0511   for (unsigned int j = 0; j < m_fullCutList[i].size(); j++) {
0512     binVals.push_back(m_funcMap[m_fullCutList[i][j]](pfCand));
0513   }
0514 
0515   return getBinNumbers(binVals, m_binList[i]);
0516 }
0517 
0518 int PFAnalyzer::getJetBin(const reco::PFJet jetCand, int i) {
0519   std::vector<double> binVals;
0520   for (unsigned int j = 0; j < m_fullJetCutList[i].size(); j++) {
0521     binVals.push_back(m_jetFuncMap[m_fullJetCutList[i][j]](jetCand));
0522   }
0523 
0524   return getBinNumbers(binVals, m_jetBinList[i]);
0525 }
0526 
0527 // ***********************************************************
0528 void PFAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0529   const edm::Handle<GenEventInfoProduct> genEventInfo = iEvent.getHandle(tok_ew_);
0530   double eventWeight = 1;
0531   if (genEventInfo.isValid()) {
0532     eventWeight = genEventInfo->weight();
0533   }
0534 
0535   weights_ = &iEvent.get(weightsToken_);
0536 
0537   // **** Get the TriggerResults container
0538   edm::Handle<edm::TriggerResults> triggerResults;
0539   iEvent.getByToken(triggerResultsToken_, triggerResults);
0540 
0541   // Hack to make it pass the lowest unprescaled HLT?
0542   Int_t JetHiPass = 0;
0543 
0544   if (triggerResults.isValid()) {
0545     const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0546 
0547     const unsigned int nTrig(triggerNames.size());
0548     for (unsigned int i = 0; i < nTrig; ++i) {
0549       if (triggerNames.triggerName(i).find(highPtJetExpr_.label()) != std::string::npos && triggerResults->accept(i)) {
0550         JetHiPass = 1;
0551       }
0552     }
0553   }
0554 
0555   //Vertex information
0556   edm::Handle<reco::VertexCollection> vertexHandle;
0557   iEvent.getByToken(vertexToken_, vertexHandle);
0558 
0559   if (!vertexHandle.isValid()) {
0560     LogDebug("") << "PFAnalyzer: Could not find vertex collection" << std::endl;
0561   }
0562   int numPV = 0;
0563 
0564   if (vertexHandle.isValid()) {
0565     reco::VertexCollection vertex = *(vertexHandle.product());
0566     for (reco::VertexCollection::const_iterator v = vertex.begin(); v != vertex.end(); ++v) {
0567       if (v->isFake())
0568         continue;
0569       if (v->ndof() < 4)
0570         continue;
0571       if (fabs(v->z()) > 24.0)
0572         continue;
0573       ++numPV;
0574     }
0575   }
0576 
0577   int npvBin = getBinNumber(numPV, m_npvBins);
0578   if (npvBin < 0)
0579     return;
0580   std::string npvString = Form("npv_%.0f_%.0f", m_npvBins[npvBin], m_npvBins[npvBin + 1]);
0581 
0582   if (!JetHiPass)
0583     return;
0584 
0585   // Retrieve the PFCs
0586   edm::Handle<reco::PFCandidateCollection> pfCollection;
0587   iEvent.getByToken(thePfCandidateCollection_, pfCollection);
0588   if (!pfCollection.isValid()) {
0589     edm::LogError("PFAnalyzer") << "invalid collection: PF candidate \n";
0590     return;
0591   }
0592 
0593   edm::Handle<reco::PFJetCollection> pfJets;
0594   iEvent.getByToken(pfJetsToken_, pfJets);
0595   if (!pfJets.isValid()) {
0596     edm::LogError("PFAnalyzer") << "invalid collection: PF jets \n";
0597     return;
0598   }
0599 
0600   // Probably we want to define a few different options for how the selection will work
0601   // Currently it is just a dummy function, and we hardcode the other cuts.
0602   if (pfJets->size() < 2)
0603     return;
0604   if (pfJets->at(0).pt() < 450)
0605     return;
0606   if (pfJets->at(0).pt() / pfJets->at(1).pt() > 2)
0607     return;
0608 
0609   if (!passesEventSelection(iEvent))
0610     return;
0611 
0612   for (reco::PFCandidateCollection::const_iterator recoPF = pfCollection->begin(); recoPF != pfCollection->end();
0613        ++recoPF) {
0614     for (unsigned int j = 0; j < m_fullCutList.size(); j++) {
0615       int binNumber = getPFBin(*recoPF, j);
0616       if (binNumber < 0)
0617         continue;
0618       if (binNumber >= int(m_allSuffixes[j].size())) {
0619         continue;
0620       }
0621       std::string binString = m_allSuffixes[j][binNumber];
0622 
0623       // Eventually, we might want the hist name to include the cuts that we are applying,
0624       // so I am keepking it as a separate string for now, even though it is redundant.
0625       // Make plots of all observables
0626       for (unsigned int i = 0; i < m_observables.size(); i++) {
0627         std::string histName = Form("%s%s_%s", m_observableNames[i].c_str(), binString.c_str(), npvString.c_str());
0628         map_of_MEs[m_directory + "/allPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0629 
0630         switch (recoPF->particleId()) {
0631           case reco::PFCandidate::ParticleType::h:
0632             map_of_MEs[m_directory + "/chargedHadPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0633                                                                          eventWeight);
0634             break;
0635           case reco::PFCandidate::ParticleType::h0:
0636             map_of_MEs[m_directory + "/neutralHadPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0637                                                                          eventWeight);
0638             break;
0639           case reco::PFCandidate::ParticleType::e:
0640             map_of_MEs[m_directory + "/electronPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0641                                                                        eventWeight);
0642             break;
0643           case reco::PFCandidate::ParticleType::mu:
0644             map_of_MEs[m_directory + "/muonPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0645                                                                    eventWeight);
0646             break;
0647           case reco::PFCandidate::ParticleType::gamma:
0648             map_of_MEs[m_directory + "/gammaPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0649                                                                     eventWeight);
0650             break;
0651           case reco::PFCandidate::ParticleType::h_HF:
0652             map_of_MEs[m_directory + "/hadHFPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0653                                                                     eventWeight);
0654             break;
0655           case reco::PFCandidate::ParticleType::egamma_HF:
0656             map_of_MEs[m_directory + "/emHFPFC_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0657                                                                    eventWeight);
0658             break;
0659           default:
0660             break;
0661         }
0662       }
0663     }
0664   }
0665 
0666   for (unsigned int i = 0; i < m_eventObservableNames.size(); i++) {
0667     std::string histName = Form("%s_%s", m_eventObservableNames[i].c_str(), npvString.c_str());
0668     map_of_MEs[m_directory + "/allPFC_" + histName]->Fill(
0669         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::X), eventWeight);
0670     map_of_MEs[m_directory + "/chargedHadPFC_" + histName]->Fill(
0671         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::h), eventWeight);
0672     map_of_MEs[m_directory + "/neutralHadPFC_" + histName]->Fill(
0673         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::h0), eventWeight);
0674     map_of_MEs[m_directory + "/electronPFC_" + histName]->Fill(
0675         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::e), eventWeight);
0676     map_of_MEs[m_directory + "/muonPFC_" + histName]->Fill(
0677         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::mu), eventWeight);
0678     map_of_MEs[m_directory + "/gammaPFC_" + histName]->Fill(
0679         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::gamma), eventWeight);
0680     map_of_MEs[m_directory + "/hadHFPFC_" + histName]->Fill(
0681         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::h_HF), eventWeight);
0682     map_of_MEs[m_directory + "/emHFPFC_" + histName]->Fill(
0683         m_eventFuncMap[m_eventObservableNames[i]](*pfCollection, reco::PFCandidate::ParticleType::egamma_HF),
0684         eventWeight);
0685   }
0686 
0687   // Plots for generic debugging
0688   map_of_MEs[m_directory + "/NPV"]->Fill(numPV, eventWeight);
0689   map_of_MEs[m_directory + Form("/jetPtLead_%s", npvString.c_str())]->Fill(pfJets->begin()->pt(), eventWeight);
0690   map_of_MEs[m_directory + Form("/jetEtaLead_%s", npvString.c_str())]->Fill(pfJets->begin()->eta(), eventWeight);
0691 
0692   // Make plots of all observables, this time for PF candidates within jets
0693   for (reco::PFJetCollection::const_iterator cjet = pfJets->begin(); cjet != pfJets->end(); ++cjet) {
0694     map_of_MEs[m_directory + Form("/jetPt_%s", npvString.c_str())]->Fill(cjet->pt(), eventWeight);
0695     map_of_MEs[m_directory + Form("/jetEta_%s", npvString.c_str())]->Fill(cjet->eta(), eventWeight);
0696 
0697     for (unsigned int k = 0; k < m_fullJetCutList.size(); k++) {
0698       int jetBinNumber = getJetBin(*cjet, k);
0699       if (jetBinNumber < 0)
0700         continue;
0701       std::string jetBinString = m_allJetSuffixes[k][jetBinNumber];
0702 
0703       std::vector<reco::PFCandidatePtr> pfConstits = cjet->getPFConstituents();
0704 
0705       for (auto recoPF : pfConstits) {
0706         for (unsigned int j = 0; j < m_fullCutList.size(); j++) {
0707           int binNumber = getPFBin(*recoPF, j);
0708           if (binNumber < 0)
0709             continue;
0710           if (binNumber >= int(m_allSuffixes[j].size())) {
0711             continue;
0712           }
0713           std::string binString = m_allSuffixes[j][binNumber];
0714 
0715           for (unsigned int i = 0; i < m_observableNames.size(); i++) {
0716             std::string histName = Form("%s%s_jetCuts%s_%s",
0717                                         m_observableNames[i].c_str(),
0718                                         binString.c_str(),
0719                                         jetBinString.c_str(),
0720                                         npvString.c_str());
0721             map_of_MEs[m_directory + "/allPFC_jetMatched_" + histName]->Fill(m_funcMap[m_observableNames[i]](*recoPF),
0722                                                                              eventWeight);
0723 
0724             switch (recoPF->particleId()) {
0725               case reco::PFCandidate::ParticleType::h:
0726                 map_of_MEs[m_directory + "/chargedHadPFC_jetMatched_" + histName]->Fill(
0727                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0728                 break;
0729               case reco::PFCandidate::ParticleType::h0:
0730                 map_of_MEs[m_directory + "/neutralHadPFC_jetMatched_" + histName]->Fill(
0731                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0732                 break;
0733               case reco::PFCandidate::ParticleType::e:
0734                 map_of_MEs[m_directory + "/electronPFC_jetMatched_" + histName]->Fill(
0735                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0736                 break;
0737               case reco::PFCandidate::ParticleType::mu:
0738                 map_of_MEs[m_directory + "/muonPFC_jetMatched_" + histName]->Fill(
0739                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0740                 break;
0741               case reco::PFCandidate::ParticleType::gamma:
0742                 map_of_MEs[m_directory + "/gammaPFC_jetMatched_" + histName]->Fill(
0743                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0744                 break;
0745               case reco::PFCandidate::ParticleType::h_HF:
0746                 map_of_MEs[m_directory + "/hadHFPFC_jetMatched_" + histName]->Fill(
0747                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0748                 break;
0749               case reco::PFCandidate::ParticleType::egamma_HF:
0750                 map_of_MEs[m_directory + "/emHFPFC_jetMatched_" + histName]->Fill(
0751                     m_funcMap[m_observableNames[i]](*recoPF), eventWeight);
0752                 break;
0753               default:
0754                 break;
0755             }
0756           }
0757 
0758           for (unsigned int i = 0; i < m_pfInJetObservableNames.size(); i++) {
0759             std::string histName = Form("%s%s_jetCuts%s_%s",
0760                                         m_pfInJetObservableNames[i].c_str(),
0761                                         binString.c_str(),
0762                                         jetBinString.c_str(),
0763                                         npvString.c_str());
0764             map_of_MEs[m_directory + "/allPFC_jetMatched_" + histName]->Fill(
0765                 m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0766 
0767             switch (recoPF->particleId()) {
0768               case reco::PFCandidate::ParticleType::h:
0769                 map_of_MEs[m_directory + "/chargedHadPFC_jetMatched_" + histName]->Fill(
0770                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0771                 break;
0772               case reco::PFCandidate::ParticleType::h0:
0773                 map_of_MEs[m_directory + "/neutralHadPFC_jetMatched_" + histName]->Fill(
0774                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0775                 break;
0776               case reco::PFCandidate::ParticleType::e:
0777                 map_of_MEs[m_directory + "/electronPFC_jetMatched_" + histName]->Fill(
0778                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0779                 break;
0780               case reco::PFCandidate::ParticleType::mu:
0781                 map_of_MEs[m_directory + "/muonPFC_jetMatched_" + histName]->Fill(
0782                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0783                 break;
0784               case reco::PFCandidate::ParticleType::gamma:
0785                 map_of_MEs[m_directory + "/gammaPFC_jetMatched_" + histName]->Fill(
0786                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0787                 break;
0788               case reco::PFCandidate::ParticleType::h_HF:
0789                 map_of_MEs[m_directory + "/hadHFPFC_jetMatched_" + histName]->Fill(
0790                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0791                 break;
0792               case reco::PFCandidate::ParticleType::egamma_HF:
0793                 map_of_MEs[m_directory + "/emHFPFC_jetMatched_" + histName]->Fill(
0794                     m_pfInJetFuncMap[m_pfInJetObservableNames[i]](*recoPF, *cjet), eventWeight);
0795                 break;
0796               default:
0797                 break;
0798             }
0799           }
0800         }
0801 
0802         for (unsigned int i = 0; i < m_eventObservableNames.size(); i++) {
0803           std::string histName =
0804               Form("%s_jetCuts%s_%s", m_eventObservableNames[i].c_str(), jetBinString.c_str(), npvString.c_str());
0805           map_of_MEs[m_directory + "/allPFC_jetMatched_" + histName]->Fill(
0806               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::X), eventWeight);
0807           map_of_MEs[m_directory + "/chargedHadPFC_jetMatched_" + histName]->Fill(
0808               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::h), eventWeight);
0809           map_of_MEs[m_directory + "/neutralHadPFC_jetMatched_" + histName]->Fill(
0810               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::h0),
0811               eventWeight);
0812           map_of_MEs[m_directory + "/electronPFC_jetMatched_" + histName]->Fill(
0813               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::e), eventWeight);
0814           map_of_MEs[m_directory + "/muonPFC_jetMatched_" + histName]->Fill(
0815               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::mu),
0816               eventWeight);
0817           map_of_MEs[m_directory + "/gammaPFC_jetMatched_" + histName]->Fill(
0818               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::gamma),
0819               eventWeight);
0820           map_of_MEs[m_directory + "/hadHFPFC_jetMatched_" + histName]->Fill(
0821               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::h_HF),
0822               eventWeight);
0823           map_of_MEs[m_directory + "/emHFPFC_jetMatched_" + histName]->Fill(
0824               m_jetWideFuncMap[m_eventObservableNames[i]](pfConstits, reco::PFCandidate::ParticleType::egamma_HF),
0825               eventWeight);
0826         }
0827       }
0828     }
0829   }
0830 }
[ 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!