Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:06:36 

0001 // -*- C++ -*-
0002 //
0003 // Package:    OccupancyPlots
0004 // Class:      OccupancyPlots
0005 //
0006 /**\class OccupancyPlots OccupancyPlots.cc myTKAnalyses/DigiInvestigator/src/OccupancyPlots.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Andrea Venturi
0015 //         Created:  Mon Oct 27 17:37:53 CET 2008
0016 // $Id: OccupancyPlots.cc,v 1.3 2013/02/27 19:49:47 wmtan Exp $
0017 //
0018 //
0019 
0020 // system include files
0021 #include <memory>
0022 
0023 // user include files
0024 
0025 #include <vector>
0026 #include <map>
0027 #include <limits>
0028 
0029 #include "TProfile.h"
0030 
0031 #include "FWCore/Framework/interface/Frameworkfwd.h"
0032 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0033 #include "FWCore/Framework/interface/Event.h"
0034 #include "FWCore/Framework/interface/Run.h"
0035 #include "FWCore/Framework/interface/MakerMacros.h"
0036 #include "FWCore/Framework/interface/ESHandle.h"
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0039 #include "FWCore/Utilities/interface/InputTag.h"
0040 #include "FWCore/Utilities/interface/transform.h"
0041 
0042 #include "DPGAnalysis/SiStripTools/interface/RunHistogramManager.h"
0043 #include "CommonTools/UtilAlgos/interface/DetIdSelector.h"
0044 #include "CalibTracker/SiPixelESProducers/interface/SiPixelDetInfoFileReader.h"
0045 #include "FWCore/ParameterSet/interface/FileInPath.h"
0046 
0047 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0048 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0049 #include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
0050 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0051 
0052 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0053 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0054 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0055 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0056 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0057 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0058 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0059 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0060 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0061 
0062 //
0063 // class decleration
0064 //
0065 
0066 class OccupancyPlots : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0067 public:
0068   explicit OccupancyPlots(const edm::ParameterSet&);
0069   ~OccupancyPlots() override;
0070 
0071 private:
0072   void beginJob() override;
0073   void analyze(const edm::Event&, const edm::EventSetup&) override;
0074   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0075   void endRun(const edm::Run&, const edm::EventSetup&) override;
0076   void endJob() override;
0077 
0078   // ----------member data ---------------------------
0079 
0080   std::vector<edm::EDGetTokenT<std::map<unsigned int, int> > > m_multiplicityMapTokens;
0081   std::vector<edm::EDGetTokenT<std::map<unsigned int, int> > > m_occupancyMapTokens;
0082   edm::FileInPath m_fp;
0083   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> m_tkGeomToken;
0084   edm::ESGetToken<SiStripQuality, SiStripQualityRcd> m_stripQualityToken;
0085   edm::ESGetToken<SiPixelQuality, SiPixelQualityRcd> m_pixelQualityToken;
0086 
0087   RunHistogramManager m_rhm;
0088   std::map<unsigned int, DetIdSelector> m_wantedsubdets;
0089 
0090   TProfile** m_avemultiplicity;
0091   TProfile** m_aveoccupancy;
0092 
0093   TH1F** m_nchannels_ideal;
0094   TH1F** m_nchannels_real;
0095 
0096   TProfile** m_averadius;
0097   TProfile** m_avez;
0098   TProfile** m_avex;
0099   TProfile** m_avey;
0100   TProfile** m_zavedr;
0101   TProfile** m_zavedz;
0102   TProfile** m_zavedrphi;
0103   TProfile** m_yavedr;
0104   TProfile** m_yavedz;
0105   TProfile** m_yavedrphi;
0106   TProfile** m_xavedr;
0107   TProfile** m_xavedz;
0108   TProfile** m_xavedrphi;
0109 };
0110 
0111 //
0112 // constants, enums and typedefs
0113 //
0114 
0115 //
0116 // static data member definitions
0117 //
0118 
0119 //
0120 // constructors and destructor
0121 //
0122 OccupancyPlots::OccupancyPlots(const edm::ParameterSet& iConfig)
0123     : m_multiplicityMapTokens(edm::vector_transform(
0124           iConfig.getParameter<std::vector<edm::InputTag> >("multiplicityMaps"),
0125           [this](edm::InputTag const& tag) { return consumes<std::map<unsigned int, int> >(tag); })),
0126       m_occupancyMapTokens(edm::vector_transform(
0127           iConfig.getParameter<std::vector<edm::InputTag> >("occupancyMaps"),
0128           [this](edm::InputTag const& tag) { return consumes<std::map<unsigned int, int> >(tag); })),
0129       m_fp(iConfig.getUntrackedParameter<edm::FileInPath>(
0130           "file", edm::FileInPath("CalibTracker/SiPixelESProducers/data/PixelSkimmedGeometry.txt"))),
0131       m_tkGeomToken(esConsumes<edm::Transition::EndRun>()),
0132       m_stripQualityToken(esConsumes<edm::Transition::EndRun>()),
0133       m_pixelQualityToken(esConsumes<edm::Transition::EndRun>()),
0134       m_rhm(consumesCollector()),
0135       m_wantedsubdets() {
0136   //now do what ever initialization is needed
0137 
0138   m_avemultiplicity = m_rhm.makeTProfile("avemult", "Average Multiplicty", 6000, 0.5, 6000.5);
0139   m_aveoccupancy = m_rhm.makeTProfile("aveoccu", "Average Occupancy", 6000, 0.5, 6000.5);
0140 
0141   m_nchannels_ideal = m_rhm.makeTH1F("nchannels_ideal", "Number of channels (ideal)", 6000, 0.5, 6000.5);
0142   m_nchannels_real = m_rhm.makeTH1F("nchannels_real", "Number of channels (real)", 6000, 0.5, 6000.5);
0143 
0144   m_averadius = m_rhm.makeTProfile("averadius", "Average Module Radius", 6000, 0.5, 6000.5);
0145   m_avez = m_rhm.makeTProfile("avez", "Average Module z coordinate", 6000, 0.5, 6000.5);
0146   m_avex = m_rhm.makeTProfile("avex", "Average Module x coordinate", 6000, 0.5, 6000.5);
0147   m_avey = m_rhm.makeTProfile("avey", "Average Module y coordinate", 6000, 0.5, 6000.5);
0148 
0149   m_zavedr = m_rhm.makeTProfile("zavedr", "Average z unit vector dr", 6000, 0.5, 6000.5);
0150   m_zavedz = m_rhm.makeTProfile("zavedz", "Average z unit vector dz", 6000, 0.5, 6000.5);
0151   m_zavedrphi = m_rhm.makeTProfile("zavedrphi", "Average z unit vector drphi", 6000, 0.5, 6000.5);
0152   m_xavedr = m_rhm.makeTProfile("xavedr", "Average x unit vector dr", 6000, 0.5, 6000.5);
0153   m_xavedz = m_rhm.makeTProfile("xavedz", "Average x unit vctor dz", 6000, 0.5, 6000.5);
0154   m_xavedrphi = m_rhm.makeTProfile("xavedrphi", "Average Module x unit vector drphi", 6000, 0.5, 6000.5);
0155   m_yavedr = m_rhm.makeTProfile("yavedr", "Average y unit vector dr", 6000, 0.5, 6000.5);
0156   m_yavedz = m_rhm.makeTProfile("yavedz", "Average y unit vector dz", 6000, 0.5, 6000.5);
0157   m_yavedrphi = m_rhm.makeTProfile("yavedrphi", "Average y unit vector drphi", 6000, 0.5, 6000.5);
0158 
0159   std::vector<edm::ParameterSet> wantedsubdets_ps =
0160       iConfig.getParameter<std::vector<edm::ParameterSet> >("wantedSubDets");
0161 
0162   for (std::vector<edm::ParameterSet>::const_iterator wsdps = wantedsubdets_ps.begin(); wsdps != wantedsubdets_ps.end();
0163        ++wsdps) {
0164     unsigned int detsel = wsdps->getParameter<unsigned int>("detSelection");
0165     std::vector<std::string> selstr = wsdps->getUntrackedParameter<std::vector<std::string> >("selection");
0166     m_wantedsubdets[detsel] = DetIdSelector(selstr);
0167   }
0168 }
0169 
0170 OccupancyPlots::~OccupancyPlots() {
0171   // do anything here that needs to be done at desctruction time
0172   // (e.g. close files, deallocate resources etc.)
0173 }
0174 
0175 //
0176 // member functions
0177 //
0178 
0179 // ------------ method called to for each event  ------------
0180 void OccupancyPlots::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0181   using namespace edm;
0182 
0183   for (std::vector<edm::EDGetTokenT<std::map<unsigned int, int> > >::const_iterator mapToken =
0184            m_multiplicityMapTokens.begin();
0185        mapToken != m_multiplicityMapTokens.end();
0186        ++mapToken) {
0187     Handle<std::map<unsigned int, int> > mults;
0188     iEvent.getByToken(*mapToken, mults);
0189 
0190     for (std::map<unsigned int, int>::const_iterator mult = mults->begin(); mult != mults->end(); mult++) {
0191       if (m_avemultiplicity && *m_avemultiplicity)
0192         (*m_avemultiplicity)->Fill(mult->first, mult->second);
0193     }
0194   }
0195 
0196   for (std::vector<edm::EDGetTokenT<std::map<unsigned int, int> > >::const_iterator mapToken =
0197            m_occupancyMapTokens.begin();
0198        mapToken != m_occupancyMapTokens.end();
0199        ++mapToken) {
0200     Handle<std::map<unsigned int, int> > occus;
0201     iEvent.getByToken(*mapToken, occus);
0202 
0203     for (std::map<unsigned int, int>::const_iterator occu = occus->begin(); occu != occus->end(); occu++) {
0204       if (m_aveoccupancy && *m_aveoccupancy)
0205         (*m_aveoccupancy)->Fill(occu->first, occu->second);
0206     }
0207   }
0208 }
0209 
0210 // ------------ method called once each job just before starting event loop  ------------
0211 void OccupancyPlots::beginJob() {}
0212 
0213 void OccupancyPlots::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) { m_rhm.beginRun(iRun); }
0214 
0215 void OccupancyPlots::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0216   const auto& trkgeo = iSetup.getData(m_tkGeomToken);
0217 
0218   // Test new TrackerGeometry features
0219   LogDebug("IsThereTest") << "Test of TrackerGeometry::isThere";
0220   LogTrace("IsThereTest") << " is there PixelBarrel: " << trkgeo.isThere(GeomDetEnumerators::PixelBarrel);
0221   LogTrace("IsThereTest") << " is there PixelEndcap: " << trkgeo.isThere(GeomDetEnumerators::PixelEndcap);
0222   LogTrace("IsThereTest") << " is there P1PXB: " << trkgeo.isThere(GeomDetEnumerators::P1PXB);
0223   LogTrace("IsThereTest") << " is there P1PXEC: " << trkgeo.isThere(GeomDetEnumerators::P1PXEC);
0224   LogTrace("IsThereTest") << " is there P2PXB: " << trkgeo.isThere(GeomDetEnumerators::P2PXB);
0225   LogTrace("IsThereTest") << " is there P2PXEC: " << trkgeo.isThere(GeomDetEnumerators::P2PXEC);
0226   LogTrace("IsThereTest") << " is there TIB: " << trkgeo.isThere(GeomDetEnumerators::TIB);
0227   LogTrace("IsThereTest") << " is there TID: " << trkgeo.isThere(GeomDetEnumerators::TID);
0228   LogTrace("IsThereTest") << " is there TOB: " << trkgeo.isThere(GeomDetEnumerators::TOB);
0229   LogTrace("IsThereTest") << " is there TEC: " << trkgeo.isThere(GeomDetEnumerators::TEC);
0230   LogTrace("IsThereTest") << " is there P2OTB: " << trkgeo.isThere(GeomDetEnumerators::P2OTB);
0231   LogTrace("IsThereTest") << " is there P2OTEC: " << trkgeo.isThere(GeomDetEnumerators::P2OTEC);
0232 
0233   const Local2DPoint center(0., 0.);
0234   const Local3DPoint locz(0., 0., 1.);
0235   const Local3DPoint locx(1., 0., 0.);
0236   const Local3DPoint locy(0., 1., 0.);
0237   const GlobalPoint origin(0., 0., 0.);
0238 
0239   TrackingGeometry::DetIdContainer detunits = trkgeo.detUnitIds();
0240 
0241   for (TrackingGeometry::DetIdContainer::const_iterator det = detunits.begin(); det != detunits.end(); ++det) {
0242     if (det->det() != DetId::Tracker)
0243       continue;
0244 
0245     edm::LogInfo("DetIdFromGeometry") << det->rawId();
0246 
0247     GlobalPoint position = trkgeo.idToDet(*det)->toGlobal(center);
0248     GlobalPoint zpos = trkgeo.idToDet(*det)->toGlobal(locz);
0249     GlobalPoint xpos = trkgeo.idToDet(*det)->toGlobal(locx);
0250     GlobalPoint ypos = trkgeo.idToDet(*det)->toGlobal(locy);
0251     GlobalVector posvect = position - origin;
0252     GlobalVector dz = zpos - position;
0253     GlobalVector dx = xpos - position;
0254     GlobalVector dy = ypos - position;
0255 
0256     double dzdr = posvect.perp() > 0 ? (dz.x() * posvect.x() + dz.y() * posvect.y()) / posvect.perp() : 0.;
0257     double dxdr = posvect.perp() > 0 ? (dx.x() * posvect.x() + dx.y() * posvect.y()) / posvect.perp() : 0.;
0258     double dydr = posvect.perp() > 0 ? (dy.x() * posvect.x() + dy.y() * posvect.y()) / posvect.perp() : 0.;
0259 
0260     double dzdrphi = posvect.perp() > 0 ? (dz.y() * posvect.x() - dz.x() * posvect.y()) / posvect.perp() : 0.;
0261     double dxdrphi = posvect.perp() > 0 ? (dx.y() * posvect.x() - dx.x() * posvect.y()) / posvect.perp() : 0.;
0262     double dydrphi = posvect.perp() > 0 ? (dy.y() * posvect.x() - dy.x() * posvect.y()) / posvect.perp() : 0.;
0263 
0264     for (std::map<unsigned int, DetIdSelector>::const_iterator sel = m_wantedsubdets.begin();
0265          sel != m_wantedsubdets.end();
0266          ++sel) {
0267       if (sel->second.isSelected(*det)) {
0268         edm::LogInfo("SelectedDetId") << sel->first;
0269         // average positions
0270         if (m_averadius && *m_averadius)
0271           (*m_averadius)->Fill(sel->first, position.perp());
0272         if (m_avez && *m_avez)
0273           (*m_avez)->Fill(sel->first, position.z());
0274         if (m_avex && *m_avex)
0275           (*m_avex)->Fill(sel->first, position.x());
0276         if (m_avey && *m_avey)
0277           (*m_avey)->Fill(sel->first, position.y());
0278         if (m_zavedr && *m_zavedr)
0279           (*m_zavedr)->Fill(sel->first, dzdr);
0280         if (m_zavedz && *m_zavedz)
0281           (*m_zavedz)->Fill(sel->first, dz.z());
0282         if (m_zavedrphi && *m_zavedrphi)
0283           (*m_zavedrphi)->Fill(sel->first, dzdrphi);
0284         if (m_xavedr && *m_xavedr)
0285           (*m_xavedr)->Fill(sel->first, dxdr);
0286         if (m_xavedz && *m_xavedz)
0287           (*m_xavedz)->Fill(sel->first, dx.z());
0288         if (m_xavedrphi && *m_xavedrphi)
0289           (*m_xavedrphi)->Fill(sel->first, dxdrphi);
0290         if (m_yavedr && *m_yavedr)
0291           (*m_yavedr)->Fill(sel->first, dydr);
0292         if (m_yavedz && *m_yavedz)
0293           (*m_yavedz)->Fill(sel->first, dy.z());
0294         if (m_yavedrphi && *m_yavedrphi)
0295           (*m_yavedrphi)->Fill(sel->first, dydrphi);
0296       }
0297     }
0298   }
0299 
0300   // counting the number of channels per module subset
0301 
0302   // the histograms have to be reset to avoid double counting if endRun is called more than once
0303 
0304   if (m_nchannels_ideal && *m_nchannels_ideal)
0305     (*m_nchannels_ideal)->Reset();
0306   if (m_nchannels_real && *m_nchannels_real)
0307     (*m_nchannels_real)->Reset();
0308 
0309   const auto& stripQuality = iSetup.getData(m_stripQualityToken);
0310 
0311   for (const auto det : trkgeo.detUnits()) {
0312     const StripGeomDetUnit* stripDet = dynamic_cast<const StripGeomDetUnit*>(det);
0313     if (stripDet != nullptr) {
0314       const DetId detid = stripDet->geographicalId();
0315 
0316       int nchannideal = stripDet->specificTopology().nstrips();
0317       //     int nchannreal = stripDet->specificTopology().nstrips();
0318       int nchannreal = 0;
0319       for (int strip = 0; strip < nchannideal; ++strip) {
0320         if (!stripQuality.IsStripBad(detid, strip))
0321           ++nchannreal;
0322       }
0323 
0324       for (std::map<unsigned int, DetIdSelector>::const_iterator sel = m_wantedsubdets.begin();
0325            sel != m_wantedsubdets.end();
0326            ++sel) {
0327         if (sel->second.isSelected(detid)) {
0328           if (m_nchannels_ideal && *m_nchannels_ideal)
0329             (*m_nchannels_ideal)->Fill(sel->first, nchannideal);
0330           if (m_nchannels_real && *m_nchannels_real)
0331             (*m_nchannels_real)->Fill(sel->first, nchannreal);
0332         }
0333       }
0334     }
0335   }
0336 
0337   const auto& pxlquality = iSetup.getData(m_pixelQualityToken);
0338 
0339   SiPixelDetInfoFileReader pxlreader(m_fp.fullPath());
0340 
0341   const std::vector<uint32_t>& pxldetids = pxlreader.getAllDetIds();
0342 
0343   for (std::vector<uint32_t>::const_iterator detid = pxldetids.begin(); detid != pxldetids.end(); ++detid) {
0344     int nchannideal = pxlreader.getDetUnitDimensions(*detid).first * pxlreader.getDetUnitDimensions(*detid).second;
0345     int nchannreal = 0;
0346     if (!pxlquality.IsModuleBad(*detid)) {
0347       nchannreal = pxlreader.getDetUnitDimensions(*detid).first * pxlreader.getDetUnitDimensions(*detid).second;
0348     }
0349     /*
0350      int nchannreal = 0;
0351      for(int strip = 0; strip < nchannideal; ++strip) {
0352        if(!stripquality.IsStripBad(*detid,strip)) ++nchannreal;
0353      }
0354      */
0355 
0356     for (std::map<unsigned int, DetIdSelector>::const_iterator sel = m_wantedsubdets.begin();
0357          sel != m_wantedsubdets.end();
0358          ++sel) {
0359       if (sel->second.isSelected(*detid)) {
0360         if (m_nchannels_ideal && *m_nchannels_ideal)
0361           (*m_nchannels_ideal)->Fill(sel->first, nchannideal);
0362         if (m_nchannels_real && *m_nchannels_real)
0363           (*m_nchannels_real)->Fill(sel->first, nchannreal);
0364       }
0365     }
0366   }
0367 }
0368 // ------------ method called once each job just after ending the event loop  ------------
0369 void OccupancyPlots::endJob() {}
0370 //define this as a plug-in
0371 DEFINE_FWK_MODULE(OccupancyPlots);
[ 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!