Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:08:25

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   //----------------------------------------------------------------------------
0912   // Residuals:
0913   //
0914   for (reco::TrackCollection::const_iterator iTrack = TracksForRes->begin(); iTrack != TracksForRes->end(); ++iTrack) {
0915     // count
0916     // Calculate minimal track pt before curling
0917     // cpt = cqRB = 0.3*R[m]*B[T] = 1.14*R[m] for B=3.8T
0918     // D = 2R = 2*pt/1.14
0919     // calo: D = 1.3 m => pt = 0.74 GeV/c
0920     double pt = iTrack->pt();
0921     if (pt < 0.75)
0922       continue;  // curls up
0923     if (abs(iTrack->dxy(vtxP)) > 5 * iTrack->dxyError())
0924       continue;  // not prompt
0925 
0926     double charge = iTrack->charge();
0927 
0928     reco::TransientTrack tTrack = theB->build(*iTrack);
0929     // get curvature of the track, needed for the residuals
0930     double kap = tTrack.initialFreeState().transverseCurvature();
0931     // needed for the TransienTrackingRecHitBuilder
0932     TrajectoryStateOnSurface initialTSOS = tTrack.innermostMeasurementState();
0933     if (iTrack->extra().isNonnull() && iTrack->extra().isAvailable()) {
0934       double x1 = 0;
0935       double y1 = 0;
0936       double z1 = 0;
0937       double x2 = 0;
0938       double y2 = 0;
0939       double z2 = 0;
0940       double x3 = 0;
0941       double y3 = 0;
0942       double z3 = 0;
0943       int n1 = 0;
0944       int n2 = 0;
0945       int n3 = 0;
0946 
0947       // for saving the pixel barrel hits
0948       vector<TransientTrackingRecHit::RecHitPointer> GoodPixBarrelHits;
0949       // looping through the RecHits of the track
0950       for (trackingRecHit_iterator irecHit = iTrack->recHitsBegin(); irecHit != iTrack->recHitsEnd(); ++irecHit) {
0951         if ((*irecHit)->isValid()) {
0952           DetId detId = (*irecHit)->geographicalId();
0953           // enum Detector { Tracker=1, Muon=2, Ecal=3, Hcal=4, Calo=5 };
0954           if (detId.det() != 1) {
0955             if (debug_) {
0956               cout << "rec hit ID = " << detId.det() << " not in tracker!?!?\n";
0957             }
0958             continue;
0959           }
0960           uint32_t subDet = detId.subdetId();
0961 
0962           // enum SubDetector{ PixelBarrel=1, PixelEndcap=2 };
0963           // enum SubDetector{ TIB=3, TID=4, TOB=5, TEC=6 };
0964 
0965           TransientTrackingRecHit::RecHitPointer trecHit = theTrackerRecHitBuilder->build(&*(*irecHit), initialTSOS);
0966 
0967           double gX = trecHit->globalPosition().x();
0968           double gY = trecHit->globalPosition().y();
0969           double gZ = trecHit->globalPosition().z();
0970 
0971           if (subDet == PixelSubdetector::PixelBarrel) {
0972             int ilay = tTopo->pxbLayer(detId);
0973 
0974             if (ilay == 1) {
0975               n1++;
0976               x1 = gX;
0977               y1 = gY;
0978               z1 = gZ;
0979 
0980               GoodPixBarrelHits.push_back((trecHit));
0981             }  // PXB1
0982             if (ilay == 2) {
0983               n2++;
0984               x2 = gX;
0985               y2 = gY;
0986               z2 = gZ;
0987 
0988               GoodPixBarrelHits.push_back((trecHit));
0989 
0990             }  // PXB2
0991             if (ilay == 3) {
0992               n3++;
0993               x3 = gX;
0994               y3 = gY;
0995               z3 = gZ;
0996               GoodPixBarrelHits.push_back((trecHit));
0997             }
0998           }  // PXB
0999 
1000         }  // valid
1001       }    // loop rechits
1002 
1003       // CS extra plots
1004 
1005       if (n1 + n2 + n3 == 3 && n1 * n2 * n3 > 0) {
1006         for (unsigned int i = 0; i < GoodPixBarrelHits.size(); i++) {
1007           if (GoodPixBarrelHits[i]->isValid()) {
1008             DetId detId = GoodPixBarrelHits[i]->geographicalId().rawId();
1009             int ilay = tTopo->pxbLayer(detId);
1010             if (pt > ptminres_) {
1011               double dca2 = 0.0, dz2 = 0.0;
1012               double ptsig = pt;
1013               if (charge < 0.)
1014                 ptsig = -pt;
1015               // Filling the histograms in modules
1016 
1017               MeasurementPoint Test;
1018               MeasurementPoint Test2;
1019               Test = MeasurementPoint(0, 0);
1020               Test2 = MeasurementPoint(0, 0);
1021               Measurement2DVector residual;
1022 
1023               if (ilay == 1) {
1024                 triplets(x2, y2, z2, x1, y1, z1, x3, y3, z3, ptsig, dca2, dz2, kap);
1025 
1026                 Test = MeasurementPoint(dca2 * 1E4, dz2 * 1E4);
1027                 residual = Test - Test2;
1028               }
1029 
1030               if (ilay == 2) {
1031                 triplets(x1, y1, z1, x2, y2, z2, x3, y3, z3, ptsig, dca2, dz2, kap);
1032 
1033                 Test = MeasurementPoint(dca2 * 1E4, dz2 * 1E4);
1034                 residual = Test - Test2;
1035               }
1036 
1037               if (ilay == 3) {
1038                 triplets(x1, y1, z1, x3, y3, z3, x2, y2, z2, ptsig, dca2, dz2, kap);
1039 
1040                 Test = MeasurementPoint(dca2 * 1E4, dz2 * 1E4);
1041                 residual = Test - Test2;
1042               }
1043               // fill the residual histograms
1044 
1045               std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd = theSiPixelStructure.find(detId);
1046               if (pxd != theSiPixelStructure.end())
1047                 (*pxd).second->fill(residual, reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1048 
1049               if (ladOn && pxd != theSiPixelStructure.end()) {
1050                 meResidualXSummedLay.at(PixelBarrelNameUpgrade((*pxd).first).layerName() - 1)->Fill(residual.x());
1051                 meResidualYSummedLay.at(PixelBarrelNameUpgrade((*pxd).first).layerName() - 1)->Fill(residual.y());
1052               }
1053 
1054             }  // three hits
1055           }    // is valid
1056         }      // rechits loop
1057       }        // pt 4
1058     }
1059 
1060   }  //-----Tracks
1061   ////////////////////////////
1062   // get tracks
1063   edm::Handle<std::vector<reco::Track>> trackCollectionHandle;
1064   // iEvent.getByLabel(tracksrc_,trackCollectionHandle);
1065   iEvent.getByToken(trackToken_, trackCollectionHandle);
1066   auto const &trackColl = *(trackCollectionHandle.product());
1067 
1068   // get clusters
1069   edm::Handle<edmNew::DetSetVector<SiPixelCluster>> clusterColl;
1070   // iEvent.getByLabel( clustersrc_, clusterColl );
1071   iEvent.getByToken(clustersrcToken_, clusterColl);
1072   auto const &clustColl = *(clusterColl.product());
1073 
1074   // get digis
1075   edm::Handle<edm::DetSetVector<PixelDigi>> digiinput;
1076   iEvent.getByToken(digisrcToken_, digiinput);
1077   edm::DetSetVector<PixelDigi> const &diginp = *(digiinput.product());
1078 
1079   std::set<SiPixelCluster> clusterSet;
1080   TrajectoryStateCombiner tsoscomb;
1081   int tracks = 0, pixeltracks = 0, bpixtracks = 0, fpixtracks = 0;
1082   int trackclusters = 0, barreltrackclusters = 0, endcaptrackclusters = 0;
1083   int otherclusters = 0, barrelotherclusters = 0, endcapotherclusters = 0;
1084 
1085   // get trajectories
1086   edm::Handle<std::vector<Trajectory>> trajCollectionHandle;
1087   iEvent.getByToken(tracksrcToken_, trajCollectionHandle);
1088 
1089   if (debug_) {
1090     edm::EDConsumerBase::Labels labels;
1091     labelsForToken(tracksrcToken_, labels);
1092     std::cout << "Trajectories for Res from " << labels.module << std::endl;
1093   }
1094 
1095   if (trajCollectionHandle.isValid()) {
1096     auto const &trajColl = *(trajCollectionHandle.product());
1097     // get the map
1098     edm::Handle<TrajTrackAssociationCollection> match;
1099     iEvent.getByToken(trackAssociationToken_, match);
1100     auto const &ttac = *(match.product());
1101 
1102     if (debug_) {
1103       std::cout << "Trajectories\t : " << trajColl.size() << std::endl;
1104       std::cout << "recoTracks  \t : " << trackColl.size() << std::endl;
1105       std::cout << "Map entries \t : " << ttac.size() << std::endl;
1106     }
1107 
1108     // Loop over map entries
1109     for (TrajTrackAssociationCollection::const_iterator it = ttac.begin(); it != ttac.end(); ++it) {
1110       const edm::Ref<std::vector<Trajectory>> traj_iterator = it->key;
1111       // Trajectory Map, extract Trajectory for this track
1112       reco::TrackRef trackref = it->val;
1113       tracks++;
1114 
1115       bool isBpixtrack = false, isFpixtrack = false, crossesPixVol = false;
1116 
1117       // find out whether track crosses pixel fiducial volume (for cosmic
1118       // tracks)
1119 
1120       double d0 = (*trackref).d0(), dz = (*trackref).dz();
1121 
1122       if (abs(d0) < 15 && abs(dz) < 50)
1123         crossesPixVol = true;
1124 
1125       const std::vector<TrajectoryMeasurement> &tmeasColl = traj_iterator->measurements();
1126       std::vector<TrajectoryMeasurement>::const_iterator tmeasIt;
1127       // loop on measurements to find out whether there are bpix and/or fpix
1128       // hits
1129       for (tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end(); tmeasIt++) {
1130         if (!tmeasIt->updatedState().isValid())
1131           continue;
1132         TransientTrackingRecHit::ConstRecHitPointer testhit = tmeasIt->recHit();
1133         if (!testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker)
1134           continue;
1135         uint testSubDetID = (testhit->geographicalId().subdetId());
1136         if (testSubDetID == PixelSubdetector::PixelBarrel)
1137           isBpixtrack = true;
1138         if (testSubDetID == PixelSubdetector::PixelEndcap)
1139           isFpixtrack = true;
1140       }  // end loop on measurements
1141       if (isBpixtrack) {
1142         bpixtracks++;
1143         if (debug_)
1144           std::cout << "bpixtrack\n";
1145       }
1146       if (isFpixtrack) {
1147         fpixtracks++;
1148         if (debug_)
1149           std::cout << "fpixtrack\n";
1150       }
1151       if (isBpixtrack || isFpixtrack) {
1152         pixeltracks++;
1153 
1154         if (crossesPixVol)
1155           meNofTracksInPixVol_->Fill(0, 1);
1156 
1157         const std::vector<TrajectoryMeasurement> &tmeasColl = traj_iterator->measurements();
1158         for (std::vector<TrajectoryMeasurement>::const_iterator tmeasIt = tmeasColl.begin(); tmeasIt != tmeasColl.end();
1159              tmeasIt++) {
1160           if (!tmeasIt->updatedState().isValid())
1161             continue;
1162 
1163           TrajectoryStateOnSurface tsos = tsoscomb(tmeasIt->forwardPredictedState(), tmeasIt->backwardPredictedState());
1164           if (!tsos.isValid())
1165             continue;  // Happens rarely, due to singular matrix or similar
1166 
1167           TransientTrackingRecHit::ConstRecHitPointer hit = tmeasIt->recHit();
1168           if (!hit->isValid() || hit->geographicalId().det() != DetId::Tracker) {
1169             continue;
1170           } else {
1171             //    //residual
1172             const DetId &hit_detId = hit->geographicalId();
1173             // uint IntRawDetID = (hit_detId.rawId());
1174             uint IntSubDetID = (hit_detId.subdetId());
1175 
1176             if (IntSubDetID == 0)
1177               continue;  // don't look at SiStrip hits!
1178 
1179             // get the enclosed persistent hit
1180             const TrackingRecHit *persistentHit = hit->hit();
1181             // check if it's not null, and if it's a valid pixel hit
1182             if ((persistentHit != nullptr) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) {
1183               // tell the C++ compiler that the hit is a pixel hit
1184               const SiPixelRecHit *pixhit = static_cast<const SiPixelRecHit *>(hit->hit());
1185               // Hit probability:
1186               float hit_prob = -1.;
1187               if (pixhit->hasFilledProb()) {
1188                 hit_prob = pixhit->clusterProbability(0);
1189                 // std::cout<<"HITPROB= "<<hit_prob<<std::endl;
1190                 if (hit_prob < pow(10., -15.))
1191                   NLowProb++;
1192                 NTotal++;
1193                 if (NTotal > 0)
1194                   meHitProbability->Fill(float(NLowProb / NTotal));
1195               }
1196 
1197               // get the edm::Ref to the cluster
1198               edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = (*pixhit).cluster();
1199 
1200               //  check if the ref is not null
1201               if (clust.isNonnull()) {
1202                 // define tracker and pixel geometry and topology
1203                 const TrackerGeometry &theTracker(*theTrackerGeometry);
1204                 const PixelGeomDetUnit *theGeomDet =
1205                     static_cast<const PixelGeomDetUnit *>(theTracker.idToDet(hit_detId));
1206 
1207                 // test if PixelGeomDetUnit exists
1208                 if (theGeomDet == nullptr) {
1209                   if (debug_)
1210                     std::cout << "NO THEGEOMDET\n";
1211                   continue;
1212                 }
1213 
1214                 const PixelTopology *topol = &(theGeomDet->specificTopology());
1215                 // fill histograms for clusters on tracks
1216                 // correct SiPixelTrackResidualModule
1217                 std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd =
1218                     theSiPixelStructure.find((*hit).geographicalId().rawId());
1219 
1220                 // CHARGE CORRECTION (for track impact angle)
1221                 // calculate alpha and beta from cluster position
1222                 LocalTrajectoryParameters ltp = tsos.localParameters();
1223                 LocalVector localDir = ltp.momentum() / ltp.momentum().mag();
1224 
1225                 float clust_alpha = atan2(localDir.z(), localDir.x());
1226                 float clust_beta = atan2(localDir.z(), localDir.y());
1227                 double corrCharge = clust->charge() *
1228                                     sqrt(1.0 / (1.0 / pow(tan(clust_alpha), 2) + 1.0 / pow(tan(clust_beta), 2) + 1.0)) /
1229                                     1000.;
1230 
1231                 if (pxd != theSiPixelStructure.end())
1232                   (*pxd).second->fill(
1233                       (*clust), true, corrCharge, reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1234 
1235                 trackclusters++;
1236                 // CORR CHARGE
1237                 meClChargeOnTrack_all->Fill(corrCharge);
1238                 meClSizeOnTrack_all->Fill((*clust).size());
1239                 meClSizeXOnTrack_all->Fill((*clust).sizeX());
1240                 meClSizeYOnTrack_all->Fill((*clust).sizeY());
1241                 clusterSet.insert(*clust);
1242 
1243                 // find cluster global position (rphi, z)
1244                 // get cluster center of gravity (of charge)
1245                 float xcenter = clust->x();
1246                 float ycenter = clust->y();
1247                 // get the cluster position in local coordinates (cm)
1248                 LocalPoint clustlp = topol->localPosition(MeasurementPoint(xcenter, ycenter));
1249                 // get the cluster position in global coordinates (cm)
1250                 GlobalPoint clustgp = theGeomDet->surface().toGlobal(clustlp);
1251 
1252                 // find location of hit (barrel or endcap, same for cluster)
1253                 bool barrel =
1254                     DetId((*hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1255                 bool endcap =
1256                     DetId((*hit).geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelEndcap);
1257                 if (barrel) {
1258                   barreltrackclusters++;
1259 
1260                   DetId detId = (*hit).geographicalId();
1261                   if (detId >= 302055684 && detId <= 352477708) {
1262                     getrococcupancy(detId, diginp, tTopo, meZeroRocLadvsModOnTrackBarrel);
1263                   }
1264 
1265                   // CORR CHARGE
1266                   meClChargeOnTrack_bpix->Fill(corrCharge);
1267                   meClSizeOnTrack_bpix->Fill((*clust).size());
1268                   meClSizeXOnTrack_bpix->Fill((*clust).sizeX());
1269                   meClSizeYOnTrack_bpix->Fill((*clust).sizeY());
1270                   int DBlayer;
1271                   DBlayer = PixelBarrelName(DetId((*hit).geographicalId()), tTopo, isUpgrade).layerName();
1272                   float phi = clustgp.phi();
1273                   float z = clustgp.z();
1274 
1275                   PixelBarrelName pbn(DetId((*hit).geographicalId()), tTopo, isUpgrade);
1276                   int ladder = pbn.ladderName();
1277                   int module = pbn.moduleName();
1278 
1279                   PixelBarrelName::Shell sh = pbn.shell();  // enum
1280                   int ladderSigned = ladder;
1281                   int moduleSigned = module;
1282                   // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
1283                   int shell = int(sh);
1284                   // change the module sign for z<0
1285                   if (shell == 1 || shell == 2) {
1286                     moduleSigned = -module;
1287                   }
1288                   // change ladeer sign for Outer )x<0)
1289                   if (shell == 1 || shell == 3) {
1290                     ladderSigned = -ladder;
1291                   }
1292 
1293                   for (int i = 0; i < noOfLayers; i++) {
1294                     if (DBlayer == i + 1) {
1295                       meClPosLayersOnTrack.at(i)->Fill(z, phi);
1296                       meClPosLayersLadVsModOnTrack.at(i)->Fill(moduleSigned, ladderSigned);
1297                       meClChargeOnTrack_layers.at(i)->Fill(corrCharge);
1298                       meClSizeOnTrack_layers.at(i)->Fill((*clust).size());
1299                       meClSizeXOnTrack_layers.at(i)->Fill((*clust).sizeX());
1300                       meClSizeYOnTrack_layers.at(i)->Fill((*clust).sizeY());
1301                       meNofClustersvsPhiOnTrack_layers.at(i)->Fill(phi);
1302                     }
1303                   }
1304                 }
1305                 if (endcap) {
1306                   endcaptrackclusters++;
1307                   // CORR CHARGE
1308                   meClChargeOnTrack_fpix->Fill(corrCharge);
1309                   meClSizeOnTrack_fpix->Fill((*clust).size());
1310                   meClSizeXOnTrack_fpix->Fill((*clust).sizeX());
1311                   meClSizeYOnTrack_fpix->Fill((*clust).sizeY());
1312                   int DBdisk = 0;
1313                   DBdisk = PixelEndcapName(DetId((*hit).geographicalId()), tTopo, isUpgrade).diskName();
1314                   float x = clustgp.x();
1315                   float y = clustgp.y();
1316                   float z = clustgp.z();
1317                   float phi = clustgp.phi();
1318 
1319                   float xclust = clust->x();
1320                   float yclust = clust->y();
1321 
1322                   getepixrococcupancyontrk(tTopo, hit, xclust, yclust, z, meRocBladevsDiskEndcapOnTrk);
1323 
1324                   if (z > 0) {
1325                     for (int i = 0; i < noOfDisks; i++) {
1326                       if (DBdisk == i + 1) {
1327                         meClPosDiskspzOnTrack.at(i)->Fill(x, y);
1328                         meClChargeOnTrack_diskps.at(i)->Fill(corrCharge);
1329                         meClSizeOnTrack_diskps.at(i)->Fill((*clust).size());
1330                         meClSizeXOnTrack_diskps.at(i)->Fill((*clust).sizeX());
1331                         meClSizeYOnTrack_diskps.at(i)->Fill((*clust).sizeY());
1332                         meNofClustersvsPhiOnTrack_diskps.at(i)->Fill(phi);
1333                       }
1334                     }
1335                   } else {
1336                     for (int i = 0; i < noOfDisks; i++) {
1337                       if (DBdisk == i + 1) {
1338                         meClPosDisksmzOnTrack.at(i)->Fill(x, y);
1339                         meClChargeOnTrack_diskms.at(i)->Fill(corrCharge);
1340                         meClSizeOnTrack_diskms.at(i)->Fill((*clust).size());
1341                         meClSizeXOnTrack_diskms.at(i)->Fill((*clust).sizeX());
1342                         meClSizeYOnTrack_diskms.at(i)->Fill((*clust).sizeY());
1343                         meNofClustersvsPhiOnTrack_diskms.at(i)->Fill(phi);
1344                       }
1345                     }
1346                   }
1347                 }
1348 
1349               }  // end if (cluster exists)
1350 
1351             }  // end if (persistent hit exists and is pixel hit)
1352 
1353           }  // end of else
1354 
1355         }  // end for (all traj measurements of pixeltrack)
1356       }    // end if (is pixeltrack)
1357       else {
1358         if (debug_)
1359           std::cout << "no pixeltrack:\n";
1360         if (crossesPixVol)
1361           meNofTracksInPixVol_->Fill(1, 1);
1362       }
1363 
1364     }  // end loop on map entries
1365 
1366   }  // end valid trajectory:
1367 
1368   // find clusters that are NOT on track
1369   // edmNew::DetSet<SiPixelCluster>::const_iterator  di;
1370   if (debug_)
1371     std::cout << "clusters not on track: (size " << clustColl.size() << ") ";
1372 
1373   for (TrackerGeometry::DetContainer::const_iterator it = TG->dets().begin(); it != TG->dets().end(); it++) {
1374     // if(dynamic_cast<PixelGeomDetUnit const *>((*it))!=0){
1375     DetId detId = (*it)->geographicalId();
1376     if (detId >= 302055684 && detId <= 352477708) {  // make sure it's a Pixel module WITHOUT using
1377                                                      // dynamic_cast!
1378       int nofclOnTrack = 0, nofclOffTrack = 0;
1379       float z = 0.;
1380       edmNew::DetSetVector<SiPixelCluster>::const_iterator isearch = clustColl.find(detId);
1381       if (isearch != clustColl.end()) {  // Not an empty iterator
1382         edmNew::DetSet<SiPixelCluster>::const_iterator di;
1383         for (di = isearch->begin(); di != isearch->end(); di++) {
1384           unsigned int temp = clusterSet.size();
1385           clusterSet.insert(*di);
1386           // check if cluster is off track
1387           if (clusterSet.size() > temp) {
1388             otherclusters++;
1389             nofclOffTrack++;
1390             // fill histograms for clusters off tracks
1391             // correct SiPixelTrackResidualModule
1392             bool barrel = DetId(detId).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
1393             if (barrel) {
1394               getrococcupancy(detId, diginp, tTopo, meZeroRocLadvsModOffTrackBarrel);
1395             }
1396 
1397             std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd =
1398                 theSiPixelStructure.find((*it)->geographicalId().rawId());
1399 
1400             if (pxd != theSiPixelStructure.end())
1401               (*pxd).second->fill((*di), false, -1., reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1402 
1403             meClSizeNotOnTrack_all->Fill((*di).size());
1404             meClSizeXNotOnTrack_all->Fill((*di).sizeX());
1405             meClSizeYNotOnTrack_all->Fill((*di).sizeY());
1406             meClChargeNotOnTrack_all->Fill((*di).charge() / 1000);
1407 
1408             /////////////////////////////////////////////////
1409             // find cluster global position (rphi, z) get cluster
1410             // define tracker and pixel geometry and topology
1411             const TrackerGeometry &theTracker(*theTrackerGeometry);
1412             const PixelGeomDetUnit *theGeomDet = static_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
1413             // test if PixelGeomDetUnit exists
1414             if (theGeomDet == nullptr) {
1415               if (debug_)
1416                 std::cout << "NO THEGEOMDET\n";
1417               continue;
1418             }
1419             const PixelTopology *topol = &(theGeomDet->specificTopology());
1420 
1421             // center of gravity (of charge)
1422             float xcenter = di->x();
1423             float ycenter = di->y();
1424             // get the cluster position in local coordinates (cm)
1425             LocalPoint clustlp = topol->localPosition(MeasurementPoint(xcenter, ycenter));
1426             // get the cluster position in global coordinates (cm)
1427             GlobalPoint clustgp = theGeomDet->surface().toGlobal(clustlp);
1428 
1429             /////////////////////////////////////////////////
1430 
1431             // barrel
1432             if (DetId(detId).subdetId() == 1) {
1433               meClSizeNotOnTrack_bpix->Fill((*di).size());
1434               meClSizeXNotOnTrack_bpix->Fill((*di).sizeX());
1435               meClSizeYNotOnTrack_bpix->Fill((*di).sizeY());
1436               meClChargeNotOnTrack_bpix->Fill((*di).charge() / 1000);
1437               barrelotherclusters++;
1438               int DBlayer = PixelBarrelName(DetId(detId), tTopo, isUpgrade).layerName();
1439               float phi = clustgp.phi();
1440               // float r = clustgp.perp();
1441               z = clustgp.z();
1442               for (int i = 0; i < noOfLayers; i++) {
1443                 if (DBlayer == i + 1) {
1444                   meClPosLayersNotOnTrack.at(i)->Fill(z, phi);
1445                   meClSizeNotOnTrack_layers.at(i)->Fill((*di).size());
1446                   meClSizeXNotOnTrack_layers.at(i)->Fill((*di).sizeX());
1447                   meClSizeYNotOnTrack_layers.at(i)->Fill((*di).sizeY());
1448                   meClChargeNotOnTrack_layers.at(i)->Fill((*di).charge() / 1000);
1449                 }
1450               }
1451             }
1452             // endcap
1453             if (DetId(detId).subdetId() == 2) {
1454               meClSizeNotOnTrack_fpix->Fill((*di).size());
1455               meClSizeXNotOnTrack_fpix->Fill((*di).sizeX());
1456               meClSizeYNotOnTrack_fpix->Fill((*di).sizeY());
1457               meClChargeNotOnTrack_fpix->Fill((*di).charge() / 1000);
1458               endcapotherclusters++;
1459               int DBdisk = PixelEndcapName(DetId(detId), tTopo, isUpgrade).diskName();
1460               float x = clustgp.x();
1461               float y = clustgp.y();
1462               z = clustgp.z();
1463 
1464               getepixrococcupancyofftrk(DetId(detId), tTopo, xcenter, ycenter, z, meRocBladevsDiskEndcapOffTrk);
1465 
1466               if (z > 0) {
1467                 for (int i = 0; i < noOfDisks; i++) {
1468                   if (DBdisk == i + 1) {
1469                     meClPosDiskspzNotOnTrack.at(i)->Fill(x, y);
1470                     meClSizeNotOnTrack_diskps.at(i)->Fill((*di).size());
1471                     meClSizeXNotOnTrack_diskps.at(i)->Fill((*di).sizeX());
1472                     meClSizeYNotOnTrack_diskps.at(i)->Fill((*di).sizeY());
1473                     meClChargeNotOnTrack_diskps.at(i)->Fill((*di).charge() / 1000);
1474                   }
1475                 }
1476               } else {
1477                 for (int i = 0; i < noOfDisks; i++) {
1478                   if (DBdisk == i + 1) {
1479                     meClPosDisksmzNotOnTrack.at(i)->Fill(x, y);
1480                     meClSizeNotOnTrack_diskms.at(i)->Fill((*di).size());
1481                     meClSizeXNotOnTrack_diskms.at(i)->Fill((*di).sizeX());
1482                     meClSizeYNotOnTrack_diskms.at(i)->Fill((*di).sizeY());
1483                     meClChargeNotOnTrack_diskms.at(i)->Fill((*di).charge() / 1000);
1484                   }
1485                 }
1486               }
1487             }
1488           }  // end "if cluster off track"
1489           else {
1490             nofclOnTrack++;
1491             if (z == 0) {
1492               // find cluster global position (rphi, z) get cluster
1493               // define tracker and pixel geometry and topology
1494               const TrackerGeometry &theTracker(*theTrackerGeometry);
1495               const PixelGeomDetUnit *theGeomDet = static_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
1496               // test if PixelGeomDetUnit exists
1497               if (theGeomDet == nullptr) {
1498                 if (debug_)
1499                   std::cout << "NO THEGEOMDET\n";
1500                 continue;
1501               }
1502               const PixelTopology *topol = &(theGeomDet->specificTopology());
1503               // center of gravity (of charge)
1504               float xcenter = di->x();
1505               float ycenter = di->y();
1506               // get the cluster position in local coordinates (cm)
1507               LocalPoint clustlp = topol->localPosition(MeasurementPoint(xcenter, ycenter));
1508               // get the cluster position in global coordinates (cm)
1509               GlobalPoint clustgp = theGeomDet->surface().toGlobal(clustlp);
1510               z = clustgp.z();
1511             }
1512           }
1513         }
1514       }
1515       //++ fill the number of clusters on a module
1516       std::map<uint32_t, SiPixelTrackResidualModule *>::iterator pxd =
1517           theSiPixelStructure.find((*it)->geographicalId().rawId());
1518       if (pxd != theSiPixelStructure.end())
1519         (*pxd).second->nfill(
1520             nofclOnTrack, nofclOffTrack, reducedSet, modOn, ladOn, layOn, phiOn, bladeOn, diskOn, ringOn);
1521       if (nofclOnTrack != 0)
1522         meNClustersOnTrack_all->Fill(nofclOnTrack);
1523       if (nofclOffTrack != 0)
1524         meNClustersNotOnTrack_all->Fill(nofclOffTrack);
1525       // barrel
1526       if (DetId(detId).subdetId() == 1) {
1527         if (nofclOnTrack != 0)
1528           meNClustersOnTrack_bpix->Fill(nofclOnTrack);
1529         if (nofclOffTrack != 0)
1530           meNClustersNotOnTrack_bpix->Fill(nofclOffTrack);
1531         int DBlayer = PixelBarrelName(DetId(detId), tTopo, isUpgrade).layerName();
1532         for (int i = 0; i < noOfLayers; i++) {
1533           if (DBlayer == i + 1) {
1534             if (nofclOnTrack != 0)
1535               meNClustersOnTrack_layers.at(i)->Fill(nofclOnTrack);
1536             if (nofclOffTrack != 0)
1537               meNClustersNotOnTrack_layers.at(i)->Fill(nofclOffTrack);
1538           }
1539         }
1540       }  // end barrel
1541       // endcap
1542       if (DetId(detId).subdetId() == 2) {
1543         int DBdisk = PixelEndcapName(DetId(detId)).diskName();
1544         // z = clustgp.z();
1545         if (nofclOnTrack != 0)
1546           meNClustersOnTrack_fpix->Fill(nofclOnTrack);
1547         if (nofclOffTrack != 0)
1548           meNClustersNotOnTrack_fpix->Fill(nofclOffTrack);
1549         if (z > 0) {
1550           for (int i = 0; i < noOfDisks; i++) {
1551             if (DBdisk == i + 1) {
1552               if (nofclOnTrack != 0)
1553                 meNClustersOnTrack_diskps.at(i)->Fill(nofclOnTrack);
1554               if (nofclOffTrack != 0)
1555                 meNClustersNotOnTrack_diskps.at(i)->Fill(nofclOffTrack);
1556             }
1557           }
1558         }
1559         if (z < 0) {
1560           for (int i = 0; i < noOfDisks; i++) {
1561             if (DBdisk == i + 1) {
1562               if (nofclOnTrack != 0)
1563                 meNClustersOnTrack_diskms.at(i)->Fill(nofclOnTrack);
1564               if (nofclOffTrack != 0)
1565                 meNClustersNotOnTrack_diskms.at(i)->Fill(nofclOffTrack);
1566             }
1567           }
1568         }
1569       }
1570 
1571     }  // end if it's a Pixel module
1572   }    // end for loop over tracker detector geometry modules
1573 
1574   if (trackclusters > 0)
1575     (meNofClustersOnTrack_)->Fill(0, trackclusters);
1576   if (barreltrackclusters > 0)
1577     (meNofClustersOnTrack_)->Fill(1, barreltrackclusters);
1578   if (endcaptrackclusters > 0)
1579     (meNofClustersOnTrack_)->Fill(2, endcaptrackclusters);
1580   if (otherclusters > 0)
1581     (meNofClustersNotOnTrack_)->Fill(0, otherclusters);
1582   if (barrelotherclusters > 0)
1583     (meNofClustersNotOnTrack_)->Fill(1, barrelotherclusters);
1584   if (endcapotherclusters > 0)
1585     (meNofClustersNotOnTrack_)->Fill(2, endcapotherclusters);
1586   if (tracks > 0)
1587     (meNofTracks_)->Fill(0, tracks);
1588   if (pixeltracks > 0)
1589     (meNofTracks_)->Fill(1, pixeltracks);
1590   if (bpixtracks > 0)
1591     (meNofTracks_)->Fill(2, bpixtracks);
1592   if (fpixtracks > 0)
1593     (meNofTracks_)->Fill(3, fpixtracks);
1594 }
1595 
1596 void SiPixelTrackResidualSource::getrococcupancy(DetId detId,
1597                                                  const edm::DetSetVector<PixelDigi> &diginp,
1598                                                  const TrackerTopology *tTopo,
1599                                                  std::vector<MonitorElement *> meinput) {
1600   edm::DetSetVector<PixelDigi>::const_iterator ipxsearch = diginp.find(detId);
1601   if (ipxsearch != diginp.end()) {
1602     // Look at digis now
1603     edm::DetSet<PixelDigi>::const_iterator pxdi;
1604     for (pxdi = ipxsearch->begin(); pxdi != ipxsearch->end(); pxdi++) {
1605       bool isHalfModule = PixelBarrelName(DetId(detId), tTopo, isUpgrade).isHalfModule();
1606       int DBlayer = PixelBarrelName(DetId(detId), tTopo, isUpgrade).layerName();
1607       int DBmodule = PixelBarrelName(DetId(detId), tTopo, isUpgrade).moduleName();
1608       int DBladder = PixelBarrelName(DetId(detId), tTopo, isUpgrade).ladderName();
1609       int DBshell = PixelBarrelName(DetId(detId), tTopo, isUpgrade).shell();
1610 
1611       // add sign to the modules
1612       if (DBshell == 1 || DBshell == 2) {
1613         DBmodule = -DBmodule;
1614       }
1615       if (DBshell == 1 || DBshell == 3) {
1616         DBladder = -DBladder;
1617       }
1618 
1619       int col = pxdi->column();
1620       int row = pxdi->row();
1621 
1622       float modsign = (float)DBmodule / (abs((float)DBmodule));
1623       float ladsign = (float)DBladder / (abs((float)DBladder));
1624       float rocx = ((float)col / (52. * 8.)) * modsign + ((float)DBmodule - (modsign)*0.5);
1625       float rocy = ((float)row / (80. * 2.)) * ladsign + ((float)DBladder - (ladsign)*0.5);
1626 
1627       // do the flip where need
1628       bool flip = false;
1629       if ((DBladder % 2 == 0) && (!isHalfModule)) {
1630         flip = true;
1631       }
1632       if ((flip) && (DBladder > 0)) {
1633         if ((((float)DBladder - (ladsign)*0.5) <= rocy) && (rocy < (float)DBladder)) {
1634           rocy = rocy + ladsign * 0.5;
1635         } else if ((((float)DBladder) <= rocy) && (rocy < ((float)DBladder + (ladsign)*0.5))) {
1636           rocy = rocy - ladsign * 0.5;
1637         }
1638       }
1639 
1640       // tweak border effect for negative modules/ladders
1641       if (modsign < 0) {
1642         rocx = rocx - 0.0001;
1643       }
1644       if (ladsign < 0) {
1645         rocy = rocy - 0.0001;
1646       } else {
1647         rocy = rocy + 0.0001;
1648       }
1649       if (abs(DBladder) == 1) {
1650         rocy = rocy + ladsign * 0.5;
1651       }  // take care of the half module
1652 
1653       if (DBlayer == 1) {
1654         meinput.at(0)->Fill(rocx, rocy);
1655       }
1656       if (DBlayer == 2) {
1657         meinput.at(1)->Fill(rocx, rocy);
1658       }
1659       if (DBlayer == 3) {
1660         meinput.at(2)->Fill(rocx, rocy);
1661       }
1662     }  // end of looping over pxdi
1663   }
1664 }
1665 
1666 void SiPixelTrackResidualSource::triplets(double x1,
1667                                           double y1,
1668                                           double z1,
1669                                           double x2,
1670                                           double y2,
1671                                           double z2,
1672                                           double x3,
1673                                           double y3,
1674                                           double z3,
1675                                           double ptsig,
1676                                           double &dca2,
1677                                           double &dz2,
1678                                           double kap) {
1679   // Define some constants
1680   using namespace std;
1681 
1682   // Curvature kap from global Track
1683 
1684   // inverse of the curvature is the radius in the transverse plane
1685   double rho = 1 / kap;
1686   // Check that the hits are in the correct layers
1687   double r1 = sqrt(x1 * x1 + y1 * y1);
1688   double r3 = sqrt(x3 * x3 + y3 * y3);
1689 
1690   if (r3 - r1 < 2.0)
1691     cout << "warn r1 = " << r1 << ", r3 = " << r3 << endl;
1692 
1693   // Calculate the centre of the helix in xy-projection with radius rho from the
1694   // track.
1695   // start with a line (sekante) connecting the two points (x1,y1) and (x3,y3)
1696   // vec_L = vec_x3-vec_x1 with L being the length of that vector.
1697   double L = sqrt((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1));
1698   // lam is the line from the middel point of vec_q towards the center of the
1699   // circle X0,Y0
1700   // we already have kap and rho = 1/kap
1701   double lam = sqrt(rho * rho - L * L / 4);
1702 
1703   // There are two solutions, the sign of kap gives the information
1704   // which of them is correct.
1705   //
1706   if (kap > 0)
1707     lam = -lam;
1708 
1709   //
1710   // ( X0, Y0 ) is the centre of the circle that describes the helix in
1711   // xy-projection.
1712   //
1713   double x0 = 0.5 * (x1 + x3) + lam / L * (-y1 + y3);
1714   double y0 = 0.5 * (y1 + y3) + lam / L * (x1 - x3);
1715 
1716   // Calculate the dipangle in z direction (needed later for z residual) :
1717   // Starting from the heliz equation whihc has to hold for both points z1,z3
1718   double num = (y3 - y0) * (x1 - x0) - (x3 - x0) * (y1 - y0);
1719   double den = (x1 - x0) * (x3 - x0) + (y1 - y0) * (y3 - y0);
1720   double tandip = kap * (z3 - z1) / atan(num / den);
1721 
1722   // angle from first hit to dca point:
1723   //
1724   double dphi = atan(((x1 - x0) * y0 - (y1 - y0) * x0) / ((x1 - x0) * x0 + (y1 - y0) * y0));
1725   // z position of the track based on the middle of the circle
1726   // track equation for the z component
1727   double uz0 = z1 + tandip * dphi * rho;
1728 
1729   /////////////////////////
1730   // RESIDUAL IN R-PHI
1731   ////////////////////////////////
1732   // Calculate distance dca2 from point (x2,y2) to the circle which is given by
1733   // the distance of the point to the middlepoint dcM = sqrt((x0-x2)^2+(y0-y2))
1734   // and rho dca = rho +- dcM
1735   if (kap > 0)
1736     dca2 = rho - sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2));
1737   else
1738     dca2 = rho + sqrt((-x0 + x2) * (-x0 + x2) + (-y0 + y2) * (-y0 + y2));
1739 
1740   /////////////////////////
1741   // RESIDUAL IN Z
1742   /////////////////////////
1743   double xx = 0;
1744   double yy = 0;
1745   // sign of kappa determines the calculation
1746   // xx and yy are the new coordinates starting from x2, y2 that are on the
1747   // track itself vec_X2+-dca2*vec(X0-X2)/|(X0-X2)|
1748   if (kap < 0) {
1749     xx = x2 + (dca2 * ((x0 - x2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1750     yy = y2 + (dca2 * ((y0 - y2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1751   } else if (kap >= 0) {
1752     xx = x2 - (dca2 * ((x0 - x2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1753     yy = y2 - (dca2 * ((y0 - y2)) / sqrt((x0 - x2) * (x0 - x2) + (y0 - y2) * (y0 - y2)));
1754   }
1755 
1756   // to get residual in z start with calculating the new uz2 position if one has
1757   // moved to xx, yy on the track. First calculate the change in phi2 with
1758   // respect to the center X0, Y0
1759   double dphi2 = atan(((xx - x0) * y0 - (yy - y0) * x0) / ((xx - x0) * x0 + (yy - y0) * y0));
1760   // Solve track equation for this new z depending on the dip angle of the track
1761   // (see above calculated based on X1, X3 and X0, use uz0 as reference point
1762   // again.
1763   double uz2 = uz0 - dphi2 * tandip * rho;
1764 
1765   // subtract new z position from the old one
1766   dz2 = z2 - uz2;
1767 
1768   // if we are interested in the arclength this is unsigned though
1769   //  double cosphi2 = (x2*xx+y2*yy)/(sqrt(x2*x2+y2*y2)*sqrt(xx*xx+yy*yy));
1770   // double arcdca2=sqrt(x2*x2+y2*y2)*acos(cosphi2);
1771 }
1772 
1773 void SiPixelTrackResidualSource::getepixrococcupancyontrk(const TrackerTopology *tTopo,
1774                                                           TransientTrackingRecHit::ConstRecHitPointer hit,
1775                                                           float xclust,
1776                                                           float yclust,
1777                                                           float z,
1778                                                           MonitorElement *meinput) {
1779   int pxfpanel = tTopo->pxfPanel((*hit).geographicalId());
1780   int pxfmodule = tTopo->pxfModule((*hit).geographicalId());
1781   int pxfdisk = tTopo->pxfDisk((*hit).geographicalId());
1782   int pxfblade_off = tTopo->pxfBlade((*hit).geographicalId());
1783 
1784   // translate to online conventions
1785   // each EPIX disk is split in 2 half-disk - each one consists of 12 blades
1786   // in offline: blades num 0->24; here translate numbering to online convetion
1787   // positive blades pointing to beam;  negative pointing away
1788   // each blade has two panels: each consisting of an array of ROC plaquettes
1789   // front (rear) pannel: 3 (4) plaquettes
1790   // number of ROCs in each plaquette depends on the position on the panel
1791   // each ROC has 80x52 pixel cells
1792   if (z < 0.) {
1793     pxfdisk = -1. * pxfdisk;
1794   }
1795   int pxfblade = -99;
1796   if (pxfblade_off <= 6 && pxfblade_off >= 1) {
1797     pxfblade = 7 - pxfblade_off;
1798   } else if (pxfblade_off <= 18 && pxfblade_off >= 7) {
1799     pxfblade = 6 - pxfblade_off;
1800   } else if (pxfblade_off <= 24 && pxfblade_off >= 19) {
1801     pxfblade = 31 - pxfblade_off;
1802   }
1803 
1804   int clu_sdpx = ((pxfdisk > 0) ? 1 : -1) * (2 * (abs(pxfdisk) - 1) + pxfpanel);
1805   int binselx = (pxfpanel == 1 && (pxfmodule == 1 || pxfmodule == 4))
1806                     ? (pxfmodule == 1)
1807                     : ((pxfpanel == 1 && xclust < 80.0) || (pxfpanel == 2 && xclust >= 80.0));
1808   int nperpan = 2 * pxfmodule + pxfpanel - 1 + binselx;
1809   int clu_roc_binx =
1810       ((pxfdisk > 0) ? nperpan : 9 - nperpan) + (clu_sdpx + 4) * 8 - 2 * ((abs(pxfdisk) == 1) ? pxfdisk : 0);
1811 
1812   int clu_roc_biny = -99.;
1813   int nrocly = pxfmodule + pxfpanel;
1814   for (int i = 0; i < nrocly; i++) {
1815     int j = (pxfdisk < 0) ? i : nrocly - 1 - i;
1816     if (yclust >= (j * 52.0) && yclust < ((j + 1) * 52.0))
1817       clu_roc_biny = 6 - nrocly + 2 * i + ((pxfblade > 0) ? pxfblade - 1 : pxfblade + 12) * 12 + 1;
1818   }
1819   if (pxfblade > 0) {
1820     clu_roc_biny = clu_roc_biny + 144;
1821   }
1822 
1823   meinput->setBinContent(clu_roc_binx, clu_roc_biny, meinput->getBinContent(clu_roc_binx, clu_roc_biny) + 1);
1824   meinput->setBinContent(clu_roc_binx, clu_roc_biny + 1, meinput->getBinContent(clu_roc_binx, clu_roc_biny + 1) + 1);
1825 }
1826 
1827 void SiPixelTrackResidualSource::getepixrococcupancyofftrk(
1828     DetId detId, const TrackerTopology *tTopo, float xclust, float yclust, float z, MonitorElement *meinput) {
1829   PXFDetId pxfid = PXFDetId(detId);
1830 
1831   // int pxfside   = PixelEndcapName(detId,tTopo,isUpgrade).halfCylinder();
1832   int pxfpanel = pxfid.panel();
1833   int pxfmodule = pxfid.module();
1834   int pxfdisk = pxfid.disk();
1835   int pxfblade_off = pxfid.blade();
1836 
1837   if (z < 0.) {
1838     pxfdisk = -1. * pxfdisk;
1839   }
1840 
1841   int pxfblade = -99;
1842   if (pxfblade_off <= 6 && pxfblade_off >= 1) {
1843     pxfblade = 7 - pxfblade_off;
1844   } else if (pxfblade_off <= 18 && pxfblade_off >= 7) {
1845     pxfblade = 6 - pxfblade_off;
1846   } else if (pxfblade_off <= 24 && pxfblade_off >= 19) {
1847     pxfblade = 31 - pxfblade_off;
1848   }
1849 
1850   int clu_sdpx = ((pxfdisk > 0) ? 1 : -1) * (2 * (abs(pxfdisk) - 1) + pxfpanel);
1851   int binselx = (pxfpanel == 1 && (pxfmodule == 1 || pxfmodule == 4))
1852                     ? (pxfmodule == 1)
1853                     : ((pxfpanel == 1 && xclust < 80.0) || (pxfpanel == 2 && xclust >= 80.0));
1854   int nperpan = 2 * pxfmodule + pxfpanel - 1 + binselx;
1855   int clu_roc_binx =
1856       ((pxfdisk > 0) ? nperpan : 9 - nperpan) + (clu_sdpx + 4) * 8 - 2 * ((abs(pxfdisk) == 1) ? pxfdisk : 0);
1857 
1858   int clu_roc_biny = -99.;
1859   int nrocly = pxfmodule + pxfpanel;
1860   for (int i = 0; i < nrocly; i++) {
1861     int j = (pxfdisk < 0) ? i : nrocly - 1 - i;
1862     if (yclust >= (j * 52.0) && yclust < ((j + 1) * 52.0))
1863       clu_roc_biny = 6 - nrocly + 2 * i + ((pxfblade > 0) ? pxfblade - 1 : pxfblade + 12) * 12 + 1;
1864   }
1865   if (pxfblade > 0) {
1866     clu_roc_biny = clu_roc_biny + 144;
1867   }
1868 
1869   meinput->setBinContent(clu_roc_binx, clu_roc_biny, meinput->getBinContent(clu_roc_binx, clu_roc_biny) + 1);
1870   meinput->setBinContent(clu_roc_binx, clu_roc_biny + 1, meinput->getBinContent(clu_roc_binx, clu_roc_biny + 1) + 1);
1871 }
1872 
1873 DEFINE_FWK_MODULE(SiPixelTrackResidualSource);  // define this as a plug-in