CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HGCalRecHitValidation.cc
Go to the documentation of this file.
1 // system include files
2 #include <cmath>
3 #include <fstream>
4 #include <iostream>
5 #include <map>
6 #include <string>
7 #include <vector>
8 
15 
17 
20 
29 
37 
38 #include "CLHEP/Units/GlobalSystemOfUnits.h"
39 #include "TVector3.h"
40 
42 
43 public:
44  struct energysum{
46  double e15, e25, e50, e100, e250, e1000;
47  };
48 
49  struct HitsInfo{
51  x=y=z=time=energy=phi=eta=0.0;
52  layer=0;
53  }
54  float x, y, z, time, energy, phi, eta ;
55  float layer;
56  };
57 
58 
60  ~HGCalRecHitValidation() override {}
61 
62  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
63  void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
64  void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
65  void analyze(const edm::Event&, const edm::EventSetup&) override;
66 
67 private:
68  template<class T1, class T2>
69  void recHitValidation(DetId & detId, int layer, const T1* geom, T2 it);
70  void fillHitsInfo();
71  void fillHitsInfo(HitsInfo& hits);
72  void fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer);
73 
74  // ----------member data ---------------------------
77  bool ifHCAL_;
79  unsigned int layers_;
80  std::map<int, int> OccupancyMap_plus;
81  std::map<int, int> OccupancyMap_minus;
82 
83  std::vector<MonitorElement*> EtaPhi_Plus_;
84  std::vector<MonitorElement*> EtaPhi_Minus_;
85  std::vector<MonitorElement*> energy_;
86  std::vector<MonitorElement*> HitOccupancy_Plus_;
87  std::vector<MonitorElement*> HitOccupancy_Minus_;
90 };
91 
93  nameDetector_(iConfig.getParameter<std::string>("DetectorName")),
94  ifHCAL_(iConfig.getParameter<bool>("ifHCAL")),
95  verbosity_(iConfig.getUntrackedParameter<int>("Verbosity",0)) {
96 
97  auto temp = iConfig.getParameter<edm::InputTag>("RecHitSource");
98  if (nameDetector_ == "HGCalEESensitive" ||
99  nameDetector_ == "HGCalHESiliconSensitive" ||
100  nameDetector_ == "HGCalHEScintillatorSensitive") {
101  recHitSource_ = consumes<HGCRecHitCollection>(temp);
102  } else if (nameDetector_ == "HCal") {
103  if (ifHCAL_) recHitSource_ = consumes<HBHERecHitCollection>(temp);
104  else recHitSource_ = consumes<HGChebRecHitCollection>(temp);
105  } else {
106  throw cms::Exception("BadHGCRecHitSource")
107  << "HGCal DetectorName given as " << nameDetector_ << " must be: "
108  << "\"HGCalHESiliconSensitive\", \"HGCalHESiliconSensitive\", "
109  << "\"HGCalHEScintillatorSensitive\", or \"HCal\"!";
110  }
111 }
112 
113 
116  desc.add<std::string>("DetectorName","HGCalEESensitive");
117  desc.add<edm::InputTag>("RecHitSource",edm::InputTag("HGCalRecHit","HGCEERecHits"));
118  desc.add<bool>("ifHCAL",false);
119  desc.addUntracked<int>("Verbosity",0);
120  descriptions.add("hgcalRecHitValidationEE",desc);
121 }
122 
124  const edm::EventSetup& iSetup) {
125  OccupancyMap_plus.clear();
126  OccupancyMap_minus.clear();
127 
128  bool ok(true);
129  unsigned int ntot(0), nused(0);
130  if (nameDetector_ == "HCal") {
132  iSetup.get<CaloGeometryRecord>().get(geom);
133  if (!geom.isValid()) {
134  edm::LogVerbatim("HGCalValidation") << "Cannot get valid HGCalGeometry "
135  << "Object for " << nameDetector_;
136  } else {
137  const CaloGeometry* geom0 = geom.product();
138 
139  if (ifHCAL_) {
141  iEvent.getByToken(recHitSource_, hbhecoll);
142  if (hbhecoll.isValid()) {
143  if (verbosity_>0)
144  edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with "
145  << hbhecoll->size()
146  << " element(s)";
147  for (const auto & it : *(hbhecoll.product())) {
148  DetId detId = it.id();
149  ntot++;
150  if (detId.subdetId() == HcalEndcap) {
151  nused++;
152  int layer = HcalDetId(detId).depth();
153  recHitValidation(detId, layer, geom0, &it);
154  }
155  }
156  } else {
157  ok = false;
158  edm::LogVerbatim("HGCalValidation") << "HBHERecHitCollection Handle does not "
159  << "exist !!!";
160  }
161  } else {
163  iEvent.getByToken(recHitSource_, hbhecoll);
164  if (hbhecoll.isValid()) {
165  if (verbosity_>0)
166  edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with "
167  << hbhecoll->size()
168  << " element(s)";
169  for (const auto & it : *(hbhecoll.product())) {
170  DetId detId = it.id();
171  ntot++; nused++;
172  int layer = HcalDetId(detId).depth();
173  recHitValidation(detId, layer, geom0, &it);
174  }
175  } else {
176  ok = false;
177  edm::LogVerbatim("HGCalValidation") << "HGChebRecHitCollection Handle does not "
178  << "exist !!!";
179  }
180  }
181  }
182  } else {
184  iSetup.get<IdealGeometryRecord>().get(nameDetector_, geom);
185  if (!geom.isValid()) {
186  edm::LogVerbatim("HGCalValidation") << "Cannot get valid HGCalGeometry Object for "
187  << nameDetector_;
188  } else {
189  const HGCalGeometry* geom0 = geom.product();
191  int geomType = (((mode == HGCalGeometryMode::Hexagon8) ||
192  (mode == HGCalGeometryMode::Hexagon8Full)) ? 1 :
193  ((mode == HGCalGeometryMode::Trapezoid) ? 2 : 0));
194 
195  edm::Handle<HGCRecHitCollection> theRecHitContainers;
196  iEvent.getByToken(recHitSource_, theRecHitContainers);
197  if (theRecHitContainers.isValid()) {
198  if (verbosity_>0)
199  edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with "
200  << theRecHitContainers->size()
201  << " element(s)";
202  for (const auto & it : *(theRecHitContainers.product())) {
203  ntot++; nused++;
204  DetId detId = it.id();
205  int layer = ((geomType == 0) ? HGCalDetId(detId).layer() :
206  ((geomType == 1) ? HGCSiliconDetId(detId).layer() :
207  HGCScintillatorDetId(detId).layer()));
208  recHitValidation(detId, layer, geom0, &it);
209  }
210  } else {
211  ok = false;
212  edm::LogVerbatim("HGCalValidation") << "HGCRecHitCollection Handle does not "
213  << "exist !!!";
214  }
215  }
216  }
217  if (ok) fillHitsInfo();
218  edm::LogVerbatim("HGCalValidation") << "Event " << iEvent.id().event()
219  << " with " << ntot << " total and "
220  << nused << " used recHits";
221 }
222 
223 template<class T1, class T2>
225  const T1* geom, T2 it) {
226 
227  const GlobalPoint& global = geom->getPosition(detId);
228  double energy = it->energy();
229 
230  float globalx = global.x();
231  float globaly = global.y();
232  float globalz = global.z();
233 
234  HitsInfo hinfo;
235  hinfo.energy = energy;
236  hinfo.x = globalx;
237  hinfo.y = globaly;
238  hinfo.z = globalz;
239  hinfo.layer = std::min((unsigned)layer,layers_);
240  hinfo.phi = global.phi();
241  hinfo.eta = global.eta();
242 
243  if (verbosity_>1)
244  edm::LogVerbatim("HGCalValidation") << "-------------------------- gx = "
245  << globalx << " gy = " << globaly
246  << " gz = " << globalz << " phi = "
247  << hinfo.phi << " eta = " << hinfo.eta << " lay = " << hinfo.layer;
248 
249  fillHitsInfo(hinfo);
250 
251  if (hinfo.eta > 0) fillOccupancyMap(OccupancyMap_plus, hinfo.layer -1);
252  else fillOccupancyMap(OccupancyMap_minus, hinfo.layer -1);
253 
254 }
255 
256 void HGCalRecHitValidation::fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer){
257  if (OccupancyMap.find(layer) != OccupancyMap.end()) OccupancyMap[layer]++;
258  else OccupancyMap[layer] = 1;
259 }
260 
262 
263  for (auto itr = OccupancyMap_plus.begin() ; itr != OccupancyMap_plus.end(); ++itr) {
264  int layer = (*itr).first;
265  int occupancy = (*itr).second;
266  HitOccupancy_Plus_.at(layer)->Fill(occupancy);
267  }
268 
269  for (auto itr = OccupancyMap_minus.begin() ; itr != OccupancyMap_minus.end(); ++itr) {
270  int layer = (*itr).first;
271  int occupancy = (*itr).second;
272  HitOccupancy_Minus_.at(layer)->Fill(occupancy);
273  }
274 
275 }
276 
278 
279  unsigned int ilayer = hits.layer -1;
280  energy_.at(ilayer)->Fill(hits.energy);
281 
282  EtaPhi_Plus_.at(ilayer)->Fill(hits.eta , hits.phi);
283  EtaPhi_Minus_.at(ilayer)->Fill(hits.eta, hits.phi);
284 
285 }
286 
288  const edm::EventSetup& iSetup) {
289 
290  if (nameDetector_ == "HCal") {
292  iSetup.get<HcalRecNumberingRecord>().get( pHRNDC );
293  const HcalDDDRecConstants *hcons = &(*pHRNDC);
294  layers_ = hcons->getMaxDepth(1);
295  } else {
297  iSetup.get<IdealGeometryRecord>().get(nameDetector_, pHGDC);
298  const HGCalDDDConstants & hgcons_ = (*pHGDC);
299  layers_ = hgcons_.layers(true);
300  }
301 }
302 
304  edm::Run const&,
305  edm::EventSetup const&) {
306 
307  iB.setCurrentFolder("HGCAL/HGCalRecHitsV/"+nameDetector_);
308  std::ostringstream histoname;
309  for (unsigned int ilayer = 0; ilayer < layers_; ilayer++ ) {
310  histoname.str(""); histoname << "HitOccupancy_Plus_layer_" << ilayer;
311  HitOccupancy_Plus_.push_back(iB.book1D( histoname.str().c_str(), "RecHitOccupancy_Plus", 100, 0, 10000));
312  histoname.str(""); histoname << "HitOccupancy_Minus_layer_" << ilayer;
313  HitOccupancy_Minus_.push_back(iB.book1D( histoname.str().c_str(), "RecHitOccupancy_Minus", 100, 0, 10000));
314 
315  histoname.str(""); histoname << "EtaPhi_Plus_" << "layer_" << ilayer;
316  EtaPhi_Plus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, 1.45, 3.0, 72, -CLHEP::pi, CLHEP::pi));
317  histoname.str(""); histoname << "EtaPhi_Minus_" << "layer_" << ilayer;
318  EtaPhi_Minus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, -3.0, -1.45, 72, -CLHEP::pi, CLHEP::pi));
319 
320  histoname.str(""); histoname << "energy_layer_" << ilayer;
321  energy_.push_back(iB.book1D(histoname.str().c_str(),"energy_",100,0,0.002));
322  }//loop over layers ends here
323 
324  histoname.str(""); histoname << "SUMOfRecHitOccupancy_Plus";
325  MeanHitOccupancy_Plus_= iB.book1D( histoname.str().c_str(), "SUMOfRecHitOccupancy_Plus", layers_, -0.5, layers_-0.5);
326  histoname.str(""); histoname << "SUMOfRecHitOccupancy_Minus";
327  MeanHitOccupancy_Minus_ = iB.book1D( histoname.str().c_str(), "SUMOfRecHitOccupancy_Minus", layers_, -0.5,layers_-0.5);
328 }
329 
330 
331 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
332 
334 
335 //define this as a plug-in
HGCalGeometryMode::GeometryMode geomMode() const
Geometry mode.
Definition: HGCalTopology.h:86
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
HGCalRecHitValidation(const edm::ParameterSet &)
std::vector< MonitorElement * > EtaPhi_Minus_
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:579
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
std::vector< MonitorElement * > EtaPhi_Plus_
std::vector< MonitorElement * > HitOccupancy_Minus_
const Double_t pi
MonitorElement * MeanHitOccupancy_Plus_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
int depth() const
get the tower depth
Definition: HcalDetId.h:162
int iEvent
Definition: GenABIO.cc:230
std::map< int, int > OccupancyMap_minus
unsigned int layers(bool reco) const
int layer() const
get the layer #
T z() const
Definition: PV3DBase.h:64
const HGCalTopology & topology() const
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:41
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:118
MonitorElement * MeanHitOccupancy_Minus_
T min(T a, T b)
Definition: MathUtil.h:58
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool isValid() const
Definition: HandleBase.h:74
void dqmBeginRun(const edm::Run &, const edm::EventSetup &) override
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: DetId.h:18
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
void fillOccupancyMap(std::map< int, int > &OccupancyMap, int layer)
void analyze(const edm::Event &, const edm::EventSetup &) override
T const * product() const
Definition: Handle.h:81
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:136
int getMaxDepth(const int &type) const
void add(std::string const &label, ParameterSetDescription const &psetDescription)
int layer() const
get the layer #
T eta() const
Definition: PV3DBase.h:76
std::map< int, int > OccupancyMap_plus
edm::EventID id() const
Definition: EventBase.h:60
size_type size() const
T get() const
Definition: EventSetup.h:63
std::vector< MonitorElement * > HitOccupancy_Plus_
bool isValid() const
Definition: ESHandle.h:47
T x() const
Definition: PV3DBase.h:62
T const * product() const
Definition: ESHandle.h:86
void recHitValidation(DetId &detId, int layer, const T1 *geom, T2 it)
std::vector< MonitorElement * > energy_
int layer() const
get the layer #
Definition: HGCalDetId.h:48
Definition: Run.h:44