Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-21 22:40:24

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