Phase2TrackerRecHits/test/RecHitsValidation.cc

0001 // system includes
0002 #include <iostream>
0003 #include <map>
0004 #include <vector>
0005 #include <algorithm>
0006
0007 // user file includes
0008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0009 #include "CommonTools/Utils/interface/TFileDirectory.h"
0010 #include "DataFormats/Common/interface/DetSetVector.h"
0011 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "DataFormats/DetId/interface/DetId.h"
0014 #include "DataFormats/GeometrySurface/interface/LocalError.h"
0015 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0016 #include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h"
0017 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0018 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0019 #include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 #include "FWCore/ServiceRegistry/interface/Service.h"
0028 #include "FWCore/Utilities/interface/InputTag.h"
0029 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0030 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0031 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0032 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0033 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0034 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0035 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0036 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0037
0038 // ROOT includes
0039 #include <TH2F.h>
0040 #include <TH1F.h>
0041
0042 struct RecHitHistos {
0043   // use TH1D instead of TH1F to avoid stauration at 2^31
0044   // above this increments with +1 don't work for float, need double
0045
0046   TH1D* numberRecHits[3];
0047   TH1D* clusterSize[3];
0048
0049   TH2F* globalPosXY[3][5];
0050   TH2F* localPosXY[3][5];
0051
0052   TH1F* deltaX[3][5];
0053   TH1F* deltaY[3][5];
0054   TH1F* deltaX_P[3][5];
0055   TH1F* deltaY_P[3][5];
0056
0057   TH1F* pullX[3][5];
0058   TH1F* pullY[3][5];
0059   TH1F* pullX_P[3][5];
0060   TH1F* pullY_P[3][5];
0061
0062   TH2F* deltaX_eta[3][5];
0063   TH2F* deltaY_eta[3][5];
0064   TH2F* deltaX_eta_P[3][5];
0065   TH2F* deltaY_eta_P[3][5];
0066
0067   TH2F* pullX_eta[3][5];
0068   TH2F* pullY_eta[3][5];
0069   TH2F* pullX_eta_P[3][5];
0070   TH2F* pullY_eta_P[3][5];
0071
0072   TH1D* primarySimHits[3];
0073   TH1D* otherSimHits[3];
0074 };
0075
0076 class Phase2TrackerRecHitsValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0077 public:
0078   typedef std::map<unsigned int, std::vector<PSimHit> > SimHitsMap;
0079   typedef std::map<unsigned int, SimTrack> SimTracksMap;
0080
0081   explicit Phase2TrackerRecHitsValidation(const edm::ParameterSet&);
0082   ~Phase2TrackerRecHitsValidation();
0083   void beginJob() override;
0084   void analyze(const edm::Event&, const edm::EventSetup&) override;
0085
0086 private:
0087   std::map<unsigned int, RecHitHistos>::iterator createLayerHistograms(unsigned int);
0088   std::vector<unsigned int> getSimTrackId(const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >&,
0089                                           const DetId&,
0090                                           unsigned int);
0091
0092   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> esTokenGeom_;
0093   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> esTokenTopo_;
0094   const edm::EDGetTokenT<Phase2TrackerRecHit1DCollectionNew> tokenRecHits_;
0095   const edm::EDGetTokenT<Phase2TrackerCluster1DCollectionNew> tokenClusters_;
0096   const edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > tokenLinks_;
0097   const edm::EDGetTokenT<edm::PSimHitContainer> tokenSimHitsB_;
0098   const edm::EDGetTokenT<edm::PSimHitContainer> tokenSimHitsE_;
0099   const edm::EDGetTokenT<edm::SimTrackContainer> tokenSimTracks_;
0100
0101   const bool catECasRings_;
0102   const double simtrackminpt_;
0103   const bool makeEtaPlots_;
0104   const double mineta_;
0105   const double maxeta_;
0106
0107   TH2F* trackerLayout_;
0108   TH2F* trackerLayoutXY_;
0109   TH2F* trackerLayoutXYBar_;
0110   TH2F* trackerLayoutXYEC_;
0111
0112   std::map<unsigned int, RecHitHistos> histograms_;
0113 };
0114
0115 Phase2TrackerRecHitsValidation::Phase2TrackerRecHitsValidation(const edm::ParameterSet& conf)
0116     : esTokenGeom_(esConsumes()),
0117       esTokenTopo_(esConsumes()),
0118       tokenRecHits_(consumes<Phase2TrackerRecHit1DCollectionNew>(conf.getParameter<edm::InputTag>("src"))),
0119       tokenClusters_(consumes<Phase2TrackerCluster1DCollectionNew>(conf.getParameter<edm::InputTag>("clusters"))),
0120       tokenLinks_(consumes<edm::DetSetVector<PixelDigiSimLink> >(conf.getParameter<edm::InputTag>("links"))),
0121       tokenSimHitsB_(consumes<edm::PSimHitContainer>(conf.getParameter<edm::InputTag>("simhitsbarrel"))),
0122       tokenSimHitsE_(consumes<edm::PSimHitContainer>(conf.getParameter<edm::InputTag>("simhitsendcap"))),
0123       tokenSimTracks_(consumes<edm::SimTrackContainer>(conf.getParameter<edm::InputTag>("simtracks"))),
0124       catECasRings_(conf.getParameter<bool>("ECasRings")),
0125       simtrackminpt_(conf.getParameter<double>("SimTrackMinPt")),
0126       makeEtaPlots_(conf.getParameter<bool>("MakeEtaPlots")),
0127       mineta_(conf.getParameter<double>("MinEta")),
0128       maxeta_(conf.getParameter<double>("MaxEta")) {
0129   usesResource(TFileService::kSharedResource);
0130 }
0131
0132 Phase2TrackerRecHitsValidation::~Phase2TrackerRecHitsValidation() = default;
0133
0134 void Phase2TrackerRecHitsValidation::beginJob() {
0135   edm::Service<TFileService> fs;
0136   fs->file().cd("/");
0137   TFileDirectory td = fs->mkdir("Common");
0138   // Create common histograms
0139   trackerLayout_ = td.make<TH2F>("RVsZ", "R vs. z position", 6000, -300.0, 300.0, 1200, 0.0, 120.0);
0140   trackerLayoutXY_ = td.make<TH2F>("XVsY", "x vs. y position", 2400, -120.0, 120.0, 2400, -120.0, 120.0);
0141   trackerLayoutXYBar_ = td.make<TH2F>("XVsYBar", "x vs. y position", 2400, -120.0, 120.0, 2400, -120.0, 120.0);
0142   trackerLayoutXYEC_ = td.make<TH2F>("XVsYEC", "x vs. y position", 2400, -120.0, 120.0, 2400, -120.0, 120.0);
0143 }
0144
0145 void Phase2TrackerRecHitsValidation::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0146   /*
0147      * Get the needed objects
0148      */
0149
0150   // Get the RecHits
0151   edm::Handle<Phase2TrackerRecHit1DCollectionNew> rechits;
0152   event.getByToken(tokenRecHits_, rechits);
0153
0154   // Get the Clusters
0155   edm::Handle<Phase2TrackerCluster1DCollectionNew> clusters;
0156   event.getByToken(tokenClusters_, clusters);
0157
0158   // Get the PixelDigiSimLinks
0159   edm::Handle<edm::DetSetVector<PixelDigiSimLink> > pixelSimLinks;
0160   event.getByToken(tokenLinks_, pixelSimLinks);
0161
0162   // Get the SimHits
0163   edm::Handle<edm::PSimHitContainer> simHitsRaw[2];
0164   event.getByToken(tokenSimHitsB_, simHitsRaw[0]);
0165   event.getByToken(tokenSimHitsE_, simHitsRaw[1]);
0166
0167   // Get the SimTracks
0168   edm::Handle<edm::SimTrackContainer> simTracksRaw;
0169   event.getByToken(tokenSimTracks_, simTracksRaw);
0170
0171   // Get the geometry
0172   const TrackerGeometry* tkGeom = &eventSetup.getData(esTokenGeom_);
0173   const TrackerTopology* tTopo = &eventSetup.getData(esTokenTopo_);
0174
0175   /*
0176      * Rearrange the simTracks
0177      */
0178
0179   // Rearrange the simTracks for ease of use <simTrackID, simTrack>
0180   SimTracksMap simTracks;
0181   for (edm::SimTrackContainer::const_iterator simTrackIt(simTracksRaw->begin()); simTrackIt != simTracksRaw->end();
0182        ++simTrackIt) {
0183     if (simTrackIt->momentum().pt() > simtrackminpt_) {
0184       simTracks.insert(std::pair<unsigned int, SimTrack>(simTrackIt->trackId(), *simTrackIt));
0185     }
0186   }
0187
0188   /*
0189      * Validation
0190      */
0191
0192   // Number of rechits
0193   std::map<unsigned int, unsigned int> nRecHits[3];
0194   std::map<unsigned int, unsigned int> nPrimarySimHits[3];
0195   std::map<unsigned int, unsigned int> nOtherSimHits[3];
0196
0197   // Loop over modules
0198   for (Phase2TrackerRecHit1DCollectionNew::const_iterator DSViter = rechits->begin(); DSViter != rechits->end();
0199        ++DSViter) {
0200     // Get the detector unit's id
0201     unsigned int rawid(DSViter->detId());
0202     DetId detId(rawid);
0203     unsigned int layer = (tTopo->side(detId) != 0) * 1000;  // don't split up endcap sides
0204     if (!layer) {
0205       layer += tTopo->layer(detId);
0206     } else {
0207       layer += (catECasRings_ ? tTopo->tidRing(detId) * 10 : tTopo->layer(detId));
0208     }
0209
0210     // determine the detector we are in
0211     TrackerGeometry::ModuleType mType = tkGeom->getDetectorType(detId);
0212     unsigned int det = 0;
0213     if (mType == TrackerGeometry::ModuleType::Ph2PSP) {
0214       det = 1;
0215     } else if (mType == TrackerGeometry::ModuleType::Ph2PSS || mType == TrackerGeometry::ModuleType::Ph2SS) {
0216       det = 2;
0217     } else {
0218       std::cout << "UNKNOWN DETECTOR TYPE!" << std::endl;
0219     }
0220
0221     // Get the geomdet
0222     const GeomDetUnit* geomDetUnit(tkGeom->idToDetUnit(detId));
0223     if (!geomDetUnit)
0224       continue;
0225
0226     // initialize the nhit counters if they don't exist for this layer
0227     auto nhitit(nRecHits[det].find(layer));
0228     if (nhitit == nRecHits[det].end()) {
0229       nRecHits[det].emplace(layer, 0);
0230       nPrimarySimHits[det].emplace(layer, 0);
0231       nOtherSimHits[det].emplace(layer, 0);
0232     }
0233
0234     // Create histograms if they do not yet exist for this layer
0235     std::map<unsigned int, RecHitHistos>::iterator histogramLayer(histograms_.find(layer));
0236     if (histogramLayer == histograms_.end())
0237       histogramLayer = createLayerHistograms(layer);
0238
0239     // Loop over the rechits in the detector unit
0240     for (edmNew::DetSet<Phase2TrackerRecHit1D>::const_iterator rechitIt = DSViter->begin(); rechitIt != DSViter->end();
0241          ++rechitIt) {
0242       // determine the position
0243       LocalPoint localPosClu = rechitIt->localPosition();
0244       Global3DPoint globalPosClu = geomDetUnit->surface().toGlobal(localPosClu);
0245
0246       // restrict eta range
0247       float eta = globalPosClu.eta();
0248       if (fabs(eta) < mineta_ || fabs(eta) > maxeta_)
0249         continue;
0250
0251       // Get the cluster from the rechit
0252       const Phase2TrackerCluster1D* clustIt = &*rechitIt->cluster();
0253
0254       // Get all the simTracks that form the cluster
0255       std::vector<unsigned int> clusterSimTrackIds;
0256       for (unsigned int i(0); i < clustIt->size(); ++i) {
0257         unsigned int channel(Phase2TrackerDigi::pixelToChannel(clustIt->firstRow() + i, clustIt->column()));
0258         std::vector<unsigned int> simTrackIds(getSimTrackId(pixelSimLinks, detId, channel));
0259         for (unsigned int i = 0; i < simTrackIds.size(); ++i) {
0260           bool add = true;
0261           for (unsigned int j = 0; j < clusterSimTrackIds.size(); ++j) {
0262             // only save simtrackids that are not present yet
0263             if (simTrackIds.at(i) == clusterSimTrackIds.at(j))
0264               add = false;
0265           }
0266           if (add)
0267             clusterSimTrackIds.push_back(simTrackIds.at(i));
0268         }
0269       }
0270
0271       // find the closest simhit
0272       // this is needed because otherwise you get cases with simhits and clusters being swapped
0273       // when there are more than 1 cluster with common simtrackids
0274       const PSimHit* simhit = 0;  // bad naming to avoid changing code below. This is the closest simhit in x
0275       float minx = 10000;
0276       for (unsigned int simhitidx = 0; simhitidx < 2; ++simhitidx) {  // loop over both barrel and endcap hits
0277         for (edm::PSimHitContainer::const_iterator simhitIt(simHitsRaw[simhitidx]->begin());
0278              simhitIt != simHitsRaw[simhitidx]->end();
0279              ++simhitIt) {
0280           if (rawid == simhitIt->detUnitId()) {
0281             //std::cout << "=== " << rawid << " " << &*simhitIt << " " << simhitIt->trackId() << " " << simhitIt->localPosition().x() << " " << simhitIt->localPosition().y() << std::endl;
0282             auto it = std::lower_bound(clusterSimTrackIds.begin(), clusterSimTrackIds.end(), simhitIt->trackId());
0283             if (it != clusterSimTrackIds.end() && *it == simhitIt->trackId()) {
0284               if (!simhit || fabs(simhitIt->localPosition().x() - localPosClu.x()) < minx) {
0285                 minx = fabs(simhitIt->localPosition().x() - localPosClu.x());
0286                 simhit = &*simhitIt;
0287               }
0288             }
0289           }
0290         }
0291       }
0292       if (!simhit)
0293         continue;
0294
0295       // only look at simhits from highpT tracks
0296       std::map<unsigned int, SimTrack>::const_iterator simTrackIt(simTracks.find(simhit->trackId()));
0297       if (simTrackIt == simTracks.end())
0298         continue;
0299
0300       /*
0301              * Rechit related variables
0302              */
0303
0304       ++(nRecHits[det].at(layer));
0305       ++(nOtherSimHits[det].at(layer));
0306
0307       // cluster size
0308       unsigned int nch = rechitIt->cluster()->size();
0309       histogramLayer->second.clusterSize[det]->Fill(nch);
0310       if (nch > 4)
0311         nch = 4;  // collapse 4 or more strips to 4
0312
0313       // Fill the position histograms
0314       trackerLayout_->Fill(globalPosClu.z(), globalPosClu.perp());
0315       trackerLayoutXY_->Fill(globalPosClu.x(), globalPosClu.y());
0316       if (layer < 1000) {
0317         trackerLayoutXYBar_->Fill(globalPosClu.x(), globalPosClu.y());
0318       } else {
0319         trackerLayoutXYEC_->Fill(globalPosClu.x(), globalPosClu.y());
0320       }
0321
0322       histogramLayer->second.localPosXY[det][0]->Fill(localPosClu.x(), localPosClu.y());
0323       histogramLayer->second.localPosXY[det][nch]->Fill(localPosClu.x(), localPosClu.y());
0324       if (layer < 1000) {
0325         histogramLayer->second.globalPosXY[det][0]->Fill(globalPosClu.z(), globalPosClu.perp());
0326         histogramLayer->second.globalPosXY[det][nch]->Fill(globalPosClu.z(), globalPosClu.perp());
0327       } else {
0328         histogramLayer->second.globalPosXY[det][0]->Fill(globalPosClu.x(), globalPosClu.y());
0329         histogramLayer->second.globalPosXY[det][nch]->Fill(globalPosClu.x(), globalPosClu.y());
0330       }
0331
0332       // now get the position of the closest hit
0333       Local3DPoint localPosHit(simhit->localPosition());
0334
0335       // and fill bias and pull histograms
0336       histogramLayer->second.deltaX[det][0]->Fill(localPosClu.x() - localPosHit.x());
0337       histogramLayer->second.deltaX[det][nch]->Fill(localPosClu.x() - localPosHit.x());
0338       histogramLayer->second.deltaY[det][0]->Fill(localPosClu.y() - localPosHit.y());
0339       histogramLayer->second.deltaY[det][nch]->Fill(localPosClu.y() - localPosHit.y());
0340       if (rechitIt->localPositionError().xx() && rechitIt->localPositionError().yy()) {
0341         histogramLayer->second.pullX[det][0]->Fill((localPosClu.x() - localPosHit.x()) /
0342                                                    sqrt(rechitIt->localPositionError().xx()));
0343         histogramLayer->second.pullX[det][nch]->Fill((localPosClu.x() - localPosHit.x()) /
0344                                                      sqrt(rechitIt->localPositionError().xx()));
0345         histogramLayer->second.pullY[det][0]->Fill((localPosClu.y() - localPosHit.y()) /
0346                                                    sqrt(rechitIt->localPositionError().yy()));
0347         histogramLayer->second.pullY[det][nch]->Fill((localPosClu.y() - localPosHit.y()) /
0348                                                      sqrt(rechitIt->localPositionError().yy()));
0349       }
0350       if (makeEtaPlots_) {
0351         histogramLayer->second.deltaX_eta[det][0]->Fill(eta, localPosClu.x() - localPosHit.x());
0352         histogramLayer->second.deltaX_eta[det][nch]->Fill(eta, localPosClu.x() - localPosHit.x());
0353         histogramLayer->second.deltaY_eta[det][0]->Fill(eta, localPosClu.y() - localPosHit.y());
0354         histogramLayer->second.deltaY_eta[det][nch]->Fill(eta, localPosClu.y() - localPosHit.y());
0355         if (rechitIt->localPositionError().xx() && rechitIt->localPositionError().yy()) {
0356           histogramLayer->second.pullX_eta[det][0]->Fill(
0357               eta, (localPosClu.x() - localPosHit.x()) / sqrt(rechitIt->localPositionError().xx()));
0358           histogramLayer->second.pullX_eta[det][nch]->Fill(
0359               eta, (localPosClu.x() - localPosHit.x()) / sqrt(rechitIt->localPositionError().xx()));
0360           histogramLayer->second.pullY_eta[det][0]->Fill(
0361               eta, (localPosClu.y() - localPosHit.y()) / sqrt(rechitIt->localPositionError().yy()));
0362           histogramLayer->second.pullY_eta[det][nch]->Fill(
0363               eta, (localPosClu.y() - localPosHit.y()) / sqrt(rechitIt->localPositionError().yy()));
0364         }
0365       }
0366
0367       // fill histos for primary particles only
0368       unsigned int procT(simhit->processType());
0369       if (simTrackIt->second.vertIndex() == 0 and
0370           (procT == 2 || procT == 7 || procT == 9 || procT == 11 || procT == 13 || procT == 15)) {
0371         ++(nPrimarySimHits[det].at(layer));
0372         --(nOtherSimHits[det].at(layer));  // avoid double counting
0373         histogramLayer->second.deltaX_P[det][0]->Fill(localPosClu.x() - localPosHit.x());
0374         histogramLayer->second.deltaX_P[det][nch]->Fill(localPosClu.x() - localPosHit.x());
0375         histogramLayer->second.deltaY_P[det][0]->Fill(localPosClu.y() - localPosHit.y());
0376         histogramLayer->second.deltaY_P[det][nch]->Fill(localPosClu.y() - localPosHit.y());
0377         if (rechitIt->localPositionError().xx() && rechitIt->localPositionError().yy()) {
0378           histogramLayer->second.pullX_P[det][0]->Fill((localPosClu.x() - localPosHit.x()) /
0379                                                        sqrt(rechitIt->localPositionError().xx()));
0380           histogramLayer->second.pullX_P[det][nch]->Fill((localPosClu.x() - localPosHit.x()) /
0381                                                          sqrt(rechitIt->localPositionError().xx()));
0382           histogramLayer->second.pullY_P[det][0]->Fill((localPosClu.y() - localPosHit.y()) /
0383                                                        sqrt(rechitIt->localPositionError().yy()));
0384           histogramLayer->second.pullY_P[det][nch]->Fill((localPosClu.y() - localPosHit.y()) /
0385                                                          sqrt(rechitIt->localPositionError().yy()));
0386         }
0387         if (makeEtaPlots_) {
0388           histogramLayer->second.deltaX_eta_P[det][0]->Fill(eta, localPosClu.x() - localPosHit.x());
0389           histogramLayer->second.deltaX_eta_P[det][nch]->Fill(eta, localPosClu.x() - localPosHit.x());
0390           histogramLayer->second.deltaY_eta_P[det][0]->Fill(eta, localPosClu.y() - localPosHit.y());
0391           histogramLayer->second.deltaY_eta_P[det][nch]->Fill(eta, localPosClu.y() - localPosHit.y());
0392           if (rechitIt->localPositionError().xx() && rechitIt->localPositionError().yy()) {
0393             histogramLayer->second.pullX_eta_P[det][0]->Fill(
0394                 eta, (localPosClu.x() - localPosHit.x()) / sqrt(rechitIt->localPositionError().xx()));
0395             histogramLayer->second.pullX_eta_P[det][nch]->Fill(
0396                 eta, (localPosClu.x() - localPosHit.x()) / sqrt(rechitIt->localPositionError().xx()));
0397             histogramLayer->second.pullY_eta_P[det][0]->Fill(
0398                 eta, (localPosClu.y() - localPosHit.y()) / sqrt(rechitIt->localPositionError().yy()));
0399             histogramLayer->second.pullY_eta_P[det][nch]->Fill(
0400                 eta, (localPosClu.y() - localPosHit.y()) / sqrt(rechitIt->localPositionError().yy()));
0401           }
0402         }
0403       }
0404     }
0405   }
0406
0407   // fill the counter histos per layer
0408   for (unsigned int det = 1; det < 3; ++det) {
0409     for (auto it : nRecHits[det]) {
0410       auto histogramLayer(histograms_.find(it.first));
0411       if (histogramLayer == histograms_.end())
0412         std::cout << "*** SL *** No histogram for an existing counter! This should not happen!" << std::endl;
0413       histogramLayer->second.numberRecHits[det]->Fill(it.second);
0414     }
0415     for (auto it : nPrimarySimHits[det]) {
0416       auto histogramLayer(histograms_.find(it.first));
0417       if (histogramLayer == histograms_.end())
0418         std::cout << "*** SL *** No histogram for an existing counter! This should not happen!" << std::endl;
0419       histogramLayer->second.primarySimHits[det]->Fill(it.second);
0420     }
0421     for (auto it : nOtherSimHits[det]) {
0422       auto histogramLayer(histograms_.find(it.first));
0423       if (histogramLayer == histograms_.end())
0424         std::cout << "*** SL *** No histogram for an existing counter! This should not happen!" << std::endl;
0425       histogramLayer->second.otherSimHits[det]->Fill(it.second);
0426     }
0427   }
0428 }
0429
0430 // Create the histograms
0431 std::map<unsigned int, RecHitHistos>::iterator Phase2TrackerRecHitsValidation::createLayerHistograms(unsigned int ival) {
0432   std::ostringstream fname1, fname2;
0433
0434   edm::Service<TFileService> fs;
0435   fs->file().cd("/");
0436
0437   std::string tag;
0438   unsigned int id;
0439   if (ival < 1000) {
0440     id = ival;
0441     fname1 << "Barrel";
0442     fname2 << "Layer_" << id;
0443     tag = "_layer_";
0444   } else {
0445     int side = ival / 1000;
0446     id = ival - side * 1000;
0447     if (ival > 10) {
0448       id /= 10;
0449       //        fname1 << "EndCap_Side_" << side;
0450       fname1 << "EndCap";
0451       fname2 << "Ring_" << id;
0452       tag = "_ring_";
0453     } else {
0454       id = ival;
0455       //        fname1 << "EndCap_Side_" << side;
0456       fname1 << "EndCap";
0457       fname2 << "Disk_" << id;
0458       tag = "_disk_";
0459     }
0460   }
0461
0462   TFileDirectory td1 = fs->mkdir(fname1.str().c_str());
0463   TFileDirectory td = td1.mkdir(fname2.str().c_str());
0464
0465   RecHitHistos local_histos;
0466
0467   std::ostringstream histoName;
0468
0469   /*
0470      * Number of rechits
0471      */
0472
0473   histoName.str("");
0474   histoName << "Number_RecHits_Pixel" << tag.c_str() << id;
0475   local_histos.numberRecHits[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 50, 0., 0.);
0476   histoName.str("");
0477   histoName << "Number_RecHits_Strip" << tag.c_str() << id;
0478   local_histos.numberRecHits[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 50, 0., 0.);
0479
0480   /*
0481      * Cluster size
0482      */
0483
0484   histoName.str("");
0485   histoName << "Cluster_Size_Pixel" << tag.c_str() << id;
0486   local_histos.clusterSize[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 21, -0.5, 20.5);
0487   histoName.str("");
0488   histoName << "Cluster_Size_Strip" << tag.c_str() << id;
0489   local_histos.clusterSize[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 21, -0.5, 20.5);
0490
0491   /*
0492      * Local and Global positions
0493      */
0494
0495   for (int cls = 0; cls < 5; ++cls) {
0496     std::string clsstr = "";
0497     if (cls > 0)
0498       clsstr = "_ClS_" + std::to_string(cls);
0499
0500     histoName.str("");
0501     histoName << "Local_Position_XY_Pixel" << tag.c_str() << id << clsstr.c_str();
0502     local_histos.localPosXY[1][cls] =
0503         td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, 0., 0., 500, 0., 0.);
0504
0505     histoName.str("");
0506     histoName << "Local_Position_XY_Strip" << tag.c_str() << id << clsstr.c_str();
0507     local_histos.localPosXY[2][cls] =
0508         td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, 0., 0., 500, 0., 0.);
0509
0510     histoName.str("");
0511     histoName << "Global_Position_XY_Pixel" << tag.c_str() << id << clsstr.c_str();
0512     local_histos.globalPosXY[1][cls] =
0513         td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, -120., 120., 500, -120., 120.);
0514
0515     histoName.str("");
0516     histoName << "Global_Position_XY_Strip" << tag.c_str() << id << clsstr.c_str();
0517     local_histos.globalPosXY[2][cls] =
0518         td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 500, -120., 120., 500, -120., 120.);
0519
0520     /*
0521        * Delta positions with SimHits
0522        */
0523
0524     histoName.str("");
0525     histoName << "Delta_X_Pixel" << tag.c_str() << id << clsstr.c_str();
0526     local_histos.deltaX[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0527
0528     histoName.str("");
0529     histoName << "Delta_X_Strip" << tag.c_str() << id << clsstr.c_str();
0530     local_histos.deltaX[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0531
0532     histoName.str("");
0533     histoName << "Delta_Y_Pixel" << tag.c_str() << id << clsstr.c_str();
0534     local_histos.deltaY[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.2, 0.2);
0535
0536     histoName.str("");
0537     histoName << "Delta_Y_Strip" << tag.c_str() << id << clsstr.c_str();
0538     local_histos.deltaY[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3, 3);
0539
0540     if (makeEtaPlots_) {
0541       histoName.str("");
0542       histoName << "Delta_X_vs_Eta_Pixel" << tag.c_str() << id << clsstr.c_str();
0543       local_histos.deltaX_eta[1][cls] =
0544           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -0.02, 0.02);
0545
0546       histoName.str("");
0547       histoName << "Delta_X_vs_Eta_Strip" << tag.c_str() << id << clsstr.c_str();
0548       local_histos.deltaX_eta[2][cls] =
0549           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -0.02, 0.02);
0550
0551       histoName.str("");
0552       histoName << "Delta_Y_vs_Eta_Pixel" << tag.c_str() << id << clsstr.c_str();
0553       local_histos.deltaY_eta[1][cls] =
0554           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -0.2, 0.2);
0555
0556       histoName.str("");
0557       histoName << "Delta_Y_vs_Eta_Strip" << tag.c_str() << id << clsstr.c_str();
0558       local_histos.deltaY_eta[2][cls] =
0559           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3, 3);
0560     }
0561
0562     /*
0563        * Pulls
0564        */
0565
0566     histoName.str("");
0567     histoName << "Pull_X_Pixel" << tag.c_str() << id << clsstr.c_str();
0568     local_histos.pullX[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0569
0570     histoName.str("");
0571     histoName << "Pull_X_Strip" << tag.c_str() << id << clsstr.c_str();
0572     local_histos.pullX[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0573
0574     histoName.str("");
0575     histoName << "Pull_Y_Pixel" << tag.c_str() << id << clsstr.c_str();
0576     local_histos.pullY[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0577
0578     histoName.str("");
0579     histoName << "Pull_Y_Strip" << tag.c_str() << id << clsstr.c_str();
0580     local_histos.pullY[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0581
0582     if (makeEtaPlots_) {
0583       histoName.str("");
0584       histoName << "Pull_X_Eta_Pixel" << tag.c_str() << id << clsstr.c_str();
0585       local_histos.pullX_eta[1][cls] =
0586           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0587
0588       histoName.str("");
0589       histoName << "Pull_X_Eta_Strip" << tag.c_str() << id << clsstr.c_str();
0590       local_histos.pullX_eta[2][cls] =
0591           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0592
0593       histoName.str("");
0594       histoName << "Pull_Y_Eta_Pixel" << tag.c_str() << id << clsstr.c_str();
0595       local_histos.pullY_eta[1][cls] =
0596           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0597
0598       histoName.str("");
0599       histoName << "Pull_Y_Eta_Strip" << tag.c_str() << id << clsstr.c_str();
0600       local_histos.pullY_eta[2][cls] =
0601           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0602     }
0603
0604     /*
0605        * Delta position with simHits for primary tracks only
0606        */
0607
0608     histoName.str("");
0609     histoName << "Delta_X_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0610     local_histos.deltaX_P[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0611
0612     histoName.str("");
0613     histoName << "Delta_X_Strip_P" << tag.c_str() << id << clsstr.c_str();
0614     local_histos.deltaX_P[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.02, 0.02);
0615
0616     histoName.str("");
0617     histoName << "Delta_Y_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0618     local_histos.deltaY_P[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -0.2, 0.2);
0619
0620     histoName.str("");
0621     histoName << "Delta_Y_Strip_P" << tag.c_str() << id << clsstr.c_str();
0622     local_histos.deltaY_P[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0623
0624     if (makeEtaPlots_) {
0625       histoName.str("");
0626       histoName << "Delta_X_vs_Eta_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0627       local_histos.deltaX_eta_P[1][cls] =
0628           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -0.02, 0.02);
0629
0630       histoName.str("");
0631       histoName << "Delta_X_vs_Eta_Strip_P" << tag.c_str() << id << clsstr.c_str();
0632       local_histos.deltaX_eta_P[2][cls] =
0633           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -0.02, 0.02);
0634
0635       histoName.str("");
0636       histoName << "Delta_Y_vs_Eta_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0637       local_histos.deltaY_eta_P[1][cls] =
0638           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -0.2, 0.2);
0639
0640       histoName.str("");
0641       histoName << "Delta_Y_vs_Eta_Strip_P" << tag.c_str() << id << clsstr.c_str();
0642       local_histos.deltaY_eta_P[2][cls] =
0643           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3, 3);
0644     }
0645
0646     /*
0647        * Pulls for primary tracks only
0648        */
0649
0650     histoName.str("");
0651     histoName << "Pull_X_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0652     local_histos.pullX_P[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0653
0654     histoName.str("");
0655     histoName << "Pull_X_Strip_P" << tag.c_str() << id << clsstr.c_str();
0656     local_histos.pullX_P[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0657
0658     histoName.str("");
0659     histoName << "Pull_Y_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0660     local_histos.pullY_P[1][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0661
0662     histoName.str("");
0663     histoName << "Pull_Y_Strip_P" << tag.c_str() << id << clsstr.c_str();
0664     local_histos.pullY_P[2][cls] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, -3., 3.);
0665
0666     if (makeEtaPlots_) {
0667       histoName.str("");
0668       histoName << "Pull_X_Eta_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0669       local_histos.pullX_eta_P[1][cls] =
0670           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0671
0672       histoName.str("");
0673       histoName << "Pull_X_Eta_Strip_P" << tag.c_str() << id << clsstr.c_str();
0674       local_histos.pullX_eta_P[2][cls] =
0675           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0676
0677       histoName.str("");
0678       histoName << "Pull_Y_Eta_Pixel_P" << tag.c_str() << id << clsstr.c_str();
0679       local_histos.pullY_eta_P[1][cls] =
0680           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0681
0682       histoName.str("");
0683       histoName << "Pull_Y_Eta_Strip_P" << tag.c_str() << id << clsstr.c_str();
0684       local_histos.pullY_eta_P[2][cls] =
0685           td.make<TH2F>(histoName.str().c_str(), histoName.str().c_str(), 50, -2.5, 2.5, 100, -3., 3.);
0686     }
0687   }
0688
0689   /*
0690      * Information on the Digis per cluster
0691      */
0692
0693   histoName.str("");
0694   histoName << "Primary_Digis_Pixel" << tag.c_str() << id;
0695   local_histos.primarySimHits[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0696   histoName.str("");
0697   histoName << "Primary_Digis_Strip" << tag.c_str() << id;
0698   local_histos.primarySimHits[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0699
0700   histoName.str("");
0701   histoName << "Other_Digis_Pixel" << tag.c_str() << id;
0702   local_histos.otherSimHits[1] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0703   histoName.str("");
0704   histoName << "Other_Digis_Strip" << tag.c_str() << id;
0705   local_histos.otherSimHits[2] = td.make<TH1D>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
0706
0707   /*
0708      * End
0709      */
0710
0711   std::pair<std::map<unsigned int, RecHitHistos>::iterator, bool> insertedIt(
0712       histograms_.insert(std::make_pair(ival, local_histos)));
0713   fs->file().cd("/");
0714
0715   return insertedIt.first;
0716 }
0717
0718 std::vector<unsigned int> Phase2TrackerRecHitsValidation::getSimTrackId(
0719     const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& pixelSimLinks, const DetId& detId, unsigned int channel) {
0720   std::vector<unsigned int> retvec;
0721   edm::DetSetVector<PixelDigiSimLink>::const_iterator DSViter(pixelSimLinks->find(detId));
0722   if (DSViter == pixelSimLinks->end())
0723     return retvec;
0724   for (edm::DetSet<PixelDigiSimLink>::const_iterator it = DSViter->data.begin(); it != DSViter->data.end(); ++it) {
0725     if (channel == it->channel()) {
0726       retvec.push_back(it->SimTrackId());
0727     }
0728   }
0729   return retvec;
0730 }
0731
0732 DEFINE_FWK_MODULE(Phase2TrackerRecHitsValidation);