File indexing completed on 2021-12-24 02:22:28
0001
0002 #include <cmath>
0003 #include <iostream>
0004 #include <fstream>
0005 #include <vector>
0006 #include <map>
0007 #include <string>
0008
0009 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0010 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0011 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0012 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0013
0014 #include "DetectorDescription/Core/interface/DDCompactView.h"
0015 #include "DetectorDescription/Core/interface/DDSpecifics.h"
0016 #include "DetectorDescription/Core/interface/DDSolid.h"
0017 #include "DetectorDescription/Core/interface/DDFilter.h"
0018 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0019 #include "DetectorDescription/DDCMS/interface/DDCompactView.h"
0020 #include "DetectorDescription/DDCMS/interface/DDFilteredView.h"
0021
0022 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0023 #include "DQMServices/Core/interface/DQMStore.h"
0024
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0031 #include "FWCore/ServiceRegistry/interface/Service.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033
0034 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0035 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0036 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0037
0038 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0039 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0040 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0041 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0042
0043 #include "CLHEP/Geometry/Point3D.h"
0044 #include "CLHEP/Geometry/Transform3D.h"
0045 #include "CLHEP/Geometry/Vector3D.h"
0046 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0047 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0048
0049 class HGCalSimHitValidation : public DQMEDAnalyzer {
0050 public:
0051 struct energysum {
0052 energysum() {
0053 etotal = 0;
0054 for (int i = 0; i < 6; ++i)
0055 eTime[i] = 0.;
0056 }
0057 double eTime[6], etotal;
0058 };
0059
0060 struct hitsinfo {
0061 hitsinfo() {
0062 x = y = z = phi = eta = 0.0;
0063 cell = cell2 = sector = sector2 = type = layer = 0;
0064 }
0065 double x, y, z, phi, eta;
0066 int cell, cell2, sector, sector2, type, layer;
0067 };
0068
0069 explicit HGCalSimHitValidation(const edm::ParameterSet&);
0070 ~HGCalSimHitValidation() override {}
0071
0072 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0073
0074 protected:
0075 void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0076 void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0077 void analyze(const edm::Event&, const edm::EventSetup&) override;
0078
0079 private:
0080 void analyzeHits(std::vector<PCaloHit>& hits);
0081 void fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer);
0082 void fillHitsInfo(std::pair<hitsinfo, energysum> hit_, unsigned int itimeslice, double esum);
0083 bool defineGeometry(const DDCompactView* ddViewH);
0084 bool defineGeometry(const cms::DDCompactView* ddViewH);
0085
0086 TH1F* createHisto(std::string histname, const int nbins, float minIndexX, float maxIndexX, bool isLogX = true);
0087 void histoSetting(TH1F*& histo,
0088 const char* xTitle,
0089 const char* yTitle = "",
0090 Color_t lineColor = kBlack,
0091 Color_t markerColor = kBlack,
0092 int linewidth = 1);
0093 void histoSetting(TH2F*& histo,
0094 const char* xTitle,
0095 const char* yTitle = "",
0096 Color_t lineColor = kBlack,
0097 Color_t markerColor = kBlack,
0098 int linewidth = 1);
0099 void fillMuonTomoHistos(int partialType, std::pair<hitsinfo, energysum> hit_);
0100
0101
0102 const std::string nameDetector_, caloHitSource_;
0103 const HGCalDDDConstants* hgcons_;
0104 const std::vector<double> times_;
0105 const int verbosity_;
0106 const bool fromDDD_;
0107 const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> tok_hgcal_;
0108 const edm::ESGetToken<DDCompactView, IdealGeometryRecord> tok_cpv_;
0109 const edm::ESGetToken<cms::DDCompactView, IdealGeometryRecord> tok_cpvc_;
0110 edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0111 edm::EDGetTokenT<edm::HepMCProduct> tok_hepMC_;
0112 bool symmDet_;
0113 unsigned int layers_;
0114 int firstLayer_;
0115 std::map<uint32_t, HepGeom::Transform3D> transMap_;
0116
0117 std::vector<MonitorElement*> HitOccupancy_Plus_, HitOccupancy_Minus_;
0118 std::vector<MonitorElement*> EtaPhi_Plus_, EtaPhi_Minus_;
0119 MonitorElement *MeanHitOccupancy_Plus_, *MeanHitOccupancy_Minus_;
0120 static const unsigned int maxTime_ = 6;
0121 std::vector<MonitorElement*> energy_[maxTime_];
0122 std::vector<MonitorElement*> energyFWF_, energyFWCN_, energyFWCK_;
0123 std::vector<MonitorElement*> energyPWF_, energyPWCN_, energyPWCK_;
0124 std::vector<MonitorElement*> hitXYFWF_, hitXYFWCN_, hitXYFWCK_, hitXYB_;
0125 unsigned int nTimes_;
0126 };
0127
0128 HGCalSimHitValidation::HGCalSimHitValidation(const edm::ParameterSet& iConfig)
0129 : nameDetector_(iConfig.getParameter<std::string>("DetectorName")),
0130 caloHitSource_(iConfig.getParameter<std::string>("CaloHitSource")),
0131 times_(iConfig.getParameter<std::vector<double> >("TimeSlices")),
0132 verbosity_(iConfig.getUntrackedParameter<int>("Verbosity", 0)),
0133 fromDDD_(iConfig.getUntrackedParameter<bool>("fromDDD", true)),
0134 tok_hgcal_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0135 edm::ESInputTag{"", nameDetector_})),
0136 tok_cpv_(esConsumes<DDCompactView, IdealGeometryRecord, edm::Transition::BeginRun>()),
0137 tok_cpvc_(esConsumes<cms::DDCompactView, IdealGeometryRecord, edm::Transition::BeginRun>()),
0138 symmDet_(true),
0139 firstLayer_(1) {
0140 tok_hepMC_ = consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"));
0141 tok_hits_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", caloHitSource_));
0142 nTimes_ = (times_.size() > maxTime_) ? maxTime_ : times_.size();
0143 }
0144
0145 void HGCalSimHitValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0146 edm::ParameterSetDescription desc;
0147 std::vector<double> times = {25.0, 1000.0};
0148 desc.add<std::string>("DetectorName", "HGCalEESensitive");
0149 desc.add<std::string>("CaloHitSource", "HGCHitsEE");
0150 desc.add<std::vector<double> >("TimeSlices", times);
0151 desc.addUntracked<int>("Verbosity", 0);
0152 desc.addUntracked<bool>("TestNumber", true);
0153 desc.addUntracked<bool>("fromDDD", true);
0154 descriptions.add("hgcalSimHitValidationEE", desc);
0155 }
0156
0157 void HGCalSimHitValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0158
0159 if (verbosity_ > 0) {
0160 edm::Handle<edm::HepMCProduct> evtMC;
0161 iEvent.getByToken(tok_hepMC_, evtMC);
0162 if (!evtMC.isValid()) {
0163 edm::LogVerbatim("HGCalValidation") << "no HepMCProduct found";
0164 } else {
0165 const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
0166 unsigned int k(0);
0167 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0168 ++p, ++k) {
0169 edm::LogVerbatim("HGCalValidation") << "Particle[" << k << "] with pt " << (*p)->momentum().perp() << " eta "
0170 << (*p)->momentum().eta() << " phi " << (*p)->momentum().phi();
0171 }
0172 }
0173 }
0174
0175
0176 edm::Handle<edm::PCaloHitContainer> theCaloHitContainers;
0177 iEvent.getByToken(tok_hits_, theCaloHitContainers);
0178 if (theCaloHitContainers.isValid()) {
0179 if (verbosity_ > 0)
0180 edm::LogVerbatim("HGCalValidation") << " PcalohitItr = " << theCaloHitContainers->size();
0181 std::vector<PCaloHit> caloHits;
0182 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0183 analyzeHits(caloHits);
0184 } else if (verbosity_ > 0) {
0185 edm::LogVerbatim("HGCalValidation") << "PCaloHitContainer does not exist!";
0186 }
0187 }
0188
0189 void HGCalSimHitValidation::analyzeHits(std::vector<PCaloHit>& hits) {
0190 std::map<int, int> OccupancyMap_plus, OccupancyMap_minus;
0191 OccupancyMap_plus.clear();
0192 OccupancyMap_minus.clear();
0193
0194 std::map<uint32_t, std::pair<hitsinfo, energysum> > map_hits;
0195 map_hits.clear();
0196
0197 if (verbosity_ > 0)
0198 edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << hits.size() << " PcaloHit elements";
0199 unsigned int nused(0);
0200 for (unsigned int i = 0; i < hits.size(); i++) {
0201 double energy = hits[i].energy();
0202 double time = hits[i].time();
0203 uint32_t id_ = hits[i].id();
0204 int cell, sector, subsector(0), layer, zside;
0205 int cell2(0), type(0);
0206 if (hgcons_->waferHexagon8()) {
0207 HGCSiliconDetId detId = HGCSiliconDetId(id_);
0208 cell = detId.cellU();
0209 cell2 = detId.cellV();
0210 sector = detId.waferU();
0211 subsector = detId.waferV();
0212 type = detId.type();
0213 layer = detId.layer();
0214 zside = detId.zside();
0215 } else if (hgcons_->tileTrapezoid()) {
0216 HGCScintillatorDetId detId = HGCScintillatorDetId(id_);
0217 sector = detId.ietaAbs();
0218 cell = detId.iphi();
0219 subsector = 1;
0220 type = detId.type();
0221 layer = detId.layer();
0222 zside = detId.zside();
0223 } else {
0224 int subdet;
0225 HGCalTestNumbering::unpackHexagonIndex(id_, subdet, zside, layer, sector, type, cell);
0226 }
0227 nused++;
0228 if (verbosity_ > 1)
0229 edm::LogVerbatim("HGCalValidation")
0230 << "Detector " << nameDetector_ << " zside = " << zside << " sector|wafer = " << sector << ":" << subsector
0231 << " type = " << type << " layer = " << layer << " cell = " << cell << ":" << cell2 << " energy = " << energy
0232 << " energyem = " << hits[i].energyEM() << " energyhad = " << hits[i].energyHad() << " time = " << time;
0233
0234 HepGeom::Point3D<float> gcoord;
0235 std::pair<float, float> xy;
0236 if (hgcons_->waferHexagon8()) {
0237 xy = hgcons_->locateCell(layer, sector, subsector, cell, cell2, false, true);
0238 } else if (hgcons_->tileTrapezoid()) {
0239 xy = hgcons_->locateCellTrap(layer, sector, cell, false);
0240 } else {
0241 xy = hgcons_->locateCell(cell, layer, sector, false);
0242 }
0243 double zp = hgcons_->waferZ(layer, false);
0244 if (zside < 0)
0245 zp = -zp;
0246 float xp = (zp < 0) ? -xy.first : xy.first;
0247 gcoord = HepGeom::Point3D<float>(xp, xy.second, zp);
0248 double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
0249 if (verbosity_ > 1)
0250 edm::LogVerbatim("HGCalValidation")
0251 << std::hex << id_ << std::dec << " global coordinate " << gcoord << " time " << time << ":" << tof;
0252 time -= tof;
0253
0254 energysum esum;
0255 hitsinfo hinfo;
0256 if (map_hits.count(id_) != 0) {
0257 hinfo = map_hits[id_].first;
0258 esum = map_hits[id_].second;
0259 } else {
0260 hinfo.x = gcoord.x();
0261 hinfo.y = gcoord.y();
0262 hinfo.z = gcoord.z();
0263 hinfo.sector = sector;
0264 hinfo.sector2 = subsector;
0265 hinfo.cell = cell;
0266 hinfo.cell2 = cell2;
0267 hinfo.type = type;
0268 hinfo.layer = layer - firstLayer_;
0269 hinfo.phi = gcoord.getPhi();
0270 hinfo.eta = gcoord.getEta();
0271 }
0272 esum.etotal += energy;
0273 for (unsigned int k = 0; k < nTimes_; ++k) {
0274 if (time > 0 && time < times_[k])
0275 esum.eTime[k] += energy;
0276 }
0277
0278 if (verbosity_ > 1)
0279 edm::LogVerbatim("HGCalValidation") << " ----------------------- gx = " << hinfo.x << " gy = " << hinfo.y
0280 << " gz = " << hinfo.z << " phi = " << hinfo.phi << " eta = " << hinfo.eta;
0281 map_hits[id_] = std::pair<hitsinfo, energysum>(hinfo, esum);
0282 }
0283 if (verbosity_ > 0)
0284 edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << map_hits.size()
0285 << " detector elements being hit";
0286
0287 std::map<uint32_t, std::pair<hitsinfo, energysum> >::iterator itr;
0288 for (itr = map_hits.begin(); itr != map_hits.end(); ++itr) {
0289 hitsinfo hinfo = (*itr).second.first;
0290 energysum esum = (*itr).second.second;
0291 int layer = hinfo.layer;
0292 double eta = hinfo.eta;
0293 int type, part, orient;
0294 int partialType = -1;
0295 if (nameDetector_ == "HGCalEESensitive" or nameDetector_ == "HGCalHESiliconSensitive") {
0296 HGCSiliconDetId detId = HGCSiliconDetId((*itr).first);
0297 std::tie(type, part, orient) = hgcons_->waferType(detId);
0298 partialType = part;
0299 }
0300
0301 for (unsigned int itimeslice = 0; itimeslice < nTimes_; itimeslice++) {
0302 fillHitsInfo((*itr).second, itimeslice, esum.eTime[itimeslice]);
0303 }
0304
0305 if (eta > 0.0)
0306 fillOccupancyMap(OccupancyMap_plus, layer);
0307 else
0308 fillOccupancyMap(OccupancyMap_minus, layer);
0309
0310 fillMuonTomoHistos(partialType, (*itr).second);
0311 }
0312 if (verbosity_ > 0)
0313 edm::LogVerbatim("HGCalValidation") << "With map:used:total " << hits.size() << "|" << nused << "|"
0314 << map_hits.size() << " hits";
0315
0316 for (auto const& itr : OccupancyMap_plus) {
0317 int layer = itr.first;
0318 int occupancy = itr.second;
0319 HitOccupancy_Plus_.at(layer)->Fill(occupancy);
0320 }
0321 for (auto const& itr : OccupancyMap_minus) {
0322 int layer = itr.first;
0323 int occupancy = itr.second;
0324 HitOccupancy_Minus_.at(layer)->Fill(occupancy);
0325 }
0326 }
0327
0328 void HGCalSimHitValidation::fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer) {
0329 if (OccupancyMap.find(layer) != OccupancyMap.end()) {
0330 ++OccupancyMap[layer];
0331 } else {
0332 OccupancyMap[layer] = 1;
0333 }
0334 }
0335
0336 void HGCalSimHitValidation::fillHitsInfo(std::pair<hitsinfo, energysum> hits, unsigned int itimeslice, double esum) {
0337 unsigned int ilayer = hits.first.layer;
0338 if (ilayer < layers_) {
0339 energy_[itimeslice].at(ilayer)->Fill(esum);
0340 if (itimeslice == 0) {
0341 EtaPhi_Plus_.at(ilayer)->Fill(hits.first.eta, hits.first.phi);
0342 EtaPhi_Minus_.at(ilayer)->Fill(hits.first.eta, hits.first.phi);
0343 }
0344 } else {
0345 if (verbosity_ > 0)
0346 edm::LogVerbatim("HGCalValidation")
0347 << "Problematic Hit for " << nameDetector_ << " at sector " << hits.first.sector << ":" << hits.first.sector2
0348 << " layer " << hits.first.layer << " cell " << hits.first.cell << ":" << hits.first.cell2 << " energy "
0349 << hits.second.etotal;
0350 }
0351 }
0352
0353 void HGCalSimHitValidation::fillMuonTomoHistos(int partialType, std::pair<hitsinfo, energysum> hits) {
0354 hitsinfo hinfo = hits.first;
0355 energysum esum = hits.second;
0356 double edep =
0357 esum.eTime[0] * CLHEP::GeV /
0358 CLHEP::keV;
0359
0360 unsigned int ilayer = hinfo.layer;
0361 double x = hinfo.x * CLHEP::mm / CLHEP::cm;
0362 double y = hinfo.y * CLHEP::mm / CLHEP::cm;
0363 if (ilayer < layers_) {
0364 if (nameDetector_ == "HGCalEESensitive" or nameDetector_ == "HGCalHESiliconSensitive") {
0365
0366 if (!TMath::AreEqualAbs(edep, 0.0, 1.e-5)) {
0367 if (hinfo.type == HGCSiliconDetId::HGCalFine) {
0368 if (partialType == 0)
0369 energyFWF_.at(ilayer)->Fill(edep);
0370 if (partialType > 0)
0371 energyPWF_.at(ilayer)->Fill(edep);
0372 }
0373 if (hinfo.type == HGCSiliconDetId::HGCalCoarseThin) {
0374 if (partialType == 0)
0375 energyFWCN_.at(ilayer)->Fill(edep);
0376 if (partialType > 0)
0377 energyPWCN_.at(ilayer)->Fill(edep);
0378 }
0379 if (hinfo.type == HGCSiliconDetId::HGCalCoarseThick) {
0380 if (partialType == 0)
0381 energyFWCK_.at(ilayer)->Fill(edep);
0382 if (partialType > 0)
0383 energyPWCK_.at(ilayer)->Fill(edep);
0384 }
0385 }
0386
0387
0388 if (hinfo.type == HGCSiliconDetId::HGCalFine)
0389 hitXYFWF_.at(ilayer)->Fill(x, y);
0390
0391 if (hinfo.type == HGCSiliconDetId::HGCalCoarseThin)
0392 hitXYFWCN_.at(ilayer)->Fill(x, y);
0393
0394 if (hinfo.type == HGCSiliconDetId::HGCalCoarseThick)
0395 hitXYFWCK_.at(ilayer)->Fill(x, y);
0396
0397 }
0398 if (nameDetector_ == "HGCalHEScintillatorSensitive") {
0399 hitXYB_.at(ilayer)->Fill(x, y);
0400 }
0401 }
0402 }
0403
0404 bool HGCalSimHitValidation::defineGeometry(const DDCompactView* ddViewH) {
0405 if (verbosity_ > 0)
0406 edm::LogVerbatim("HGCalValidation") << "Initialize HGCalDDDConstants (DDD) for " << nameDetector_ << " : "
0407 << hgcons_;
0408
0409 if (hgcons_->geomMode() == HGCalGeometryMode::Square) {
0410 const DDCompactView& cview = *ddViewH;
0411 std::string attribute = "Volume";
0412 std::string value = nameDetector_;
0413
0414 DDSpecificsMatchesValueFilter filter{DDValue(attribute, value, 0)};
0415 DDFilteredView fv(cview, filter);
0416 bool dodet = fv.firstChild();
0417
0418 while (dodet) {
0419 const DDSolid& sol = fv.logicalPart().solid();
0420 const std::string& name = sol.name().fullname();
0421 int isd = (name.find(nameDetector_) == std::string::npos) ? -1 : 1;
0422 if (isd > 0) {
0423 std::vector<int> copy = fv.copyNumbers();
0424 int nsiz = static_cast<int>(copy.size());
0425 int lay = (nsiz > 0) ? copy[nsiz - 1] : -1;
0426 int sec = (nsiz > 1) ? copy[nsiz - 2] : -1;
0427 int zp = (nsiz > 3) ? copy[nsiz - 4] : -1;
0428 if (zp != 1)
0429 zp = -1;
0430 const DDTrap& trp = static_cast<DDTrap>(sol);
0431 int subs = (trp.alpha1() > 0 ? 1 : 0);
0432 symmDet_ = (trp.alpha1() == 0 ? true : false);
0433 uint32_t id = HGCalTestNumbering::packSquareIndex(zp, lay, sec, subs, 0);
0434 DD3Vector x, y, z;
0435 fv.rotation().GetComponents(x, y, z);
0436 const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0437 const CLHEP::HepRotation hr(rotation);
0438 const CLHEP::Hep3Vector h3v(fv.translation().X(), fv.translation().Y(), fv.translation().Z());
0439 const HepGeom::Transform3D ht3d(hr, h3v);
0440 transMap_.insert(std::make_pair(id, ht3d));
0441 if (verbosity_ > 2)
0442 edm::LogVerbatim("HGCalValidation") << HGCalDetId(id) << " Transform using " << h3v << " and " << hr;
0443 }
0444 dodet = fv.next();
0445 }
0446 if (verbosity_ > 0)
0447 edm::LogVerbatim("HGCalValidation") << "Finds " << transMap_.size() << " elements and SymmDet_ = " << symmDet_;
0448 }
0449 return true;
0450 }
0451
0452 bool HGCalSimHitValidation::defineGeometry(const cms::DDCompactView* ddViewH) {
0453 if (verbosity_ > 0)
0454 edm::LogVerbatim("HGCalValidation") << "Initialize HGCalDDDConstants (DD4hep) for " << nameDetector_ << " : "
0455 << hgcons_;
0456
0457 if (hgcons_->geomMode() == HGCalGeometryMode::Square) {
0458 const cms::DDCompactView& cview = *ddViewH;
0459 const cms::DDFilter filter("Volume", nameDetector_);
0460 cms::DDFilteredView fv(cview, filter);
0461
0462 while (fv.firstChild()) {
0463 const auto& name = fv.name();
0464 int isd = (name.find(nameDetector_) == std::string::npos) ? -1 : 1;
0465 if (isd > 0) {
0466 const auto& copy = fv.copyNos();
0467 int nsiz = static_cast<int>(copy.size());
0468 int lay = (nsiz > 0) ? copy[0] : -1;
0469 int sec = (nsiz > 1) ? copy[1] : -1;
0470 int zp = (nsiz > 3) ? copy[3] : -1;
0471 if (zp != 1)
0472 zp = -1;
0473 const auto& pars = fv.parameters();
0474 int subs = (pars[6] > 0 ? 1 : 0);
0475 symmDet_ = (pars[6] == 0 ? true : false);
0476 uint32_t id = HGCalTestNumbering::packSquareIndex(zp, lay, sec, subs, 0);
0477 DD3Vector x, y, z;
0478 fv.rotation().GetComponents(x, y, z);
0479 const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0480 const CLHEP::HepRotation hr(rotation);
0481 const CLHEP::Hep3Vector h3v(fv.translation().X(), fv.translation().Y(), fv.translation().Z());
0482 const HepGeom::Transform3D ht3d(hr, h3v);
0483 transMap_.insert(std::make_pair(id, ht3d));
0484 if (verbosity_ > 2)
0485 edm::LogVerbatim("HGCalValidation") << HGCalDetId(id) << " Transform using " << h3v << " and " << hr;
0486 }
0487 }
0488 if (verbosity_ > 0)
0489 edm::LogVerbatim("HGCalValidation") << "Finds " << transMap_.size() << " elements and SymmDet_ = " << symmDet_;
0490 }
0491 return true;
0492 }
0493
0494
0495 void HGCalSimHitValidation::dqmBeginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0496 hgcons_ = &iSetup.getData(tok_hgcal_);
0497 layers_ = hgcons_->layers(false);
0498 firstLayer_ = hgcons_->firstLayer();
0499 if (fromDDD_) {
0500 const DDCompactView* pDD = &iSetup.getData(tok_cpv_);
0501 defineGeometry(pDD);
0502 } else {
0503 const cms::DDCompactView* pDD = &iSetup.getData(tok_cpvc_);
0504 defineGeometry(pDD);
0505 }
0506 if (verbosity_ > 0)
0507 edm::LogVerbatim("HGCalValidation") << nameDetector_ << " defined with " << layers_ << " Layers with first at "
0508 << firstLayer_;
0509 }
0510
0511 void HGCalSimHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const&, edm::EventSetup const&) {
0512 iB.setCurrentFolder("HGCAL/HGCalSimHitsV/" + nameDetector_);
0513
0514 std::ostringstream histoname;
0515 for (unsigned int il = 0; il < layers_; ++il) {
0516 int ilayer = firstLayer_ + static_cast<int>(il);
0517 auto istr1 = std::to_string(ilayer);
0518 while (istr1.size() < 2) {
0519 istr1.insert(0, "0");
0520 }
0521 histoname.str("");
0522 histoname << "HitOccupancy_Plus_layer_" << istr1;
0523 HitOccupancy_Plus_.push_back(iB.book1D(histoname.str().c_str(), "HitOccupancy_Plus", 501, -0.5, 500.5));
0524 histoname.str("");
0525 histoname << "HitOccupancy_Minus_layer_" << istr1;
0526 HitOccupancy_Minus_.push_back(iB.book1D(histoname.str().c_str(), "HitOccupancy_Minus", 501, -0.5, 500.5));
0527
0528 histoname.str("");
0529 histoname << "EtaPhi_Plus_"
0530 << "layer_" << istr1;
0531 EtaPhi_Plus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, 1.45, 3.0, 72, -CLHEP::pi, CLHEP::pi));
0532 histoname.str("");
0533 histoname << "EtaPhi_Minus_"
0534 << "layer_" << istr1;
0535 EtaPhi_Minus_.push_back(
0536 iB.book2D(histoname.str().c_str(), "Occupancy", 31, -3.0, -1.45, 72, -CLHEP::pi, CLHEP::pi));
0537
0538 for (unsigned int itimeslice = 0; itimeslice < nTimes_; itimeslice++) {
0539 histoname.str("");
0540 histoname << "energy_time_" << itimeslice << "_layer_" << istr1;
0541 energy_[itimeslice].push_back(iB.book1D(histoname.str().c_str(), "energy_", 100, 0, 0.1));
0542 }
0543
0544
0545 if (nameDetector_ == "HGCalEESensitive" or nameDetector_ == "HGCalHESiliconSensitive") {
0546 histoname.str("");
0547 histoname << "energy_FullWafer_Fine_layer_" << istr1;
0548 TH1F* hEdepFWF = createHisto(histoname.str(), 100, 0., 400., false);
0549 histoSetting(hEdepFWF, "Eloss (keV)", "", kRed, kRed, 2);
0550 energyFWF_.push_back(iB.book1D(histoname.str().c_str(), hEdepFWF));
0551 hEdepFWF->Delete();
0552
0553 histoname.str("");
0554 histoname << "energy_FullWafer_CoarseThin_layer_" << istr1;
0555 TH1F* hEdepFWCN = createHisto(histoname.str(), 100, 0., 400., false);
0556 histoSetting(hEdepFWCN, "Eloss (keV)", "", kGreen + 1, kGreen + 1, 2);
0557 energyFWCN_.push_back(iB.book1D(histoname.str().c_str(), hEdepFWCN));
0558 hEdepFWCN->Delete();
0559
0560 histoname.str("");
0561 histoname << "energy_FullWafer_CoarseThick_layer_" << istr1;
0562 TH1F* hEdepFWCK = createHisto(histoname.str(), 100, 0., 400., false);
0563 histoSetting(hEdepFWCK, "Eloss (keV)", "", kMagenta, kMagenta, 2);
0564 energyFWCK_.push_back(iB.book1D(histoname.str().c_str(), hEdepFWCK));
0565 hEdepFWCK->Delete();
0566 }
0567
0568
0569
0570
0571 if (nameDetector_ == "HGCalEESensitive" or nameDetector_ == "HGCalHESiliconSensitive") {
0572 histoname.str("");
0573 histoname << "energy_PartialWafer_Fine_layer_" << istr1;
0574 TH1F* hEdepPWF = createHisto(histoname.str(), 100, 0., 400., false);
0575 histoSetting(hEdepPWF, "Eloss (keV)", "", kRed, kRed, 2);
0576 energyPWF_.push_back(iB.book1D(histoname.str().c_str(), hEdepPWF));
0577 hEdepPWF->Delete();
0578
0579 histoname.str("");
0580 histoname << "energy_PartialWafer_CoarseThin_layer_" << istr1;
0581 TH1F* hEdepPWCN = createHisto(histoname.str(), 100, 0., 400., false);
0582 histoSetting(hEdepPWCN, "Eloss (keV)", "", kGreen + 1, kGreen + 1, 2);
0583 energyPWCN_.push_back(iB.book1D(histoname.str().c_str(), hEdepPWCN));
0584 hEdepPWCN->Delete();
0585
0586 histoname.str("");
0587 histoname << "energy_PartialWafer_CoarseThick_layer_" << istr1;
0588 TH1F* hEdepPWCK = createHisto(histoname.str(), 100, 0., 400., false);
0589 histoSetting(hEdepPWCK, "Eloss (keV)", "", kMagenta, kMagenta, 2);
0590 energyPWCK_.push_back(iB.book1D(histoname.str().c_str(), hEdepPWCK));
0591 hEdepPWCK->Delete();
0592 }
0593
0594
0595
0596 if (nameDetector_ == "HGCalEESensitive" or nameDetector_ == "HGCalHESiliconSensitive") {
0597 histoname.str("");
0598 histoname << "hitXY_FullWafer_Fine_layer_" << istr1;
0599 TH2F* hitXYFWF = new TH2F(
0600 Form("hitXYFWF_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0601 histoSetting(hitXYFWF, "x (cm)", "y (cm)", kRed, kRed);
0602 hitXYFWF_.push_back(iB.book2D(histoname.str().c_str(), hitXYFWF));
0603 hitXYFWF->Delete();
0604
0605 histoname.str("");
0606 histoname << "hitXY_FullWafer_CoarseThin_layer_" << istr1;
0607 TH2F* hitXYFWCN = new TH2F(
0608 Form("hitXYFWCN_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0609 histoSetting(hitXYFWCN, "x (cm)", "y (cm)", kGreen + 1, kGreen + 1);
0610 hitXYFWCN_.push_back(iB.book2D(histoname.str().c_str(), hitXYFWCN));
0611 hitXYFWCN->Delete();
0612
0613 histoname.str("");
0614 histoname << "hitXY_FullWafer_CoarseThick_layer_" << istr1;
0615 TH2F* hitXYFWCK = new TH2F(
0616 Form("hitXYFWCK_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0617 histoSetting(hitXYFWCK, "x (cm)", "y (cm)", kMagenta, kMagenta);
0618 hitXYFWCK_.push_back(iB.book2D(histoname.str().c_str(), hitXYFWCK));
0619 hitXYFWCK->Delete();
0620 }
0621
0622 if (nameDetector_ == "HGCalHEScintillatorSensitive") {
0623 histoname.str("");
0624 histoname << "hitXY_Scintillator_layer_" << istr1;
0625 TH2F* hitXYB = new TH2F(
0626 Form("hitXYB_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0627 histoSetting(hitXYB, "x (cm)", "y (cm)", kBlue, kBlue);
0628 hitXYB_.push_back(iB.book2D(histoname.str().c_str(), hitXYB));
0629 hitXYB->Delete();
0630 }
0631
0632 }
0633 MeanHitOccupancy_Plus_ = iB.book1D("MeanHitOccupancy_Plus", "MeanHitOccupancy_Plus", layers_, 0.5, layers_ + 0.5);
0634 MeanHitOccupancy_Minus_ = iB.book1D("MeanHitOccupancy_Minus", "MeanHitOccupancy_Minus", layers_, 0.5, layers_ + 0.5);
0635 }
0636
0637 TH1F* HGCalSimHitValidation::createHisto(
0638 std::string histname, const int nbins, float minIndexX, float maxIndexX, bool isLogX) {
0639 TH1F* hist = nullptr;
0640 if (isLogX) {
0641 Double_t xbins[nbins + 1];
0642 double dx = (maxIndexX - minIndexX) / nbins;
0643 for (int i = 0; i <= nbins; i++) {
0644 xbins[i] = TMath::Power(10, (minIndexX + i * dx));
0645 }
0646 hist = new TH1F(Form("hEdep_%s", histname.c_str()), histname.c_str(), nbins, xbins);
0647 } else {
0648 hist = new TH1F(Form("hEdep_%s", histname.c_str()), histname.c_str(), nbins, minIndexX, maxIndexX);
0649 }
0650 return hist;
0651 }
0652
0653 void HGCalSimHitValidation::histoSetting(
0654 TH1F*& histo, const char* xTitle, const char* yTitle, Color_t lineColor, Color_t markerColor, int lineWidth) {
0655 histo->SetStats();
0656 histo->SetLineColor(lineColor);
0657 histo->SetLineWidth(lineWidth);
0658 histo->SetMarkerColor(markerColor);
0659 histo->GetXaxis()->SetTitle(xTitle);
0660 histo->GetYaxis()->SetTitle(yTitle);
0661 }
0662
0663 void HGCalSimHitValidation::histoSetting(
0664 TH2F*& histo, const char* xTitle, const char* yTitle, Color_t lineColor, Color_t markerColor, int lineWidth) {
0665 histo->SetStats();
0666 histo->SetLineColor(lineColor);
0667 histo->SetLineWidth(lineWidth);
0668 histo->SetMarkerColor(markerColor);
0669 histo->SetMarkerStyle(kFullCircle);
0670 histo->SetMarkerSize(0.6);
0671 histo->GetXaxis()->SetTitle(xTitle);
0672 histo->GetYaxis()->SetTitle(yTitle);
0673 }
0674 #include "FWCore/Framework/interface/MakerMacros.h"
0675
0676 DEFINE_FWK_MODULE(HGCalSimHitValidation);
