Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:39 

0001 // system include files
0002 #include <cmath>
0003 #include <fstream>
0004 #include <iostream>
0005 #include <map>
0006 #include <string>
0007 #include <vector>
0008 
0009 #include "TH2D.h"
0010 #include "TH1D.h"
0011 
0012 #include "DataFormats/DetId/interface/DetId.h"
0013 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0014 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0015 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0016 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0017 
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0028 #include "FWCore/Utilities/interface/InputTag.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0031 
0032 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0033 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0034 
0035 class HGCalDigiStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0036 public:
0037   explicit HGCalDigiStudy(const edm::ParameterSet&);
0038   ~HGCalDigiStudy() override = default;
0039 
0040   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0041 
0042 protected:
0043   void beginJob() override;
0044   void endJob() override {}
0045   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0046   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0047   void analyze(edm::Event const&, edm::EventSetup const&) override;
0048 
0049 private:
0050   struct digiInfo {
0051     digiInfo() {
0052       phi = eta = r = z = charge = 0.0;
0053       layer = adc = 0;
0054     }
0055     double phi, eta, r, z, charge;
0056     int layer, adc;
0057   };
0058 
0059   template <class T1, class T2>
0060   void digiValidation(const T1& detId, const T2* geom, int layer, uint16_t adc, double charge);
0061 
0062   // ----------member data ---------------------------
0063   const std::string nameDetector_;
0064   const edm::InputTag source_;
0065   const bool ifNose_, ifLayer_;
0066   const int verbosity_, SampleIndx_, nbinR_, nbinZ_, nbinEta_, nLayers_;
0067   const double rmin_, rmax_, zmin_, zmax_, etamin_, etamax_;
0068   const edm::EDGetTokenT<HGCalDigiCollection> digiSource_;
0069   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> tok_hgcGeom_;
0070   const HGCalGeometry* hgcGeom_;
0071   int layers_, layerFront_, geomType_;
0072   TH1D *h_Charge_, *h_ADC_, *h_LayZp_, *h_LayZm_;
0073   TH2D *h_RZ_, *h_EtaPhi_;
0074   std::vector<TH2D*> h_XY_;
0075   TH1D *h_W1_, *h_W2_, *h_C1_, *h_C2_, *h_Ly_;
0076 };
0077 
0078 HGCalDigiStudy::HGCalDigiStudy(const edm::ParameterSet& iConfig)
0079     : nameDetector_(iConfig.getParameter<std::string>("detectorName")),
0080       source_(iConfig.getParameter<edm::InputTag>("digiSource")),
0081       ifNose_(iConfig.getUntrackedParameter<bool>("ifNose", false)),
0082       ifLayer_(iConfig.getUntrackedParameter<bool>("ifLayer", false)),
0083       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0084       SampleIndx_(iConfig.getUntrackedParameter<int>("sampleIndex", 5)),
0085       nbinR_(iConfig.getUntrackedParameter<int>("nBinR", 300)),
0086       nbinZ_(iConfig.getUntrackedParameter<int>("nBinZ", 300)),
0087       nbinEta_(iConfig.getUntrackedParameter<int>("nBinEta", 200)),
0088       nLayers_(iConfig.getUntrackedParameter<int>("layers", 28)),
0089       rmin_(iConfig.getUntrackedParameter<double>("rMin", 0.0)),
0090       rmax_(iConfig.getUntrackedParameter<double>("rMax", 300.0)),
0091       zmin_(iConfig.getUntrackedParameter<double>("zMin", 300.0)),
0092       zmax_(iConfig.getUntrackedParameter<double>("zMax", 600.0)),
0093       etamin_(iConfig.getUntrackedParameter<double>("etaMin", 1.0)),
0094       etamax_(iConfig.getUntrackedParameter<double>("etaMax", 3.0)),
0095       digiSource_(consumes<HGCalDigiCollection>(source_)),
0096       tok_hgcGeom_(esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
0097           edm::ESInputTag{"", nameDetector_})) {
0098   usesResource(TFileService::kSharedResource);
0099   edm::LogVerbatim("HGCalValidation") << "HGCalDigiStudy: request for Digi collection " << source_ << " for "
0100                                       << nameDetector_;
0101 }
0102 
0103 void HGCalDigiStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0104   edm::ParameterSetDescription desc;
0105   desc.add<std::string>("detectorName", "HGCalEESensitive");
0106   desc.add<edm::InputTag>("digiSource", edm::InputTag("simHGCalUnsuppressedDigis", "EE"));
0107   desc.addUntracked<bool>("ifNose", false);
0108   desc.addUntracked<bool>("ifLayer", false);
0109   desc.addUntracked<int>("verbosity", 0);
0110   desc.addUntracked<int>("sampleIndex", 0);
0111   desc.addUntracked<double>("rMin", 0.0);
0112   desc.addUntracked<double>("rMax", 300.0);
0113   desc.addUntracked<double>("zMin", 300.0);
0114   desc.addUntracked<double>("zMax", 600.0);
0115   desc.addUntracked<double>("etaMin", 1.0);
0116   desc.addUntracked<double>("etaMax", 3.0);
0117   desc.addUntracked<int>("nBinR", 300);
0118   desc.addUntracked<int>("nBinZ", 300);
0119   desc.addUntracked<int>("nBinEta", 200);
0120   desc.addUntracked<int>("layers", 28);
0121   descriptions.add("hgcalDigiStudyEE", desc);
0122 }
0123 
0124 void HGCalDigiStudy::beginJob() {
0125   edm::Service<TFileService> fs;
0126   std::ostringstream hname, title;
0127   hname.str("");
0128   title.str("");
0129   hname << "RZ_" << nameDetector_;
0130   title << "R vs Z for " << nameDetector_;
0131   h_RZ_ = fs->make<TH2D>(hname.str().c_str(), title.str().c_str(), nbinZ_, zmin_, zmax_, nbinR_, rmin_, rmax_);
0132   if (ifLayer_) {
0133     for (int ly = 0; ly < nLayers_; ++ly) {
0134       hname.str("");
0135       title.str("");
0136       hname << "XY_L" << (ly + 1);
0137       title << "Y vs X at Layer " << (ly + 1);
0138       h_XY_.emplace_back(
0139           fs->make<TH2D>(hname.str().c_str(), title.str().c_str(), nbinR_, -rmax_, rmax_, nbinR_, -rmax_, rmax_));
0140     }
0141   } else {
0142     hname.str("");
0143     title.str("");
0144     hname << "EtaPhi_" << nameDetector_;
0145     title << "#phi vs #eta for " << nameDetector_;
0146     h_EtaPhi_ = fs->make<TH2D>(hname.str().c_str(), title.str().c_str(), nbinEta_, etamin_, etamax_, 200, -M_PI, M_PI);
0147   }
0148   hname.str("");
0149   title.str("");
0150   hname << "Charge_" << nameDetector_;
0151   title << "Charge for " << nameDetector_;
0152   h_Charge_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 100, -25, 25);
0153   hname.str("");
0154   title.str("");
0155   hname << "ADC_" << nameDetector_;
0156   title << "ADC for " << nameDetector_;
0157   h_ADC_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 200, 0, 50);
0158   hname.str("");
0159   title.str("");
0160   hname << "LayerZp_" << nameDetector_;
0161   title << "Charge vs Layer (+z) for " << nameDetector_;
0162   h_LayZp_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 60, 0.0, 60.0);
0163   hname.str("");
0164   title.str("");
0165   hname << "LayerZm_" << nameDetector_;
0166   title << "Charge vs Layer (-z) for " << nameDetector_;
0167   h_LayZm_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 60, 0.0, 60.0);
0168   hname.str("");
0169   title.str("");
0170   hname << "LY_" << nameDetector_;
0171   title << "Layer number for " << nameDetector_;
0172   h_Ly_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 200, 0, 100);
0173   if (nameDetector_ == "HGCalHEScintillatorSensitive") {
0174     hname.str("");
0175     title.str("");
0176     hname << "IR_" << nameDetector_;
0177     title << "Radius index for " << nameDetector_;
0178     h_W1_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 200, -50, 50);
0179     hname.str("");
0180     title.str("");
0181     hname << "FI_" << nameDetector_;
0182     title << "#phi index for " << nameDetector_;
0183     h_C1_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 720, 0, 360);
0184   } else {
0185     hname.str("");
0186     title.str("");
0187     hname << "WU_" << nameDetector_;
0188     title << "u index of wafers for " << nameDetector_;
0189     h_W1_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 200, -50, 50);
0190     hname.str("");
0191     title.str("");
0192     hname << "WV_" << nameDetector_;
0193     title << "v index of wafers for " << nameDetector_;
0194     h_W2_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 100, -50, 50);
0195     hname.str("");
0196     title.str("");
0197     hname << "CU_" << nameDetector_;
0198     title << "u index of cells for " << nameDetector_;
0199     h_C1_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 100, 0, 50);
0200     hname.str("");
0201     title.str("");
0202     hname << "CV_" << nameDetector_;
0203     title << "v index of cells for " << nameDetector_;
0204     h_C2_ = fs->make<TH1D>(hname.str().c_str(), title.str().c_str(), 100, 0, 50);
0205   }
0206 }
0207 
0208 void HGCalDigiStudy::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0209   edm::ESHandle<HGCalGeometry> geom = iSetup.getHandle(tok_hgcGeom_);
0210   if (!geom.isValid())
0211     edm::LogVerbatim("HGCalValidation") << "HGCalDigiStudy: Cannot get "
0212                                         << "valid Geometry Object for " << nameDetector_;
0213   hgcGeom_ = geom.product();
0214   layerFront_ = hgcGeom_->topology().dddConstants().firstLayer();
0215   layers_ = hgcGeom_->topology().dddConstants().layers(true);
0216   if (hgcGeom_->topology().waferHexagon8())
0217     geomType_ = 1;
0218   else
0219     geomType_ = 2;
0220   if (nameDetector_ == "HGCalHFNoseSensitive")
0221     geomType_ = 3;
0222   edm::LogVerbatim("HGCalValidation") << "HGCalDigiStudy: gets Geometry for " << nameDetector_ << " of type "
0223                                       << geomType_ << " with " << layers_ << " layers and front layer " << layerFront_;
0224 }
0225 
0226 void HGCalDigiStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0227   unsigned int ntot(0), nused(0);
0228 
0229   if ((nameDetector_ == "HGCalEESensitive") || (nameDetector_ == "HGCalHFNoseSensitive")) {
0230     //HGCalEE
0231     const edm::Handle<HGCalDigiCollection>& theHGCEEDigiContainer = iEvent.getHandle(digiSource_);
0232     if (theHGCEEDigiContainer.isValid()) {
0233       if (verbosity_ > 0)
0234         edm::LogVerbatim("HGCalValidation")
0235             << nameDetector_ << " with " << theHGCEEDigiContainer->size() << " element(s)";
0236       for (const auto& it : *(theHGCEEDigiContainer.product())) {
0237         ntot++;
0238         nused++;
0239         DetId detId = it.id();
0240         int layer =
0241             ((geomType_ == 1) ? HGCSiliconDetId(detId).layer()
0242                               : ((geomType_ == 2) ? HGCScintillatorDetId(detId).layer() : HFNoseDetId(detId).layer()));
0243         const HGCSample& hgcSample = it.sample(SampleIndx_);
0244         uint16_t adc = hgcSample.data();
0245         double charge = adc;
0246         //      uint16_t   gain      = hgcSample.toa();
0247         //      double     charge    = adc*gain;
0248         digiValidation(detId, hgcGeom_, layer, adc, charge);
0249         if (geomType_ == 1) {
0250           HGCSiliconDetId id = HGCSiliconDetId(detId);
0251           h_Ly_->Fill(id.layer());
0252           h_W1_->Fill(id.waferU());
0253           h_W2_->Fill(id.waferV());
0254           h_C1_->Fill(id.cellU());
0255           h_C2_->Fill(id.cellV());
0256         } else if (geomType_ == 2) {
0257           HGCScintillatorDetId id = HGCScintillatorDetId(detId);
0258           h_Ly_->Fill(id.layer());
0259           h_W1_->Fill(id.ieta());
0260           h_C1_->Fill(id.iphi());
0261         } else {
0262           HFNoseDetId id = HFNoseDetId(detId);
0263           h_Ly_->Fill(id.layer());
0264           h_W1_->Fill(id.waferU());
0265           h_W2_->Fill(id.waferV());
0266           h_C1_->Fill(id.cellU());
0267           h_C2_->Fill(id.cellV());
0268         }
0269       }
0270     } else {
0271       edm::LogVerbatim("HGCalValidation")
0272           << "DigiCollection handle " << source_ << " does not exist for " << nameDetector_ << " !!!";
0273     }
0274   } else if ((nameDetector_ == "HGCalHESiliconSensitive") || (nameDetector_ == "HGCalHEScintillatorSensitive")) {
0275     //HGCalHE
0276     const edm::Handle<HGCalDigiCollection>& theHGCHEDigiContainer = iEvent.getHandle(digiSource_);
0277     if (theHGCHEDigiContainer.isValid()) {
0278       if (verbosity_ > 0)
0279         edm::LogVerbatim("HGCalValidation")
0280             << nameDetector_ << " with " << theHGCHEDigiContainer->size() << " element(s)";
0281       for (const auto& it : *(theHGCHEDigiContainer.product())) {
0282         ntot++;
0283         nused++;
0284         DetId detId = it.id();
0285         int layer =
0286             ((geomType_ == 1) ? HGCSiliconDetId(detId).layer()
0287                               : ((geomType_ == 2) ? HGCScintillatorDetId(detId).layer() : HFNoseDetId(detId).layer()));
0288         const HGCSample& hgcSample = it.sample(SampleIndx_);
0289         uint16_t adc = hgcSample.data();
0290         double charge = adc;
0291         //      uint16_t   gain      = hgcSample.toa();
0292         //      double     charge    = adc*gain;
0293         digiValidation(detId, hgcGeom_, layer, adc, charge);
0294         if (geomType_ == 1) {
0295           HGCSiliconDetId id = HGCSiliconDetId(detId);
0296           h_Ly_->Fill(id.layer());
0297           h_W1_->Fill(id.waferU());
0298           h_W2_->Fill(id.waferV());
0299           h_C1_->Fill(id.cellU());
0300           h_C2_->Fill(id.cellV());
0301         } else if (geomType_ == 2) {
0302           HGCScintillatorDetId id = HGCScintillatorDetId(detId);
0303           h_Ly_->Fill(id.layer());
0304           h_W1_->Fill(id.ieta());
0305           h_C1_->Fill(id.iphi());
0306         } else {
0307           HFNoseDetId id = HFNoseDetId(detId);
0308           h_Ly_->Fill(id.layer());
0309           h_W1_->Fill(id.waferU());
0310           h_W2_->Fill(id.waferV());
0311           h_C1_->Fill(id.cellU());
0312           h_C2_->Fill(id.cellV());
0313         }
0314       }
0315     } else {
0316       edm::LogVerbatim("HGCalValidation")
0317           << "DigiCollection handle " << source_ << " does not exist for " << nameDetector_ << " !!!";
0318     }
0319   } else {
0320     edm::LogWarning("HGCalValidation") << "invalid detector name !! " << nameDetector_ << " source " << source_;
0321   }
0322   edm::LogVerbatim("HGCalValidation") << "Event " << iEvent.id().event() << ":" << nameDetector_ << " with " << ntot
0323                                       << " total and " << nused << " used digis";
0324 }
0325 
0326 template <class T1, class T2>
0327 void HGCalDigiStudy::digiValidation(const T1& detId, const T2* geom, int layer, uint16_t adc, double charge) {
0328   if (verbosity_ > 1)
0329     edm::LogVerbatim("HGCalValidation") << std::hex << detId.rawId() << std::dec << " adc = " << adc
0330                                         << " charge = " << charge;
0331 
0332   DetId id1 = DetId(detId.rawId());
0333   const GlobalPoint& gcoord = geom->getPosition(id1);
0334 
0335   digiInfo hinfo;
0336   hinfo.r = gcoord.perp();
0337   hinfo.z = gcoord.z();
0338   hinfo.eta = std::abs(gcoord.eta());
0339   hinfo.phi = gcoord.phi();
0340   hinfo.adc = adc;
0341   hinfo.charge = charge;
0342   hinfo.layer = layer + layerFront_;
0343   if (verbosity_ > 1)
0344     edm::LogVerbatim("HGCalValidation") << "R =  " << hinfo.r << " z = " << hinfo.z << " eta = " << hinfo.eta
0345                                         << " phi = " << hinfo.phi;
0346 
0347   h_Charge_->Fill(hinfo.charge);
0348   h_ADC_->Fill(hinfo.adc);
0349   h_RZ_->Fill(std::abs(hinfo.z), hinfo.r);
0350   if (ifLayer_) {
0351     if (layer <= static_cast<int>(h_XY_.size()))
0352       h_XY_[layer - 1]->Fill(gcoord.x(), gcoord.y());
0353   } else {
0354     h_EtaPhi_->Fill(hinfo.eta, hinfo.phi);
0355   }
0356   if (hinfo.z > 0)
0357     h_LayZp_->Fill(hinfo.layer, hinfo.charge);
0358   else
0359     h_LayZm_->Fill(hinfo.layer, hinfo.charge);
0360 }
0361 
0362 //define this as a plug-in
0363 DEFINE_FWK_MODULE(HGCalDigiStudy);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!