DTRecHits/plugins/DTSegment2DQuality.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/DTRecSegment2DCollection.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/DTGeometry.h"
0017 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0018 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0019 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
0020
0021 #include "DTSegment2DQuality.h"
0022 #include "Histograms.h"
0023
0024 using namespace std;
0025 using namespace edm;
0026
0027 namespace dtsegment2d {
0028   struct Histograms {
0029     std::unique_ptr<HRes2DHit> h2DHitRPhi;
0030     std::unique_ptr<HRes2DHit> h2DHitRZ;
0031     std::unique_ptr<HRes2DHit> h2DHitRZ_W0;
0032     std::unique_ptr<HRes2DHit> h2DHitRZ_W1;
0033     std::unique_ptr<HRes2DHit> h2DHitRZ_W2;
0034
0035     std::unique_ptr<HEff2DHit> h2DHitEff_RPhi;
0036     std::unique_ptr<HEff2DHit> h2DHitEff_RZ;
0037     std::unique_ptr<HEff2DHit> h2DHitEff_RZ_W0;
0038     std::unique_ptr<HEff2DHit> h2DHitEff_RZ_W1;
0039     std::unique_ptr<HEff2DHit> h2DHitEff_RZ_W2;
0040   };
0041 }  // namespace dtsegment2d
0042
0043 using namespace dtsegment2d;
0044
0045 // Constructor
0046 DTSegment2DQuality::DTSegment2DQuality(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0047   // get the debug parameter for verbose output
0048   debug_ = pset.getUntrackedParameter<bool>("debug");
0049   DTHitQualityUtils::debug = debug_;
0050   // the name of the simhit collection
0051   simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
0052   simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
0053   // the name of the 2D rec hit collection
0054   segment2DLabel_ = pset.getUntrackedParameter<InputTag>("segment2DLabel");
0055   segment2DToken_ = consumes<DTRecSegment2DCollection>(pset.getUntrackedParameter<InputTag>("segment2DLabel"));
0056
0057   // sigma resolution on position
0058   sigmaResPos_ = pset.getParameter<double>("sigmaResPos");
0059   // sigma resolution on angle
0060   sigmaResAngle_ = pset.getParameter<double>("sigmaResAngle");
0061
0062   if (debug_) {
0063     cout << "[DTSegment2DQuality] Constructor called " << endl;
0064   }
0065 }
0066
0067 void DTSegment2DQuality::bookHistograms(DQMStore::IBooker &booker,
0068                                         edm::Run const &run,
0069                                         edm::EventSetup const &setup,
0070                                         Histograms &histograms) const {
0071   histograms.h2DHitRPhi = std::make_unique<HRes2DHit>("RPhi", booker, true, true);
0072   histograms.h2DHitRZ = std::make_unique<HRes2DHit>("RZ", booker, true, true);
0073   histograms.h2DHitRZ_W0 = std::make_unique<HRes2DHit>("RZ_W0", booker, true, true);
0074   histograms.h2DHitRZ_W1 = std::make_unique<HRes2DHit>("RZ_W1", booker, true, true);
0075   histograms.h2DHitRZ_W2 = std::make_unique<HRes2DHit>("RZ_W2", booker, true, true);
0076
0077   histograms.h2DHitEff_RPhi = std::make_unique<HEff2DHit>("RPhi", booker);
0078   histograms.h2DHitEff_RZ = std::make_unique<HEff2DHit>("RZ", booker);
0079   histograms.h2DHitEff_RZ_W0 = std::make_unique<HEff2DHit>("RZ_W0", booker);
0080   histograms.h2DHitEff_RZ_W1 = std::make_unique<HEff2DHit>("RZ_W1", booker);
0081   histograms.h2DHitEff_RZ_W2 = std::make_unique<HEff2DHit>("RZ_W2", booker);
0082   if (debug_) {
0083     cout << "[DTSegment2DQuality] hitsos created " << endl;
0084   }
0085 }
0086
0087 // The real analysis
0088 void DTSegment2DQuality::dqmAnalyze(edm::Event const &event,
0089                                     edm::EventSetup const &setup,
0090                                     Histograms const &histograms) const {
0091   // Get the DT Geometry
0092   const DTGeometry &dtGeom = setup.getData(muonGeomToken_);
0093
0094   // Get the SimHit collection from the event
0095   edm::Handle<PSimHitContainer> simHits;
0096   event.getByToken(simHitToken_, simHits);  // FIXME: second string to be removed
0097
0098   // Map simHits by sl
0099   map<DTSuperLayerId, PSimHitContainer> simHitsPerSl;
0100   for (const auto &simHit : *simHits) {
0101     // Create the id of the sl (the simHits in the DT known their wireId)
0102     DTSuperLayerId slId = ((DTWireId(simHit.detUnitId())).layerId()).superlayerId();
0103     // Fill the map
0104     simHitsPerSl[slId].push_back(simHit);
0105   }
0106
0107   // Get the 2D rechits from the event
0108   Handle<DTRecSegment2DCollection> segment2Ds;
0109   event.getByToken(segment2DToken_, segment2Ds);
0110
0111   if (not segment2Ds.isValid()) {
0112     if (debug_) {
0113       cout << "[DTSegment2DQuality]**Warning: no 2DSegments with label: " << segment2DLabel_
0114            << " in this event, skipping !" << endl;
0115     }
0116     return;
0117   }
0118
0119   // Loop over all superlayers containing a segment
0120   DTRecSegment2DCollection::id_iterator slId;
0121   for (slId = segment2Ds->id_begin(); slId != segment2Ds->id_end(); ++slId) {
0122     //------------------------- simHits ---------------------------//
0123     // Get simHits of each superlayer
0124     PSimHitContainer simHits = simHitsPerSl[(*slId)];
0125
0126     // Map simhits per wire
0127     map<DTWireId, PSimHitContainer> simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
0128     map<DTWireId, const PSimHit *> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
0129     int nMuSimHit = muSimHitPerWire.size();
0130     if (nMuSimHit == 0 or nMuSimHit == 1) {
0131       if (debug_ and nMuSimHit == 1) {
0132         cout << "[DTSegment2DQuality] Only " << nMuSimHit << " mu SimHit in this SL, skipping !" << endl;
0133       }
0134       continue;  // If no or only one mu SimHit is found skip this SL
0135     }
0136     if (debug_) {
0137       cout << "=== SL " << (*slId) << " has " << nMuSimHit << " SimHits" << endl;
0138     }
0139
0140     // Find outer and inner mu SimHit to build a segment
0141     pair<const PSimHit *, const PSimHit *> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
0142     // Check that outermost and innermost SimHit are not the same
0143     if (inAndOutSimHit.first == inAndOutSimHit.second) {
0144       cout << "[DTHitQualityUtils]***Warning: outermost and innermost SimHit "
0145               "are the same !"
0146            << endl;
0147       continue;
0148     }
0149
0150     // Find direction and position of the sim Segment in SL RF
0151     pair<LocalVector, LocalPoint> dirAndPosSimSegm =
0152         DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit, (*slId), dtGeom);
0153
0154     LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
0155     LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
0156     if (debug_) {
0157       cout << "  Simulated segment:  local direction " << simSegmLocalDir << endl
0158            << "                      local position  " << simSegmLocalPos << endl;
0159     }
0160     const DTSuperLayer *superLayer = dtGeom.superLayer(*slId);
0161     GlobalPoint simSegmGlobalPos = superLayer->toGlobal(simSegmLocalPos);
0162
0163     // Atan(x/z) angle and x position in SL RF
0164     float angleSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
0165     float posSimSeg = simSegmLocalPos.x();
0166     // Position (in eta, phi coordinates) in the global RF
0167     float etaSimSeg = simSegmGlobalPos.eta();
0168     float phiSimSeg = simSegmGlobalPos.phi();
0169
0170     //---------------------------- recHits --------------------------//
0171     // Get the range of rechit for the corresponding slId
0172     bool recHitFound = false;
0173     DTRecSegment2DCollection::range range = segment2Ds->get(*slId);
0174     int nsegm = distance(range.first, range.second);
0175     if (debug_) {
0176       cout << "   Sl: " << *slId << " has " << nsegm << " 2D segments" << endl;
0177     }
0178
0179     if (nsegm != 0) {
0180       // Find the best RecHit: look for the 2D RecHit with the angle closest
0181       // to that of segment made of SimHits.
0182       // RecHits must have delta alpha and delta position within 5 sigma of
0183       // the residual distribution (we are looking for residuals of segments
0184       // usefull to the track fit) for efficency purpose
0185       const DTRecSegment2D *bestRecHit = nullptr;
0186       bool bestRecHitFound = false;
0187       double deltaAlpha = 99999;
0188
0189       // Loop over the recHits of this slId
0190       for (DTRecSegment2DCollection::const_iterator segment2D = range.first; segment2D != range.second; ++segment2D) {
0191         // Check the dimension
0192         if ((*segment2D).dimension() != 2) {
0193           if (debug_) {
0194             cout << "[DTSegment2DQuality]***Error: This is not 2D segment !!!" << endl;
0195           }
0196           abort();
0197         }
0198         // Segment Local Direction and position (in SL RF)
0199         LocalVector recSegDirection = (*segment2D).localDirection();
0200         LocalPoint recSegPosition = (*segment2D).localPosition();
0201
0202         float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
0203         if (debug_) {
0204           cout << "  RecSegment direction: " << recSegDirection << endl
0205                << "             position : " << recSegPosition << endl
0206                << "             alpha    : " << recSegAlpha << endl;
0207         }
0208
0209         if (fabs(recSegAlpha - angleSimSeg) < deltaAlpha) {
0210           deltaAlpha = fabs(recSegAlpha - angleSimSeg);
0211           bestRecHit = &(*segment2D);
0212           bestRecHitFound = true;
0213         }
0214       }  // End of Loop over all 2D RecHits
0215
0216       if (bestRecHitFound) {
0217         // Best rechit direction and position in SL RF
0218         LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
0219         LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
0220
0221         LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
0222         LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
0223
0224         float angleBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
0225
0226         if (fabs(angleBestRHit - angleSimSeg) < 5 * sigmaResAngle_ and
0227             fabs(bestRecHitLocalPos.x() - posSimSeg) < 5 * sigmaResPos_) {
0228           recHitFound = true;
0229         }
0230
0231         // Fill Residual histos
0232         HRes2DHit *hRes = nullptr;
0233         if ((*slId).superlayer() == 1 or (*slId).superlayer() == 3) {  // RPhi SL
0234           hRes = histograms.h2DHitRPhi.get();
0235         } else if ((*slId).superlayer() == 2) {  // RZ SL
0236           histograms.h2DHitRZ->fill(angleSimSeg,
0237                                     angleBestRHit,
0238                                     posSimSeg,
0239                                     bestRecHitLocalPos.x(),
0240                                     etaSimSeg,
0241                                     phiSimSeg,
0242                                     sqrt(bestRecHitLocalPosErr.xx()),
0243                                     sqrt(bestRecHitLocalDirErr.xx()));
0244           if (abs((*slId).wheel()) == 0) {
0245             hRes = histograms.h2DHitRZ_W0.get();
0246           } else if (abs((*slId).wheel()) == 1) {
0247             hRes = histograms.h2DHitRZ_W1.get();
0248           } else if (abs((*slId).wheel()) == 2) {
0249             hRes = histograms.h2DHitRZ_W2.get();
0250           }
0251         }
0252         hRes->fill(angleSimSeg,
0253                    angleBestRHit,
0254                    posSimSeg,
0255                    bestRecHitLocalPos.x(),
0256                    etaSimSeg,
0257                    phiSimSeg,
0258                    sqrt(bestRecHitLocalPosErr.xx()),
0259                    sqrt(bestRecHitLocalDirErr.xx()));
0260       }
0261     }  // end of if (nsegm != 0)
0262
0263     // Fill Efficiency plot
0264     HEff2DHit *hEff = nullptr;
0265     if ((*slId).superlayer() == 1 or (*slId).superlayer() == 3) {  // RPhi SL
0266       hEff = histograms.h2DHitEff_RPhi.get();
0267     } else if ((*slId).superlayer() == 2) {  // RZ SL
0268       histograms.h2DHitEff_RZ->fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
0269       if (abs((*slId).wheel()) == 0) {
0270         hEff = histograms.h2DHitEff_RZ_W0.get();
0271       } else if (abs((*slId).wheel()) == 1) {
0272         hEff = histograms.h2DHitEff_RZ_W1.get();
0273       } else if (abs((*slId).wheel()) == 2) {
0274         hEff = histograms.h2DHitEff_RZ_W2.get();
0275       }
0276     }
0277     hEff->fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
0278   }  // End of loop over superlayers
0279 }
0280
0281 // declare this as a framework plugin
0282 #include "FWCore/Framework/interface/MakerMacros.h"
0283 DEFINE_FWK_MODULE(DTSegment2DQuality);