Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-24 02:22:28

0001 // system include files
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   // ----------member data ---------------------------
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   //Generator input
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   //Now the hits
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;  //index 0 and 1 corresponds to 25 ns and 1000 ns, respectively. In addititon, chaging energy loss unit to keV.
0359 
0360   unsigned int ilayer = hinfo.layer;
0361   double x = hinfo.x * CLHEP::mm / CLHEP::cm;  // chaging length unit to cm.
0362   double y = hinfo.y * CLHEP::mm / CLHEP::cm;
0363   if (ilayer < layers_) {
0364     if (nameDetector_ == "HGCalEESensitive" or nameDetector_ == "HGCalHESiliconSensitive") {
0365       // Fill the energy loss histograms for MIP
0366       if (!TMath::AreEqualAbs(edep, 0.0, 1.e-5)) {  //to avoid peak at zero due Eloss less than 10 mili eV.
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       // Fill the XY distribution of detector hits
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     }  //is Silicon
0398     if (nameDetector_ == "HGCalHEScintillatorSensitive") {
0399       hitXYB_.at(ilayer)->Fill(x, y);
0400     }  //is Scintillator
0401   }    //layer condition
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 // ------------ method called when starting to processes a run  ------------
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     ///////////// Histograms for Energy loss in full wafers////////////
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     ///////////// Histograms for Energy loss in partial wafers////////////
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     // ///////////// Histograms for the XY distribution of fired cells/scintillator tiles ///////////////
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 //define this as a plug-in
0676 DEFINE_FWK_MODULE(HGCalSimHitValidation);