CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
PFJetDQMPostProcessor Class Reference
Inheritance diagram for PFJetDQMPostProcessor:
DQMEDHarvester edm::one::EDProducer< edm::EndLuminosityBlockProducer, edm::EndRunProducer, edm::EndProcessBlockProducer, edm::one::WatchLuminosityBlocks, edm::one::WatchRuns, edm::one::SharedResources, edm::Accumulator > edm::one::EDProducerBase edm::ProducerBase edm::EDConsumerBase edm::ProductRegistryHelper

Public Member Functions

 PFJetDQMPostProcessor (const edm::ParameterSet &)
 
 ~PFJetDQMPostProcessor () override
 
- Public Member Functions inherited from DQMEDHarvester
void accumulate (edm::Event const &ev, edm::EventSetup const &es) final
 
void beginJob () override
 
void beginLuminosityBlock (edm::LuminosityBlock const &lumi, edm::EventSetup const &) final
 
void beginRun (edm::Run const &run, edm::EventSetup const &) override
 
virtual void dqmAnalyze (DQMStore::IBooker &, DQMStore::IGetter &, edm::Event const &, edm::EventSetup const &)
 
 DQMEDHarvester (edm::ParameterSet const &iConfig)
 
 DQMEDHarvester ()
 
virtual void dqmEndLuminosityBlock (DQMStore::IBooker &, DQMStore::IGetter &, edm::LuminosityBlock const &, edm::EventSetup const &)
 
virtual void dqmEndRun (DQMStore::IBooker &, DQMStore::IGetter &, edm::Run const &, edm::EventSetup const &)
 
void endLuminosityBlock (edm::LuminosityBlock const &, edm::EventSetup const &) final
 
void endLuminosityBlockProduce (edm::LuminosityBlock &lumi, edm::EventSetup const &es) final
 
void endProcessBlockProduce (edm::ProcessBlock &) final
 
void endRun (edm::Run const &, edm::EventSetup const &) override
 
void endRunProduce (edm::Run &run, edm::EventSetup const &es) final
 
 ~DQMEDHarvester () override=default
 
- Public Member Functions inherited from edm::one::EDProducer< edm::EndLuminosityBlockProducer, edm::EndRunProducer, edm::EndProcessBlockProducer, edm::one::WatchLuminosityBlocks, edm::one::WatchRuns, edm::one::SharedResources, edm::Accumulator >
 EDProducer ()=default
 
 EDProducer (const EDProducer &)=delete
 
SerialTaskQueueglobalLuminosityBlocksQueue () final
 
SerialTaskQueueglobalRunsQueue () final
 
bool hasAbilityToProduceInBeginLumis () const final
 
bool hasAbilityToProduceInBeginProcessBlocks () const final
 
bool hasAbilityToProduceInBeginRuns () const final
 
bool hasAbilityToProduceInEndLumis () const final
 
bool hasAbilityToProduceInEndProcessBlocks () const final
 
bool hasAbilityToProduceInEndRuns () const final
 
const EDProduceroperator= (const EDProducer &)=delete
 
bool wantsGlobalLuminosityBlocks () const final
 
bool wantsGlobalRuns () const final
 
bool wantsInputProcessBlocks () const final
 
bool wantsProcessBlocks () const final
 
- Public Member Functions inherited from edm::one::EDProducerBase
 EDProducerBase ()
 
ModuleDescription const & moduleDescription () const
 
bool wantsStreamLuminosityBlocks () const
 
bool wantsStreamRuns () const
 
 ~EDProducerBase () override
 
- Public Member Functions inherited from edm::ProducerBase
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
std::vector< edm::ProductResolverIndex > const & indiciesForPutProducts (BranchType iBranchType) const
 
 ProducerBase ()
 
std::vector< edm::ProductResolverIndex > const & putTokenIndexToProductResolverIndex () const
 
std::vector< bool > const & recordProvenanceList () const
 
void registerProducts (ProducerBase *, ProductRegistry *, ModuleDescription const &)
 
std::function< void(BranchDescription const &)> registrationCallback () const
 used by the fwk to register list of products More...
 
void resolvePutIndicies (BranchType iBranchType, ModuleToResolverIndicies const &iIndicies, std::string const &moduleLabel)
 
TypeLabelList const & typeLabelList () const
 used by the fwk to register the list of products of this module More...
 
 ~ProducerBase () noexcept(false) override
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfoconsumesInfo () const
 
void convertCurrentProcessAlias (std::string const &processName)
 Convert "@currentProcess" in InputTag process names to the actual current process name. More...
 
 EDConsumerBase ()
 
 EDConsumerBase (EDConsumerBase const &)=delete
 
 EDConsumerBase (EDConsumerBase &&)=default
 
ESResolverIndex const * esGetTokenIndices (edm::Transition iTrans) const
 
std::vector< ESResolverIndex > const & esGetTokenIndicesVector (edm::Transition iTrans) const
 
std::vector< ESRecordIndex > const & esGetTokenRecordIndicesVector (edm::Transition iTrans) const
 
ProductResolverIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
std::vector< ProductResolverIndexAndSkipBit > const & itemsToGetFrom (BranchType iType) const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesWhoseProductsAreConsumed (std::array< std::vector< ModuleDescription const *> *, NumBranchTypes > &modulesAll, std::vector< ModuleProcessName > &modulesInPreviousProcesses, ProductRegistry const &preg, std::map< std::string, ModuleDescription const *> const &labelsToDesc, std::string const &processName) const
 
EDConsumerBase const & operator= (EDConsumerBase const &)=delete
 
EDConsumerBaseoperator= (EDConsumerBase &&)=default
 
bool registeredToConsume (ProductResolverIndex, bool, BranchType) const
 
void selectInputProcessBlocks (ProductRegistry const &productRegistry, ProcessBlockHelperBase const &processBlockHelperBase)
 
ProductResolverIndexAndSkipBit uncheckedIndexFrom (EDGetToken) const
 
void updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
 
void updateLookup (eventsetup::ESRecordsToProductResolverIndices const &)
 
virtual ~EDConsumerBase () noexcept(false)
 

Private Member Functions

void dqmEndJob (DQMStore::IBooker &, DQMStore::IGetter &) override
 
void fitResponse (TH1F *hreso, TH1F *h_genjet_pt, int ptbinlow, int ietahigh, double recoptcut, double &resp, double &resp_err, double &reso, double &reso_err)
 
double getRespUnc (double width, double width_err, double mean, double mean_err)
 
std::string seta (double eta)
 
std::string spt (double ptmin, double ptmax)
 

Private Attributes

bool debug = false
 
std::vector< double > etaBins
 
std::string genjetDir
 
std::vector< std::string > jetResponseDir
 
std::string offsetDir
 
std::vector< double > ptBins
 
double recoptcut
 

Additional Inherited Members

- Public Types inherited from DQMEDHarvester
typedef dqm::harvesting::DQMStore DQMStore
 
typedef dqm::harvesting::MonitorElement MonitorElement
 
- Public Types inherited from edm::one::EDProducerBase
typedef EDProducerBase ModuleType
 
- Public Types inherited from edm::ProducerBase
template<typename T >
using BranchAliasSetterT = ProductRegistryHelper::BranchAliasSetterT< T >
 
using ModuleToResolverIndicies = std::unordered_multimap< std::string, std::tuple< edm::TypeID const *, const char *, edm::ProductResolverIndex > >
 
typedef ProductRegistryHelper::TypeLabelList TypeLabelList
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::one::EDProducerBase
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &descriptions)
 
- Protected Member Functions inherited from edm::ProducerBase
template<Transition Tr = Transition::Event>
auto produces (std::string instanceName) noexcept
 declare what type of product will make and with which optional label More...
 
template<Transition B>
BranchAliasSetter produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)
 
template<BranchType B>
BranchAliasSetter produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)
 
BranchAliasSetter produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)
 
template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType > produces (std::string instanceName)
 
template<class ProductType >
BranchAliasSetterT< ProductType > produces ()
 
template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType > produces (std::string instanceName)
 
template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType > produces ()
 
template<class ProductType >
BranchAliasSetterT< ProductType > produces (std::string instanceName)
 
template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType > produces ()
 
template<Transition Tr = Transition::Event>
auto produces () noexcept
 
ProducesCollector producesCollector ()
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
template<BranchType B = InEvent>
EDConsumerBaseAdaptor< Bconsumes (edm::InputTag tag) noexcept
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto esConsumes ()
 
template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto esConsumes (ESInputTag const &tag)
 
template<Transition Tr = Transition::Event>
constexpr auto esConsumes ()
 
template<Transition Tr = Transition::Event>
auto esConsumes (ESInputTag tag)
 
template<Transition Tr = Transition::Event>
ESGetTokenGeneric esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
 Used with EventSetupRecord::doGet. More...
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
void resetItemsToGetFrom (BranchType iType)
 
- Protected Attributes inherited from DQMEDHarvester
DQMStoredqmstore_
 
edm::GetterOfProducts< DQMTokenjobmegetter_
 
edm::EDPutTokenT< DQMTokenjobToken_
 
edm::GetterOfProducts< DQMTokenlumimegetter_
 
edm::EDPutTokenT< DQMTokenlumiToken_
 
edm::GetterOfProducts< DQMTokenrunmegetter_
 
edm::EDPutTokenT< DQMTokenrunToken_
 

Detailed Description

Definition at line 27 of file PFJetDQMPostProcessor.cc.

Constructor & Destructor Documentation

◆ PFJetDQMPostProcessor()

PFJetDQMPostProcessor::PFJetDQMPostProcessor ( const edm::ParameterSet iConfig)
explicit

Definition at line 72 of file PFJetDQMPostProcessor.cc.

References etaBins, genjetDir, edm::ParameterSet::getParameter(), jetResponseDir, offsetDir, ptBins, recoptcut, and AlCaHLTBitMon_QueryRunRegistry::string.

72  {
73  jetResponseDir = iConfig.getParameter<std::vector<std::string>>("jetResponseDir");
74  genjetDir = iConfig.getParameter<std::string>("genjetDir");
75  offsetDir = iConfig.getParameter<std::string>("offsetDir");
76  ptBins = iConfig.getParameter<std::vector<double>>("ptBins");
77  etaBins = iConfig.getParameter<std::vector<double>>("etaBins");
78  recoptcut = iConfig.getParameter<double>("recoPtCut");
79 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
std::vector< std::string > jetResponseDir
std::vector< double > etaBins
std::vector< double > ptBins

◆ ~PFJetDQMPostProcessor()

PFJetDQMPostProcessor::~PFJetDQMPostProcessor ( )
override

Definition at line 81 of file PFJetDQMPostProcessor.cc.

81 {}

Member Function Documentation

◆ dqmEndJob()

void PFJetDQMPostProcessor::dqmEndJob ( DQMStore::IBooker ibook_,
DQMStore::IGetter iget_ 
)
overrideprivatevirtual

Implements DQMEDHarvester.

Definition at line 84 of file PFJetDQMPostProcessor.cc.

References dqm::implementation::IBooker::book1D(), filterCSVwithJSON::copy, debug, submitPVResolutionJobs::err, etaBins, spr::find(), fitResponse(), genjetDir, dqm::implementation::IGetter::get(), dqm::implementation::IGetter::getMEs(), getRespUnc(), dqm::legacy::MonitorElement::getTH1F(), mps_fire::i, hcalRecHitTable_cff::ieta, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, jetResponseDir, hlt_dqm_clientPB-live_cfg::me, SiStripPI::mean, median(), createIOVlist::nEvents, HLTTauDQMOffline_cfi::nPtBins, offsetDir, ptBins, recoptcut, alignCSCRings::s, seta(), dqm::implementation::NavigatorBase::setCurrentFolder(), spt(), AlCaHLTBitMon_QueryRunRegistry::string, x, testProducerWithPsetDescEmpty_cfi::x1, y, and testProducerWithPsetDescEmpty_cfi::y1.

84  {
85  for (unsigned int idir = 0; idir < jetResponseDir.size(); idir++) {
87  std::vector<std::string> sME_genjets = iget_.getMEs();
88  std::for_each(sME_genjets.begin(), sME_genjets.end(), [&](auto& s) { s.insert(0, genjetDir); });
89 
91  std::vector<std::string> sME_offset = iget_.getMEs();
92  std::for_each(sME_offset.begin(), sME_offset.end(), [&](auto& s) { s.insert(0, offsetDir); });
93 
94  iget_.setCurrentFolder(jetResponseDir[idir]);
95  std::vector<std::string> sME_response = iget_.getMEs();
96  std::for_each(sME_response.begin(), sME_response.end(), [&](auto& s) { s.insert(0, jetResponseDir[idir]); });
97 
98  iget_.setCurrentFolder(jetResponseDir[idir]);
99 
100  double ptBinsArray[ptBins.size()];
101  unsigned int nPtBins = ptBins.size() - 1;
102  std::copy(ptBins.begin(), ptBins.end(), ptBinsArray);
103  //for(unsigned int ipt = 0; ipt < ptBins.size(); ++ipt) std::cout << ptBins[ipt] << std::endl;
104 
105  std::string stitle;
106  char ctitle[50];
107  std::vector<MonitorElement*> vME_presponse;
108  std::vector<MonitorElement*> vME_presponse_mean;
109  std::vector<MonitorElement*> vME_presponse_median;
110  std::vector<MonitorElement*> vME_preso;
111  std::vector<MonitorElement*> vME_preso_rms;
112  std::vector<MonitorElement*> vME_efficiency;
113  std::vector<MonitorElement*> vME_purity;
114  std::vector<MonitorElement*> vME_ratePUJet;
115 
117  TH1F* h_resp;
118  TH1F *h_genjet_pt, *h_genjet_matched_pt = nullptr;
119  TH1F *h_recojet_pt = nullptr, *h_recojet_matched_pt = nullptr, *h_recojet_unmatched_pt = nullptr;
120 
121  stitle = offsetDir + "mu";
122  std::vector<std::string>::const_iterator it = std::find(sME_offset.begin(), sME_offset.end(), stitle);
123  if (it == sME_offset.end())
124  continue;
125  me = iget_.get(stitle);
126  int nEvents = ((TH1F*)me->getTH1F())->GetEntries();
127  if (nEvents == 0)
128  continue;
129  iget_.setCurrentFolder(jetResponseDir[idir]);
130 
131  bool isNoJEC = (jetResponseDir[idir].find("noJEC") != std::string::npos);
132  bool isJEC = false;
133  if (!isNoJEC)
134  isJEC = (jetResponseDir[idir].find("JEC") != std::string::npos);
135 
136  //
137  // Response distributions
138  //
139  for (unsigned int ieta = 1; ieta < etaBins.size(); ++ieta) {
140  stitle = genjetDir + "genjet_pt" + "_eta" + seta(etaBins[ieta]);
141 
142  std::vector<std::string>::const_iterator it = std::find(sME_genjets.begin(), sME_genjets.end(), stitle);
143  if (it == sME_genjets.end())
144  continue;
145  me = iget_.get(stitle);
146  h_genjet_pt = (TH1F*)me->getTH1F();
147 
148  if (isNoJEC) {
149  // getting the histogram for matched gen jets
150  stitle = jetResponseDir[idir] + "genjet_pt" + "_eta" + seta(etaBins[ieta]) + "_matched";
151  me = iget_.get(stitle);
152  h_genjet_matched_pt = (TH1F*)me->getTH1F();
153 
154  /*// getting the histogram for unmatched gen jets
155  stitle = jetResponseDir[idir] + "genjet_pt" + "_eta" + seta(etaBins[ieta]) + "_unmatched";
156  me = iget_.get(stitle);
157  h_genjet_unmatched_pt = (TH1F*)me->getTH1F();*/
158  }
159  if (isJEC) {
160  // getting the histogram for reco jets
161  stitle = jetResponseDir[idir] + "recojet_pt" + "_eta" + seta(etaBins[ieta]);
162  me = iget_.get(stitle);
163  h_recojet_pt = (TH1F*)me->getTH1F();
164 
165  // getting the histogram for matched reco jets
166  stitle = jetResponseDir[idir] + "recojet_pt" + "_eta" + seta(etaBins[ieta]) + "_matched";
167  me = iget_.get(stitle);
168  h_recojet_matched_pt = (TH1F*)me->getTH1F();
169 
170  // getting the histogram for unmatched reco jets
171  stitle = jetResponseDir[idir] + "recojet_pt" + "_eta" + seta(etaBins[ieta]) + "_unmatched";
172  me = iget_.get(stitle);
173  h_recojet_unmatched_pt = (TH1F*)me->getTH1F();
174  }
175 
176  stitle = "presponse_eta" + seta(etaBins[ieta]);
177  // adding "Raw" to the title of raw jet response histograms
178  if (isNoJEC) {
179  sprintf(ctitle, "Raw Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
180  } else {
181  sprintf(ctitle, "Jet pT response, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
182  }
183  TH1F* h_presponse = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
184 
185  stitle = "presponse_eta" + seta(etaBins[ieta]) + "_mean";
186  // adding "Raw" to the title of raw jet response histograms
187  if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
188  sprintf(ctitle, "Raw Jet pT response using Mean, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
189  } else {
190  sprintf(ctitle, "Jet pT response using Mean, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
191  }
192  TH1F* h_presponse_mean = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
193 
194  stitle = "presponse_eta" + seta(etaBins[ieta]) + "_median";
195  // adding "Raw" to the title of raw jet response histograms
196  if (isNoJEC) {
197  sprintf(ctitle, "Raw Jet pT response using Med., %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
198  } else {
199  sprintf(ctitle, "Jet pT response using Med., %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
200  }
201  TH1F* h_presponse_median = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
202 
203  stitle = "preso_eta" + seta(etaBins[ieta]);
204  if (isNoJEC) {
205  sprintf(ctitle, "Raw Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
206  } else {
207  sprintf(ctitle, "Jet pT resolution, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
208  }
209  TH1F* h_preso = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
210 
211  stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
212  if (jetResponseDir[idir].find("noJEC") != std::string::npos) {
213  sprintf(ctitle, "Raw Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
214  } else {
215  sprintf(ctitle, "Jet pT resolution using RMS, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
216  }
217  TH1F* h_preso_rms = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
218 
219  //Booking histogram for Jet Efficiency vs pT
220  stitle = "efficiency_eta" + seta(etaBins[ieta]);
221  sprintf(ctitle, "Efficiency, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
222  TH1F* h_efficiency = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
223 
224  //Booking histogram for Jet Purity vs pT
225  stitle = "purity_eta" + seta(etaBins[ieta]);
226  sprintf(ctitle, "Purity, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
227  TH1F* h_purity = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
228 
229  //Booking histogram for #PU jets vs pT
230  stitle = "ratePUJet_eta" + seta(etaBins[ieta]);
231  sprintf(ctitle, "PU Jet Rate, %4.1f<|#eta|<%4.1f", etaBins[ieta - 1], etaBins[ieta]);
232  TH1F* h_ratePUJet = new TH1F(stitle.c_str(), ctitle, nPtBins, ptBinsArray);
233 
234  if (isNoJEC) {
235  h_efficiency->Divide(h_genjet_matched_pt, h_genjet_pt, 1, 1, "B");
236  }
237  if (isJEC) {
238  h_purity->Divide(h_recojet_matched_pt, h_recojet_pt, 1, 1, "B");
239  h_ratePUJet = (TH1F*)h_recojet_unmatched_pt->Clone();
240  h_ratePUJet->SetName("h_ratePUJet");
241  h_ratePUJet->Scale(1. / double(nEvents));
242  }
243 
244  for (unsigned int ipt = 0; ipt < ptBins.size() - 1; ++ipt) {
245  stitle = jetResponseDir[idir] + "reso_dist_" + spt(ptBins[ipt], ptBins[ipt + 1]) + "_eta" + seta(etaBins[ieta]);
246  std::vector<std::string>::const_iterator it = std::find(sME_response.begin(), sME_response.end(), stitle);
247  if (it == sME_response.end())
248  continue;
249  me = iget_.get(stitle);
250  h_resp = (TH1F*)me->getTH1F();
251 
252  // Fit-based
253  double resp = 1.0, resp_err = 0.0, reso = 0.0, reso_err = 0.0;
254  fitResponse(h_resp, h_genjet_pt, ipt, ieta, recoptcut, resp, resp_err, reso, reso_err);
255 
256  h_presponse->SetBinContent(ipt + 1, resp);
257  h_presponse->SetBinError(ipt + 1, resp_err);
258  h_preso->SetBinContent(ipt + 1, reso);
259  h_preso->SetBinError(ipt + 1, reso_err);
260 
261  // RMS-based
262  double std = h_resp->GetStdDev();
263  double std_error = h_resp->GetStdDevError();
264 
265  // Scale each bin with mean response
266  double mean = 1.0;
267  double mean_error = 0.0;
268  double err = 0.0;
269  if (h_resp->GetMean() > 0) {
270  mean = h_resp->GetMean();
271  mean_error = h_resp->GetMeanError();
272 
273  // Scale resolution by response.
274  std /= mean;
275  std_error /= mean;
276 
277  err = getRespUnc(std, std_error, mean, mean_error);
278  }
279 
280  h_preso_rms->SetBinContent(ipt + 1, std);
281  h_preso_rms->SetBinError(ipt + 1, err);
282 
283  // Mean-based
284  h_presponse_mean->SetBinContent(ipt + 1, h_resp->GetMean());
285  h_presponse_mean->SetBinError(ipt + 1, h_resp->GetMeanError());
286 
287  // Median-based
288  if (h_resp->GetEntries() > 0) {
289  int numBins = h_resp->GetXaxis()->GetNbins();
290  Double_t x1[numBins];
291  Double_t y1[numBins];
292  for (int i = 0; i < numBins; i++) {
293  x1[i] = h_resp->GetBinCenter(i + 1);
294  y1[i] = h_resp->GetBinContent(i + 1) > 0 ? h_resp->GetBinContent(i + 1) : 0.0;
295  }
296  const auto x = x1, y = y1;
297  double median = TMath::Median(numBins, x, y);
298  h_presponse_median->SetBinContent(ipt + 1, median);
299  h_presponse_median->SetBinError(ipt + 1, 1.2533 * (h_resp->GetMeanError()));
300  }
301 
302  } // ipt
303 
304  if (isNoJEC) {
305  stitle = "efficiency_eta" + seta(etaBins[ieta]);
306  me = ibook_.book1D(stitle.c_str(), h_efficiency);
307  vME_efficiency.push_back(me);
308  }
309 
310  if (isJEC) {
311  stitle = "purity_eta" + seta(etaBins[ieta]);
312  me = ibook_.book1D(stitle.c_str(), h_purity);
313  vME_purity.push_back(me);
314 
315  stitle = "ratePUJet_eta" + seta(etaBins[ieta]);
316  me = ibook_.book1D(stitle.c_str(), h_ratePUJet);
317  vME_ratePUJet.push_back(me);
318  }
319 
320  stitle = "presponse_eta" + seta(etaBins[ieta]);
321  me = ibook_.book1D(stitle.c_str(), h_presponse);
322  vME_presponse.push_back(me);
323 
324  stitle = "presponse_eta" + seta(etaBins[ieta]) + "_mean";
325  me = ibook_.book1D(stitle.c_str(), h_presponse_mean);
326  vME_presponse_mean.push_back(me);
327 
328  stitle = "presponse_eta" + seta(etaBins[ieta]) + "_median";
329  me = ibook_.book1D(stitle.c_str(), h_presponse_median);
330  vME_presponse_median.push_back(me);
331 
332  stitle = "preso_eta" + seta(etaBins[ieta]);
333  me = ibook_.book1D(stitle.c_str(), h_preso);
334  vME_preso.push_back(me);
335 
336  stitle = "preso_eta" + seta(etaBins[ieta]) + "_rms";
337  me = ibook_.book1D(stitle.c_str(), h_preso_rms);
338  vME_preso_rms.push_back(me);
339 
340  } // ieta
341 
342  //
343  // Checks
344  //
345  if (debug) {
346  if (isNoJEC) {
347  for (std::vector<MonitorElement*>::const_iterator i = vME_efficiency.begin(); i != vME_efficiency.end(); ++i)
348  (*i)->getTH1F()->Print();
349  }
350  if (isJEC) {
351  for (std::vector<MonitorElement*>::const_iterator i = vME_purity.begin(); i != vME_purity.end(); ++i)
352  (*i)->getTH1F()->Print();
353  for (std::vector<MonitorElement*>::const_iterator i = vME_ratePUJet.begin(); i != vME_ratePUJet.end(); ++i)
354  (*i)->getTH1F()->Print();
355  }
356 
357  for (std::vector<MonitorElement*>::const_iterator i = vME_presponse.begin(); i != vME_presponse.end(); ++i)
358  (*i)->getTH1F()->Print();
359  for (std::vector<MonitorElement*>::const_iterator i = vME_presponse_mean.begin(); i != vME_presponse_mean.end();
360  ++i)
361  (*i)->getTH1F()->Print();
362  for (std::vector<MonitorElement*>::const_iterator i = vME_presponse_median.begin();
363  i != vME_presponse_median.end();
364  ++i)
365  (*i)->getTH1F()->Print();
366  for (std::vector<MonitorElement*>::const_iterator i = vME_preso.begin(); i != vME_preso.end(); ++i)
367  (*i)->getTH1F()->Print();
368  for (std::vector<MonitorElement*>::const_iterator i = vME_preso_rms.begin(); i != vME_preso_rms.end(); ++i)
369  (*i)->getTH1F()->Print();
370  }
371  }
372 }
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
std::vector< std::string > jetResponseDir
virtual std::vector< std::string > getMEs() const
Definition: DQMStore.cc:759
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
std::string spt(double ptmin, double ptmax)
std::vector< double > etaBins
virtual TH1F * getTH1F() const
void fitResponse(TH1F *hreso, TH1F *h_genjet_pt, int ptbinlow, int ietahigh, double recoptcut, double &resp, double &resp_err, double &reso, double &reso_err)
virtual MonitorElement * get(std::string const &fullpath) const
Definition: DQMStore.cc:712
T median(std::vector< T > values)
Definition: median.h:16
std::vector< double > ptBins
std::string seta(double eta)
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
double getRespUnc(double width, double width_err, double mean, double mean_err)

◆ fitResponse()

void PFJetDQMPostProcessor::fitResponse ( TH1F *  hreso,
TH1F *  h_genjet_pt,
int  ptbinlow,
int  ietahigh,
double  recoptcut,
double &  resp,
double &  resp_err,
double &  reso,
double &  reso_err 
)
private

Definition at line 374 of file PFJetDQMPostProcessor.cc.

References etaBins, getRespUnc(), METSkim_cff::Max, ptBins, recoptcut, and seta().

Referenced by dqmEndJob().

382  {
383  // This 'smartfitter' is converted from the original Python smart_fit() -function
384  // implemented in test/helperFunctions.py. See that file for more commentary.
385  // Juska 23 May 2019
386 
387  // Only do plots if needed for debugging
388  // NOTE a directory called 'debug' must exist in the working directory
389  //
390  bool doPlots = false;
391 
392  double ptlow = ptBins[ptbinlow];
393  double pthigh = ptBins[ptbinlow + 1];
394 
395  // Take range by Mikko's advice: -1.5 and + 1.5 * RMS width
396 
397  double rmswidth = hreso->GetStdDev();
398  double rmsmean = hreso->GetMean();
399  double fitlow = rmsmean - 1.5 * rmswidth;
400  fitlow = TMath::Max(recoptcut / ptlow, fitlow);
401  double fithigh = rmsmean + 1.5 * rmswidth;
402 
403  TF1* fg = new TF1("mygaus", "gaus", fitlow, fithigh);
404  TF1* fg2 = new TF1("fg2", "TMath::Gaus(x,[0],[1],true)*[2]", fitlow, fithigh);
405 
406  hreso->Fit("mygaus", "RQN");
407 
408  fg2->SetParameter(0, fg->GetParameter(1));
409  fg2->SetParameter(1, fg->GetParameter(2));
410 
411  // Extract ngenjet in the current pT bin from the genjet histo
412  float ngenjet = h_genjet_pt->GetBinContent(ptbinlow + 1);
413 
414  // Here the fit is forced to take the area of ngenjets.
415  // The area is further normalized for the response histogram x-axis length
416  // (3) and number of bins (100)
417  fg2->FixParameter(2, ngenjet * 3. / 100.);
418 
419  hreso->Fit("fg2", "RQN");
420 
421  fitlow = fg2->GetParameter(0) - 1.5 * fg2->GetParameter(1);
422  fitlow = TMath::Max(recoptcut / ptlow, fitlow);
423  fithigh = fg2->GetParameter(0) + 1.5 * fg2->GetParameter(1);
424 
425  fg2->SetRange(fitlow, fithigh);
426 
427  hreso->Fit("fg2", "RQ");
428 
429  fg->SetRange(0, 3);
430  fg2->SetRange(0, 3);
431  fg->SetLineWidth(2);
432  fg2->SetLineColor(kGreen + 2);
433 
434  hreso->GetXaxis()->SetRangeUser(0, 2);
435 
436  // Save plots to a subdirectory if asked (debug-directory must exist!)
437  if (doPlots & (hreso->GetEntries() > 0)) {
438  TCanvas* cfit = new TCanvas(Form("respofit_%i", int(ptlow)), "respofit", 600, 600);
439  hreso->Draw("ehist");
440  fg->Draw("same");
441  fg2->Draw("same");
442  cfit->SaveAs(
443  Form("debug/respo_smartfit_%04d_%i_eta%s.pdf", (int)ptlow, (int)pthigh, seta(etaBins[ietahigh]).c_str()));
444  }
445 
446  resp = fg2->GetParameter(0);
447  resp_err = fg2->GetParError(0);
448 
449  // Scale resolution by response. Avoid division by zero.
450  if (0 == resp)
451  reso = 0;
452  else
453  reso = fg2->GetParameter(1) / resp;
454  reso_err = fg2->GetParError(1) / resp;
455 
456  reso_err = getRespUnc(reso, reso_err, resp, resp_err);
457 }
std::vector< double > etaBins
std::vector< double > ptBins
std::string seta(double eta)
double getRespUnc(double width, double width_err, double mean, double mean_err)

◆ getRespUnc()

double PFJetDQMPostProcessor::getRespUnc ( double  width,
double  width_err,
double  mean,
double  mean_err 
)
private

Definition at line 460 of file PFJetDQMPostProcessor.cc.

References SiStripPI::mean, funct::pow(), and ApeEstimator_cff::width.

Referenced by dqmEndJob(), and fitResponse().

460  {
461  if (0 == width || 0 == mean)
462  return 0;
463  return TMath::Sqrt(pow(width_err, 2) / pow(width, 2) + pow(mean_err, 2) / pow(mean, 2)) * width;
464 }
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ seta()

std::string PFJetDQMPostProcessor::seta ( double  eta)
inlineprivate

Definition at line 55 of file PFJetDQMPostProcessor.cc.

References PVValHelper::eta, AlCaHLTBitMon_QueryRunRegistry::string, and to_string().

Referenced by dqmEndJob(), and fitResponse().

55  {
56  std::string seta = std::to_string(int(eta * 10.));
57  if (seta.length() < 2)
58  seta = "0" + seta;
59  return seta;
60  }
static std::string to_string(const XMLCh *ch)
std::string seta(double eta)

◆ spt()

std::string PFJetDQMPostProcessor::spt ( double  ptmin,
double  ptmax 
)
inlineprivate

Definition at line 62 of file PFJetDQMPostProcessor.cc.

References muonTiming_cfi::ptmax, ptmin, AlCaHLTBitMon_QueryRunRegistry::string, and to_string().

Referenced by dqmEndJob().

62  {
64  return spt;
65  }
static std::string to_string(const XMLCh *ch)
std::string spt(double ptmin, double ptmax)
double ptmin
Definition: HydjetWrapper.h:86

Member Data Documentation

◆ debug

bool PFJetDQMPostProcessor::debug = false
private

Definition at line 53 of file PFJetDQMPostProcessor.cc.

Referenced by dqmEndJob().

◆ etaBins

std::vector<double> PFJetDQMPostProcessor::etaBins
private

Definition at line 49 of file PFJetDQMPostProcessor.cc.

Referenced by dqmEndJob(), fitResponse(), and PFJetDQMPostProcessor().

◆ genjetDir

std::string PFJetDQMPostProcessor::genjetDir
private

Definition at line 46 of file PFJetDQMPostProcessor.cc.

Referenced by dqmEndJob(), and PFJetDQMPostProcessor().

◆ jetResponseDir

std::vector<std::string> PFJetDQMPostProcessor::jetResponseDir
private

Definition at line 45 of file PFJetDQMPostProcessor.cc.

Referenced by dqmEndJob(), and PFJetDQMPostProcessor().

◆ offsetDir

std::string PFJetDQMPostProcessor::offsetDir
private

Definition at line 47 of file PFJetDQMPostProcessor.cc.

Referenced by dqmEndJob(), and PFJetDQMPostProcessor().

◆ ptBins

std::vector<double> PFJetDQMPostProcessor::ptBins
private

Definition at line 48 of file PFJetDQMPostProcessor.cc.

Referenced by dqmEndJob(), fitResponse(), and PFJetDQMPostProcessor().

◆ recoptcut

double PFJetDQMPostProcessor::recoptcut
private

Definition at line 51 of file PFJetDQMPostProcessor.cc.

Referenced by dqmEndJob(), fitResponse(), and PFJetDQMPostProcessor().