CSCRecHits/src/CSCSegmentValidation.cc

0001 #include "DQMServices/Core/interface/DQMStore.h"
0002 #include "Validation/CSCRecHits/interface/CSCSegmentValidation.h"
0003 #include <algorithm>
0004
0005 CSCSegmentValidation::CSCSegmentValidation(const edm::ParameterSet &ps, edm::ConsumesCollector &&iC)
0006     : CSCBaseValidation(ps), theLayerHitsPerChamber(), theChamberSegmentMap(), theShowerThreshold(10) {
0007   const auto &pset = ps.getParameterSet("cscSegment");
0008   inputTag_ = pset.getParameter<edm::InputTag>("inputTag");
0009   segments_Token_ = iC.consumes<CSCSegmentCollection>(inputTag_);
0010 }
0011
0012 CSCSegmentValidation::~CSCSegmentValidation() {}
0013
0014 void CSCSegmentValidation::bookHistograms(DQMStore::IBooker &iBooker) {
0015   theNPerEventPlot = iBooker.book1D("CSCSegmentsPerEvent", "Number of CSC segments per event", 100, 0, 50);
0016   theNRecHitsPlot = iBooker.book1D("CSCRecHitsPerSegment", "Number of CSC rec hits per segment", 8, 0, 7);
0017   theNPerChamberTypePlot =
0018       iBooker.book1D("CSCSegmentsPerChamberType", "Number of CSC segments per chamber type", 11, 0, 10);
0019   theTypePlot4HitsNoShower = iBooker.book1D("CSCSegments4HitsNoShower", "", 100, 0, 10);
0020   theTypePlot4HitsNoShowerSeg = iBooker.book1D("CSCSegments4HitsNoShowerSeg", "", 100, 0, 10);
0021   theTypePlot4HitsShower = iBooker.book1D("CSCSegments4HitsShower", "", 100, 0, 10);
0022   theTypePlot4HitsShowerSeg = iBooker.book1D("CSCSegments4HitsShowerSeg", "", 100, 0, 10);
0023   theTypePlot5HitsNoShower = iBooker.book1D("CSCSegments5HitsNoShower", "", 100, 0, 10);
0024   theTypePlot5HitsNoShowerSeg = iBooker.book1D("CSCSegments5HitsNoShowerSeg", "", 100, 0, 10);
0025   theTypePlot5HitsShower = iBooker.book1D("CSCSegments5HitsShower", "", 100, 0, 10);
0026   theTypePlot5HitsShowerSeg = iBooker.book1D("CSCSegments5HitsShowerSeg", "", 100, 0, 10);
0027   theTypePlot6HitsNoShower = iBooker.book1D("CSCSegments6HitsNoShower", "", 100, 0, 10);
0028   theTypePlot6HitsNoShowerSeg = iBooker.book1D("CSCSegments6HitsNoShowerSeg", "", 100, 0, 10);
0029   theTypePlot6HitsShower = iBooker.book1D("CSCSegments6HitsShower", "", 100, 0, 10);
0030   theTypePlot6HitsShowerSeg = iBooker.book1D("CSCSegments6HitsShowerSeg", "", 100, 0, 10);
0031
0032   for (int i = 1; i <= 10; ++i) {
0033     const std::string cn(CSCDetId::chamberName(i));
0034
0035     const std::string t1("CSCSegmentRdPhiResolution_" + cn);
0036     const std::string t2("CSCSegmentRdPhiPull_" + cn);
0037     const std::string t3("CSCSegmentThetaResolution_" + cn);
0038     const std::string t5("CSCSegmentdXdZResolution_" + cn);
0039     const std::string t6("CSCSegmentdXdZPull_" + cn);
0040     const std::string t7("CSCSegmentdYdZResolution_" + cn);
0041     const std::string t8("CSCSegmentdYdZPull_" + cn);
0042
0043     theRdPhiResolutionPlots[i - 1] = iBooker.book1D(t1, t1, 100, -0.4, 0.4);
0044     theRdPhiPullPlots[i - 1] = iBooker.book1D(t2, t2, 100, -5, 5);
0045     theThetaResolutionPlots[i - 1] = iBooker.book1D(t3, t3, 100, -1, 1);
0046     thedXdZResolutionPlots[i - 1] = iBooker.book1D(t5, t5, 100, -1, 1);
0047     thedXdZPullPlots[i - 1] = iBooker.book1D(t6, t6, 100, -5, 5);
0048     thedYdZResolutionPlots[i - 1] = iBooker.book1D(t7, t7, 100, -1, 1);
0049     thedYdZPullPlots[i - 1] = iBooker.book1D(t8, t8, 100, -5, 5);
0050   }
0051 }
0052
0053 void CSCSegmentValidation::analyze(const edm::Event &e, const edm::EventSetup &eventSetup) {
0054   // get the collection of CSCRecHsegmentItrD
0055   edm::Handle<CSCSegmentCollection> hRecHits;
0056   e.getByToken(segments_Token_, hRecHits);
0057   const CSCSegmentCollection *cscRecHits = hRecHits.product();
0058
0059   theChamberSegmentMap.clear();
0060   unsigned nPerEvent = 0;
0061   for (auto segmentItr = cscRecHits->begin(); segmentItr != cscRecHits->end(); segmentItr++) {
0062     ++nPerEvent;
0063     int detId = segmentItr->geographicalId().rawId();
0064     int chamberType = segmentItr->cscDetId().iChamberType();
0065
0066     theNRecHitsPlot->Fill(segmentItr->nRecHits());
0067     theNPerChamberTypePlot->Fill(chamberType);
0068     theChamberSegmentMap[detId].push_back(*segmentItr);
0069
0070     // do the resolution plots
0071     const PSimHit *hit = keyHit(detId);
0072     if (hit != nullptr) {
0073       const CSCLayer *layer = findLayer(hit->detUnitId());
0074       plotResolution(*hit, *segmentItr, layer, chamberType);
0075     }
0076   }
0077
0078   theNPerEventPlot->Fill(nPerEvent);
0079
0080   fillLayerHitsPerChamber();
0081   fillEfficiencyPlots();
0082 }
0083
0084 void CSCSegmentValidation::fillEfficiencyPlots() {
0085   // now plot efficiency by looping over all chambers with hits
0086   for (auto mapItr = theLayerHitsPerChamber.begin(), mapEnd = theLayerHitsPerChamber.end(); mapItr != mapEnd;
0087        ++mapItr) {
0088     int chamberId = mapItr->first;
0089     int nHitsInChamber = mapItr->second.size();
0090     bool isShower = (nHitsInChamber > theShowerThreshold);
0091     bool hasSeg = hasSegment(chamberId);
0092     int chamberType = CSCDetId(chamberId).iChamberType();
0093     // find how many layers were hit in this chamber
0094     std::vector<int> v = mapItr->second;
0095     std::sort(v.begin(), v.end());
0096     // maybe can just count
0097     v.erase(std::unique(v.begin(), v.end()), v.end());
0098     int nLayersHit = v.size();
0099
0100     if (nLayersHit == 4) {
0101       if (isShower)
0102         theTypePlot4HitsShower->Fill(chamberType);
0103       else
0104         theTypePlot4HitsNoShower->Fill(chamberType);
0105
0106       if (hasSeg) {
0107         if (isShower)
0108           theTypePlot4HitsShowerSeg->Fill(chamberType);
0109         else
0110           theTypePlot4HitsNoShowerSeg->Fill(chamberType);
0111       }
0112     }
0113
0114     if (nLayersHit == 5) {
0115       if (isShower)
0116         theTypePlot5HitsShower->Fill(chamberType);
0117       else
0118         theTypePlot5HitsNoShower->Fill(chamberType);
0119
0120       if (hasSeg) {
0121         if (isShower)
0122           theTypePlot5HitsShowerSeg->Fill(chamberType);
0123         else
0124           theTypePlot5HitsNoShowerSeg->Fill(chamberType);
0125       }
0126     }
0127
0128     if (nLayersHit == 6) {
0129       if (isShower)
0130         theTypePlot6HitsShower->Fill(chamberType);
0131       else
0132         theTypePlot6HitsNoShower->Fill(chamberType);
0133
0134       if (hasSeg) {
0135         if (isShower)
0136           theTypePlot6HitsShowerSeg->Fill(chamberType);
0137         else
0138           theTypePlot6HitsNoShowerSeg->Fill(chamberType);
0139       }
0140     }
0141   }
0142 }
0143
0144 bool CSCSegmentValidation::hasSegment(int chamberId) const {
0145   return (theChamberSegmentMap.find(chamberId) != theChamberSegmentMap.end());
0146 }
0147
0148 void CSCSegmentValidation::plotResolution(const PSimHit &simHit,
0149                                           const CSCSegment &segment,
0150                                           const CSCLayer *layer,
0151                                           int chamberType) {
0152   GlobalPoint simHitPos = layer->toGlobal(simHit.localPosition());
0153   GlobalPoint segmentPos = layer->toGlobal(segment.localPosition());
0154   LocalVector simHitDir = simHit.localDirection();
0155   LocalVector segmentDir = segment.localDirection();
0156
0157   double dphi = segmentPos.phi() - simHitPos.phi();
0158   double rdphi = segmentPos.perp() * dphi;
0159   double dtheta = segmentPos.theta() - simHitPos.theta();
0160
0161   double sigmax = sqrt(segment.localPositionError().xx());
0162
0163   double ddxdz = segmentDir.x() / segmentDir.z() - simHitDir.x() / simHitDir.z();
0164   double ddydz = segmentDir.y() / segmentDir.z() - simHitDir.y() / simHitDir.z();
0165   double sigmadxdz = sqrt(segment.localDirectionError().xx());
0166   double sigmadydz = sqrt(segment.localDirectionError().yy());
0167
0168   theRdPhiResolutionPlots[chamberType - 1]->Fill(rdphi);
0169   theRdPhiPullPlots[chamberType - 1]->Fill(rdphi / sigmax);
0170   theThetaResolutionPlots[chamberType - 1]->Fill(dtheta);
0171
0172   thedXdZResolutionPlots[chamberType - 1]->Fill(ddxdz);
0173   thedXdZPullPlots[chamberType - 1]->Fill(ddxdz / sigmadxdz);
0174   thedYdZResolutionPlots[chamberType - 1]->Fill(ddydz);
0175   thedYdZPullPlots[chamberType - 1]->Fill(ddydz / sigmadydz);
0176 }
0177
0178 void CSCSegmentValidation::fillLayerHitsPerChamber() {
0179   theLayerHitsPerChamber.clear();
0180   std::vector<int> layersHit = theSimHitMap->detsWithHits();
0181   for (auto layerItr = layersHit.begin(), layersHitEnd = layersHit.end(); layerItr != layersHitEnd; ++layerItr) {
0182     CSCDetId layerId(*layerItr);
0183     CSCDetId chamberId = layerId.chamberId();
0184     int nhits = theSimHitMap->hits(*layerItr).size();
0185     // multiple entries, so we can see showers
0186     for (int i = 0; i < nhits; ++i) {
0187       theLayerHitsPerChamber[chamberId.rawId()].push_back(*layerItr);
0188     }
0189   }
0190 }
0191
0192 const PSimHit *CSCSegmentValidation::keyHit(int chamberId) const {
0193   auto SimHitPabsLessThan = [](const PSimHit &p1, const PSimHit &p2) -> bool { return p1.pabs() < p2.pabs(); };
0194
0195   const PSimHit *result = nullptr;
0196   int layerId = chamberId + 3;
0197   const auto &layerHits = theSimHitMap->hits(layerId);
0198
0199   if (!layerHits.empty()) {
0200     // pick the hit with maximum energy
0201     auto hitItr = std::max_element(layerHits.begin(), layerHits.end(), SimHitPabsLessThan);
0202     result = &(*hitItr);
0203   }
0204   return result;
0205 }