CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HGCalHitCalibration.cc
Go to the documentation of this file.
1 // user include files
2 
6 
10 
12 
14 
19 
26 
31 
37 
39 
40 #include <map>
41 #include <array>
42 #include <string>
43 #include <numeric>
44 
46  public:
47  explicit HGCalHitCalibration(const edm::ParameterSet&);
48  ~HGCalHitCalibration() override;
49 
50  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
51 
52  private:
54  edm::EventSetup const&) override;
55  void analyze(const edm::Event&, const edm::EventSetup&) override;
56 
57  void fillWithRecHits(std::map<DetId, const HGCRecHit*>&, DetId, unsigned int,
58  float, int&, float&);
59 
67 
68  int algo_;
70  int debug_;
72 
73  std::map<int, MonitorElement*> h_EoP_CPene_calib_fraction_;
74  std::map<int, MonitorElement*> hgcal_EoP_CPene_calib_fraction_;
75  std::map<int, MonitorElement*> hgcal_ele_EoP_CPene_calib_fraction_;
76  std::map<int, MonitorElement*> hgcal_photon_EoP_CPene_calib_fraction_;
78 
79  static constexpr int layers_ = 60;
80  std::array<float, layers_> Energy_layer_calib_;
81  std::array<float, layers_> Energy_layer_calib_fraction_;
82 };
83 
85  : rawRecHits_(iConfig.getParameter<bool>("rawRecHits")),
86  debug_(iConfig.getParameter<int>("debug")) {
87  auto detector = iConfig.getParameter<std::string>("detector");
88  auto recHitsEE = iConfig.getParameter<edm::InputTag>("recHitsEE");
89  auto recHitsFH = iConfig.getParameter<edm::InputTag>("recHitsFH");
90  auto recHitsBH = iConfig.getParameter<edm::InputTag>("recHitsBH");
91  auto caloParticles = iConfig.getParameter<edm::InputTag>("caloParticles");
92  auto hgcalMultiClusters = iConfig.getParameter<edm::InputTag>("hgcalMultiClusters");
93  auto electrons = iConfig.getParameter<edm::InputTag>("electrons");
94  auto photons = iConfig.getParameter<edm::InputTag>("photons");
95  if (detector == "all") {
96  recHitsEE_ = consumes<HGCRecHitCollection>(recHitsEE);
97  recHitsFH_ = consumes<HGCRecHitCollection>(recHitsFH);
98  recHitsBH_ = consumes<HGCRecHitCollection>(recHitsBH);
99  algo_ = 1;
100  } else if (detector == "EM") {
101  recHitsEE_ = consumes<HGCRecHitCollection>(recHitsEE);
102  algo_ = 2;
103  } else if (detector == "HAD") {
104  recHitsFH_ = consumes<HGCRecHitCollection>(recHitsFH);
105  recHitsBH_ = consumes<HGCRecHitCollection>(recHitsBH);
106  algo_ = 3;
107  }
108  caloParticles_ = consumes<std::vector<CaloParticle> >(caloParticles);
109  hgcalMultiClusters_ = consumes<std::vector<reco::PFCluster> >(hgcalMultiClusters);
110  electrons_ = consumes<std::vector<reco::GsfElectron> >(electrons);
111  photons_ = consumes<std::vector<reco::Photon> >(photons);
112 
113  // size should be HGC layers 52 is enough
114  Energy_layer_calib_.fill(0.);
116 }
117 
119  // do anything here that needs to be done at desctruction time
120  // (e.g. close files, deallocate resources etc.)
121 }
122 
124  edm::Run const& iRun,
125  edm::EventSetup const& /* iSetup */) {
126  ibooker.cd();
127  ibooker.setCurrentFolder("HGCalHitCalibration");
129  ibooker.book1D("h_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
131  ibooker.book1D("h_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
133  ibooker.book1D("h_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
135  ibooker.book1D("hgcal_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
137  ibooker.book1D("hgcal_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
139  ibooker.book1D("hgcal_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
141  ibooker.book1D("hgcal_ele_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
143  ibooker.book1D("hgcal_ele_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
145  ibooker.book1D("hgcal_ele_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
147  ibooker.book1D("hgcal_photon_EoP_CPene_100_calib_fraction", "", 1000, -0.5, 2.5);
149  ibooker.book1D("hgcal_photon_EoP_CPene_200_calib_fraction", "", 1000, -0.5, 2.5);
151  ibooker.book1D("hgcal_photon_EoP_CPene_300_calib_fraction", "", 1000, -0.5, 2.5);
152  LayerOccupancy_ = ibooker.book1D("LayerOccupancy", "", layers_, 0., (float)layers_);
153 }
154 
156  std::map<DetId, const HGCRecHit*>& hitmap, const DetId hitid,
157  const unsigned int hitlayer, const float fraction, int& seedDet,
158  float& seedEnergy) {
159  if (hitmap.find(hitid) == hitmap.end()) {
160  // Hit was not reconstructed
161  IfLogTrace(debug_ > 0, "HGCalHitCalibration")
162  << ">>> Failed to find detid " << std::hex << hitid.rawId() << std::dec
163  << " Det " << hitid.det()
164  << " Subdet " << hitid.subdetId() << std::endl;
165  return;
166  }
167  if (hitid.det() != DetId::Forward) {
168  IfLogTrace(debug_ > 0, "HGCalHitCalibration")
169  << ">>> Wrong type of detid " << std::hex << hitid.rawId() << std::dec
170  << " Det " << hitid.det()
171  << " Subdet " << hitid.subdetId() << std::endl;
172  return;
173  }
174  unsigned int layer = recHitTools_.getLayerWithOffset(hitid);
175  assert(hitlayer == layer);
176  Energy_layer_calib_fraction_[layer] += hitmap[hitid]->energy() * fraction;
177  LayerOccupancy_->Fill(layer);
178  if (seedEnergy < hitmap[hitid]->energy()) {
179  seedEnergy = hitmap[hitid]->energy();
180  seedDet = recHitTools_.getSiThickness(hitid);
181  }
182 }
183 
185  const edm::EventSetup& iSetup) {
186  using namespace edm;
187 
188  recHitTools_.getEventSetup(iSetup);
189 
190  Handle<HGCRecHitCollection> recHitHandleEE;
191  Handle<HGCRecHitCollection> recHitHandleFH;
192  Handle<HGCRecHitCollection> recHitHandleBH;
193 
194  Handle<std::vector<CaloParticle> > caloParticleHandle;
195  iEvent.getByToken(caloParticles_, caloParticleHandle);
196  const std::vector<CaloParticle>& caloParticles = *caloParticleHandle;
197 
198  Handle<std::vector<reco::PFCluster> > hgcalMultiClustersHandle;
199  iEvent.getByToken(hgcalMultiClusters_, hgcalMultiClustersHandle);
200 
201  Handle<std::vector<reco::GsfElectron> > PFElectronHandle;
202  iEvent.getByToken(electrons_, PFElectronHandle);
203 
204  Handle<std::vector<reco::Photon> > PFPhotonHandle;
205  iEvent.getByToken(photons_, PFPhotonHandle);
206 
207  // make a map detid-rechit
208  std::map<DetId, const HGCRecHit*> hitmap;
209  switch (algo_) {
210  case 1: {
211  iEvent.getByToken(recHitsEE_, recHitHandleEE);
212  iEvent.getByToken(recHitsFH_, recHitHandleFH);
213  iEvent.getByToken(recHitsBH_, recHitHandleBH);
214  const auto& rechitsEE = *recHitHandleEE;
215  const auto& rechitsFH = *recHitHandleFH;
216  const auto& rechitsBH = *recHitHandleBH;
217  for (unsigned int i = 0; i < rechitsEE.size(); ++i) {
218  hitmap[rechitsEE[i].detid()] = &rechitsEE[i];
219  }
220  for (unsigned int i = 0; i < rechitsFH.size(); ++i) {
221  hitmap[rechitsFH[i].detid()] = &rechitsFH[i];
222  }
223  for (unsigned int i = 0; i < rechitsBH.size(); ++i) {
224  hitmap[rechitsBH[i].detid()] = &rechitsBH[i];
225  }
226  break;
227  }
228  case 2: {
229  iEvent.getByToken(recHitsEE_, recHitHandleEE);
230  const HGCRecHitCollection& rechitsEE = *recHitHandleEE;
231  for (unsigned int i = 0; i < rechitsEE.size(); i++) {
232  hitmap[rechitsEE[i].detid()] = &rechitsEE[i];
233  }
234  break;
235  }
236  case 3: {
237  iEvent.getByToken(recHitsFH_, recHitHandleFH);
238  iEvent.getByToken(recHitsBH_, recHitHandleBH);
239  const auto& rechitsFH = *recHitHandleFH;
240  const auto& rechitsBH = *recHitHandleBH;
241  for (unsigned int i = 0; i < rechitsFH.size(); i++) {
242  hitmap[rechitsFH[i].detid()] = &rechitsFH[i];
243  }
244  for (unsigned int i = 0; i < rechitsBH.size(); i++) {
245  hitmap[rechitsBH[i].detid()] = &rechitsBH[i];
246  }
247  break;
248  }
249  default:
250  assert(false);
251  break;
252  }
253 
254  // loop over caloParticles
255  int seedDet = 0;
256  float seedEnergy = 0.;
257  IfLogTrace(debug_ > 0, "HGCalHitCalibration")
258  << "Number of caloParticles: " << caloParticles.size() << std::endl;
259  for (const auto& it_caloPart : caloParticles) {
260  const SimClusterRefVector& simClusterRefVector = it_caloPart.simClusters();
261  Energy_layer_calib_.fill(0.);
263 
264  seedDet = 0;
265  seedEnergy = 0.;
266  for (const auto& it_sc : simClusterRefVector) {
267  const SimCluster& simCluster = (*(it_sc));
268  IfLogTrace(debug_ > 1, "HGCalHitCalibration")
269  << ">>> SC.energy(): " << simCluster.energy()
270  << " SC.simEnergy(): " << simCluster.simEnergy()
271  << std::endl;
272  const std::vector<std::pair<uint32_t, float> >& hits_and_fractions =
273  simCluster.hits_and_fractions();
274 
275  // loop over hits
276  for (const auto& it_haf : hits_and_fractions) {
277  unsigned int hitlayer = recHitTools_.getLayerWithOffset(it_haf.first);
278  DetId hitid = (it_haf.first);
279  // dump raw RecHits and match
280  if (rawRecHits_) {
281  if ((hitid.det() == DetId::Forward &&
282  (hitid.subdetId() == HGCEE or hitid.subdetId() == HGCHEF or
283  hitid.subdetId() == HGCHEB)) ||
284  (hitid.det() == DetId::Hcal && hitid.subdetId() == HcalEndcap))
285  fillWithRecHits(hitmap, hitid, hitlayer, it_haf.second, seedDet,
286  seedEnergy);
287  }
288  } // end simHit
289  } // end simCluster
290 
291  auto sumCalibRecHitCalib_fraction = std::accumulate(Energy_layer_calib_fraction_.begin(),
292  Energy_layer_calib_fraction_.end(), 0.);
293  IfLogTrace(debug_ > 0, "HGCalHitCalibration")
294  << ">>> MC Energy: " << it_caloPart.energy()
295  << " reco energy: " << sumCalibRecHitCalib_fraction
296  << std::endl;
297  if (h_EoP_CPene_calib_fraction_.count(seedDet))
298  h_EoP_CPene_calib_fraction_[seedDet]->Fill(sumCalibRecHitCalib_fraction /
299  it_caloPart.energy());
300 
301  // Loop on reconstructed SC.
302  const auto& clusters = *hgcalMultiClustersHandle;
303  float total_energy = 0.;
304  float max_dR2 = 0.0025;
305  auto closest = std::min_element(clusters.begin(),
306  clusters.end(),
307  [&](const reco::PFCluster &a,
308  const reco::PFCluster &b) {
309  auto dR2_a = reco::deltaR2(it_caloPart, a);
310  auto dR2_b = reco::deltaR2(it_caloPart, b);
311  auto ERatio_a = a.correctedEnergy()/it_caloPart.energy();
312  auto ERatio_b = b.correctedEnergy()/it_caloPart.energy();
313  // If both clusters are within 0.0025, mark as first (min) the
314  // element with the highest ratio against the SimCluster
315  if (dR2_a < max_dR2 && dR2_b < max_dR2)
316  return ERatio_a > ERatio_b;
317  return dR2_a < dR2_b;
318  });
319  if (closest != clusters.end()
320  && reco::deltaR2(*closest, it_caloPart) < 0.01) {
321  total_energy = closest->correctedEnergy();
322  seedDet = recHitTools_.getSiThickness(closest->seed());
323  if (hgcal_EoP_CPene_calib_fraction_.count(seedDet)) {
324  hgcal_EoP_CPene_calib_fraction_[seedDet]->Fill(total_energy /
325  it_caloPart.energy());
326  }
327  }
328 
329  auto closest_fcn = [&](auto const &a, auto const&b){
330  auto dR2_a = reco::deltaR2(it_caloPart, a);
331  auto dR2_b = reco::deltaR2(it_caloPart, b);
332  auto ERatio_a = a.energy()/it_caloPart.energy();
333  auto ERatio_b = b.energy()/it_caloPart.energy();
334  // If both clusters are within 0.0025, mark as first (min) the
335  // element with the highest ratio against the SimCluster
336  if (dR2_a < max_dR2 && dR2_b < max_dR2)
337  return ERatio_a > ERatio_b;
338  return dR2_a < dR2_b;
339  };
340  // ELECTRONS in HGCAL
341  if (PFElectronHandle.isValid()) {
342  auto const & ele = (*PFElectronHandle);
343  auto closest = std::min_element(ele.begin(),
344  ele.end(),
345  closest_fcn);
346  if (closest != ele.end()
347  && closest->superCluster()->seed()->seed().det() == DetId::Forward
348  && reco::deltaR2(*closest, it_caloPart) < 0.01) {
349  seedDet = recHitTools_.getSiThickness(closest->superCluster()->seed()->seed());
350  if (hgcal_ele_EoP_CPene_calib_fraction_.count(seedDet)) {
351  hgcal_ele_EoP_CPene_calib_fraction_[seedDet]->Fill(closest->energy() /
352  it_caloPart.energy());
353  }
354  }
355  }
356 
357  // PHOTONS in HGCAL
358  if (PFPhotonHandle.isValid()) {
359  auto const & photon = (*PFPhotonHandle);
360  auto closest = std::min_element(photon.begin(),
361  photon.end(),
362  closest_fcn);
363  if (closest != photon.end()
364  && closest->superCluster()->seed()->seed().det() == DetId::Forward
365  && reco::deltaR2(*closest, it_caloPart) < 0.01) {
366  seedDet = recHitTools_.getSiThickness(closest->superCluster()->seed()->seed());
367  if (hgcal_photon_EoP_CPene_calib_fraction_.count(seedDet)) {
368  hgcal_photon_EoP_CPene_calib_fraction_[seedDet]->Fill(closest->energy() /
369  it_caloPart.energy());
370  }
371  }
372  }
373  } // end caloparticle
374 }
375 
376 // ------------ method fills 'descriptions' with the allowed parameters for the
377 // module ------------
379  edm::ConfigurationDescriptions& descriptions) {
381  desc.add<int>("debug", 0);
382  desc.add<bool>("rawRecHits", true);
383  desc.add<std::string>("detector", "all");
384  desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
385  desc.add<edm::InputTag>("recHitsEE", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
386  desc.add<edm::InputTag>("recHitsFH", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
387  desc.add<edm::InputTag>("recHitsBH", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
388  desc.add<edm::InputTag>("hgcalMultiClusters", edm::InputTag("particleFlowClusterHGCalFromMultiCl"));
389  desc.add<edm::InputTag>("electrons", edm::InputTag("ecalDrivenGsfElectronsFromMultiCl"));
390  desc.add<edm::InputTag>("photons", edm::InputTag("photonsFromMultiCl"));
391  descriptions.add("hgcalHitCalibration", desc);
392 }
393 
394 // define this as a plug-in
T getParameter(std::string const &) const
edm::EDGetTokenT< std::vector< CaloParticle > > caloParticles_
edm::EDGetTokenT< std::vector< reco::Photon > > photons_
std::array< float, layers_ > Energy_layer_calib_fraction_
edm::EDGetTokenT< std::vector< reco::PFCluster > > hgcalMultiClusters_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
std::map< int, MonitorElement * > hgcal_EoP_CPene_calib_fraction_
double correctedEnergy() const
Definition: CaloCluster.h:125
edm::EDGetTokenT< HGCRecHitCollection > recHitsFH_
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:47
edm::EDGetTokenT< HGCRecHitCollection > recHitsEE_
std::vector< std::pair< uint32_t, float > > hits_and_fractions() const
Returns list of rechit IDs and fractions for this SimCluster.
Definition: SimCluster.h:210
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:519
std::map< int, MonitorElement * > hgcal_ele_EoP_CPene_calib_fraction_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
#define IfLogTrace(cond, cat)
HGCalHitCalibration(const edm::ParameterSet &)
std::map< int, MonitorElement * > hgcal_photon_EoP_CPene_calib_fraction_
#define constexpr
void getEventSetup(const edm::EventSetup &)
Definition: RecHitTools.cc:61
void Fill(long long x)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
MonitorElement * LayerOccupancy_
int iEvent
Definition: GenABIO.cc:230
edm::EDGetTokenT< HGCRecHitCollection > recHitsBH_
Monte Carlo truth information used for tracking validation.
Definition: SimCluster.h:28
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:118
float energy() const
Energy. Note this is taken from the first SimTrack only.
Definition: SimCluster.h:111
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool isValid() const
Definition: HandleBase.h:74
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
double energy() const
cluster energy
Definition: PFCluster.h:82
hgcal::RecHitTools recHitTools_
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
float simEnergy() const
returns the accumulated sim energy in the cluster
Definition: SimCluster.h:225
std::float_t getSiThickness(const DetId &) const
Definition: RecHitTools.cc:103
Definition: DetId.h:18
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:279
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:179
edm::EDGetTokenT< std::vector< reco::GsfElectron > > electrons_
std::map< int, MonitorElement * > h_EoP_CPene_calib_fraction_
double b
Definition: hdecay.h:120
void add(std::string const &label, ParameterSetDescription const &psetDescription)
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
Definition: deltaR.h:36
void analyze(const edm::Event &, const edm::EventSetup &) override
HLT enums.
void fillWithRecHits(std::map< DetId, const HGCRecHit * > &, DetId, unsigned int, float, int &, float &)
size_type size() const
double a
Definition: hdecay.h:121
std::array< float, layers_ > Energy_layer_calib_
Detector det() const
get the detector field from this detid
Definition: DetId.h:35
Definition: Run.h:43