DTRecHits/plugins/DTSegment2DSLPhiQuality.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author S. Bolognesi and G. Cerminara - INFN Torino
0005  */
0006
0007 #include <iostream>
0008 #include <map>
0009
0010 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0011 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0012 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "Geometry/DTGeometry/interface/DTChamber.h"
0017 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0018 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0019 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
0020
0021 #include "DTSegment2DSLPhiQuality.h"
0022 #include "Histograms.h"
0023
0024 using namespace std;
0025 using namespace edm;
0026
0027 namespace dtsegment2dsl {
0028   struct Histograms {
0029     std::unique_ptr<HRes2DHit> h2DHitSuperPhi;
0030     std::unique_ptr<HEff2DHit> h2DHitEff_SuperPhi;
0031   };
0032 }  // namespace dtsegment2dsl
0033
0034 using namespace dtsegment2dsl;
0035
0036 // Constructor
0037 DTSegment2DSLPhiQuality::DTSegment2DSLPhiQuality(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0038   // Get the debug parameter for verbose output
0039   debug_ = pset.getUntrackedParameter<bool>("debug");
0040   DTHitQualityUtils::debug = debug_;
0041
0042   // the name of the simhit collection
0043   simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
0044   simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
0045   // the name of the 2D rec hit collection
0046   segment4DLabel_ = pset.getUntrackedParameter<InputTag>("segment4DLabel");
0047   segment4DToken_ = consumes<DTRecSegment4DCollection>(pset.getUntrackedParameter<InputTag>("segment4DLabel"));
0048
0049   // sigma resolution on position
0050   sigmaResPos_ = pset.getParameter<double>("sigmaResPos");
0051   // sigma resolution on angle
0052   sigmaResAngle_ = pset.getParameter<double>("sigmaResAngle");
0053   doall_ = pset.getUntrackedParameter<bool>("doall", false);
0054   local_ = pset.getUntrackedParameter<bool>("local", false);
0055 }
0056
0057 void DTSegment2DSLPhiQuality::bookHistograms(DQMStore::IBooker &booker,
0058                                              edm::Run const &run,
0059                                              edm::EventSetup const &setup,
0060                                              Histograms &histograms) const {
0061   // Book the histos
0062   histograms.h2DHitSuperPhi = std::make_unique<HRes2DHit>("SuperPhi", booker, doall_, local_);
0063   if (doall_) {
0064     histograms.h2DHitEff_SuperPhi = std::make_unique<HEff2DHit>("SuperPhi", booker);
0065   }
0066 }
0067
0068 // The real analysis
0069 void DTSegment2DSLPhiQuality::dqmAnalyze(edm::Event const &event,
0070                                          edm::EventSetup const &setup,
0071                                          Histograms const &histograms) const {
0072   // Get the DT Geometry
0073   const DTGeometry &dtGeom = setup.getData(muonGeomToken_);
0074
0075   // Get the SimHit collection from the event
0076   edm::Handle<PSimHitContainer> simHits;
0077   event.getByToken(simHitToken_, simHits);  // FIXME: second string to be removed
0078
0079   // Map simHits by chamber
0080   map<DTChamberId, PSimHitContainer> simHitsPerCh;
0081   for (const auto &simHit : *simHits) {
0082     // Create the id of the chamber (the simHits in the DT known their wireId)
0083     DTChamberId chamberId = (((DTWireId(simHit.detUnitId())).layerId()).superlayerId()).chamberId();
0084     // Fill the map
0085     simHitsPerCh[chamberId].push_back(simHit);
0086   }
0087
0088   // Get the 4D rechits from the event
0089   Handle<DTRecSegment4DCollection> segment4Ds;
0090   event.getByToken(segment4DToken_, segment4Ds);
0091
0092   if (!segment4Ds.isValid()) {
0093     if (debug_) {
0094       cout << "[DTSegment2DSLPhiQuality]**Warning: no 4D Segments with label: " << segment4DLabel_
0095            << " in this event, skipping!" << endl;
0096     }
0097     return;
0098   }
0099
0100   // Loop over all chambers containing a segment
0101   DTRecSegment4DCollection::id_iterator chamberId;
0102   for (chamberId = segment4Ds->id_begin(); chamberId != segment4Ds->id_end(); ++chamberId) {
0103     //------------------------- simHits ---------------------------//
0104     // Get simHits of each chamber
0105     PSimHitContainer simHits = simHitsPerCh[(*chamberId)];
0106
0107     // Map simhits per wire
0108     map<DTWireId, PSimHitContainer> simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
0109     map<DTWireId, const PSimHit *> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
0110     int nMuSimHit = muSimHitPerWire.size();
0111     if (nMuSimHit == 0 || nMuSimHit == 1) {
0112       if (debug_ && nMuSimHit == 1) {
0113         cout << "[DTSegment2DSLPhiQuality] Only " << nMuSimHit << " mu SimHit in this chamber, skipping!" << endl;
0114       }
0115       continue;  // If no or only one mu SimHit is found skip this chamber
0116     }
0117     if (debug_) {
0118       cout << "=== Chamber " << (*chamberId) << " has " << nMuSimHit << " SimHits" << endl;
0119     }
0120
0121     // Find outer and inner mu SimHit to build a segment
0122     pair<const PSimHit *, const PSimHit *> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
0123
0124     // Find direction and position of the sim Segment in Chamber RF
0125     pair<LocalVector, LocalPoint> dirAndPosSimSegm =
0126         DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit, (*chamberId), dtGeom);
0127
0128     LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
0129     LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
0130     const DTChamber *chamber = dtGeom.chamber(*chamberId);
0131     GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos);
0132
0133     // Atan(x/z) angle and x position in SL RF
0134     float angleSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
0135     float posSimSeg = simSegmLocalPos.x();
0136     // Position (in eta, phi coordinates) in lobal RF
0137     float etaSimSeg = simSegmGlobalPos.eta();
0138     float phiSimSeg = simSegmGlobalPos.phi();
0139
0140     if (debug_) {
0141       cout << "  Simulated segment:  local direction " << simSegmLocalDir << endl
0142            << "                      local position  " << simSegmLocalPos << endl
0143            << "                      angle           " << angleSimSeg << endl;
0144     }
0145
0146     //---------------------------- recHits --------------------------//
0147     // Get the range of rechit for the corresponding chamberId
0148     bool recHitFound = false;
0149     DTRecSegment4DCollection::range range = segment4Ds->get(*chamberId);
0150     int nsegm = distance(range.first, range.second);
0151     if (debug_) {
0152       cout << "   Chamber: " << *chamberId << " has " << nsegm << " 4D segments" << endl;
0153     }
0154
0155     if (nsegm != 0) {
0156       // Find the best RecHit: look for the 4D RecHit with the phi angle closest
0157       // to that of segment made of SimHits.
0158       // RecHits must have delta alpha and delta position within 5 sigma of
0159       // the residual distribution (we are looking for residuals of segments
0160       // usefull to the track fit) for efficency purpose
0161       const DTRecSegment2D *bestRecHit = nullptr;
0162       bool bestRecHitFound = false;
0163       double deltaAlpha = 99999;
0164
0165       // Loop over the recHits of this chamberId
0166       for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0167         // Check the dimension
0168         if ((*segment4D).dimension() != 4) {
0169           if (debug_) {
0170             cout << "[DTSegment2DSLPhiQuality]***Error: This is not 4D "
0171                     "segment!!!"
0172                  << endl;
0173           }
0174           continue;
0175         }
0176
0177         // Get 2D superPhi segments from 4D segments
0178         const DTChamberRecSegment2D *phiSegment2D = (*segment4D).phiSegment();
0179         if ((*phiSegment2D).dimension() != 2) {
0180           if (debug_) {
0181             cout << "[DTSegment2DQuality]***Error: This is not 2D segment!!!" << endl;
0182           }
0183           abort();
0184         }
0185
0186         // Segment Local Direction and position (in Chamber RF)
0187         LocalVector recSegDirection = (*phiSegment2D).localDirection();
0188
0189         float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
0190         if (debug_) {
0191           cout << "  RecSegment direction: " << recSegDirection << endl
0192                << "             position : " << (*phiSegment2D).localPosition() << endl
0193                << "             alpha    : " << recSegAlpha << endl;
0194         }
0195
0196         if (fabs(recSegAlpha - angleSimSeg) < deltaAlpha) {
0197           deltaAlpha = fabs(recSegAlpha - angleSimSeg);
0198           bestRecHit = &(*phiSegment2D);
0199           bestRecHitFound = true;
0200         }
0201       }  // End of Loop over all 4D RecHits of this chambers
0202
0203       if (bestRecHitFound) {
0204         // Best rechit direction and position in Chamber RF
0205         LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
0206         LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
0207
0208         LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
0209         LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
0210
0211         float angleBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
0212         if (fabs(angleBestRHit - angleSimSeg) < 5 * sigmaResAngle_ &&
0213             fabs(bestRecHitLocalPos.x() - posSimSeg) < 5 * sigmaResPos_) {
0214           recHitFound = true;
0215         }
0216
0217         // Fill Residual histos
0218         histograms.h2DHitSuperPhi->fill(angleSimSeg,
0219                                         angleBestRHit,
0220                                         posSimSeg,
0221                                         bestRecHitLocalPos.x(),
0222                                         etaSimSeg,
0223                                         phiSimSeg,
0224                                         sqrt(bestRecHitLocalPosErr.xx()),
0225                                         sqrt(bestRecHitLocalDirErr.xx()));
0226       }
0227     }  // end of if (nsegm!= 0)
0228
0229     // Fill Efficiency plot
0230     if (doall_) {
0231       histograms.h2DHitEff_SuperPhi->fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
0232     }
0233   }  // End of loop over chambers
0234 }
0235
0236 // declare this as a framework plugin
0237 #include "FWCore/Framework/interface/MakerMacros.h"
0238 DEFINE_FWK_MODULE(DTSegment2DSLPhiQuality);