Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:11:23

0001 // Package:    SiPixelMonitorTrack
0002 // Class:      SiPixelTrackResidualSource
0003 //
0004 // class SiPixelTrackResidualSource SiPixelTrackResidualSource.cc
0005 //       DQM/SiPixelMonitorTrack/src/SiPixelTrackResidualSource.cc
0006 //
0007 // Description: SiPixel hit-to-track residual data quality monitoring modules
0008 // Implementation: prototype -> improved -> never final - end of the 1st step
0009 //
0010 // Original Author: Shan-Huei Chuang
0011 //         Created: Fri Mar 23 18:41:42 CET 2007
0012 //         Updated by Lukas Wehrli (plots for clusters on/off track added)
0013 
0014 #include <iostream>
0015 #include <map>
0016 #include <string>
0017 #include <utility>
0018 #include <vector>
0019 
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 
0023 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0024 
0025 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0026 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0027 #include "DataFormats/SiPixelDetId/interface/PixelBarrelNameUpgrade.h"
0028 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0029 
0030 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0031 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0032 
0033 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0034 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0035 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0036 
0037 #include "DQM/SiPixelCommon/interface/SiPixelFolderOrganizer.h"
0038 #include "DQM/SiPixelMonitorTrack/interface/SiPixelTrackResidualSource.h"
0039 
0040 // Claudia new libraries
0041 #include "DataFormats/VertexReco/interface/Vertex.h"
0042 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0043 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0044 
0045 using namespace std;
0046 using namespace edm;
0047 
0048 SiPixelTrackResidualSource::SiPixelTrackResidualSource(const edm::ParameterSet &pSet)
0049     : pSet_(pSet),
0050       modOn(pSet.getUntrackedParameter<bool>("modOn", true)),
0051       reducedSet(pSet.getUntrackedParameter<bool>("reducedSet", true)),
0052       ladOn(pSet.getUntrackedParameter<bool>("ladOn", false)),
0053       layOn(pSet.getUntrackedParameter<bool>("layOn", false)),
0054       phiOn(pSet.getUntrackedParameter<bool>("phiOn", false)),
0055       ringOn(pSet.getUntrackedParameter<bool>("ringOn", false)),
0056       bladeOn(pSet.getUntrackedParameter<bool>("bladeOn", false)),
0057       diskOn(pSet.getUntrackedParameter<bool>("diskOn", false)),
0058       isUpgrade(pSet.getUntrackedParameter<bool>("isUpgrade", false)),
0059       noOfLayers(0),
0060       noOfDisks(0) {
0061   pSet_ = pSet;
0062   debug_ = pSet_.getUntrackedParameter<bool>("debug", false);
0063   src_ = pSet_.getParameter<edm::InputTag>("src");
0064   clustersrc_ = pSet_.getParameter<edm::InputTag>("clustersrc");
0065   tracksrc_ = pSet_.getParameter<edm::InputTag>("trajectoryInput");
0066   ttrhbuilder_ = pSet_.getParameter<std::string>("TTRHBuilder");
0067   ptminres_ = pSet.getUntrackedParameter<double>("PtMinRes", 4.0);
0068   beamSpotToken_ = consumes<reco::BeamSpot>(std::string("offlineBeamSpot"));
0069   vtxsrc_ = pSet_.getUntrackedParameter<std::string>("vtxsrc", "offlinePrimaryVertices");
0070   offlinePrimaryVerticesToken_ =
0071       consumes<reco::VertexCollection>(vtxsrc_);  // consumes<reco::VertexCollection>(std::string("hiSelectedVertex"));
0072                                                   // //"offlinePrimaryVertices"));
0073   generalTracksToken_ = consumes<reco::TrackCollection>(pSet_.getParameter<edm::InputTag>("tracksrc"));
0074   tracksrcToken_ = consumes<std::vector<Trajectory>>(pSet_.getParameter<edm::InputTag>("trajectoryInput"));
0075   trackToken_ = consumes<std::vector<reco::Track>>(pSet_.getParameter<edm::InputTag>("trajectoryInput"));
0076   trackAssociationToken_ =
0077       consumes<TrajTrackAssociationCollection>(pSet_.getParameter<edm::InputTag>("trajectoryInput"));
0078   clustersrcToken_ = consumes<edmNew::DetSetVector<SiPixelCluster>>(pSet_.getParameter<edm::InputTag>("clustersrc"));
0079   digisrc_ = pSet_.getParameter<edm::InputTag>("digisrc");
0080   digisrcToken_ = consumes<edm::DetSetVector<PixelDigi>>(pSet_.getParameter<edm::InputTag>("digisrc"));
0081 
0082   trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0083   trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0084   transientTrackBuilderToken_ =
0085       esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0086   transientTrackingRecHitBuilderToken_ =
0087       esConsumes<TransientTrackingRecHitBuilder, TransientRecHitRecord>(edm::ESInputTag("", ttrhbuilder_));
0088   trackerTopoTokenBeginRun_ = esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>();
0089   trackerGeomTokenBeginRun_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>();
0090 
0091   LogInfo("PixelDQM") << "SiPixelTrackResidualSource constructor" << endl;
0092   LogInfo("PixelDQM") << "Mod/Lad/Lay/Phi " << modOn << "/" << ladOn << "/" << layOn << "/" << phiOn << std::endl;
0093   LogInfo("PixelDQM") << "Blade/Disk/Ring" << bladeOn << "/" << diskOn << "/" << ringOn << std::endl;
0094 
0095   topFolderName_ = pSet_.getParameter<std::string>("TopFolderName");
0096 
0097   firstRun = true;
0098   NTotal = 0;
0099   NLowProb = 0;
0100 }
0101 
0102 SiPixelTrackResidualSource::~SiPixelTrackResidualSource() {
0103   LogInfo("PixelDQM") << "SiPixelTrackResidualSource destructor" << endl;
0104 
0105   std::map<uint32_t, SiPixelTrackResidualModule *>::iterator struct_iter;
0106   for (struct_iter = theSiPixelStructure.begin(); struct_iter != theSiPixelStructure.end(); struct_iter++) {
0107     delete struct_iter->second;
0108     struct_iter->second = nullptr;
0109   }
0110 }
0111 
0112 void SiPixelTrackResidualSource::dqmBeginRun(const edm::Run &r, edm::EventSetup const &iSetup) {
0113   LogInfo("PixelDQM") << "SiPixelTrackResidualSource beginRun()" << endl;
0114   // retrieve TrackerGeometry for pixel dets
0115   edm::ESHandle<TrackerGeometry> TG = iSetup.getHandle(trackerGeomTokenBeginRun_);
0116   if (debug_)
0117     LogVerbatim("PixelDQM") << "TrackerGeometry " << &(*TG) << " size is " << TG->dets().size() << endl;
0118   edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(trackerTopoTokenBeginRun_);
0119   const TrackerTopology *pTT = tTopoHandle.product();
0120 
0121   // build theSiPixelStructure with the pixel barrel and endcap dets from
0122   // TrackerGeometry
0123   for (TrackerGeometry::DetContainer::const_iterator pxb = TG->detsPXB().begin(); pxb != TG->detsPXB().end(); pxb++) {
0124     if (dynamic_cast<PixelGeomDetUnit const *>((*pxb)) != nullptr) {
0125       SiPixelTrackResidualModule *module = new SiPixelTrackResidualModule((*pxb)->geographicalId().rawId());
0126       theSiPixelStructure.insert(
0127           pair<uint32_t, SiPixelTrackResidualModule *>((*pxb)->geographicalId().rawId(), module));
0128       // int DBlayer = PixelBarrelNameWrapper(pSet_,
0129       // DetId((*pxb)->geographicalId())).layerName();
0130       int DBlayer = PixelBarrelName(DetId((*pxb)->geographicalId()), pTT, isUpgrade).layerName();
0131       if (noOfLayers < DBlayer)
0132         noOfLayers = DBlayer;
0133     }
0134   }
0135   for (TrackerGeometry::DetContainer::const_iterator pxf = TG->detsPXF().begin(); pxf != TG->detsPXF().end(); pxf++) {
0136     if (dynamic_cast<PixelGeomDetUnit const *>((*pxf)) != nullptr) {
0137       SiPixelTrackResidualModule *module = new SiPixelTrackResidualModule((*pxf)->geographicalId().rawId());
0138       theSiPixelStructure.insert(
0139           pair<uint32_t, SiPixelTrackResidualModule *>((*pxf)->geographicalId().rawId(), module));
0140       int DBdisk;
0141       DBdisk = PixelEndcapName(DetId((*pxf)->geographicalId()), pTT, isUpgrade).diskName();
0142       if (noOfDisks < DBdisk)
0143         noOfDisks = DBdisk;
0144     }
0145   }
0146   LogInfo("PixelDQM") << "SiPixelStructure size is " << theSiPixelStructure.size() << endl;
0147 }
0148 
0149 void SiPixelTrackResidualSource::bookHistograms(DQMStore::IBooker &iBooker,
0150                                                 edm::Run const &,
0151                                                 edm::EventSetup const &iSetup) {
0152   // book residual histograms in theSiPixelFolder - one (x,y) pair of histograms
0153   // per det
0154   SiPixelFolderOrganizer theSiPixelFolder(false);
0155 
0156   edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(trackerTopoTokenBeginRun_);
0157   const TrackerTopology *pTT = tTopoHandle.product();
0158 
0159   std::stringstream nameX, titleX, nameY, titleY;
0160 
0161   if (ladOn) {
0162     iBooker.setCurrentFolder(topFolderName_ + "/Barrel");
0163     for (int i = 1; i <= noOfLayers; i++) {
0164       nameX.str(std::string());
0165       nameX << "siPixelTrackResidualsX_SummedLayer_" << i;
0166       titleX.str(std::string());
0167       titleX << "Layer" << i << "Hit-to-Track Residual in r-phi";
0168       meResidualXSummedLay.push_back(iBooker.book1D(nameX.str(), titleX.str(), 100, -150, 150));
0169       meResidualXSummedLay.at(i - 1)->setAxisTitle("hit-to-track residual in r-phi (um)", 1);
0170       nameY.str(std::string());
0171       nameY << "siPixelTrackResidualsY_SummedLayer_" << i;
0172       titleY.str(std::string());
0173       titleY << "Layer" << i << "Hit-to-Track Residual in Z";
0174       meResidualYSummedLay.push_back(iBooker.book1D(nameY.str(), titleY.str(), 100, -300, 300));
0175       meResidualYSummedLay.at(i - 1)->setAxisTitle("hit-to-track residual in z (um)", 1);
0176     }
0177   }
0178 
0179   for (std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd = theSiPixelStructure.begin();
0180        pxd != theSiPixelStructure.end();
0181        pxd++) {
0182     if (modOn) {
0183       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 0, isUpgrade))
0184         (*pxd).second->book(pSet_, pTT, iBooker, reducedSet, 0, isUpgrade);
0185       else
0186         throw cms::Exception("LogicError") << "SiPixelTrackResidualSource Folder Creation Failed! ";
0187     }
0188     if (ladOn) {
0189       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 1, isUpgrade)) {
0190         (*pxd).second->book(pSet_, pTT, iBooker, reducedSet, 1, isUpgrade);
0191       } else
0192         throw cms::Exception("LogicError") << "SiPixelTrackResidualSource ladder Folder Creation Failed! ";
0193     }
0194     if (layOn) {
0195       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 2, isUpgrade))
0196         (*pxd).second->book(pSet_, pTT, iBooker, reducedSet, 2, isUpgrade);
0197       else
0198         throw cms::Exception("LogicError") << "SiPixelTrackResidualSource layer Folder Creation Failed! ";
0199     }
0200     if (phiOn) {
0201       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 3, isUpgrade))
0202         (*pxd).second->book(pSet_, pTT, iBooker, reducedSet, 3, isUpgrade);
0203       else
0204         throw cms::Exception("LogicError") << "SiPixelTrackResidualSource phi Folder Creation Failed! ";
0205     }
0206     if (bladeOn) {
0207       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 4, isUpgrade))
0208         (*pxd).second->book(pSet_, pTT, iBooker, reducedSet, 4, isUpgrade);
0209       else
0210         throw cms::Exception("LogicError") << "SiPixelTrackResidualSource Blade Folder Creation Failed! ";
0211     }
0212     if (diskOn) {
0213       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 5, isUpgrade))
0214         (*pxd).second->book(pSet_, pTT, iBooker, reducedSet, 5, isUpgrade);
0215       else
0216         throw cms::Exception("LogicError") << "SiPixelTrackResidualSource Disk Folder Creation Failed! ";
0217     }
0218     if (ringOn) {
0219       if (theSiPixelFolder.setModuleFolder(iBooker, (*pxd).first, 6, isUpgrade))
0220         (*pxd).second->book(pSet_, pTT, iBooker, reducedSet, 6, isUpgrade);
0221       else
0222         throw cms::Exception("LogicError") << "SiPixelTrackResidualSource Ring Folder Creation Failed! ";
0223     }
0224   }
0225 
0226   //   edm::InputTag tracksrc =
0227   //   pSet_.getParameter<edm::InputTag>("trajectoryInput"); edm::InputTag
0228   //   clustersrc = pSet_.getParameter<edm::InputTag>("clustersrc");
0229 
0230   // number of tracks
0231   iBooker.setCurrentFolder(topFolderName_ + "/Tracks");
0232   meNofTracks_ = iBooker.book1D("ntracks_" + tracksrc_.label(), "Number of Tracks", 4, 0, 4);
0233   meNofTracks_->setAxisTitle("Number of Tracks", 1);
0234   meNofTracks_->setBinLabel(1, "All");
0235   meNofTracks_->setBinLabel(2, "Pixel");
0236   meNofTracks_->setBinLabel(3, "BPix");
0237   meNofTracks_->setBinLabel(4, "FPix");
0238 
0239   // number of tracks in pixel fiducial volume
0240   iBooker.setCurrentFolder(topFolderName_ + "/Tracks");
0241   meNofTracksInPixVol_ = iBooker.book1D(
0242       "ntracksInPixVol_" + tracksrc_.label(), "Number of Tracks crossing Pixel fiducial Volume", 2, 0, 2);
0243   meNofTracksInPixVol_->setAxisTitle("Number of Tracks", 1);
0244   meNofTracksInPixVol_->setBinLabel(1, "With Hits");
0245   meNofTracksInPixVol_->setBinLabel(2, "Without Hits");
0246 
0247   // number of clusters (associated to track / not associated)
0248   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OnTrack");
0249   meNofClustersOnTrack_ =
0250       iBooker.book1D("nclusters_" + clustersrc_.label() + "_tot", "Number of Clusters (on track)", 3, 0, 3);
0251   meNofClustersOnTrack_->setAxisTitle("Number of Clusters on Track", 1);
0252   meNofClustersOnTrack_->setBinLabel(1, "All");
0253   meNofClustersOnTrack_->setBinLabel(2, "BPix");
0254   meNofClustersOnTrack_->setBinLabel(3, "FPix");
0255   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OffTrack");
0256   meNofClustersNotOnTrack_ =
0257       iBooker.book1D("nclusters_" + clustersrc_.label() + "_tot", "Number of Clusters (off track)", 3, 0, 3);
0258   meNofClustersNotOnTrack_->setAxisTitle("Number of Clusters off Track", 1);
0259   meNofClustersNotOnTrack_->setBinLabel(1, "All");
0260   meNofClustersNotOnTrack_->setBinLabel(2, "BPix");
0261   meNofClustersNotOnTrack_->setBinLabel(3, "FPix");
0262 
0263   // cluster charge and size
0264   // charge
0265   // on track
0266   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OnTrack");
0267   std::stringstream ss1, ss2;
0268   meClChargeOnTrack_all = iBooker.book1D("charge_" + clustersrc_.label(), "Charge (on track)", 500, 0., 500.);
0269   meClChargeOnTrack_all->setAxisTitle("Charge size (in ke)", 1);
0270   meClChargeOnTrack_bpix =
0271       iBooker.book1D("charge_" + clustersrc_.label() + "_Barrel", "Charge (on track, barrel)", 500, 0., 500.);
0272   meClChargeOnTrack_bpix->setAxisTitle("Charge size (in ke)", 1);
0273   meClChargeOnTrack_fpix =
0274       iBooker.book1D("charge_" + clustersrc_.label() + "_Endcap", "Charge (on track, endcap)", 500, 0., 500.);
0275   meClChargeOnTrack_fpix->setAxisTitle("Charge size (in ke)", 1);
0276   for (int i = 1; i <= noOfLayers; i++) {
0277     ss1.str(std::string());
0278     ss1 << "charge_" + clustersrc_.label() + "_Layer_" << i;
0279     ss2.str(std::string());
0280     ss2 << "Charge (on track, layer" << i << ")";
0281     meClChargeOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0282     meClChargeOnTrack_layers.at(i - 1)->setAxisTitle("Charge size (in ke)", 1);
0283   }
0284   for (int i = 1; i <= noOfDisks; i++) {
0285     ss1.str(std::string());
0286     ss1 << "charge_" + clustersrc_.label() + "_Disk_p" << i;
0287     ss2.str(std::string());
0288     ss2 << "Charge (on track, diskp" << i << ")";
0289     meClChargeOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0290     meClChargeOnTrack_diskps.at(i - 1)->setAxisTitle("Charge size (in ke)", 1);
0291   }
0292   for (int i = 1; i <= noOfDisks; i++) {
0293     ss1.str(std::string());
0294     ss1 << "charge_" + clustersrc_.label() + "_Disk_m" << i;
0295     ss2.str(std::string());
0296     ss2 << "Charge (on track, diskm" << i << ")";
0297     meClChargeOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0298     meClChargeOnTrack_diskms.at(i - 1)->setAxisTitle("Charge size (in ke)", 1);
0299   }
0300   // off track
0301   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OffTrack");
0302   meClChargeNotOnTrack_all = iBooker.book1D("charge_" + clustersrc_.label(), "Charge (off track)", 500, 0., 500.);
0303   meClChargeNotOnTrack_all->setAxisTitle("Charge size (in ke)", 1);
0304   meClChargeNotOnTrack_bpix =
0305       iBooker.book1D("charge_" + clustersrc_.label() + "_Barrel", "Charge (off track, barrel)", 500, 0., 500.);
0306   meClChargeNotOnTrack_bpix->setAxisTitle("Charge size (in ke)", 1);
0307   meClChargeNotOnTrack_fpix =
0308       iBooker.book1D("charge_" + clustersrc_.label() + "_Endcap", "Charge (off track, endcap)", 500, 0., 500.);
0309   meClChargeNotOnTrack_fpix->setAxisTitle("Charge size (in ke)", 1);
0310   for (int i = 1; i <= noOfLayers; i++) {
0311     ss1.str(std::string());
0312     ss1 << "charge_" + clustersrc_.label() + "_Layer_" << i;
0313     ss2.str(std::string());
0314     ss2 << "Charge (off track, layer" << i << ")";
0315     meClChargeNotOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0316     meClChargeNotOnTrack_layers.at(i - 1)->setAxisTitle("Charge size (in ke)", 1);
0317   }
0318   for (int i = 1; i <= noOfDisks; i++) {
0319     ss1.str(std::string());
0320     ss1 << "charge_" + clustersrc_.label() + "_Disk_p" << i;
0321     ss2.str(std::string());
0322     ss2 << "Charge (off track, diskp" << i << ")";
0323     meClChargeNotOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0324     meClChargeNotOnTrack_diskps.at(i - 1)->setAxisTitle("Charge size (in ke)", 1);
0325   }
0326   for (int i = 1; i <= noOfDisks; i++) {
0327     ss1.str(std::string());
0328     ss1 << "charge_" + clustersrc_.label() + "_Disk_m" << i;
0329     ss2.str(std::string());
0330     ss2 << "Charge (off track, diskm" << i << ")";
0331     meClChargeNotOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0332     meClChargeNotOnTrack_diskms.at(i - 1)->setAxisTitle("Charge size (in ke)", 1);
0333   }
0334 
0335   // size
0336   // on track
0337   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OnTrack");
0338   meClSizeOnTrack_all = iBooker.book1D("size_" + clustersrc_.label(), "Size (on track)", 100, 0., 100.);
0339   meClSizeOnTrack_all->setAxisTitle("Cluster size (in pixels)", 1);
0340   meClSizeOnTrack_bpix =
0341       iBooker.book1D("size_" + clustersrc_.label() + "_Barrel", "Size (on track, barrel)", 100, 0., 100.);
0342   meClSizeOnTrack_bpix->setAxisTitle("Cluster size (in pixels)", 1);
0343   meClSizeOnTrack_fpix =
0344       iBooker.book1D("size_" + clustersrc_.label() + "_Endcap", "Size (on track, endcap)", 100, 0., 100.);
0345   meClSizeOnTrack_fpix->setAxisTitle("Cluster size (in pixels)", 1);
0346   for (int i = 1; i <= noOfLayers; i++) {
0347     ss1.str(std::string());
0348     ss1 << "size_" + clustersrc_.label() + "_Layer_" << i;
0349     ss2.str(std::string());
0350     ss2 << "Size (on track, layer" << i << ")";
0351     meClSizeOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0352     meClSizeOnTrack_layers.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0353   }
0354   for (int i = 1; i <= noOfDisks; i++) {
0355     ss1.str(std::string());
0356     ss1 << "size_" + clustersrc_.label() + "_Disk_p" << i;
0357     ss2.str(std::string());
0358     ss2 << "Size (on track, diskp" << i << ")";
0359     meClSizeOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0360     meClSizeOnTrack_diskps.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0361   }
0362   for (int i = 1; i <= noOfDisks; i++) {
0363     ss1.str(std::string());
0364     ss1 << "size_" + clustersrc_.label() + "_Disk_m1" << i;
0365     ss2.str(std::string());
0366     ss2 << "Size (on track, diskm" << i << ")";
0367     meClSizeOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0368     meClSizeOnTrack_diskms.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0369   }
0370   meClSizeXOnTrack_all = iBooker.book1D("sizeX_" + clustersrc_.label(), "SizeX (on track)", 100, 0., 100.);
0371   meClSizeXOnTrack_all->setAxisTitle("Cluster sizeX (in pixels)", 1);
0372   meClSizeXOnTrack_bpix =
0373       iBooker.book1D("sizeX_" + clustersrc_.label() + "_Barrel", "SizeX (on track, barrel)", 100, 0., 100.);
0374   meClSizeXOnTrack_bpix->setAxisTitle("Cluster sizeX (in pixels)", 1);
0375   meClSizeXOnTrack_fpix =
0376       iBooker.book1D("sizeX_" + clustersrc_.label() + "_Endcap", "SizeX (on track, endcap)", 100, 0., 100.);
0377   meClSizeXOnTrack_fpix->setAxisTitle("Cluster sizeX (in pixels)", 1);
0378   for (int i = 1; i <= noOfLayers; i++) {
0379     ss1.str(std::string());
0380     ss1 << "sizeX_" + clustersrc_.label() + "_Layer_" << i;
0381     ss2.str(std::string());
0382     ss2 << "SizeX (on track, layer" << i << ")";
0383     meClSizeXOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0384     meClSizeXOnTrack_layers.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0385   }
0386   for (int i = 1; i <= noOfDisks; i++) {
0387     ss1.str(std::string());
0388     ss1 << "sizeX_" + clustersrc_.label() + "_Disk_p" << i;
0389     ss2.str(std::string());
0390     ss2 << "SizeX (on track, diskp" << i << ")";
0391     meClSizeXOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0392     meClSizeXOnTrack_diskps.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0393   }
0394   for (int i = 1; i <= noOfDisks; i++) {
0395     ss1.str(std::string());
0396     ss1 << "sizeX_" + clustersrc_.label() + "_Disk_m" << i;
0397     ss2.str(std::string());
0398     ss2 << "SizeX (on track, diskm" << i << ")";
0399     meClSizeXOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0400     meClSizeXOnTrack_diskms.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0401   }
0402   meClSizeYOnTrack_all = iBooker.book1D("sizeY_" + clustersrc_.label(), "SizeY (on track)", 100, 0., 100.);
0403   meClSizeYOnTrack_all->setAxisTitle("Cluster sizeY (in pixels)", 1);
0404   meClSizeYOnTrack_bpix =
0405       iBooker.book1D("sizeY_" + clustersrc_.label() + "_Barrel", "SizeY (on track, barrel)", 100, 0., 100.);
0406   meClSizeYOnTrack_bpix->setAxisTitle("Cluster sizeY (in pixels)", 1);
0407   meClSizeYOnTrack_fpix =
0408       iBooker.book1D("sizeY_" + clustersrc_.label() + "_Endcap", "SizeY (on track, endcap)", 100, 0., 100.);
0409   meClSizeYOnTrack_fpix->setAxisTitle("Cluster sizeY (in pixels)", 1);
0410   for (int i = 1; i <= noOfLayers; i++) {
0411     ss1.str(std::string());
0412     ss1 << "sizeY_" + clustersrc_.label() + "_Layer_" << i;
0413     ss2.str(std::string());
0414     ss2 << "SizeY (on track, layer" << i << ")";
0415     meClSizeYOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0416     meClSizeYOnTrack_layers.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0417   }
0418   for (int i = 1; i <= noOfDisks; i++) {
0419     ss1.str(std::string());
0420     ss1 << "sizeY_" + clustersrc_.label() + "_Disk_p" << i;
0421     ss2.str(std::string());
0422     ss2 << "SizeY (on track, diskp" << i << ")";
0423     meClSizeYOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0424     meClSizeYOnTrack_diskps.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0425   }
0426   for (int i = 1; i <= noOfDisks; i++) {
0427     ss1.str(std::string());
0428     ss1 << "sizeY_" + clustersrc_.label() + "_Disk_m" << i;
0429     ss2.str(std::string());
0430     ss2 << "SizeY (on track, diskm" << i << ")";
0431     meClSizeYOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0432     meClSizeYOnTrack_diskms.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0433   }
0434   // off track
0435   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OffTrack");
0436   meClSizeNotOnTrack_all = iBooker.book1D("size_" + clustersrc_.label(), "Size (off track)", 100, 0., 100.);
0437   meClSizeNotOnTrack_all->setAxisTitle("Cluster size (in pixels)", 1);
0438   meClSizeNotOnTrack_bpix =
0439       iBooker.book1D("size_" + clustersrc_.label() + "_Barrel", "Size (off track, barrel)", 100, 0., 100.);
0440   meClSizeNotOnTrack_bpix->setAxisTitle("Cluster size (in pixels)", 1);
0441   meClSizeNotOnTrack_fpix =
0442       iBooker.book1D("size_" + clustersrc_.label() + "_Endcap", "Size (off track, endcap)", 100, 0., 100.);
0443   meClSizeNotOnTrack_fpix->setAxisTitle("Cluster size (in pixels)", 1);
0444   for (int i = 1; i <= noOfLayers; i++) {
0445     ss1.str(std::string());
0446     ss1 << "size_" + clustersrc_.label() + "_Layer_" << i;
0447     ss2.str(std::string());
0448     ss2 << "Size (off track, layer" << i << ")";
0449     meClSizeNotOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0450     meClSizeNotOnTrack_layers.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0451   }
0452 
0453   for (int i = 1; i <= noOfLayers; i++) {
0454     int ybins = -1;
0455     float ymin = 0.;
0456     float ymax = 0.;
0457     if (i == 1) {
0458       ybins = 42;
0459       ymin = -10.5;
0460       ymax = 10.5;
0461     }
0462     if (i == 2) {
0463       ybins = 66;
0464       ymin = -16.5;
0465       ymax = 16.5;
0466     }
0467     if (i == 3) {
0468       ybins = 90;
0469       ymin = -22.5;
0470       ymax = 22.5;
0471     }
0472     ss1.str(std::string());
0473     ss1 << "pix_bar Occ_roc_offtrack" + digisrc_.label() + "_layer_" << i;
0474     ss2.str(std::string());
0475     ss2 << "Pixel Barrel Occupancy, ROC level (Off Track): Layer " << i;
0476     meZeroRocLadvsModOffTrackBarrel.push_back(iBooker.book2D(ss1.str(), ss2.str(), 72, -4.5, 4.5, ybins, ymin, ymax));
0477     meZeroRocLadvsModOffTrackBarrel.at(i - 1)->setAxisTitle("ROC / Module", 1);
0478     meZeroRocLadvsModOffTrackBarrel.at(i - 1)->setAxisTitle("ROC / Ladder", 2);
0479   }
0480 
0481   for (int i = 1; i <= noOfDisks; i++) {
0482     ss1.str(std::string());
0483     ss1 << "size_" + clustersrc_.label() + "_Disk_p" << i;
0484     ss2.str(std::string());
0485     ss2 << "Size (off track, diskp" << i << ")";
0486     meClSizeNotOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0487     meClSizeNotOnTrack_diskps.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0488   }
0489   for (int i = 1; i <= noOfDisks; i++) {
0490     ss1.str(std::string());
0491     ss1 << "size_" + clustersrc_.label() + "_Disk_m" << i;
0492     ss2.str(std::string());
0493     ss2 << "Size (off track, diskm" << i << ")";
0494     meClSizeNotOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0495     meClSizeNotOnTrack_diskms.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0496   }
0497   meClSizeXNotOnTrack_all = iBooker.book1D("sizeX_" + clustersrc_.label(), "SizeX (off track)", 100, 0., 100.);
0498   meClSizeXNotOnTrack_all->setAxisTitle("Cluster sizeX (in pixels)", 1);
0499   meClSizeXNotOnTrack_bpix =
0500       iBooker.book1D("sizeX_" + clustersrc_.label() + "_Barrel", "SizeX (off track, barrel)", 100, 0., 100.);
0501   meClSizeXNotOnTrack_bpix->setAxisTitle("Cluster sizeX (in pixels)", 1);
0502   meClSizeXNotOnTrack_fpix =
0503       iBooker.book1D("sizeX_" + clustersrc_.label() + "_Endcap", "SizeX (off track, endcap)", 100, 0., 100.);
0504   meClSizeXNotOnTrack_fpix->setAxisTitle("Cluster sizeX (in pixels)", 1);
0505   for (int i = 1; i <= noOfLayers; i++) {
0506     ss1.str(std::string());
0507     ss1 << "sizeX_" + clustersrc_.label() + "_Layer_" << i;
0508     ss2.str(std::string());
0509     ss2 << "SizeX (off track, layer" << i << ")";
0510     meClSizeXNotOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0511     meClSizeXNotOnTrack_layers.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0512   }
0513   for (int i = 1; i <= noOfDisks; i++) {
0514     ss1.str(std::string());
0515     ss1 << "sizeX_" + clustersrc_.label() + "_Disk_p" << i;
0516     ss2.str(std::string());
0517     ss2 << "SizeX (off track, diskp" << i << ")";
0518     meClSizeXNotOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0519     meClSizeXNotOnTrack_diskps.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0520   }
0521   for (int i = 1; i <= noOfDisks; i++) {
0522     ss1.str(std::string());
0523     ss1 << "sizeX_" + clustersrc_.label() + "_Disk_m" << i;
0524     ss2.str(std::string());
0525     ss2 << "SizeX (off track, diskm" << i << ")";
0526     meClSizeXNotOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0527     meClSizeXNotOnTrack_diskms.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0528   }
0529   meClSizeYNotOnTrack_all = iBooker.book1D("sizeY_" + clustersrc_.label(), "SizeY (off track)", 100, 0., 100.);
0530   meClSizeYNotOnTrack_all->setAxisTitle("Cluster sizeY (in pixels)", 1);
0531   meClSizeYNotOnTrack_bpix =
0532       iBooker.book1D("sizeY_" + clustersrc_.label() + "_Barrel", "SizeY (off track, barrel)", 100, 0., 100.);
0533   meClSizeYNotOnTrack_bpix->setAxisTitle("Cluster sizeY (in pixels)", 1);
0534   meClSizeYNotOnTrack_fpix =
0535       iBooker.book1D("sizeY_" + clustersrc_.label() + "_Endcap", "SizeY (off track, endcap)", 100, 0., 100.);
0536   meClSizeYNotOnTrack_fpix->setAxisTitle("Cluster sizeY (in pixels)", 1);
0537   for (int i = 1; i <= noOfLayers; i++) {
0538     ss1.str(std::string());
0539     ss1 << "sizeY_" + clustersrc_.label() + "_Layer_" << i;
0540     ss2.str(std::string());
0541     ss2 << "SizeY (off track, layer" << i << ")";
0542     meClSizeYNotOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0543     meClSizeYNotOnTrack_layers.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0544   }
0545   for (int i = 1; i <= noOfDisks; i++) {
0546     ss1.str(std::string());
0547     ss1 << "sizeY_" + clustersrc_.label() + "_Disk_p" << i;
0548     ss2.str(std::string());
0549     ss2 << "SizeY (off track, diskp" << i << ")";
0550     meClSizeYNotOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0551     meClSizeYNotOnTrack_diskps.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0552   }
0553   for (int i = 1; i <= noOfDisks; i++) {
0554     ss1.str(std::string());
0555     ss1 << "sizeY_" + clustersrc_.label() + "_Disk_m" << i;
0556     ss2.str(std::string());
0557     ss2 << "SizeY (off track, diskm" << i << ")";
0558     meClSizeYNotOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0559     meClSizeYNotOnTrack_diskms.at(i - 1)->setAxisTitle("Cluster size (in pixels)", 1);
0560   }
0561 
0562   // cluster global position
0563   // on track
0564   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OnTrack");
0565   // bpix
0566   for (int i = 1; i <= noOfLayers; i++) {
0567     ss1.str(std::string());
0568     ss1 << "position_" + clustersrc_.label() + "_Layer_" << i;
0569     ss2.str(std::string());
0570     ss2 << "Clusters Layer" << i << " (on track)";
0571     meClPosLayersOnTrack.push_back(iBooker.book2D(ss1.str(), ss2.str(), 200, -30., 30., 128, -3.2, 3.2));
0572     meClPosLayersOnTrack.at(i - 1)->setAxisTitle("Global Z (cm)", 1);
0573     meClPosLayersOnTrack.at(i - 1)->setAxisTitle("Global #phi", 2);
0574 
0575     int ybins = -1;
0576     float ymin = 0.;
0577     float ymax = 0.;
0578     if (i == 1) {
0579       ybins = 23;
0580       ymin = -11.5;
0581       ymax = 11.5;
0582     }
0583     if (i == 2) {
0584       ybins = 33;
0585       ymin = -17.5;
0586       ymax = 17.5;
0587     }
0588     if (i == 3) {
0589       ybins = 45;
0590       ymin = -24.5;
0591       ymax = 24.5;
0592     }
0593     ss1.str(std::string());
0594     ss1 << "position_" + clustersrc_.label() + "_LadvsMod_Layer_" << i;
0595     ss2.str(std::string());
0596     ss2 << "Clusters Layer" << i << "_LadvsMod (on track)";
0597     meClPosLayersLadVsModOnTrack.push_back(iBooker.book2D(ss1.str(), ss2.str(), 11, -5.5, 5.5, ybins, ymin, ymax));
0598     meClPosLayersLadVsModOnTrack.at(i - 1)->setAxisTitle("z-module", 1);
0599     meClPosLayersLadVsModOnTrack.at(i - 1)->setAxisTitle("Ladder", 2);
0600   }
0601 
0602   for (int i = 1; i <= noOfLayers; i++) {
0603     int ybins = -1;
0604     float ymin = 0.;
0605     float ymax = 0.;
0606     if (i == 1) {
0607       ybins = 42;
0608       ymin = -10.5;
0609       ymax = 10.5;
0610     }
0611     if (i == 2) {
0612       ybins = 66;
0613       ymin = -16.5;
0614       ymax = 16.5;
0615     }
0616     if (i == 3) {
0617       ybins = 90;
0618       ymin = -22.5;
0619       ymax = 22.5;
0620     }
0621     ss1.str(std::string());
0622     ss1 << "pix_bar Occ_roc_ontrack" + digisrc_.label() + "_layer_" << i;
0623     ss2.str(std::string());
0624     ss2 << "Pixel Barrel Occupancy, ROC level (On Track): Layer " << i;
0625     meZeroRocLadvsModOnTrackBarrel.push_back(iBooker.book2D(ss1.str(), ss2.str(), 72, -4.5, 4.5, ybins, ymin, ymax));
0626     meZeroRocLadvsModOnTrackBarrel.at(i - 1)->setAxisTitle("ROC / Module", 1);
0627     meZeroRocLadvsModOnTrackBarrel.at(i - 1)->setAxisTitle("ROC / Ladder", 2);
0628   }
0629 
0630   for (int i = 1; i <= noOfLayers; i++) {
0631     ss1.str(std::string());
0632     ss1 << "nclustersvsPhi_" + clustersrc_.label() + "_Layer_" << i;
0633     ss2.str(std::string());
0634     ss2 << "nclusters (on track, layer" << i << ")";
0635     meNofClustersvsPhiOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 1400., -3.5, 3.5));
0636     meNofClustersvsPhiOnTrack_layers.at(i - 1)->setAxisTitle("Global #Phi", 1);
0637     meNofClustersvsPhiOnTrack_layers.at(i - 1)->setAxisTitle("Number of Clusters/Layer on Track", 2);
0638   }
0639 
0640   // fpix
0641   for (int i = 1; i <= noOfDisks; i++) {
0642     ss1.str(std::string());
0643     ss1 << "position_" + clustersrc_.label() + "_pz_Disk_" << i;
0644     ss2.str(std::string());
0645     ss2 << "Clusters +Z Disk" << i << " (on track)";
0646     meClPosDiskspzOnTrack.push_back(iBooker.book2D(ss1.str(), ss2.str(), 80, -20., 20., 80, -20., 20.));
0647     meClPosDiskspzOnTrack.at(i - 1)->setAxisTitle("Global X (cm)", 1);
0648     meClPosDiskspzOnTrack.at(i - 1)->setAxisTitle("Global Y (cm)", 2);
0649   }
0650   for (int i = 1; i <= noOfDisks; i++) {
0651     ss1.str(std::string());
0652     ss1 << "position_" + clustersrc_.label() + "_mz_Disk_" << i;
0653     ss2.str(std::string());
0654     ss2 << "Clusters -Z Disk" << i << " (on track)";
0655     meClPosDisksmzOnTrack.push_back(iBooker.book2D(ss1.str(), ss2.str(), 80, -20., 20., 80, -20., 20.));
0656     meClPosDisksmzOnTrack.at(i - 1)->setAxisTitle("Global X (cm)", 1);
0657     meClPosDisksmzOnTrack.at(i - 1)->setAxisTitle("Global Y (cm)", 2);
0658   }
0659   meNClustersOnTrack_all =
0660       iBooker.book1D("nclusters_" + clustersrc_.label(), "Number of Clusters (on Track)", 50, 0., 50.);
0661   meNClustersOnTrack_all->setAxisTitle("Number of Clusters", 1);
0662   meNClustersOnTrack_bpix = iBooker.book1D(
0663       "nclusters_" + clustersrc_.label() + "_Barrel", "Number of Clusters (on track, barrel)", 50, 0., 50.);
0664   meNClustersOnTrack_bpix->setAxisTitle("Number of Clusters", 1);
0665   meNClustersOnTrack_fpix = iBooker.book1D(
0666       "nclusters_" + clustersrc_.label() + "_Endcap", "Number of Clusters (on track, endcap)", 50, 0., 50.);
0667   meNClustersOnTrack_fpix->setAxisTitle("Number of Clusters", 1);
0668   for (int i = 1; i <= noOfLayers; i++) {
0669     ss1.str(std::string());
0670     ss1 << "nclusters_" + clustersrc_.label() + "_Layer_" << i;
0671     ss2.str(std::string());
0672     ss2 << "Number of Clusters (on track, layer" << i << ")";
0673     meNClustersOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0674     meNClustersOnTrack_layers.at(i - 1)->setAxisTitle("Number of Clusters", 1);
0675   }
0676   for (int i = 1; i <= noOfDisks; i++) {
0677     ss1.str(std::string());
0678     ss1 << "nclusters_" + clustersrc_.label() + "_Disk_p" << i;
0679     ss2.str(std::string());
0680     ss2 << "Number of Clusters (on track, diskp" << i << ")";
0681     meNClustersOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 50, 0., 50.));
0682     meNClustersOnTrack_diskps.at(i - 1)->setAxisTitle("Number of Clusters", 1);
0683 
0684     ss1.str(std::string());
0685     ss1 << "nclustersvsPhi_" + clustersrc_.label() + "_Disk_p" << i;
0686     ss2.str(std::string());
0687     ss2 << "nclusters (on track, diskp" << i << ")";
0688     meNofClustersvsPhiOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 1400., -3.5, 3.5));
0689     meNofClustersvsPhiOnTrack_diskps.at(i - 1)->setAxisTitle("Global #Phi", 1);
0690     meNofClustersvsPhiOnTrack_diskps.at(i - 1)->setAxisTitle("Number of Clusters/Disk on Track", 2);
0691   }
0692   for (int i = 1; i <= noOfDisks; i++) {
0693     ss1.str(std::string());
0694     ss1 << "nclusters_" + clustersrc_.label() + "_Disk_m" << i;
0695     ss2.str(std::string());
0696     ss2 << "Number of Clusters (on track, diskm" << i << ")";
0697     meNClustersOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0698     meNClustersOnTrack_diskms.at(i - 1)->setAxisTitle("Number of Clusters", 1);
0699 
0700     ss1.str(std::string());
0701     ss1 << "nclustersvsPhi_" + clustersrc_.label() + "_Disk_m" << i;
0702     ss2.str(std::string());
0703     ss2 << "nclusters (on track, diskm" << i << ")";
0704     meNofClustersvsPhiOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 1400., -3.5, 3.5));
0705     meNofClustersvsPhiOnTrack_diskms.at(i - 1)->setAxisTitle("Global #Phi", 1);
0706     meNofClustersvsPhiOnTrack_diskms.at(i - 1)->setAxisTitle("Number of Clusters/Disk on Track", 2);
0707   }
0708 
0709   meRocBladevsDiskEndcapOnTrk = iBooker.book2D(
0710       "ROC_endcap_occupancy_ontrk", "Pixel Endcap Occupancy, ROC level (On Track)", 72, -4.5, 4.5, 288, -12.5, 12.5);
0711   meRocBladevsDiskEndcapOnTrk->setBinLabel(1, "Disk-2 Pnl2", 1);
0712   meRocBladevsDiskEndcapOnTrk->setBinLabel(9, "Disk-2 Pnl1", 1);
0713   meRocBladevsDiskEndcapOnTrk->setBinLabel(19, "Disk-1 Pnl2", 1);
0714   meRocBladevsDiskEndcapOnTrk->setBinLabel(27, "Disk-1 Pnl1", 1);
0715   meRocBladevsDiskEndcapOnTrk->setBinLabel(41, "Disk+1 Pnl1", 1);
0716   meRocBladevsDiskEndcapOnTrk->setBinLabel(49, "Disk+1 Pnl2", 1);
0717   meRocBladevsDiskEndcapOnTrk->setBinLabel(59, "Disk+2 Pnl1", 1);
0718   meRocBladevsDiskEndcapOnTrk->setBinLabel(67, "Disk+2 Pnl2", 1);
0719   meRocBladevsDiskEndcapOnTrk->setAxisTitle("Blades in Inner (>0) / Outer(<) Halves", 2);
0720   meRocBladevsDiskEndcapOnTrk->setAxisTitle("ROC occupancy", 3);
0721 
0722   // not on track
0723   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OffTrack");
0724   // bpix
0725   for (int i = 1; i <= noOfLayers; i++) {
0726     ss1.str(std::string());
0727     ss1 << "position_" + clustersrc_.label() + "_Layer_" << i;
0728     ss2.str(std::string());
0729     ss2 << "Clusters Layer" << i << " (off track)";
0730     meClPosLayersNotOnTrack.push_back(iBooker.book2D(ss1.str(), ss2.str(), 200, -30., 30., 128, -3.2, 3.2));
0731     meClPosLayersNotOnTrack.at(i - 1)->setAxisTitle("Global Z (cm)", 1);
0732     meClPosLayersNotOnTrack.at(i - 1)->setAxisTitle("Global #phi", 2);
0733   }
0734   // fpix
0735   for (int i = 1; i <= noOfDisks; i++) {
0736     ss1.str(std::string());
0737     ss1 << "position_" + clustersrc_.label() + "_pz_Disk_" << i;
0738     ss2.str(std::string());
0739     ss2 << "Clusters +Z Disk" << i << " (off track)";
0740     meClPosDiskspzNotOnTrack.push_back(iBooker.book2D(ss1.str(), ss2.str(), 80, -20., 20., 80, -20., 20.));
0741     meClPosDiskspzNotOnTrack.at(i - 1)->setAxisTitle("Global X (cm)", 1);
0742     meClPosDiskspzNotOnTrack.at(i - 1)->setAxisTitle("Global Y (cm)", 2);
0743   }
0744   for (int i = 1; i <= noOfDisks; i++) {
0745     ss1.str(std::string());
0746     ss1 << "position_" + clustersrc_.label() + "_mz_Disk_" << i;
0747     ss2.str(std::string());
0748     ss2 << "Clusters -Z Disk" << i << " (off track)";
0749     meClPosDisksmzNotOnTrack.push_back(iBooker.book2D(ss1.str(), ss2.str(), 80, -20., 20., 80, -20., 20.));
0750     meClPosDisksmzNotOnTrack.at(i - 1)->setAxisTitle("Global X (cm)", 1);
0751     meClPosDisksmzNotOnTrack.at(i - 1)->setAxisTitle("Global Y (cm)", 2);
0752   }
0753   meNClustersNotOnTrack_all =
0754       iBooker.book1D("nclusters_" + clustersrc_.label(), "Number of Clusters (off Track)", 50, 0., 50.);
0755   meNClustersNotOnTrack_all->setAxisTitle("Number of Clusters", 1);
0756   meNClustersNotOnTrack_bpix = iBooker.book1D(
0757       "nclusters_" + clustersrc_.label() + "_Barrel", "Number of Clusters (off track, barrel)", 50, 0., 50.);
0758   meNClustersNotOnTrack_bpix->setAxisTitle("Number of Clusters", 1);
0759   meNClustersNotOnTrack_fpix = iBooker.book1D(
0760       "nclusters_" + clustersrc_.label() + "_Endcap", "Number of Clusters (off track, endcap)", 50, 0., 50.);
0761   meNClustersNotOnTrack_fpix->setAxisTitle("Number of Clusters", 1);
0762   for (int i = 1; i <= noOfLayers; i++) {
0763     ss1.str(std::string());
0764     ss1 << "nclusters_" + clustersrc_.label() + "_Layer_" << i;
0765     ss2.str(std::string());
0766     ss2 << "Number of Clusters (off track, layer" << i << ")";
0767     meNClustersNotOnTrack_layers.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0768     meNClustersNotOnTrack_layers.at(i - 1)->setAxisTitle("Number of Clusters", 1);
0769   }
0770   for (int i = 1; i <= noOfDisks; i++) {
0771     ss1.str(std::string());
0772     ss1 << "nclusters_" + clustersrc_.label() + "_Disk_p" << i;
0773     ss2.str(std::string());
0774     ss2 << "Number of Clusters (off track, diskp" << i << ")";
0775     meNClustersNotOnTrack_diskps.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0776     meNClustersNotOnTrack_diskps.at(i - 1)->setAxisTitle("Number of Clusters", 1);
0777   }
0778   for (int i = 1; i <= noOfDisks; i++) {
0779     ss1.str(std::string());
0780     ss1 << "nclusters_" + clustersrc_.label() + "_Disk_m" << i;
0781     ss2.str(std::string());
0782     ss2 << "Number of Clusters (off track, diskm" << i << ")";
0783     meNClustersNotOnTrack_diskms.push_back(iBooker.book1D(ss1.str(), ss2.str(), 500, 0., 500.));
0784     meNClustersNotOnTrack_diskms.at(i - 1)->setAxisTitle("Number of Clusters", 1);
0785   }
0786 
0787   meRocBladevsDiskEndcapOffTrk = iBooker.book2D(
0788       "ROC_endcap_occupancy_offtrk", "Pixel Endcap Occupancy, ROC level (Off Track)", 72, -4.5, 4.5, 288, -12.5, 12.5);
0789   meRocBladevsDiskEndcapOffTrk->setBinLabel(1, "Disk-2 Pnl2", 1);
0790   meRocBladevsDiskEndcapOffTrk->setBinLabel(9, "Disk-2 Pnl1", 1);
0791   meRocBladevsDiskEndcapOffTrk->setBinLabel(19, "Disk-1 Pnl2", 1);
0792   meRocBladevsDiskEndcapOffTrk->setBinLabel(27, "Disk-1 Pnl1", 1);
0793   meRocBladevsDiskEndcapOffTrk->setBinLabel(41, "Disk+1 Pnl1", 1);
0794   meRocBladevsDiskEndcapOffTrk->setBinLabel(49, "Disk+1 Pnl2", 1);
0795   meRocBladevsDiskEndcapOffTrk->setBinLabel(59, "Disk+2 Pnl1", 1);
0796   meRocBladevsDiskEndcapOffTrk->setBinLabel(67, "Disk+2 Pnl2", 1);
0797   meRocBladevsDiskEndcapOffTrk->setAxisTitle("Blades in Inner (>0) / Outer(<) Halves", 2);
0798   meRocBladevsDiskEndcapOffTrk->setAxisTitle("ROC occupancy", 3);
0799 
0800   // HitProbability
0801   // on track
0802   iBooker.setCurrentFolder(topFolderName_ + "/Clusters/OnTrack");
0803   meHitProbability = iBooker.book1D(
0804       "FractionLowProb", "Fraction of hits with low probability;FractionLowProb;#HitsOnTrack", 100, 0., 1.);
0805 
0806   if (debug_) {
0807     // book summary residual histograms in a debugging folder - one (x,y) pair
0808     // of histograms per subdetector
0809     iBooker.setCurrentFolder("debugging");
0810     char hisID[80];
0811     for (int s = 0; s < 3; s++) {
0812       sprintf(hisID, "residual_x_subdet_%i", s);
0813       meSubdetResidualX[s] = iBooker.book1D(hisID, "Pixel Hit-to-Track Residual in X", 500, -5., 5.);
0814 
0815       sprintf(hisID, "residual_y_subdet_%i", s);
0816       meSubdetResidualY[s] = iBooker.book1D(hisID, "Pixel Hit-to-Track Residual in Y", 500, -5., 5.);
0817     }
0818   }
0819 
0820   firstRun = false;
0821 }
0822 
0823 void SiPixelTrackResidualSource::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0824   // Retrieve tracker topology from geometry
0825   edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(trackerTopoToken_);
0826   const TrackerTopology *tTopo = tTopoHandle.product();
0827 
0828   // retrieve TrackerGeometry again and MagneticField for use in transforming
0829   // a TrackCandidate's P(ersistent)TrajectoryStateoOnDet (PTSoD) to a
0830   // TrajectoryStateOnSurface (TSoS)
0831   edm::ESHandle<TrackerGeometry> TG = iSetup.getHandle(trackerGeomToken_);
0832   const TrackerGeometry *theTrackerGeometry = TG.product();
0833 
0834   // analytic triplet method to calculate the track residuals in the pixe barrel
0835   // detector
0836 
0837   //--------------------------------------------------------------------
0838   // beam spot:
0839   //
0840   edm::Handle<reco::BeamSpot> rbs;
0841   // iEvent.getByLabel( "offlineBeamSpot", rbs );
0842   iEvent.getByToken(beamSpotToken_, rbs);
0843   math::XYZPoint bsP = math::XYZPoint(0, 0, 0);
0844   if (!rbs.failedToGet() && rbs.isValid()) {
0845     bsP = math::XYZPoint(rbs->x0(), rbs->y0(), rbs->z0());
0846   }
0847 
0848   //--------------------------------------------------------------------
0849   // primary vertices:
0850   //
0851   edm::Handle<reco::VertexCollection> vertices;
0852   // iEvent.getByLabel("offlinePrimaryVertices", vertices );
0853   iEvent.getByToken(offlinePrimaryVerticesToken_, vertices);
0854 
0855   if (vertices.failedToGet())
0856     return;
0857   if (!vertices.isValid())
0858     return;
0859 
0860   math::XYZPoint vtxN = math::XYZPoint(0, 0, 0);
0861   math::XYZPoint vtxP = math::XYZPoint(0, 0, 0);
0862 
0863   double bestNdof = 0.0;
0864   double maxSumPt = 0.0;
0865   reco::Vertex bestPvx;
0866   for (reco::VertexCollection::const_iterator iVertex = vertices->begin(); iVertex != vertices->end(); ++iVertex) {
0867     if (iVertex->ndof() > bestNdof) {
0868       bestNdof = iVertex->ndof();
0869       vtxN = math::XYZPoint(iVertex->x(), iVertex->y(), iVertex->z());
0870     }  // ndof
0871     if (iVertex->p4().pt() > maxSumPt) {
0872       maxSumPt = iVertex->p4().pt();
0873       vtxP = math::XYZPoint(iVertex->x(), iVertex->y(), iVertex->z());
0874       bestPvx = *iVertex;
0875     }  // sumpt
0876 
0877   }  // vertex
0878 
0879   if (maxSumPt < 1.0)
0880     vtxP = vtxN;
0881 
0882   //---------------------------------------------
0883   // get Tracks
0884   //
0885   edm::Handle<reco::TrackCollection> TracksForRes;
0886   // iEvent.getByLabel( "generalTracks", TracksForRes );
0887   iEvent.getByToken(generalTracksToken_, TracksForRes);
0888 
0889   if (debug_) {
0890     edm::EDConsumerBase::Labels labels;
0891     labelsForToken(generalTracksToken_, labels);
0892     std::cout << "Track for Res from " << labels.module << std::endl;
0893   }
0894 
0895   //
0896   // transient track builder, needs B-field from data base (global tag in .py)
0897   //
0898   edm::ESHandle<TransientTrackBuilder> theB = iSetup.getHandle(transientTrackBuilderToken_);
0899 
0900   // get the TransienTrackingRecHitBuilder needed for extracting the global
0901   // position of the hits in the pixel
0902   edm::ESHandle<TransientTrackingRecHitBuilder> theTrackerRecHitBuilder =
0903       iSetup.getHandle(transientTrackingRecHitBuilderToken_);
0904 
0905   // check that tracks are valid
0906   if (TracksForRes.failedToGet())
0907     return;
0908   if (!TracksForRes.isValid())
0909     return;
0910 
0911   int kk = -1;
0912   //----------------------------------------------------------------------------
0913   // Residuals:
0914   //
0915   for (reco::TrackCollection::const_iterator iTrack = TracksForRes->begin(); iTrack != TracksForRes->end(); ++iTrack) {
0916     // count
0917     kk++;
0918     // Calculate minimal track pt before curling
0919     // cpt = cqRB = 0.3*R[m]*B[T] = 1.14*R[m] for B=3.8T
0920     // D = 2R = 2*pt/1.14
0921     // calo: D = 1.3 m => pt = 0.74 GeV/c
0922     double pt = iTrack->pt();
0923     if (pt < 0.75)
0924       continue;  // curls up
0925     if (abs(iTrack->dxy(vtxP)) > 5 * iTrack->dxyError())
0926       continue;  // not prompt
0927 
0928     double charge = iTrack->charge();
0929 
0930     reco::TransientTrack tTrack = theB->build(*iTrack);
0931     // get curvature of the track, needed for the residuals
0932     double kap = tTrack.initialFreeState().transverseCurvature();
0933     // needed for the TransienTrackingRecHitBuilder
0934     TrajectoryStateOnSurface initialTSOS = tTrack.innermostMeasurementState();
0935     if (iTrack->extra().isNonnull() && iTrack->extra().isAvailable()) {
0936       double x1 = 0;
0937       double y1 = 0;
0938       double z1 = 0;
0939       double x2 = 0;
0940       double y2 = 0;
0941       double z2 = 0;
0942       double x3 = 0;
0943       double y3 = 0;
0944       double z3 = 0;
0945       int n1 = 0;
0946       int n2 = 0;
0947       int n3 = 0;
0948 
0949       // for saving the pixel barrel hits
0950       vector<TransientTrackingRecHit::RecHitPointer> GoodPixBarrelHits;
0951       // looping through the RecHits of the track
0952       for (trackingRecHit_iterator irecHit = iTrack->recHitsBegin(); irecHit != iTrack->recHitsEnd(); ++irecHit) {
0953         if ((*irecHit)->isValid()) {
0954           DetId detId = (*irecHit)->geographicalId();
0955           // enum Detector { Tracker=1, Muon=2, Ecal=3, Hcal=4, Calo=5 };
0956           if (detId.det() != 1) {
0957             if (debug_) {
0958               cout << "rec hit ID = " << detId.det() << " not in tracker!?!?\n";
0959             }
0960             continue;
0961           }
0962           uint32_t subDet = detId.subdetId();
0963 
0964           // enum SubDetector{ PixelBarrel=1, PixelEndcap=2 };
0965           // enum SubDetector{ TIB=3, TID=4, TOB=5, TEC=6 };
0966 
0967           TransientTrackingRecHit::RecHitPointer trecHit = theTrackerRecHitBuilder->build(&*(*irecHit), initialTSOS);
0968 
0969           double gX = trecHit->globalPosition().x();
0970           double gY = trecHit->globalPosition().y();
0971           double gZ = trecHit->globalPosition().z();
0972 
0973           if (subDet == PixelSubdetector::PixelBarrel) {
0974             int ilay = tTopo->pxbLayer(detId);
0975 
0976             if (ilay == 1) {
0977               n1++;
0978               x1 = gX;
0979               y1 = gY;
0980               z1 = gZ;
0981 
0982               GoodPixBarrelHits.push_back((trecHit));
0983             }  // PXB1
0984             if (ilay == 2) {
0985               n2++;
0986               x2 = gX;
0987               y2 = gY;
0988               z2 = gZ;
0989 
0990               GoodPixBarrelHits.push_back((trecHit));
0991 
0992             }  // PXB2
0993             if (ilay == 3) {
0994               n3++;
0995               x3 = gX;
0996               y3 = gY;
0997               z3 = gZ;
0998               GoodPixBarrelHits.push_back((trecHit));
0999             }
1000           }  // PXB
1001 
1002         }  // valid
1003       }    // loop rechits
1004 
1005       // CS extra plots
1006 
1007       if (n1 + n2 + n3 == 3 && n1 * n2 * n3 > 0) {
1008         for (unsigned int i = 0; i < GoodPixBarrelHits.size(); i++) {
1009           if (GoodPixBarrelHits[i]->isValid()) {
1010             DetId detId = GoodPixBarrelHits[i]->geographicalId().rawId();
1011             int ilay = tTopo->pxbLayer(detId);
1012             if (pt > ptminres_) {
1013               double dca2 = 0.0, dz2 = 0.0;
1014               double ptsig = pt;
1015               if (charge < 0.)
1016                 ptsig = -pt;
1017               // Filling the histograms in modules
1018 
1019               MeasurementPoint Test;
1020               MeasurementPoint Test2;
1021               Test = MeasurementPoint(0, 0);
1022               Test2 = MeasurementPoint(0, 0);
1023               Measurement2DVector residual;
1024 
1025               if (ilay == 1) {
1026                 triplets(x2, y2, z2, x1, y1, z1, x3, y3, z3, ptsig, dca2, dz2, kap);
1027 
1028                 Test = MeasurementPoint(dca2 * 1E4, dz2 * 1E4);
1029                 residual = Test - Test2;
1030               }
1031 
1032               if (ilay == 2) {
1033                 triplets(x1, y1, z1, x2, y2, z2, x3, y3, z3, ptsig, dca2, dz2, kap);
1034 
1035                 Test = MeasurementPoint(dca2 * 1E4, dz2 * 1E4);
1036                 residual = Test - Test2;
1037               }
1038 
1039               if (ilay == 3) {
1040                 triplets(x1, y1, z1, x3, y3, z3, x2, y2, z2, ptsig, dca2, dz2, kap);
1041 
1042                 Test = MeasurementPoint(dca2 * 1E4, dz2 * 1E4);
1043                 residual = Test - Test2;
1044               }
1045               // fill the residual histograms
1046 
1047               std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd = theSiPixelStructure.find(detId);
1048               if (pxd != theSiPixelStructure.end())
1049                 (*pxd).second->fill(residual, reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1050 
1051               if (ladOn && pxd != theSiPixelStructure.end()) {
1052                 meResidualXSummedLay.at(PixelBarrelNameUpgrade((*pxd).first).layerName() - 1)->Fill(residual.x());
1053                 meResidualYSummedLay.at(PixelBarrelNameUpgrade((*pxd).first).layerName() - 1)->Fill(residual.y());
1054               }
1055 
1056             }  // three hits
1057           }    // is valid
1058         }      // rechits loop
1059       }        // pt 4
1060     }
1061 
1062   }  //-----Tracks
1063   ////////////////////////////
1064   // get tracks
1065   edm::Handle<std::vector<reco::Track>> trackCollectionHandle;
1066   // iEvent.getByLabel(tracksrc_,trackCollectionHandle);
1067   iEvent.getByToken(trackToken_, trackCollectionHandle);
1068   auto const &trackColl = *(trackCollectionHandle.product());
1069 
1070   // get clusters
1071   edm::Handle<edmNew::DetSetVector<SiPixelCluster>> clusterColl;
1072   // iEvent.getByLabel( clustersrc_, clusterColl );
1073   iEvent.getByToken(clustersrcToken_, clusterColl);
1074   auto const &clustColl = *(clusterColl.product());
1075 
1076   // get digis
1077   edm::Handle<edm::DetSetVector<PixelDigi>> digiinput;
1078   iEvent.getByToken(digisrcToken_, digiinput);
1079   edm::DetSetVector<PixelDigi> const &diginp = *(digiinput.product());
1080 
1081   std::set<SiPixelCluster> clusterSet;
1082   TrajectoryStateCombiner tsoscomb;
1083   int tracks = 0, pixeltracks = 0, bpixtracks = 0, fpixtracks = 0;
1084   int trackclusters = 0, barreltrackclusters = 0, endcaptrackclusters = 0;
1085   int otherclusters = 0, barrelotherclusters = 0, endcapotherclusters = 0;
1086 
1087   // get trajectories
1088   edm::Handle<std::vector<Trajectory>> trajCollectionHandle;
1089   iEvent.getByToken(tracksrcToken_, trajCollectionHandle);
1090 
1091   if (debug_) {
1092     edm::EDConsumerBase::Labels labels;
1093     labelsForToken(tracksrcToken_, labels);
1094     std::cout << "Trajectories for Res from " << labels.module << std::endl;
1095   }
1096 
1097   if (trajCollectionHandle.isValid()) {
1098     auto const &trajColl = *(trajCollectionHandle.product());
1099     // get the map
1100     edm::Handle<TrajTrackAssociationCollection> match;
1101     iEvent.getByToken(trackAssociationToken_, match);
1102     auto const &ttac = *(match.product());
1103 
1104     if (debug_) {
1105       std::cout << "Trajectories\t : " << trajColl.size() << std::endl;
1106       std::cout << "recoTracks  \t : " << trackColl.size() << std::endl;
1107       std::cout << "Map entries \t : " << ttac.size() << std::endl;
1108     }
1109 
1110     // Loop over map entries
1111     for (TrajTrackAssociationCollection::const_iterator it = ttac.begin(); it != ttac.end(); ++it) {
1112       const edm::Ref<std::vector<Trajectory>> traj_iterator = it->key;
1113       // Trajectory Map, extract Trajectory for this track
1114       reco::TrackRef trackref = it->val;
1115       tracks++;
1116 
1117       bool isBpixtrack = false, isFpixtrack = false, crossesPixVol = false;
1118 
1119       // find out whether track crosses pixel fiducial volume (for cosmic
1120       // tracks)
1121 
1122       double d0 = (*trackref).d0(), dz = (*trackref).dz();
1123 
1124       if (abs(d0) < 15 && abs(dz) < 50)
1125         crossesPixVol = true;
1126 
1127       const std::vector<TrajectoryMeasurement> &tmeasColl = traj_iterator->measurements();
1128       std::vector<TrajectoryMeasurement>::const_iterator tmeasIt;
1129       // loop on measurements to find out whether there are bpix and/or fpix
1130       // hits
1131       for (tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end(); tmeasIt++) {
1132         if (!tmeasIt->updatedState().isValid())
1133           continue;
1134         TransientTrackingRecHit::ConstRecHitPointer testhit = tmeasIt->recHit();
1135         if (!testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker)
1136           continue;
1137         uint testSubDetID = (testhit->geographicalId().subdetId());
1138         if (testSubDetID == PixelSubdetector::PixelBarrel)
1139           isBpixtrack = true;
1140         if (testSubDetID == PixelSubdetector::PixelEndcap)
1141           isFpixtrack = true;
1142       }  // end loop on measurements
1143       if (isBpixtrack) {
1144         bpixtracks++;
1145         if (debug_)
1146           std::cout << "bpixtrack\n";
1147       }
1148       if (isFpixtrack) {
1149         fpixtracks++;
1150         if (debug_)
1151           std::cout << "fpixtrack\n";
1152       }
1153       if (isBpixtrack || isFpixtrack) {
1154         pixeltracks++;
1155 
1156         if (crossesPixVol)
1157           meNofTracksInPixVol_->Fill(0, 1);
1158 
1159         const std::vector<TrajectoryMeasurement> &tmeasColl = traj_iterator->measurements();
1160         for (std::vector<TrajectoryMeasurement>::const_iterator tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end();
1161              tmeasIt++) {
1162           if (!tmeasIt->updatedState().isValid())
1163             continue;
1164 
1165           TrajectoryStateOnSurface tsos = tsoscomb(tmeasIt->forwardPredictedState(), tmeasIt->backwardPredictedState());
1166           if (!tsos.isValid())
1167             continue;  // Happens rarely, due to singular matrix or similar
1168 
1169           TransientTrackingRecHit::ConstRecHitPointer hit = tmeasIt->recHit();
1170           if (!hit->isValid() || hit->geographicalId().det() != DetId::Tracker) {
1171             continue;
1172           } else {
1173             //    //residual
1174             const DetId &hit_detId = hit->geographicalId();
1175             // uint IntRawDetID = (hit_detId.rawId());
1176             uint IntSubDetID = (hit_detId.subdetId());
1177 
1178             if (IntSubDetID == 0)
1179               continue;  // don't look at SiStrip hits!
1180 
1181             // get the enclosed persistent hit
1182             const TrackingRecHit *persistentHit = hit->hit();
1183             // check if it's not null, and if it's a valid pixel hit
1184             if ((persistentHit != nullptr) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) {
1185               // tell the C++ compiler that the hit is a pixel hit
1186               const SiPixelRecHit *pixhit = static_cast<const SiPixelRecHit *>(hit->hit());
1187               // Hit probability:
1188               float hit_prob = -1.;
1189               if (pixhit->hasFilledProb()) {
1190                 hit_prob = pixhit->clusterProbability(0);
1191                 // std::cout<<"HITPROB= "<<hit_prob<<std::endl;
1192                 if (hit_prob < pow(10., -15.))
1193                   NLowProb++;
1194                 NTotal++;
1195                 if (NTotal > 0)
1196                   meHitProbability->Fill(float(NLowProb / NTotal));
1197               }
1198 
1199               // get the edm::Ref to the cluster
1200               edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = (*pixhit).cluster();
1201 
1202               //  check if the ref is not null
1203               if (clust.isNonnull()) {
1204                 // define tracker and pixel geometry and topology
1205                 const TrackerGeometry &theTracker(*theTrackerGeometry);
1206                 const PixelGeomDetUnit *theGeomDet =
1207                     static_cast<const PixelGeomDetUnit *>(theTracker.idToDet(hit_detId));
1208 
1209                 // test if PixelGeomDetUnit exists
1210                 if (theGeomDet == nullptr) {
1211                   if (debug_)
1212                     std::cout << "NO THEGEOMDET\n";
1213                   continue;
1214                 }
1215 
1216                 const PixelTopology *topol = &(theGeomDet->specificTopology());
1217                 // fill histograms for clusters on tracks
1218                 // correct SiPixelTrackResidualModule
1219                 std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd =
1220                     theSiPixelStructure.find((*hit).geographicalId().rawId());
1221 
1222                 // CHARGE CORRECTION (for track impact angle)
1223                 // calculate alpha and beta from cluster position
1224                 LocalTrajectoryParameters ltp = tsos.localParameters();
1225                 LocalVector localDir = ltp.momentum() / ltp.momentum().mag();
1226 
1227                 float clust_alpha = atan2(localDir.z(), localDir.x());
1228                 float clust_beta = atan2(localDir.z(), localDir.y());
1229                 double corrCharge = clust->charge() *
1230                                     sqrt(1.0 / (1.0 / pow(tan(clust_alpha), 2) + 1.0 / pow(tan(clust_beta), 2) + 1.0)) /
1231                                     1000.;
1232 
1233                 if (pxd != theSiPixelStructure.end())
1234                   (*pxd).second->fill(
1235                       (*clust), true, corrCharge, reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1236 
1237                 trackclusters++;
1238                 // CORR CHARGE
1239                 meClChargeOnTrack_all->Fill(corrCharge);
1240                 meClSizeOnTrack_all->Fill((*clust).size());
1241                 meClSizeXOnTrack_all->Fill((*clust).sizeX());
1242                 meClSizeYOnTrack_all->Fill((*clust).sizeY());
1243                 clusterSet.insert(*clust);
1244 
1245                 // find cluster global position (rphi, z)
1246                 // get cluster center of gravity (of charge)
1247                 float xcenter = clust->x();
1248                 float ycenter = clust->y();
1249                 // get the cluster position in local coordinates (cm)
1250                 LocalPoint clustlp = topol->localPosition(MeasurementPoint(xcenter, ycenter));
1251                 // get the cluster position in global coordinates (cm)
1252                 GlobalPoint clustgp = theGeomDet->surface().toGlobal(clustlp);
1253 
1254                 // find location of hit (barrel or endcap, same for cluster)
1255                 bool barrel =
1256                     DetId((*hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1257                 bool endcap =
1258                     DetId((*hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1259                 if (barrel) {
1260                   barreltrackclusters++;
1261 
1262                   DetId detId = (*hit).geographicalId();
1263                   if (detId >= 302055684 && detId <= 352477708) {
1264                     getrococcupancy(detId, diginp, tTopo, meZeroRocLadvsModOnTrackBarrel);
1265                   }
1266 
1267                   // CORR CHARGE
1268                   meClChargeOnTrack_bpix->Fill(corrCharge);
1269                   meClSizeOnTrack_bpix->Fill((*clust).size());
1270                   meClSizeXOnTrack_bpix->Fill((*clust).sizeX());
1271                   meClSizeYOnTrack_bpix->Fill((*clust).sizeY());
1272                   int DBlayer;
1273                   DBlayer = PixelBarrelName(DetId((*hit).geographicalId()), tTopo, isUpgrade).layerName();
1274                   float phi = clustgp.phi();
1275                   float z = clustgp.z();
1276 
1277                   PixelBarrelName pbn(DetId((*hit).geographicalId()), tTopo, isUpgrade);
1278                   int ladder = pbn.ladderName();
1279                   int module = pbn.moduleName();
1280 
1281                   PixelBarrelName::Shell sh = pbn.shell();  // enum
1282                   int ladderSigned = ladder;
1283                   int moduleSigned = module;
1284                   // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
1285                   int shell = int(sh);
1286                   // change the module sign for z<0
1287                   if (shell == 1 || shell == 2) {
1288                     moduleSigned = -module;
1289                   }
1290                   // change ladeer sign for Outer )x<0)
1291                   if (shell == 1 || shell == 3) {
1292                     ladderSigned = -ladder;
1293                   }
1294 
1295                   for (int i = 0; i < noOfLayers; i++) {
1296                     if (DBlayer == i + 1) {
1297                       meClPosLayersOnTrack.at(i)->Fill(z, phi);
1298                       meClPosLayersLadVsModOnTrack.at(i)->Fill(moduleSigned, ladderSigned);
1299                       meClChargeOnTrack_layers.at(i)->Fill(corrCharge);
1300                       meClSizeOnTrack_layers.at(i)->Fill((*clust).size());
1301                       meClSizeXOnTrack_layers.at(i)->Fill((*clust).sizeX());
1302                       meClSizeYOnTrack_layers.at(i)->Fill((*clust).sizeY());
1303                       meNofClustersvsPhiOnTrack_layers.at(i)->Fill(phi);
1304                     }
1305                   }
1306                 }
1307                 if (endcap) {
1308                   endcaptrackclusters++;
1309                   // CORR CHARGE
1310                   meClChargeOnTrack_fpix->Fill(corrCharge);
1311                   meClSizeOnTrack_fpix->Fill((*clust).size());
1312                   meClSizeXOnTrack_fpix->Fill((*clust).sizeX());
1313                   meClSizeYOnTrack_fpix->Fill((*clust).sizeY());
1314                   int DBdisk = 0;
1315                   DBdisk = PixelEndcapName(DetId((*hit).geographicalId()), tTopo, isUpgrade).diskName();
1316                   float x = clustgp.x();
1317                   float y = clustgp.y();
1318                   float z = clustgp.z();
1319                   float phi = clustgp.phi();
1320 
1321                   float xclust = clust->x();
1322                   float yclust = clust->y();
1323 
1324                   getepixrococcupancyontrk(tTopo, hit, xclust, yclust, z, meRocBladevsDiskEndcapOnTrk);
1325 
1326                   if (z > 0) {
1327                     for (int i = 0; i < noOfDisks; i++) {
1328                       if (DBdisk == i + 1) {
1329                         meClPosDiskspzOnTrack.at(i)->Fill(x, y);
1330                         meClChargeOnTrack_diskps.at(i)->Fill(corrCharge);
1331                         meClSizeOnTrack_diskps.at(i)->Fill((*clust).size());
1332                         meClSizeXOnTrack_diskps.at(i)->Fill((*clust).sizeX());
1333                         meClSizeYOnTrack_diskps.at(i)->Fill((*clust).sizeY());
1334                         meNofClustersvsPhiOnTrack_diskps.at(i)->Fill(phi);
1335                       }
1336                     }
1337                   } else {
1338                     for (int i = 0; i < noOfDisks; i++) {
1339                       if (DBdisk == i + 1) {
1340                         meClPosDisksmzOnTrack.at(i)->Fill(x, y);
1341                         meClChargeOnTrack_diskms.at(i)->Fill(corrCharge);
1342                         meClSizeOnTrack_diskms.at(i)->Fill((*clust).size());
1343                         meClSizeXOnTrack_diskms.at(i)->Fill((*clust).sizeX());
1344                         meClSizeYOnTrack_diskms.at(i)->Fill((*clust).sizeY());
1345                         meNofClustersvsPhiOnTrack_diskms.at(i)->Fill(phi);
1346                       }
1347                     }
1348                   }
1349                 }
1350 
1351               }  // end if (cluster exists)
1352 
1353             }  // end if (persistent hit exists and is pixel hit)
1354 
1355           }  // end of else
1356 
1357         }  // end for (all traj measurements of pixeltrack)
1358       }    // end if (is pixeltrack)
1359       else {
1360         if (debug_)
1361           std::cout << "no pixeltrack:\n";
1362         if (crossesPixVol)
1363           meNofTracksInPixVol_->Fill(1, 1);
1364       }
1365 
1366     }  // end loop on map entries
1367 
1368   }  // end valid trajectory:
1369 
1370   // find clusters that are NOT on track
1371   // edmNew::DetSet<SiPixelCluster>::const_iterator  di;
1372   if (debug_)
1373     std::cout << "clusters not on track: (size " << clustColl.size() << ") ";
1374 
1375   for (TrackerGeometry::DetContainer::const_iterator it = TG->dets().begin(); it != TG->dets().end(); it++) {
1376     // if(dynamic_cast<PixelGeomDetUnit const *>((*it))!=0){
1377     DetId detId = (*it)->geographicalId();
1378     if (detId >= 302055684 && detId <= 352477708) {  // make sure it's a Pixel module WITHOUT using
1379                                                      // dynamic_cast!
1380       int nofclOnTrack = 0, nofclOffTrack = 0;
1381       float z = 0.;
1382       edmNew::DetSetVector<SiPixelCluster>::const_iterator isearch = clustColl.find(detId);
1383       if (isearch != clustColl.end()) {  // Not an empty iterator
1384         edmNew::DetSet<SiPixelCluster>::const_iterator di;
1385         for (di = isearch->begin(); di != isearch->end(); di++) {
1386           unsigned int temp = clusterSet.size();
1387           clusterSet.insert(*di);
1388           // check if cluster is off track
1389           if (clusterSet.size() > temp) {
1390             otherclusters++;
1391             nofclOffTrack++;
1392             // fill histograms for clusters off tracks
1393             // correct SiPixelTrackResidualModule
1394             bool barrel = DetId(detId).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1395             if (barrel) {
1396               getrococcupancy(detId, diginp, tTopo, meZeroRocLadvsModOffTrackBarrel);
1397             }
1398 
1399             std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd =
1400                 theSiPixelStructure.find((*it)->geographicalId().rawId());
1401 
1402             if (pxd != theSiPixelStructure.end())
1403               (*pxd).second->fill((*di), false, -1., reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1404 
1405             meClSizeNotOnTrack_all->Fill((*di).size());
1406             meClSizeXNotOnTrack_all->Fill((*di).sizeX());
1407             meClSizeYNotOnTrack_all->Fill((*di).sizeY());
1408             meClChargeNotOnTrack_all->Fill((*di).charge() / 1000);
1409 
1410             /////////////////////////////////////////////////
1411             // find cluster global position (rphi, z) get cluster
1412             // define tracker and pixel geometry and topology
1413             const TrackerGeometry &theTracker(*theTrackerGeometry);
1414             const PixelGeomDetUnit *theGeomDet = static_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
1415             // test if PixelGeomDetUnit exists
1416             if (theGeomDet == nullptr) {
1417               if (debug_)
1418                 std::cout << "NO THEGEOMDET\n";
1419               continue;
1420             }
1421             const PixelTopology *topol = &(theGeomDet->specificTopology());
1422 
1423             // center of gravity (of charge)
1424             float xcenter = di->x();
1425             float ycenter = di->y();
1426             // get the cluster position in local coordinates (cm)
1427             LocalPoint clustlp = topol->localPosition(MeasurementPoint(xcenter, ycenter));
1428             // get the cluster position in global coordinates (cm)
1429             GlobalPoint clustgp = theGeomDet->surface().toGlobal(clustlp);
1430 
1431             /////////////////////////////////////////////////
1432 
1433             // barrel
1434             if (DetId(detId).subdetId() == 1) {
1435               meClSizeNotOnTrack_bpix->Fill((*di).size());
1436               meClSizeXNotOnTrack_bpix->Fill((*di).sizeX());
1437               meClSizeYNotOnTrack_bpix->Fill((*di).sizeY());
1438               meClChargeNotOnTrack_bpix->Fill((*di).charge() / 1000);
1439               barrelotherclusters++;
1440               int DBlayer = PixelBarrelName(DetId(detId), tTopo, isUpgrade).layerName();
1441               float phi = clustgp.phi();
1442               // float r = clustgp.perp();
1443               z = clustgp.z();
1444               for (int i = 0; i < noOfLayers; i++) {
1445                 if (DBlayer == i + 1) {
1446                   meClPosLayersNotOnTrack.at(i)->Fill(z, phi);
1447                   meClSizeNotOnTrack_layers.at(i)->Fill((*di).size());
1448                   meClSizeXNotOnTrack_layers.at(i)->Fill((*di).sizeX());
1449                   meClSizeYNotOnTrack_layers.at(i)->Fill((*di).sizeY());
1450                   meClChargeNotOnTrack_layers.at(i)->Fill((*di).charge() / 1000);
1451                 }
1452               }
1453             }
1454             // endcap
1455             if (DetId(detId).subdetId() == 2) {
1456               meClSizeNotOnTrack_fpix->Fill((*di).size());
1457               meClSizeXNotOnTrack_fpix->Fill((*di).sizeX());
1458               meClSizeYNotOnTrack_fpix->Fill((*di).sizeY());
1459               meClChargeNotOnTrack_fpix->Fill((*di).charge() / 1000);
1460               endcapotherclusters++;
1461               int DBdisk = PixelEndcapName(DetId(detId), tTopo, isUpgrade).diskName();
1462               float x = clustgp.x();
1463               float y = clustgp.y();
1464               z = clustgp.z();
1465 
1466               getepixrococcupancyofftrk(DetId(detId), tTopo, xcenter, ycenter, z, meRocBladevsDiskEndcapOffTrk);
1467 
1468               if (z > 0) {
1469                 for (int i = 0; i < noOfDisks; i++) {
1470                   if (DBdisk == i + 1) {
1471                     meClPosDiskspzNotOnTrack.at(i)->Fill(x, y);
1472                     meClSizeNotOnTrack_diskps.at(i)->Fill((*di).size());
1473                     meClSizeXNotOnTrack_diskps.at(i)->Fill((*di).sizeX());
1474                     meClSizeYNotOnTrack_diskps.at(i)->Fill((*di).sizeY());
1475                     meClChargeNotOnTrack_diskps.at(i)->Fill((*di).charge() / 1000);
1476                   }
1477                 }
1478               } else {
1479                 for (int i = 0; i < noOfDisks; i++) {
1480                   if (DBdisk == i + 1) {
1481                     meClPosDisksmzNotOnTrack.at(i)->Fill(x, y);
1482                     meClSizeNotOnTrack_diskms.at(i)->Fill((*di).size());
1483                     meClSizeXNotOnTrack_diskms.at(i)->Fill((*di).sizeX());
1484                     meClSizeYNotOnTrack_diskms.at(i)->Fill((*di).sizeY());
1485                     meClChargeNotOnTrack_diskms.at(i)->Fill((*di).charge() / 1000);
1486                   }
1487                 }
1488               }
1489             }
1490           }  // end "if cluster off track"
1491           else {
1492             nofclOnTrack++;
1493             if (z == 0) {
1494               // find cluster global position (rphi, z) get cluster
1495               // define tracker and pixel geometry and topology
1496               const TrackerGeometry &theTracker(*theTrackerGeometry);
1497               const PixelGeomDetUnit *theGeomDet = static_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
1498               // test if PixelGeomDetUnit exists
1499               if (theGeomDet == nullptr) {
1500                 if (debug_)
1501                   std::cout << "NO THEGEOMDET\n";
1502                 continue;
1503               }
1504               const PixelTopology *topol = &(theGeomDet->specificTopology());
1505               // center of gravity (of charge)
1506               float xcenter = di->x();
1507               float ycenter = di->y();
1508               // get the cluster position in local coordinates (cm)
1509               LocalPoint clustlp = topol->localPosition(MeasurementPoint(xcenter, ycenter));
1510               // get the cluster position in global coordinates (cm)
1511               GlobalPoint clustgp = theGeomDet->surface().toGlobal(clustlp);
1512               z = clustgp.z();
1513             }
1514           }
1515         }
1516       }
1517       //++ fill the number of clusters on a module
1518       std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd =
1519           theSiPixelStructure.find((*it)->geographicalId().rawId());
1520       if (pxd != theSiPixelStructure.end())
1521         (*pxd).second->nfill(
1522             nofclOnTrack, nofclOffTrack, reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1523       if (nofclOnTrack != 0)
1524         meNClustersOnTrack_all->Fill(nofclOnTrack);
1525       if (nofclOffTrack != 0)
1526         meNClustersNotOnTrack_all->Fill(nofclOffTrack);
1527       // barrel
1528       if (DetId(detId).subdetId() == 1) {
1529         if (nofclOnTrack != 0)
1530           meNClustersOnTrack_bpix->Fill(nofclOnTrack);
1531         if (nofclOffTrack != 0)
1532           meNClustersNotOnTrack_bpix->Fill(nofclOffTrack);
1533         int DBlayer = PixelBarrelName(DetId(detId), tTopo, isUpgrade).layerName();
1534         for (int i = 0; i < noOfLayers; i++) {
1535           if (DBlayer == i + 1) {
1536             if (nofclOnTrack != 0)
1537               meNClustersOnTrack_layers.at(i)->Fill(nofclOnTrack);
1538             if (nofclOffTrack != 0)
1539               meNClustersNotOnTrack_layers.at(i)->Fill(nofclOffTrack);
1540           }
1541         }
1542       }  // end barrel
1543       // endcap
1544       if (DetId(detId).subdetId() == 2) {
1545         int DBdisk = PixelEndcapName(DetId(detId)).diskName();
1546         // z = clustgp.z();
1547         if (nofclOnTrack != 0)
1548           meNClustersOnTrack_fpix->Fill(nofclOnTrack);
1549         if (nofclOffTrack != 0)
1550           meNClustersNotOnTrack_fpix->Fill(nofclOffTrack);
1551         if (z > 0) {
1552           for (int i = 0; i < noOfDisks; i++) {
1553             if (DBdisk == i + 1) {
1554               if (nofclOnTrack != 0)
1555                 meNClustersOnTrack_diskps.at(i)->Fill(nofclOnTrack);
1556               if (nofclOffTrack != 0)
1557                 meNClustersNotOnTrack_diskps.at(i)->Fill(nofclOffTrack);
1558             }
1559           }
1560         }
1561         if (z < 0) {
1562           for (int i = 0; i < noOfDisks; i++) {
1563             if (DBdisk == i + 1) {
1564               if (nofclOnTrack != 0)
1565                 meNClustersOnTrack_diskms.at(i)->Fill(nofclOnTrack);
1566               if (nofclOffTrack != 0)
1567                 meNClustersNotOnTrack_diskms.at(i)->Fill(nofclOffTrack);
1568             }
1569           }
1570         }
1571       }
1572 
1573     }  // end if it's a Pixel module
1574   }    // end for loop over tracker detector geometry modules
1575 
1576   if (trackclusters > 0)
1577     (meNofClustersOnTrack_)->Fill(0, trackclusters);
1578   if (barreltrackclusters > 0)
1579     (meNofClustersOnTrack_)->Fill(1, barreltrackclusters);
1580   if (endcaptrackclusters > 0)
1581     (meNofClustersOnTrack_)->Fill(2, endcaptrackclusters);
1582   if (otherclusters > 0)
1583     (meNofClustersNotOnTrack_)->Fill(0, otherclusters);
1584   if (barrelotherclusters > 0)
1585     (meNofClustersNotOnTrack_)->Fill(1, barrelotherclusters);
1586   if (endcapotherclusters > 0)
1587     (meNofClustersNotOnTrack_)->Fill(2, endcapotherclusters);
1588   if (tracks > 0)
1589     (meNofTracks_)->Fill(0, tracks);
1590   if (pixeltracks > 0)
1591     (meNofTracks_)->Fill(1, pixeltracks);
1592   if (bpixtracks > 0)
1593     (meNofTracks_)->Fill(2, bpixtracks);
1594   if (fpixtracks > 0)
1595     (meNofTracks_)->Fill(3, fpixtracks);
1596 }
1597 
1598 void SiPixelTrackResidualSource::getrococcupancy(DetId detId,
1599                                                  const edm::DetSetVector<PixelDigi> &diginp,
1600                                                  const TrackerTopology *tTopo,
1601                                                  std::vector<MonitorElement *> meinput) {
1602   edm::DetSetVector<PixelDigi>::const_iterator ipxsearch = diginp.find(detId);
1603   if (ipxsearch != diginp.end()) {
1604     // Look at digis now
1605     edm::DetSet<PixelDigi>::const_iterator pxdi;
1606     for (pxdi = ipxsearch->begin(); pxdi != ipxsearch->end(); pxdi++) {
1607       bool isHalfModule = PixelBarrelName(DetId(detId), tTopo, isUpgrade).isHalfModule();
1608       int DBlayer = PixelBarrelName(DetId(detId), tTopo, isUpgrade).layerName();
1609       int DBmodule = PixelBarrelName(DetId(detId), tTopo, isUpgrade).moduleName();
1610       int DBladder = PixelBarrelName(DetId(detId), tTopo, isUpgrade).ladderName();
1611       int DBshell = PixelBarrelName(DetId(detId), tTopo, isUpgrade).shell();
1612 
1613       // add sign to the modules
1614       if (DBshell == 1 || DBshell == 2) {
1615         DBmodule = -DBmodule;
1616       }
1617       if (DBshell == 1 || DBshell == 3) {
1618         DBladder = -DBladder;
1619       }
1620 
1621       int col = pxdi->column();
1622       int row = pxdi->row();
1623 
1624       float modsign = (float)DBmodule / (abs((float)DBmodule));
1625       float ladsign = (float)DBladder / (abs((float)DBladder));
1626       float rocx = ((float)col / (52. * 8.)) * modsign + ((float)DBmodule - (modsign)*0.5);
1627       float rocy = ((float)row / (80. * 2.)) * ladsign + ((float)DBladder - (ladsign)*0.5);
1628 
1629       // do the flip where need
1630       bool flip = false;
1631       if ((DBladder % 2 == 0) && (!isHalfModule)) {
1632         flip = true;
1633       }
1634       if ((flip) && (DBladder > 0)) {
1635         if ((((float)DBladder - (ladsign)*0.5) <= rocy) && (rocy < (float)DBladder)) {
1636           rocy = rocy + ladsign * 0.5;
1637         } else if ((((float)DBladder) <= rocy) && (rocy < ((float)DBladder + (ladsign)*0.5))) {
1638           rocy = rocy - ladsign * 0.5;
1639         }
1640       }
1641 
1642       // tweak border effect for negative modules/ladders
1643       if (modsign < 0) {
1644         rocx = rocx - 0.0001;
1645       }
1646       if (ladsign < 0) {
1647         rocy = rocy - 0.0001;
1648       } else {
1649         rocy = rocy + 0.0001;
1650       }
1651       if (abs(DBladder) == 1) {
1652         rocy = rocy + ladsign * 0.5;
1653       }  // take care of the half module
1654 
1655       if (DBlayer == 1) {
1656         meinput.at(0)->Fill(rocx, rocy);
1657       }
1658       if (DBlayer == 2) {
1659         meinput.at(1)->Fill(rocx, rocy);
1660       }
1661       if (DBlayer == 3) {
1662         meinput.at(2)->Fill(rocx, rocy);
1663       }
1664     }  // end of looping over pxdi
1665   }
1666 }
1667 
1668 void SiPixelTrackResidualSource::triplets(double x1,
1669                                           double y1,
1670                                           double z1,
1671                                           double x2,
1672                                           double y2,
1673                                           double z2,
1674                                           double x3,
1675                                           double y3,
1676                                           double z3,
1677                                           double ptsig,
1678                                           double &dca2,
1679                                           double &dz2,
1680                                           double kap) {
1681   // Define some constants
1682   using namespace std;
1683 
1684   // Curvature kap from global Track
1685 
1686   // inverse of the curvature is the radius in the transverse plane
1687   double rho = 1 / kap;
1688   // Check that the hits are in the correct layers
1689   double r1 = sqrt(x1 * x1 + y1 * y1);
1690   double r3 = sqrt(x3 * x3 + y3 * y3);
1691 
1692   if (r3 - r1 < 2.0)
1693     cout << "warn r1 = " << r1 << ", r3 = " << r3 << endl;
1694 
1695   // Calculate the centre of the helix in xy-projection with radius rho from the
1696   // track.
1697   // start with a line (sekante) connecting the two points (x1,y1) and (x3,y3)
1698   // vec_L = vec_x3-vec_x1 with L being the length of that vector.
1699   double L = sqrt((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1));
1700   // lam is the line from the middel point of vec_q towards the center of the
1701   // circle X0,Y0
1702   // we already have kap and rho = 1/kap
1703   double lam = sqrt(rho * rho - L * L / 4);
1704 
1705   // There are two solutions, the sign of kap gives the information
1706   // which of them is correct.
1707   //
1708   if (kap > 0)
1709     lam = -lam;
1710 
1711   //
1712   // ( X0, Y0 ) is the centre of the circle that describes the helix in
1713   // xy-projection.
1714   //
1715   double x0 = 0.5 * (x1 + x3) + lam / L * (-y1 + y3);
1716   double y0 = 0.5 * (y1 + y3) + lam / L * (x1 - x3);
1717 
1718   // Calculate the dipangle in z direction (needed later for z residual) :
1719   // Starting from the heliz equation whihc has to hold for both points z1,z3
1720   double num = (y3 - y0) * (x1 - x0) - (x3 - x0) * (y1 - y0);
1721   double den = (x1 - x0) * (x3 - x0) + (y1 - y0) * (y3 - y0);
1722   double tandip = kap * (z3 - z1) / atan(num / den);
1723 
1724   // angle from first hit to dca point:
1725   //
1726   double dphi = atan(((x1 - x0) * y0 - (y1 - y0) * x0) / ((x1 - x0) * x0 + (y1 - y0) * y0));
1727   // z position of the track based on the middle of the circle
1728   // track equation for the z component
1729   double uz0 = z1 + tandip * dphi * rho;
1730 
1731   /////////////////////////
1732   // RESIDUAL IN R-PHI
1733   ////////////////////////////////
1734   // Calculate distance dca2 from point (x2,y2) to the circle which is given by
1735   // the distance of the point to the middlepoint dcM = sqrt((x0-x2)^2+(y0-y2))
1736   // and rho dca = rho +- dcM
1737   if (kap > 0)
1738     dca2 = rho - sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2));
1739   else
1740     dca2 = rho + sqrt((-x0 + x2) * (-x0 + x2) + (-y0 + y2) * (-y0 + y2));
1741 
1742   /////////////////////////
1743   // RESIDUAL IN Z
1744   /////////////////////////
1745   double xx = 0;
1746   double yy = 0;
1747   // sign of kappa determines the calculation
1748   // xx and yy are the new coordinates starting from x2, y2 that are on the
1749   // track itself vec_X2+-dca2*vec(X0-X2)/|(X0-X2)|
1750   if (kap < 0) {
1751     xx = x2 + (dca2 * ((x0 - x2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1752     yy = y2 + (dca2 * ((y0 - y2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1753   } else if (kap >= 0) {
1754     xx = x2 - (dca2 * ((x0 - x2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1755     yy = y2 - (dca2 * ((y0 - y2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1756   }
1757 
1758   // to get residual in z start with calculating the new uz2 position if one has
1759   // moved to xx, yy on the track. First calculate the change in phi2 with
1760   // respect to the center X0, Y0
1761   double dphi2 = atan(((xx - x0) * y0 - (yy - y0) * x0) / ((xx - x0) * x0 + (yy - y0) * y0));
1762   // Solve track equation for this new z depending on the dip angle of the track
1763   // (see above calculated based on X1, X3 and X0, use uz0 as reference point
1764   // again.
1765   double uz2 = uz0 - dphi2 * tandip * rho;
1766 
1767   // subtract new z position from the old one
1768   dz2 = z2 - uz2;
1769 
1770   // if we are interested in the arclength this is unsigned though
1771   //  double cosphi2 = (x2*xx+y2*yy)/(sqrt(x2*x2+y2*y2)*sqrt(xx*xx+yy*yy));
1772   // double arcdca2=sqrt(x2*x2+y2*y2)*acos(cosphi2);
1773 }
1774 
1775 void SiPixelTrackResidualSource::getepixrococcupancyontrk(const TrackerTopology *tTopo,
1776                                                           TransientTrackingRecHit::ConstRecHitPointer hit,
1777                                                           float xclust,
1778                                                           float yclust,
1779                                                           float z,
1780                                                           MonitorElement *meinput) {
1781   int pxfpanel = tTopo->pxfPanel((*hit).geographicalId());
1782   int pxfmodule = tTopo->pxfModule((*hit).geographicalId());
1783   int pxfdisk = tTopo->pxfDisk((*hit).geographicalId());
1784   int pxfblade_off = tTopo->pxfBlade((*hit).geographicalId());
1785 
1786   // translate to online conventions
1787   // each EPIX disk is split in 2 half-disk - each one consists of 12 blades
1788   // in offline: blades num 0->24; here translate numbering to online convetion
1789   // positive blades pointing to beam;  negative pointing away
1790   // each blade has two panels: each consisting of an array of ROC plaquettes
1791   // front (rear) pannel: 3 (4) plaquettes
1792   // number of ROCs in each plaquette depends on the position on the panel
1793   // each ROC has 80x52 pixel cells
1794   if (z < 0.) {
1795     pxfdisk = -1. * pxfdisk;
1796   }
1797   int pxfblade = -99;
1798   if (pxfblade_off <= 6 && pxfblade_off >= 1) {
1799     pxfblade = 7 - pxfblade_off;
1800   } else if (pxfblade_off <= 18 && pxfblade_off >= 7) {
1801     pxfblade = 6 - pxfblade_off;
1802   } else if (pxfblade_off <= 24 && pxfblade_off >= 19) {
1803     pxfblade = 31 - pxfblade_off;
1804   }
1805 
1806   int clu_sdpx = ((pxfdisk > 0) ? 1 : -1) * (2 * (abs(pxfdisk) - 1) + pxfpanel);
1807   int binselx = (pxfpanel == 1 && (pxfmodule == 1 || pxfmodule == 4))
1808                     ? (pxfmodule == 1)
1809                     : ((pxfpanel == 1 && xclust < 80.0) || (pxfpanel == 2 && xclust >= 80.0));
1810   int nperpan = 2 * pxfmodule + pxfpanel - 1 + binselx;
1811   int clu_roc_binx =
1812       ((pxfdisk > 0) ? nperpan : 9 - nperpan) + (clu_sdpx + 4) * 8 - 2 * ((abs(pxfdisk) == 1) ? pxfdisk : 0);
1813 
1814   int clu_roc_biny = -99.;
1815   int nrocly = pxfmodule + pxfpanel;
1816   for (int i = 0; i < nrocly; i++) {
1817     int j = (pxfdisk < 0) ? i : nrocly - 1 - i;
1818     if (yclust >= (j * 52.0) && yclust < ((j + 1) * 52.0))
1819       clu_roc_biny = 6 - nrocly + 2 * i + ((pxfblade > 0) ? pxfblade - 1 : pxfblade + 12) * 12 + 1;
1820   }
1821   if (pxfblade > 0) {
1822     clu_roc_biny = clu_roc_biny + 144;
1823   }
1824 
1825   meinput->setBinContent(clu_roc_binx, clu_roc_biny, meinput->getBinContent(clu_roc_binx, clu_roc_biny) + 1);
1826   meinput->setBinContent(clu_roc_binx, clu_roc_biny + 1, meinput->getBinContent(clu_roc_binx, clu_roc_biny + 1) + 1);
1827 }
1828 
1829 void SiPixelTrackResidualSource::getepixrococcupancyofftrk(
1830     DetId detId, const TrackerTopology *tTopo, float xclust, float yclust, float z, MonitorElement *meinput) {
1831   PXFDetId pxfid = PXFDetId(detId);
1832 
1833   // int pxfside   = PixelEndcapName(detId,tTopo,isUpgrade).halfCylinder();
1834   int pxfpanel = pxfid.panel();
1835   int pxfmodule = pxfid.module();
1836   int pxfdisk = pxfid.disk();
1837   int pxfblade_off = pxfid.blade();
1838 
1839   if (z < 0.) {
1840     pxfdisk = -1. * pxfdisk;
1841   }
1842 
1843   int pxfblade = -99;
1844   if (pxfblade_off <= 6 && pxfblade_off >= 1) {
1845     pxfblade = 7 - pxfblade_off;
1846   } else if (pxfblade_off <= 18 && pxfblade_off >= 7) {
1847     pxfblade = 6 - pxfblade_off;
1848   } else if (pxfblade_off <= 24 && pxfblade_off >= 19) {
1849     pxfblade = 31 - pxfblade_off;
1850   }
1851 
1852   int clu_sdpx = ((pxfdisk > 0) ? 1 : -1) * (2 * (abs(pxfdisk) - 1) + pxfpanel);
1853   int binselx = (pxfpanel == 1 && (pxfmodule == 1 || pxfmodule == 4))
1854                     ? (pxfmodule == 1)
1855                     : ((pxfpanel == 1 && xclust < 80.0) || (pxfpanel == 2 && xclust >= 80.0));
1856   int nperpan = 2 * pxfmodule + pxfpanel - 1 + binselx;
1857   int clu_roc_binx =
1858       ((pxfdisk > 0) ? nperpan : 9 - nperpan) + (clu_sdpx + 4) * 8 - 2 * ((abs(pxfdisk) == 1) ? pxfdisk : 0);
1859 
1860   int clu_roc_biny = -99.;
1861   int nrocly = pxfmodule + pxfpanel;
1862   for (int i = 0; i < nrocly; i++) {
1863     int j = (pxfdisk < 0) ? i : nrocly - 1 - i;
1864     if (yclust >= (j * 52.0) && yclust < ((j + 1) * 52.0))
1865       clu_roc_biny = 6 - nrocly + 2 * i + ((pxfblade > 0) ? pxfblade - 1 : pxfblade + 12) * 12 + 1;
1866   }
1867   if (pxfblade > 0) {
1868     clu_roc_biny = clu_roc_biny + 144;
1869   }
1870 
1871   meinput->setBinContent(clu_roc_binx, clu_roc_biny, meinput->getBinContent(clu_roc_binx, clu_roc_biny) + 1);
1872   meinput->setBinContent(clu_roc_binx, clu_roc_biny + 1, meinput->getBinContent(clu_roc_binx, clu_roc_biny + 1) + 1);
1873 }
1874 
1875 DEFINE_FWK_MODULE(SiPixelTrackResidualSource);  // define this as a plug-in