Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:56 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TBeamTest
0004 // Class:      TBeamTest
0005 //
0006 /**\class TBeamTest TBeamTest.cc 
0007 
0008  Description: Access Digi collection and fill a few histograms to compare with TestBeam data
0009 
0010 */
0011 //
0012 // Author:  Suchandra Dutta, Suvankar RoyChoudhury
0013 // Created:  July 2015
0014 //
0015 //
0016 // system include files
0017 #include <memory>
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "FWCore/Framework/interface/ESWatcher.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0026 #include "DataFormats/Common/interface/Handle.h"
0027 #include "DataFormats/DetId/interface/DetId.h"
0028 #include "DataFormats/Common/interface/DetSetVector.h"
0029 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0030 #include "DataFormats/Common/interface/DetSetVector.h"
0031 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
0032 #include "DataFormats/Math/interface/deltaPhi.h"
0033 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0034 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0035 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0036 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0037 
0038 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0039 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0040 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0041 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0042 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0043 
0044 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0045 
0046 // DQM Histograming
0047 #include "DQMServices/Core/interface/DQMStore.h"
0048 #include <cmath>
0049 
0050 using Phase2TrackerGeomDetUnit = PixelGeomDetUnit;
0051 
0052 class TBeamTest : public DQMEDAnalyzer {
0053 public:
0054   explicit TBeamTest(const edm::ParameterSet&);
0055   ~TBeamTest() override;
0056   void bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) override;
0057   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0058 
0059   struct DigiMEs {
0060     MonitorElement* NumberOfDigis;
0061     MonitorElement* PositionOfDigis;
0062     MonitorElement* NumberOfClusters;
0063     std::vector<MonitorElement*> ClusterWidths;
0064     MonitorElement* ClusterPosition;
0065   };
0066 
0067 private:
0068   int matchedSimTrackIndex(edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& linkHandle,
0069                            edm::Handle<edm::SimTrackContainer>& simTkHandle,
0070                            DetId detId,
0071                            unsigned int& channel);
0072 
0073   void fillClusterWidth(DigiMEs& mes, float dphi, float width);
0074   edm::ParameterSet config_;
0075   std::map<std::string, DigiMEs> detMEs;
0076   const edm::InputTag otDigiSrc_;
0077   const edm::InputTag digiSimLinkSrc_;
0078   const edm::InputTag simTrackSrc_;
0079   const std::string geomType_;
0080 
0081   const std::vector<double> phiValues;
0082   const edm::EDGetTokenT<edm::DetSetVector<Phase2TrackerDigi> > otDigiToken_;
0083   const edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > otDigiSimLinkToken_;
0084   const edm::EDGetTokenT<edm::SimTrackContainer> simTrackToken_;
0085   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0086   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0087 };
0088 //
0089 // constructors
0090 //
0091 TBeamTest::TBeamTest(const edm::ParameterSet& iConfig)
0092     : config_(iConfig),
0093       otDigiSrc_(iConfig.getParameter<edm::InputTag>("OuterTrackerDigiSource")),
0094       digiSimLinkSrc_(iConfig.getParameter<edm::InputTag>("OuterTrackerDigiSimSource")),
0095       simTrackSrc_(iConfig.getParameter<edm::InputTag>("SimTrackSource")),
0096       geomType_(iConfig.getParameter<std::string>("GeometryType")),
0097       phiValues(iConfig.getParameter<std::vector<double> >("PhiAngles")),
0098       otDigiToken_(consumes<edm::DetSetVector<Phase2TrackerDigi> >(otDigiSrc_)),
0099       otDigiSimLinkToken_(consumes<edm::DetSetVector<PixelDigiSimLink> >(digiSimLinkSrc_)),
0100       simTrackToken_(consumes<edm::SimTrackContainer>(simTrackSrc_)),
0101       topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0102       geomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>(edm::ESInputTag{"", geomType_})) {
0103   edm::LogInfo("TBeamTest") << ">>> Construct TBeamTest ";
0104 }
0105 
0106 //
0107 // destructor
0108 //
0109 TBeamTest::~TBeamTest() {
0110   // do anything here that needs to be done at desctruction time
0111   // (e.g. close files, deallocate resources etc.)
0112   edm::LogInfo("TBeamTest") << ">>> Destroy TBeamTest ";
0113 }
0114 //
0115 // -- Analyze
0116 //
0117 void TBeamTest::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0118   using namespace edm;
0119 
0120   // Get digis
0121 
0122   edm::Handle<edm::DetSetVector<PixelDigiSimLink> > digiSimLinkHandle;
0123   iEvent.getByToken(otDigiSimLinkToken_, digiSimLinkHandle);
0124 
0125   edm::Handle<edm::DetSetVector<Phase2TrackerDigi> > otDigiHandle;
0126   iEvent.getByToken(otDigiToken_, otDigiHandle);
0127   const DetSetVector<Phase2TrackerDigi>* digis = otDigiHandle.product();
0128 
0129   // Get SimTrack
0130   edm::Handle<edm::SimTrackContainer> simTrackHandle;
0131   iEvent.getByToken(simTrackToken_, simTrackHandle);
0132 
0133   const TrackerTopology* tTopo = &iSetup.getData(topoToken_);
0134 
0135   edm::ESWatcher<TrackerDigiGeometryRecord> theTkDigiGeomWatcher;
0136   if (theTkDigiGeomWatcher.check(iSetup)) {
0137     const TrackerGeometry* tkGeom = &iSetup.getData(geomToken_);
0138 
0139     edm::DetSetVector<Phase2TrackerDigi>::const_iterator DSViter;
0140     std::string moduleType;
0141     for (DSViter = digis->begin(); DSViter != digis->end(); DSViter++) {
0142       unsigned int rawid = DSViter->id;
0143       DetId detId(rawid);
0144       if (detId.det() != DetId::Detector::Tracker)
0145         continue;
0146 
0147       if (detId.subdetId() != StripSubdetector::TOB)
0148         continue;
0149       int layer = tTopo->getOTLayerNumber(rawid);
0150       if (layer != 4)
0151         continue;
0152 
0153       const GeomDetUnit* geomDetUnit = tkGeom->idToDetUnit(detId);
0154 
0155       const Phase2TrackerGeomDetUnit* tkDetUnit = dynamic_cast<const Phase2TrackerGeomDetUnit*>(geomDetUnit);
0156       int nColumns = tkDetUnit->specificTopology().ncolumns();
0157 
0158       edm::LogInfo("TBeamTest") << " Det Id = " << rawid;
0159 
0160       if (layer <= 3) {
0161         if (nColumns > 2)
0162           moduleType = "PSP_Modules";
0163         else
0164           moduleType = "PSS_Modules";
0165       } else
0166         moduleType = "2S_Modules";
0167 
0168       std::map<std::string, DigiMEs>::iterator pos = detMEs.find(moduleType);
0169       if (pos != detMEs.end()) {
0170         DigiMEs local_mes = pos->second;
0171         int nDigi = 0;
0172         int row_last = -1;
0173         int col_last = -1;
0174         int nclus = 0;
0175         int width = 1;
0176         int position = 0;
0177         float dPhi = 9999.9;
0178         for (DetSet<Phase2TrackerDigi>::const_iterator di = DSViter->begin(); di != DSViter->end(); di++) {
0179           int col = di->column();  // column
0180           int row = di->row();     // row
0181           MeasurementPoint mp(row + 0.5, col + 0.5);
0182           unsigned int channel = Phase2TrackerDigi::pixelToChannel(row, col);
0183           int tkIndx = matchedSimTrackIndex(digiSimLinkHandle, simTrackHandle, detId, channel);
0184 
0185           if (geomDetUnit && tkIndx != -1)
0186             dPhi = reco::deltaPhi((*simTrackHandle)[tkIndx].momentum().phi(), geomDetUnit->position().phi());
0187 
0188           nDigi++;
0189           edm::LogInfo("TBeamTest") << "  column " << col << " row " << row << std::endl;
0190           local_mes.PositionOfDigis->Fill(row + 1);
0191 
0192           if (row_last == -1) {
0193             width = 1;
0194             position = row + 1;
0195             nclus++;
0196           } else {
0197             if (abs(row - row_last) == 1 && col == col_last) {
0198               position += row + 1;
0199               width++;
0200             } else {
0201               position /= width;
0202               fillClusterWidth(local_mes, dPhi, width);
0203               local_mes.ClusterPosition->Fill(position);
0204               width = 1;
0205               position = row + 1;
0206               nclus++;
0207             }
0208           }
0209           edm::LogInfo("TBeamTest") << " row " << row << " col " << col << " row_last " << row_last << " col_last "
0210                                     << col_last << " width " << width;
0211           row_last = row;
0212           col_last = col;
0213         }
0214         position /= width;
0215         fillClusterWidth(local_mes, dPhi, width);
0216         local_mes.ClusterPosition->Fill(position);
0217         local_mes.NumberOfClusters->Fill(nclus);
0218         local_mes.NumberOfDigis->Fill(nDigi);
0219       }
0220     }
0221   }
0222 }
0223 //
0224 // -- Book Histograms
0225 //
0226 void TBeamTest::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const& iSetup) {
0227   std::string top_folder = config_.getParameter<std::string>("TopFolderName");
0228 
0229   std::vector<std::string> types;
0230   types.push_back("2S_Modules");
0231   types.push_back("PSP_Modules");
0232   types.push_back("PSS_Modules");
0233   ibooker.cd();
0234 
0235   for (const auto& itype : types) {
0236     std::stringstream folder_name;
0237 
0238     folder_name << top_folder << "/" << itype;
0239 
0240     edm::LogInfo("TBeamTest") << " Booking Histograms in : " << folder_name.str();
0241     ibooker.setCurrentFolder(folder_name.str());
0242 
0243     std::ostringstream HistoName;
0244 
0245     DigiMEs local_mes;
0246     edm::ParameterSet Parameters = config_.getParameter<edm::ParameterSet>("NumberOfDigisH");
0247     HistoName.str("");
0248     HistoName << "numberOfHits";
0249     local_mes.NumberOfDigis = ibooker.book1D(HistoName.str(),
0250                                              HistoName.str(),
0251                                              Parameters.getParameter<int32_t>("Nbins"),
0252                                              Parameters.getParameter<double>("xmin"),
0253                                              Parameters.getParameter<double>("xmax"));
0254 
0255     Parameters = config_.getParameter<edm::ParameterSet>("PositionOfDigisH");
0256     HistoName.str("");
0257     HistoName << "hitPositions";
0258     local_mes.PositionOfDigis = ibooker.book1D(HistoName.str(),
0259                                                HistoName.str(),
0260                                                Parameters.getParameter<int32_t>("Nxbins"),
0261                                                Parameters.getParameter<double>("xmin"),
0262                                                Parameters.getParameter<double>("xmax"));
0263     Parameters = config_.getParameter<edm::ParameterSet>("NumberOfClustersH");
0264     HistoName.str("");
0265     HistoName << "numberOfCluetsrs";
0266     local_mes.NumberOfClusters = ibooker.book1D(HistoName.str(),
0267                                                 HistoName.str(),
0268                                                 Parameters.getParameter<int32_t>("Nbins"),
0269                                                 Parameters.getParameter<double>("xmin"),
0270                                                 Parameters.getParameter<double>("xmax"));
0271     Parameters = config_.getParameter<edm::ParameterSet>("ClusterWidthH");
0272 
0273     for (unsigned int i = 0; i < phiValues.size(); i++) {
0274       HistoName.str("");
0275       HistoName << "clusterWidth_";
0276       HistoName << i;
0277 
0278       local_mes.ClusterWidths.push_back(ibooker.book1D(HistoName.str(),
0279                                                        HistoName.str(),
0280                                                        Parameters.getParameter<int32_t>("Nbins"),
0281                                                        Parameters.getParameter<double>("xmin"),
0282                                                        Parameters.getParameter<double>("xmax")));
0283     }
0284     Parameters = config_.getParameter<edm::ParameterSet>("ClusterPositionH");
0285     HistoName.str("");
0286     HistoName << "clusterPositions";
0287     local_mes.ClusterPosition = ibooker.book1D(HistoName.str(),
0288                                                HistoName.str(),
0289                                                Parameters.getParameter<int32_t>("Nbins"),
0290                                                Parameters.getParameter<double>("xmin"),
0291                                                Parameters.getParameter<double>("xmax"));
0292     detMEs.insert(std::make_pair(itype, local_mes));
0293   }
0294 }
0295 int TBeamTest::matchedSimTrackIndex(edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& linkHandle,
0296                                     edm::Handle<edm::SimTrackContainer>& simTkHandle,
0297                                     DetId detId,
0298                                     unsigned int& channel) {
0299   int simTrkIndx = -1;
0300   unsigned int simTrkId = 0;
0301   edm::DetSetVector<PixelDigiSimLink>::const_iterator isearch = linkHandle->find(detId);
0302 
0303   if (isearch == linkHandle->end())
0304     return simTrkIndx;
0305 
0306   edm::DetSet<PixelDigiSimLink> link_detset = (*linkHandle)[detId];
0307   // Loop over DigiSimLink in this det unit
0308   for (edm::DetSet<PixelDigiSimLink>::const_iterator it = link_detset.data.begin(); it != link_detset.data.end();
0309        it++) {
0310     if (channel == it->channel()) {
0311       simTrkId = it->SimTrackId();
0312       break;
0313     }
0314   }
0315   if (simTrkId == 0)
0316     return simTrkIndx;
0317   edm::SimTrackContainer sim_tracks = (*simTkHandle.product());
0318   for (unsigned int itk = 0; itk < sim_tracks.size(); itk++) {
0319     if (sim_tracks[itk].trackId() == simTrkId) {
0320       simTrkIndx = itk;
0321       break;
0322     }
0323   }
0324   return simTrkIndx;
0325 }
0326 void TBeamTest::fillClusterWidth(DigiMEs& mes, float dphi, float width) {
0327   for (unsigned int i = 0; i < phiValues.size(); i++) {
0328     float angle_min = (phiValues[i] - 0.1) * std::acos(-1.0) / 180.0;
0329     float angle_max = (phiValues[i] + 0.1) * std::acos(-1.0) / 180.0;
0330     if (std::fabs(dphi) > angle_min && std::fabs(dphi) < angle_max) {
0331       mes.ClusterWidths[i]->Fill(width);
0332       break;
0333     }
0334   }
0335 }
0336 //define this as a plug-in
0337 DEFINE_FWK_MODULE(TBeamTest);
[ 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!