Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-05 03:16:19

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 "CommonTools/UtilAlgos/interface/TFileService.h"
0010 
0011 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0012 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0013 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0014 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0015 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0016 
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0024 #include "FWCore/ServiceRegistry/interface/Service.h"
0025 #include "FWCore/Utilities/interface/InputTag.h"
0026 
0027 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0028 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0029 
0030 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0031 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0032 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0033 
0034 #include "CLHEP/Geometry/Point3D.h"
0035 #include "CLHEP/Geometry/Transform3D.h"
0036 #include "CLHEP/Geometry/Vector3D.h"
0037 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0038 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0039 
0040 #include "TH1D.h"
0041 #include "TH2D.h"
0042 
0043 class HGCalSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0044 public:
0045   struct hitsinfo {
0046     hitsinfo() {
0047       phi = eta = energy = time = 0.0;
0048       layer = 0;
0049     }
0050     double phi, eta, energy, time;
0051     int layer;
0052   };
0053 
0054   explicit HGCalSimHitStudy(const edm::ParameterSet&);
0055   ~HGCalSimHitStudy() override {}
0056   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0057 
0058 protected:
0059   void beginJob() override;
0060   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0061   void analyze(edm::Event const&, edm::EventSetup const&) override;
0062   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0063 
0064 private:
0065   void analyzeHits(int, const std::string&, const std::vector<PCaloHit>&);
0066 
0067   // ----------member data ---------------------------
0068   const std::vector<std::string> nameDetectors_, caloHitSources_;
0069   const double rmin_, rmax_, zmin_, zmax_;
0070   const double etamin_, etamax_;
0071   const int nbinR_, nbinZ_, nbinEta_, nLayers_, verbosity_;
0072   const bool ifNose_, ifLayer_;
0073   std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> > tok_hgcGeom_;
0074   std::vector<const HGCalDDDConstants*> hgcons_;
0075   std::vector<edm::EDGetTokenT<edm::PCaloHitContainer> > tok_hits_;
0076   std::vector<int> layers_, layerFront_;
0077 
0078   //histogram related stuff
0079   std::vector<TH2D*> h_RZ_, h_EtaPhi_, h_EtFiZp_, h_EtFiZm_, h_XY_;
0080   std::vector<TH1D*> h_E_, h_T_, h_LayerZp_, h_LayerZm_;
0081   std::vector<TH1D*> h_W1_, h_W2_, h_C1_, h_C2_, h_Ly_;
0082 };
0083 
0084 HGCalSimHitStudy::HGCalSimHitStudy(const edm::ParameterSet& iConfig)
0085     : nameDetectors_(iConfig.getParameter<std::vector<std::string> >("detectorNames")),
0086       caloHitSources_(iConfig.getParameter<std::vector<std::string> >("caloHitSources")),
0087       rmin_(iConfig.getUntrackedParameter<double>("rMin", 0.0)),
0088       rmax_(iConfig.getUntrackedParameter<double>("rMax", 3000.0)),
0089       zmin_(iConfig.getUntrackedParameter<double>("zMin", 3000.0)),
0090       zmax_(iConfig.getUntrackedParameter<double>("zMax", 6000.0)),
0091       etamin_(iConfig.getUntrackedParameter<double>("etaMin", 1.0)),
0092       etamax_(iConfig.getUntrackedParameter<double>("etaMax", 3.0)),
0093       nbinR_(iConfig.getUntrackedParameter<int>("nBinR", 300)),
0094       nbinZ_(iConfig.getUntrackedParameter<int>("nBinZ", 300)),
0095       nbinEta_(iConfig.getUntrackedParameter<int>("nBinEta", 200)),
0096       nLayers_(iConfig.getUntrackedParameter<int>("layers", 50)),
0097       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0098       ifNose_(iConfig.getUntrackedParameter<bool>("ifNose", false)),
0099       ifLayer_(iConfig.getUntrackedParameter<bool>("ifLayer", false)) {
0100   usesResource(TFileService::kSharedResource);
0101 
0102   for (auto const& name : nameDetectors_)
0103     tok_hgcGeom_.emplace_back(
0104         esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name}));
0105   for (auto const& source : caloHitSources_)
0106     tok_hits_.emplace_back(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", source)));
0107 }
0108 
0109 void HGCalSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0110   edm::ParameterSetDescription desc;
0111   std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0112   std::vector<std::string> sources = {"HGCHitsEE", "HGCHitsHEfront", "HGCHitsHEback"};
0113   desc.add<std::vector<std::string> >("detectorNames", names);
0114   desc.add<std::vector<std::string> >("caloHitSources", sources);
0115   desc.addUntracked<double>("rMin", 0.0);
0116   desc.addUntracked<double>("rMax", 3000.0);
0117   desc.addUntracked<double>("zMin", 3000.0);
0118   desc.addUntracked<double>("zMax", 6000.0);
0119   desc.addUntracked<double>("etaMin", 1.0);
0120   desc.addUntracked<double>("etaMax", 3.0);
0121   desc.addUntracked<int>("nBinR", 300);
0122   desc.addUntracked<int>("nBinZ", 300);
0123   desc.addUntracked<int>("nBinEta", 200);
0124   desc.addUntracked<int>("layers", 50);
0125   desc.addUntracked<int>("verbosity", 0);
0126   desc.addUntracked<bool>("ifNose", false);
0127   desc.addUntracked<bool>("ifLayer", false);
0128   descriptions.add("hgcalSimHitStudy", desc);
0129 }
0130 
0131 void HGCalSimHitStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0132   //Now the hits
0133   for (unsigned int k = 0; k < tok_hits_.size(); ++k) {
0134     edm::Handle<edm::PCaloHitContainer> theCaloHitContainers;
0135     iEvent.getByToken(tok_hits_[k], theCaloHitContainers);
0136     if (theCaloHitContainers.isValid()) {
0137       if (verbosity_ > 0)
0138         edm::LogVerbatim("HGCalValidation") << " PcalohitItr = " << theCaloHitContainers->size();
0139       std::vector<PCaloHit> caloHits;
0140       caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0141       analyzeHits(k, nameDetectors_[k], caloHits);
0142     } else if (verbosity_ > 0) {
0143       edm::LogVerbatim("HGCalValidation") << "PCaloHitContainer does not "
0144                                           << "exist for " << nameDetectors_[k];
0145     }
0146   }
0147 }
0148 
0149 void HGCalSimHitStudy::analyzeHits(int ih, std::string const& name, std::vector<PCaloHit> const& hits) {
0150   if (verbosity_ > 0)
0151     edm::LogVerbatim("HGCalValidation") << name << " with " << hits.size() << " PcaloHit elements";
0152 
0153   std::map<uint32_t, hitsinfo> map_hits;
0154   map_hits.clear();
0155 
0156   unsigned int nused(0);
0157   for (auto const& hit : hits) {
0158     double energy = hit.energy();
0159     double time = hit.time();
0160     uint32_t id = hit.id();
0161     int cell, sector, sector2(0), layer, zside;
0162     int subdet(0), cell2(0), type(0);
0163     HepGeom::Point3D<float> gcoord;
0164     std::pair<float, float> xy;
0165     if (ifNose_) {
0166       HFNoseDetId detId = HFNoseDetId(id);
0167       subdet = detId.subdetId();
0168       cell = detId.cellU();
0169       cell2 = detId.cellV();
0170       sector = detId.waferU();
0171       sector2 = detId.waferV();
0172       type = detId.type();
0173       layer = detId.layer();
0174       zside = detId.zside();
0175       xy = hgcons_[ih]->locateCell(layer, sector, sector2, cell, cell2, false, true);
0176       h_W2_[ih]->Fill(sector2);
0177       h_C2_[ih]->Fill(cell2);
0178     } else if (hgcons_[ih]->waferHexagon8()) {
0179       HGCSiliconDetId detId = HGCSiliconDetId(id);
0180       subdet = static_cast<int>(detId.det());
0181       cell = detId.cellU();
0182       cell2 = detId.cellV();
0183       sector = detId.waferU();
0184       sector2 = detId.waferV();
0185       type = detId.type();
0186       layer = detId.layer();
0187       zside = detId.zside();
0188       xy = hgcons_[ih]->locateCell(layer, sector, sector2, cell, cell2, false, true);
0189       h_W2_[ih]->Fill(sector2);
0190       h_C2_[ih]->Fill(cell2);
0191     } else if (hgcons_[ih]->tileTrapezoid()) {
0192       HGCScintillatorDetId detId = HGCScintillatorDetId(id);
0193       subdet = static_cast<int>(detId.det());
0194       sector = detId.ieta();
0195       cell = detId.iphi();
0196       type = detId.type();
0197       layer = detId.layer();
0198       zside = detId.zside();
0199       xy = hgcons_[ih]->locateCellTrap(layer, sector, cell, false);
0200     } else {
0201       HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, type, cell);
0202       xy = hgcons_[ih]->locateCell(cell, layer, sector, false);
0203     }
0204     double zp = hgcons_[ih]->waferZ(layer, false);
0205     if (zside < 0)
0206       zp = -zp;
0207     double xp = (zp < 0) ? -xy.first : xy.first;
0208     gcoord = HepGeom::Point3D<float>(xp, xy.second, zp);
0209     if (verbosity_ > 2)
0210       edm::LogVerbatim("HGCalValidation")
0211           << "i/p " << subdet << ":" << zside << ":" << layer << ":" << sector << ":" << sector2 << ":" << cell << ":"
0212           << cell2 << " o/p " << xy.first << ":" << xy.second << ":" << zp;
0213     nused++;
0214     double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
0215     if (verbosity_ > 1)
0216       edm::LogVerbatim("HGCalValidation")
0217           << "Detector " << name << " zside = " << zside << " layer = " << layer << " type = " << type
0218           << " wafer = " << sector << ":" << sector2 << " cell = " << cell << ":" << cell2 << " positon = " << gcoord
0219           << " energy = " << energy << " time = " << time << ":" << tof;
0220     time -= tof;
0221     if (time < 0)
0222       time = 0;
0223     hitsinfo hinfo;
0224     if (map_hits.count(id) != 0) {
0225       hinfo = map_hits[id];
0226     } else {
0227       hinfo.layer = layer + layerFront_[ih];
0228       hinfo.phi = gcoord.getPhi();
0229       hinfo.eta = gcoord.getEta();
0230       hinfo.time = time;
0231     }
0232     hinfo.energy += energy;
0233     map_hits[id] = hinfo;
0234 
0235     //Fill in histograms
0236     h_RZ_[0]->Fill(std::abs(gcoord.z()), gcoord.rho());
0237     h_RZ_[ih + 1]->Fill(std::abs(gcoord.z()), gcoord.rho());
0238     if (ifLayer_) {
0239       if (hinfo.layer <= static_cast<int>(h_XY_.size()))
0240         h_XY_[hinfo.layer - 1]->Fill(gcoord.x(), gcoord.y());
0241     } else {
0242       h_EtaPhi_[0]->Fill(std::abs(hinfo.eta), hinfo.phi);
0243       h_EtaPhi_[ih + 1]->Fill(std::abs(hinfo.eta), hinfo.phi);
0244     }
0245     h_Ly_[ih]->Fill(layer);
0246     h_W1_[ih]->Fill(sector);
0247     h_C1_[ih]->Fill(cell);
0248   }
0249   if (verbosity_ > 0)
0250     edm::LogVerbatim("HGCalValidation") << name << " with " << map_hits.size() << ":" << nused << " detector elements"
0251                                         << " being hit";
0252 
0253   for (auto const& hit : map_hits) {
0254     hitsinfo hinfo = hit.second;
0255     if (verbosity_ > 1)
0256       edm::LogVerbatim("HGCalValidation")
0257           << " ----------------------   eta = " << hinfo.eta << " phi = " << hinfo.phi << " layer = " << hinfo.layer
0258           << " E = " << hinfo.energy << " T = " << hinfo.time;
0259     h_E_[0]->Fill(hinfo.energy);
0260     h_E_[ih + 1]->Fill(hinfo.energy);
0261     h_T_[0]->Fill(hinfo.time);
0262     h_T_[ih + 1]->Fill(hinfo.time);
0263     if (hinfo.eta > 0) {
0264       if (!ifLayer_) {
0265         h_EtFiZp_[0]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0266         h_EtFiZp_[ih + 1]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0267       }
0268       h_LayerZp_[0]->Fill(hinfo.layer, hinfo.energy);
0269       h_LayerZp_[ih + 1]->Fill(hinfo.layer, hinfo.energy);
0270     } else {
0271       if (!ifLayer_) {
0272         h_EtFiZm_[0]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0273         h_EtFiZm_[ih + 1]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0274       }
0275       h_LayerZm_[0]->Fill(hinfo.layer, hinfo.energy);
0276       h_LayerZm_[ih + 1]->Fill(hinfo.layer, hinfo.energy);
0277     }
0278   }
0279 }
0280 
0281 // ------------ method called when starting to processes a run  ------------
0282 void HGCalSimHitStudy::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0283   for (unsigned int k = 0; k < nameDetectors_.size(); ++k) {
0284     hgcons_.emplace_back(&iSetup.getData(tok_hgcGeom_[k]));
0285     layers_.emplace_back(hgcons_.back()->layers(false));
0286     layerFront_.emplace_back(hgcons_.back()->firstLayer());
0287     if (verbosity_ > 0)
0288       edm::LogVerbatim("HGCalValidation") << nameDetectors_[k] << " defined with " << layers_.back() << " Layers with "
0289                                           << layerFront_.back() << " in front";
0290   }
0291 }
0292 
0293 void HGCalSimHitStudy::beginJob() {
0294   edm::Service<TFileService> fs;
0295 
0296   std::ostringstream name, title;
0297   for (unsigned int ih = 0; ih <= nameDetectors_.size(); ++ih) {
0298     name.str("");
0299     title.str("");
0300     if (ih == 0) {
0301       name << "RZ_AllDetectors";
0302       title << "R vs Z for All Detectors";
0303     } else {
0304       name << "RZ_" << nameDetectors_[ih - 1];
0305       title << "R vs Z for " << nameDetectors_[ih - 1];
0306     }
0307     h_RZ_.emplace_back(
0308         fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinZ_, zmin_, zmax_, nbinR_, rmin_, rmax_));
0309     if (ifLayer_) {
0310       if (ih == 0) {
0311         for (int ly = 0; ly < nLayers_; ++ly) {
0312           name.str("");
0313           title.str("");
0314           name << "XY_L" << (ly + 1);
0315           title << "Y vs X at Layer " << (ly + 1);
0316           h_XY_.emplace_back(
0317               fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinR_, -rmax_, rmax_, nbinR_, -rmax_, rmax_));
0318         }
0319       }
0320     } else {
0321       name.str("");
0322       title.str("");
0323       if (ih == 0) {
0324         name << "EtaPhi_AllDetectors";
0325         title << "#phi vs #eta for All Detectors";
0326       } else {
0327         name << "EtaPhi_" << nameDetectors_[ih - 1];
0328         title << "#phi vs #eta for " << nameDetectors_[ih - 1];
0329       }
0330       h_EtaPhi_.emplace_back(
0331           fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinEta_, etamin_, etamax_, 200, -M_PI, M_PI));
0332       name.str("");
0333       title.str("");
0334       if (ih == 0) {
0335         name << "EtFiZp_AllDetectors";
0336         title << "#phi vs #eta (+z) for All Detectors";
0337       } else {
0338         name << "EtFiZp_" << nameDetectors_[ih - 1];
0339         title << "#phi vs #eta (+z) for " << nameDetectors_[ih - 1];
0340       }
0341       h_EtFiZp_.emplace_back(
0342           fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinEta_, etamin_, etamax_, 200, -M_PI, M_PI));
0343       name.str("");
0344       title.str("");
0345       if (ih == 0) {
0346         name << "EtFiZm_AllDetectors";
0347         title << "#phi vs #eta (-z) for All Detectors";
0348       } else {
0349         name << "EtFiZm_" << nameDetectors_[ih - 1];
0350         title << "#phi vs #eta (-z) for " << nameDetectors_[ih - 1];
0351       }
0352       h_EtFiZm_.emplace_back(
0353           fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinEta_, etamin_, etamax_, 200, -M_PI, M_PI));
0354     }
0355     name.str("");
0356     title.str("");
0357     if (ih == 0) {
0358       name << "LayerZp_AllDetectors";
0359       title << "Energy vs Layer (+z) for All Detectors";
0360     } else {
0361       name << "LayerZp_" << nameDetectors_[ih - 1];
0362       title << "Energy vs Layer (+z) for " << nameDetectors_[ih - 1];
0363     }
0364     h_LayerZp_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 60, 0.0, 60.0));
0365     name.str("");
0366     title.str("");
0367     if (ih == 0) {
0368       name << "LayerZm_AllDetectors";
0369       title << "Energy vs Layer (-z) for All Detectors";
0370     } else {
0371       name << "LayerZm_" << nameDetectors_[ih - 1];
0372       title << "Energy vs Layer (-z) for " << nameDetectors_[ih - 1];
0373     }
0374     h_LayerZm_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 60, 0.0, 60.0));
0375 
0376     name.str("");
0377     title.str("");
0378     if (ih == 0) {
0379       name << "E_AllDetectors";
0380       title << "Energy Deposit for All Detectors";
0381     } else {
0382       name << "E_" << nameDetectors_[ih - 1];
0383       title << "Energy Deposit for " << nameDetectors_[ih - 1];
0384     }
0385     h_E_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 1000, 0.0, 1.0));
0386 
0387     name.str("");
0388     title.str("");
0389     if (ih == 0) {
0390       name << "T_AllDetectors";
0391       title << "Time for All Detectors";
0392     } else {
0393       name << "T_" << nameDetectors_[ih - 1];
0394       title << "Time for " << nameDetectors_[ih - 1];
0395     }
0396     h_T_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 1000, 0.0, 200.0));
0397   }
0398 
0399   for (unsigned int ih = 0; ih < nameDetectors_.size(); ++ih) {
0400     name.str("");
0401     title.str("");
0402     name << "LY_" << nameDetectors_[ih];
0403     title << "Layer number for " << nameDetectors_[ih];
0404     h_Ly_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 200, 0, 100));
0405     if (nameDetectors_[ih] == "HGCalHEScintillatorSensitive") {
0406       name.str("");
0407       title.str("");
0408       name << "IR_" << nameDetectors_[ih];
0409       title << "Radius index for " << nameDetectors_[ih];
0410       h_W1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 200, -50, 50));
0411       name.str("");
0412       title.str("");
0413       name << "FI_" << nameDetectors_[ih];
0414       title << "#phi index for " << nameDetectors_[ih];
0415       h_C1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 720, 0, 360));
0416     } else {
0417       name.str("");
0418       title.str("");
0419       name << "WU_" << nameDetectors_[ih];
0420       title << "u index of wafers for " << nameDetectors_[ih];
0421       h_W1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 200, -50, 50));
0422       name.str("");
0423       title.str("");
0424       name << "WV_" << nameDetectors_[ih];
0425       title << "v index of wafers for " << nameDetectors_[ih];
0426       h_W2_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 100, -50, 50));
0427       name.str("");
0428       title.str("");
0429       name << "CU_" << nameDetectors_[ih];
0430       title << "u index of cells for " << nameDetectors_[ih];
0431       h_C1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 100, 0, 50));
0432       name.str("");
0433       title.str("");
0434       name << "CV_" << nameDetectors_[ih];
0435       title << "v index of cells for " << nameDetectors_[ih];
0436       h_C2_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 100, 0, 50));
0437     }
0438   }
0439 }
0440 
0441 #include "FWCore/Framework/interface/MakerMacros.h"
0442 //define this as a plug-in
0443 DEFINE_FWK_MODULE(HGCalSimHitStudy);