Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:27:04 

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Cerminara - INFN Torino
0005  */
0006 
0007 #include <iostream>
0008 #include <map>
0009 
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0014 #include "Geometry/DTGeometry/interface/DTLayer.h"
0015 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
0016 
0017 #include "DTRecHitQuality.h"
0018 #include "Histograms.h"
0019 
0020 using namespace std;
0021 using namespace edm;
0022 
0023 namespace dtrechit {
0024   struct Histograms {
0025     std::unique_ptr<HRes1DHit> hRes_S1RPhi;          // RecHits, 1. step, RPh
0026     std::unique_ptr<HRes1DHit> hRes_S2RPhi;          // RecHits, 2. step, RPhi
0027     std::unique_ptr<HRes1DHit> hRes_S3RPhi;          // RecHits, 3. step, RPhi
0028     std::unique_ptr<HRes1DHit> hRes_S1RZ;            // RecHits, 1. step, RZ
0029     std::unique_ptr<HRes1DHit> hRes_S2RZ;            // RecHits, 2. step, RZ
0030     std::unique_ptr<HRes1DHit> hRes_S3RZ;            // RecHits, 3. step, RZ
0031     std::unique_ptr<HRes1DHit> hRes_S1RZ_W0;         // RecHits, 1. step, RZ, wheel 0
0032     std::unique_ptr<HRes1DHit> hRes_S2RZ_W0;         // RecHits, 2. step, RZ, wheel 0
0033     std::unique_ptr<HRes1DHit> hRes_S3RZ_W0;         // RecHits, 3. step, RZ, wheel 0
0034     std::unique_ptr<HRes1DHit> hRes_S1RZ_W1;         // RecHits, 1. step, RZ, wheel +-1
0035     std::unique_ptr<HRes1DHit> hRes_S2RZ_W1;         // RecHits, 2. step, RZ, wheel +-1
0036     std::unique_ptr<HRes1DHit> hRes_S3RZ_W1;         // RecHits, 3. step, RZ, wheel +-1
0037     std::unique_ptr<HRes1DHit> hRes_S1RZ_W2;         // RecHits, 1. step, RZ, wheel +-2
0038     std::unique_ptr<HRes1DHit> hRes_S2RZ_W2;         // RecHits, 2. step, RZ, wheel +-2
0039     std::unique_ptr<HRes1DHit> hRes_S3RZ_W2;         // RecHits, 3. step, RZ, wheel +-2
0040     std::unique_ptr<HRes1DHit> hRes_S1RPhi_W0;       // RecHits, 1. step, RPhi, wheel 0
0041     std::unique_ptr<HRes1DHit> hRes_S2RPhi_W0;       // RecHits, 2. step, RPhi, wheel 0
0042     std::unique_ptr<HRes1DHit> hRes_S3RPhi_W0;       // RecHits, 3. step, RPhi, wheel 0
0043     std::unique_ptr<HRes1DHit> hRes_S1RPhi_W1;       // RecHits, 1. step, RPhi, wheel +-1
0044     std::unique_ptr<HRes1DHit> hRes_S2RPhi_W1;       // RecHits, 2. step, RPhi, wheel +-1
0045     std::unique_ptr<HRes1DHit> hRes_S3RPhi_W1;       // RecHits, 3. step, RPhi, wheel +-1
0046     std::unique_ptr<HRes1DHit> hRes_S1RPhi_W2;       // RecHits, 1. step, RPhi, wheel +-2
0047     std::unique_ptr<HRes1DHit> hRes_S2RPhi_W2;       // RecHits, 2. step, RPhi, wheel +-2
0048     std::unique_ptr<HRes1DHit> hRes_S3RPhi_W2;       // RecHits, 3. step, RPhi, wheel +-2
0049     std::unique_ptr<HRes1DHit> hRes_S3RPhiWS[3][4];  // RecHits, 3. step, by wheel/station
0050     std::unique_ptr<HRes1DHit> hRes_S3RZWS[3][4];    // RecHits, 3. step, by wheel/station
0051 
0052     std::unique_ptr<HEff1DHit> hEff_S1RPhi;          // RecHits, 1. step, RPhi
0053     std::unique_ptr<HEff1DHit> hEff_S2RPhi;          // RecHits, 2. step, RPhi
0054     std::unique_ptr<HEff1DHit> hEff_S3RPhi;          // RecHits, 3. step, RPhi
0055     std::unique_ptr<HEff1DHit> hEff_S1RZ;            // RecHits, 1. step, RZ
0056     std::unique_ptr<HEff1DHit> hEff_S2RZ;            // RecHits, 2. step, RZ
0057     std::unique_ptr<HEff1DHit> hEff_S3RZ;            // RecHits, 3. step, RZ
0058     std::unique_ptr<HEff1DHit> hEff_S1RZ_W0;         // RecHits, 1. step, RZ, wheel 0
0059     std::unique_ptr<HEff1DHit> hEff_S2RZ_W0;         // RecHits, 2. step, RZ, wheel 0
0060     std::unique_ptr<HEff1DHit> hEff_S3RZ_W0;         // RecHits, 3. step, RZ, wheel 0
0061     std::unique_ptr<HEff1DHit> hEff_S1RZ_W1;         // RecHits, 1. step, RZ, wheel +-1
0062     std::unique_ptr<HEff1DHit> hEff_S2RZ_W1;         // RecHits, 2. step, RZ, wheel +-1
0063     std::unique_ptr<HEff1DHit> hEff_S3RZ_W1;         // RecHits, 3. step, RZ, wheel +-1
0064     std::unique_ptr<HEff1DHit> hEff_S1RZ_W2;         // RecHits, 1. step, RZ, wheel +-2
0065     std::unique_ptr<HEff1DHit> hEff_S2RZ_W2;         // RecHits, 2. step, RZ, wheel +-2
0066     std::unique_ptr<HEff1DHit> hEff_S3RZ_W2;         // RecHits, 3. step, RZ, wheel +-2
0067     std::unique_ptr<HEff1DHit> hEff_S1RPhiWS[3][4];  // RecHits, 3. step, by wheel/station
0068     std::unique_ptr<HEff1DHit> hEff_S3RPhiWS[3][4];  // RecHits, 3. step, by wheel/station
0069     std::unique_ptr<HEff1DHit> hEff_S1RZWS[3][4];    // RecHits, 3. step, by wheel/station
0070     std::unique_ptr<HEff1DHit> hEff_S3RZWS[3][4];    // RecHits, 3. step, by wheel/station
0071   };
0072 }  // namespace dtrechit
0073 
0074 using namespace dtrechit;
0075 
0076 // In phi SLs, The dependency on X and angle is specular in positive
0077 // and negative wheels. Since positive and negative wheels are filled
0078 // together into the same plots, it is useful to mirror negative wheels
0079 // so that the actual dependency can be observerd instead of an artificially
0080 // simmetrized one.
0081 // Set mirrorMinusWheels to avoid this.
0082 namespace {
0083   constexpr bool mirrorMinusWheels = true;
0084 }
0085 
0086 // Constructor
0087 DTRecHitQuality::DTRecHitQuality(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0088   // Get the debug parameter for verbose output
0089   debug_ = pset.getUntrackedParameter<bool>("debug");
0090   // the name of the simhit collection
0091   simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
0092   simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
0093   // the name of the 1D rec hit collection
0094   recHitLabel_ = pset.getUntrackedParameter<InputTag>("recHitLabel");
0095   recHitToken_ = consumes<DTRecHitCollection>(pset.getUntrackedParameter<InputTag>("recHitLabel"));
0096   // the name of the 2D rec hit collection
0097   segment2DLabel_ = pset.getUntrackedParameter<InputTag>("segment2DLabel");
0098   segment2DToken_ = consumes<DTRecSegment2DCollection>(pset.getUntrackedParameter<InputTag>("segment2DLabel"));
0099   // the name of the 4D rec hit collection
0100   segment4DLabel_ = pset.getUntrackedParameter<InputTag>("segment4DLabel");
0101   segment4DToken_ = consumes<DTRecSegment4DCollection>(pset.getUntrackedParameter<InputTag>("segment4DLabel"));
0102   // Switches for analysis at various steps
0103   doStep1_ = pset.getUntrackedParameter<bool>("doStep1", false);
0104   doStep2_ = pset.getUntrackedParameter<bool>("doStep2", false);
0105   doStep3_ = pset.getUntrackedParameter<bool>("doStep3", false);
0106   doall_ = pset.getUntrackedParameter<bool>("doall", false);
0107   local_ = pset.getUntrackedParameter<bool>("local", true);
0108 }
0109 
0110 void DTRecHitQuality::bookHistograms(DQMStore::IBooker &booker,
0111                                      edm::Run const &run,
0112                                      edm::EventSetup const &setup,
0113                                      Histograms &histograms) const {
0114   if (doall_ && doStep1_) {
0115     histograms.hRes_S1RPhi = std::make_unique<HRes1DHit>("S1RPhi", booker, true, local_);  // RecHits, 1. step, RPhi
0116     histograms.hRes_S1RPhi_W0 =
0117         std::make_unique<HRes1DHit>("S1RPhi_W0", booker, true, local_);  // RecHits, 1. step, RZ, wheel 0
0118     histograms.hRes_S1RPhi_W1 =
0119         std::make_unique<HRes1DHit>("S1RPhi_W1", booker, true, local_);  // RecHits, 1. step, RZ, wheel +-1
0120     histograms.hRes_S1RPhi_W2 =
0121         std::make_unique<HRes1DHit>("S1RPhi_W2", booker, true, local_);  // RecHits, 1. step, RZ, wheel +-2
0122     histograms.hRes_S1RZ = std::make_unique<HRes1DHit>("S1RZ", booker, true, local_);  // RecHits, 1. step, RZ
0123     histograms.hRes_S1RZ_W0 =
0124         std::make_unique<HRes1DHit>("S1RZ_W0", booker, true, local_);  // RecHits, 1. step, RZ, wheel 0
0125     histograms.hRes_S1RZ_W1 =
0126         std::make_unique<HRes1DHit>("S1RZ_W1", booker, true, local_);  // RecHits, 1. step, RZ, wheel +-1
0127     histograms.hRes_S1RZ_W2 =
0128         std::make_unique<HRes1DHit>("S1RZ_W2", booker, true, local_);          // RecHits, 1. step, RZ, wheel +-2
0129     histograms.hEff_S1RPhi = std::make_unique<HEff1DHit>("S1RPhi", booker);    // RecHits, 1. step, RPhi
0130     histograms.hEff_S1RZ = std::make_unique<HEff1DHit>("S1RZ", booker);        // RecHits, 1. step, RZ
0131     histograms.hEff_S1RZ_W0 = std::make_unique<HEff1DHit>("S1RZ_W0", booker);  // RecHits, 1. step, RZ, wheel 0
0132     histograms.hEff_S1RZ_W1 = std::make_unique<HEff1DHit>("S1RZ_W1", booker);  // RecHits, 1. step, RZ, wheel +-1
0133     histograms.hEff_S1RZ_W2 = std::make_unique<HEff1DHit>("S1RZ_W2", booker);  // RecHits, 1. step, RZ, wheel +-2
0134   }
0135   if (doall_ && doStep2_) {
0136     histograms.hRes_S2RPhi = std::make_unique<HRes1DHit>("S2RPhi", booker, true, local_);  // RecHits, 2. step, RPhi
0137     histograms.hRes_S2RPhi_W0 =
0138         std::make_unique<HRes1DHit>("S2RPhi_W0", booker, true, local_);  // RecHits, 2. step, RPhi, wheel 0
0139     histograms.hRes_S2RPhi_W1 =
0140         std::make_unique<HRes1DHit>("S2RPhi_W1", booker, true, local_);  // RecHits, 2. step, RPhi, wheel +-1
0141     histograms.hRes_S2RPhi_W2 =
0142         std::make_unique<HRes1DHit>("S2RPhi_W2", booker, true, local_);  // RecHits, 2. step, RPhi, wheel +-2
0143     histograms.hRes_S2RZ = std::make_unique<HRes1DHit>("S2RZ", booker, true, local_);  // RecHits, 2. step, RZ
0144     histograms.hRes_S2RZ_W0 =
0145         std::make_unique<HRes1DHit>("S2RZ_W0", booker, true, local_);  // RecHits, 2. step, RZ, wheel 0
0146     histograms.hRes_S2RZ_W1 =
0147         std::make_unique<HRes1DHit>("S2RZ_W1", booker, true, local_);  // RecHits, 2. step, RZ, wheel +-1
0148     histograms.hRes_S2RZ_W2 =
0149         std::make_unique<HRes1DHit>("S2RZ_W2", booker, true, local_);          // RecHits, 2. step, RZ, wheel +-2
0150     histograms.hEff_S2RPhi = std::make_unique<HEff1DHit>("S2RPhi", booker);    // RecHits, 2. step, RPhi
0151     histograms.hEff_S2RZ_W0 = std::make_unique<HEff1DHit>("S2RZ_W0", booker);  // RecHits, 2. step, RZ, wheel 0
0152     histograms.hEff_S2RZ_W1 = std::make_unique<HEff1DHit>("S2RZ_W1", booker);  // RecHits, 2. step, RZ, wheel +-1
0153     histograms.hEff_S2RZ_W2 = std::make_unique<HEff1DHit>("S2RZ_W2", booker);  // RecHits, 2. step, RZ, wheel +-2
0154     histograms.hEff_S2RZ = std::make_unique<HEff1DHit>("S2RZ", booker);        // RecHits, 2. step, RZ
0155   }
0156   if (doStep3_) {
0157     histograms.hRes_S3RPhi = std::make_unique<HRes1DHit>("S3RPhi", booker, doall_, local_);  // RecHits, 3. step, RPhi
0158     histograms.hRes_S3RPhi_W0 =
0159         std::make_unique<HRes1DHit>("S3RPhi_W0", booker, doall_, local_);  // RecHits, 3. step, RPhi, wheel 0
0160     histograms.hRes_S3RPhi_W1 = std::make_unique<HRes1DHit>("S3RPhi_W1",
0161                                                             booker,
0162                                                             doall_,
0163                                                             local_);  // RecHits, 3. step, RPhi, wheel +-1
0164     histograms.hRes_S3RPhi_W2 = std::make_unique<HRes1DHit>("S3RPhi_W2",
0165                                                             booker,
0166                                                             doall_,
0167                                                             local_);  // RecHits, 3. step, RPhi, wheel +-2
0168     histograms.hRes_S3RZ = std::make_unique<HRes1DHit>("S3RZ", booker, doall_, local_);  // RecHits, 3. step, RZ
0169     histograms.hRes_S3RZ_W0 =
0170         std::make_unique<HRes1DHit>("S3RZ_W0", booker, doall_, local_);  // RecHits, 3. step, RZ, wheel 0
0171     histograms.hRes_S3RZ_W1 =
0172         std::make_unique<HRes1DHit>("S3RZ_W1", booker, doall_, local_);  // RecHits, 3. step, RZ, wheel +-1
0173     histograms.hRes_S3RZ_W2 =
0174         std::make_unique<HRes1DHit>("S3RZ_W2", booker, doall_, local_);  // RecHits, 3. step, RZ, wheel +-2
0175 
0176     if (local_) {
0177       // Plots with finer granularity, not to be included in DQM
0178       TString name1 = "RPhi_W";
0179       TString name2 = "RZ_W";
0180       for (long w = 0; w <= 2; ++w) {
0181         for (long s = 1; s <= 4; ++s) {
0182           histograms.hRes_S3RPhiWS[w][s - 1] =
0183               std::make_unique<HRes1DHit>(("S3" + name1 + w + "_St" + s).Data(), booker, doall_, local_);
0184           histograms.hEff_S1RPhiWS[w][s - 1] =
0185               std::make_unique<HEff1DHit>(("S1" + name1 + w + "_St" + s).Data(), booker);
0186           histograms.hEff_S3RPhiWS[w][s - 1] =
0187               std::make_unique<HEff1DHit>(("S3" + name1 + w + "_St" + s).Data(), booker);
0188           if (s != 4) {
0189             histograms.hRes_S3RZWS[w][s - 1] =
0190                 std::make_unique<HRes1DHit>(("S3" + name2 + w + "_St" + s).Data(), booker, doall_, local_);
0191             histograms.hEff_S1RZWS[w][s - 1] =
0192                 std::make_unique<HEff1DHit>(("S1" + name2 + w + "_St" + s).Data(), booker);
0193             histograms.hEff_S3RZWS[w][s - 1] =
0194                 std::make_unique<HEff1DHit>(("S3" + name2 + w + "_St" + s).Data(), booker);
0195           }
0196         }
0197       }
0198     }
0199 
0200     if (doall_) {
0201       histograms.hEff_S3RPhi = std::make_unique<HEff1DHit>("S3RPhi", booker);    // RecHits, 3. step, RPhi
0202       histograms.hEff_S3RZ = std::make_unique<HEff1DHit>("S3RZ", booker);        // RecHits, 3. step, RZ
0203       histograms.hEff_S3RZ_W0 = std::make_unique<HEff1DHit>("S3RZ_W0", booker);  // RecHits, 3. step, RZ, wheel 0
0204       histograms.hEff_S3RZ_W1 = std::make_unique<HEff1DHit>("S3RZ_W1", booker);  // RecHits, 3. step, RZ, wheel +-1
0205       histograms.hEff_S3RZ_W2 = std::make_unique<HEff1DHit>("S3RZ_W2", booker);  // RecHits, 3. step, RZ, wheel +-2
0206     }
0207   }
0208 }
0209 
0210 // The real analysis
0211 void DTRecHitQuality::dqmAnalyze(edm::Event const &event,
0212                                  edm::EventSetup const &setup,
0213                                  Histograms const &histograms) const {
0214   if (debug_) {
0215     cout << "--- [DTRecHitQuality] Analysing Event: #Run: " << event.id().run() << " #Event: " << event.id().event()
0216          << endl;
0217   }
0218 
0219   // Get the DT Geometry
0220   const DTGeometry &dtGeom = setup.getData(muonGeomToken_);
0221 
0222   // Get the SimHit collection from the event
0223   Handle<PSimHitContainer> simHits;
0224   event.getByToken(simHitToken_, simHits);
0225 
0226   // Map simhits per wire
0227   map<DTWireId, PSimHitContainer> simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(*(simHits.product()));
0228 
0229   //=======================================================================================
0230   // RecHit analysis at Step 1
0231   if (doStep1_ && doall_) {
0232     if (debug_) {
0233       cout << "  -- DTRecHit S1: begin analysis:" << endl;
0234     }
0235     // Get the rechit collection from the event
0236     Handle<DTRecHitCollection> dtRecHits;
0237     event.getByToken(recHitToken_, dtRecHits);
0238 
0239     if (!dtRecHits.isValid()) {
0240       if (debug_) {
0241         cout << "[DTRecHitQuality]**Warning: no 1DRechits with label: " << recHitLabel_ << " in this event, skipping!"
0242              << endl;
0243       }
0244       return;
0245     }
0246 
0247     // Map rechits per wire
0248     auto const &recHitsPerWire = map1DRecHitsPerWire(dtRecHits.product());
0249     compute(dtGeom, simHitsPerWire, recHitsPerWire, histograms, 1);
0250   }
0251 
0252   //=======================================================================================
0253   // RecHit analysis at Step 2
0254   if (doStep2_ && doall_) {
0255     if (debug_) {
0256       cout << "  -- DTRecHit S2: begin analysis:" << endl;
0257     }
0258 
0259     // Get the 2D rechits from the event
0260     Handle<DTRecSegment2DCollection> segment2Ds;
0261     event.getByToken(segment2DToken_, segment2Ds);
0262 
0263     if (!segment2Ds.isValid()) {
0264       if (debug_) {
0265         cout << "[DTRecHitQuality]**Warning: no 2DSegments with label: " << segment2DLabel_
0266              << " in this event, skipping!" << endl;
0267       }
0268 
0269     } else {
0270       // Map rechits per wire
0271       auto const &recHitsPerWire = map1DRecHitsPerWire(segment2Ds.product());
0272       compute(dtGeom, simHitsPerWire, recHitsPerWire, histograms, 2);
0273     }
0274   }
0275 
0276   //=======================================================================================
0277   // RecHit analysis at Step 3
0278   if (doStep3_) {
0279     if (debug_) {
0280       cout << "  -- DTRecHit S3: begin analysis:" << endl;
0281     }
0282 
0283     // Get the 4D rechits from the event
0284     Handle<DTRecSegment4DCollection> segment4Ds;
0285     event.getByToken(segment4DToken_, segment4Ds);
0286 
0287     if (!segment4Ds.isValid()) {
0288       if (debug_) {
0289         cout << "[DTRecHitQuality]**Warning: no 4D Segments with label: " << segment4DLabel_
0290              << " in this event, skipping!" << endl;
0291       }
0292       return;
0293     }
0294 
0295     // Map rechits per wire
0296     auto const &recHitsPerWire = map1DRecHitsPerWire(segment4Ds.product());
0297     compute(dtGeom, simHitsPerWire, recHitsPerWire, histograms, 3);
0298   }
0299 }
0300 
0301 // Return a map between DTRecHit1DPair and wireId
0302 map<DTWireId, vector<DTRecHit1DPair>> DTRecHitQuality::map1DRecHitsPerWire(
0303     const DTRecHitCollection *dt1DRecHitPairs) const {
0304   map<DTWireId, vector<DTRecHit1DPair>> ret;
0305 
0306   for (const auto &dt1DRecHitPair : *dt1DRecHitPairs) {
0307     ret[dt1DRecHitPair.wireId()].push_back(dt1DRecHitPair);
0308   }
0309 
0310   return ret;
0311 }
0312 
0313 // Return a map between DTRecHit1D at S2 and wireId
0314 map<DTWireId, vector<DTRecHit1D>> DTRecHitQuality::map1DRecHitsPerWire(
0315     const DTRecSegment2DCollection *segment2Ds) const {
0316   map<DTWireId, vector<DTRecHit1D>> ret;
0317 
0318   // Loop over all 2D segments
0319   for (const auto &segment2D : *segment2Ds) {
0320     vector<DTRecHit1D> component1DHits = segment2D.specificRecHits();
0321     // Loop over all component 1D hits
0322     for (auto &component1DHit : component1DHits) {
0323       ret[component1DHit.wireId()].push_back(component1DHit);
0324     }
0325   }
0326   return ret;
0327 }
0328 
0329 // Return a map between DTRecHit1D at S3 and wireId
0330 map<DTWireId, std::vector<DTRecHit1D>> DTRecHitQuality::map1DRecHitsPerWire(
0331     const DTRecSegment4DCollection *segment4Ds) const {
0332   map<DTWireId, vector<DTRecHit1D>> ret;
0333   // Loop over all 4D segments
0334   for (const auto &segment4D : *segment4Ds) {
0335     // Get component 2D segments
0336     vector<const TrackingRecHit *> segment2Ds = segment4D.recHits();
0337     // Loop over 2D segments:
0338     for (auto &segment2D : segment2Ds) {
0339       // Get 1D component rechits
0340       vector<const TrackingRecHit *> hits = segment2D->recHits();
0341       // Loop over them
0342       for (auto &hit : hits) {
0343         const auto *hit1D = dynamic_cast<const DTRecHit1D *>(hit);
0344         ret[hit1D->wireId()].push_back(*hit1D);
0345       }
0346     }
0347   }
0348 
0349   return ret;
0350 }
0351 
0352 // Compute SimHit distance from wire (cm)
0353 float DTRecHitQuality::simHitDistFromWire(const DTLayer *layer, const DTWireId &wireId, const PSimHit &hit) const {
0354   float xwire = layer->specificTopology().wirePosition(wireId.wire());
0355   LocalPoint entryP = hit.entryPoint();
0356   LocalPoint exitP = hit.exitPoint();
0357   float xEntry = entryP.x() - xwire;
0358   float xExit = exitP.x() - xwire;
0359 
0360   return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));  // FIXME: check...
0361 }
0362 
0363 // Compute SimHit impact angle (in direction perp to wire), in the SL RF
0364 float DTRecHitQuality::simHitImpactAngle(const DTLayer *layer, const DTWireId &wireId, const PSimHit &hit) const {
0365   LocalPoint entryP = hit.entryPoint();
0366   LocalPoint exitP = hit.exitPoint();
0367   float theta = (exitP.x() - entryP.x()) / (exitP.z() - entryP.z());
0368   return atan(theta);
0369 }
0370 
0371 // Compute SimHit distance from FrontEnd
0372 float DTRecHitQuality::simHitDistFromFE(const DTLayer *layer, const DTWireId &wireId, const PSimHit &hit) const {
0373   LocalPoint entryP = hit.entryPoint();
0374   LocalPoint exitP = hit.exitPoint();
0375   float wireLenght = layer->specificTopology().cellLenght();
0376   // FIXME: should take only wireLenght/2.;
0377   // moreover, pos+cellLenght/2. is shorter than the distance from FE.
0378   // In fact it would make more sense to make plots vs y.
0379   return (entryP.y() + exitP.y()) / 2. + wireLenght;
0380 }
0381 
0382 // Find the RecHit closest to the muon SimHit
0383 template <typename type>
0384 const type *DTRecHitQuality::findBestRecHit(const DTLayer *layer,
0385                                             const DTWireId &wireId,
0386                                             const vector<type> &recHits,
0387                                             const float simHitDist) const {
0388   float res = 99999;
0389   const type *theBestRecHit = nullptr;
0390   // Loop over RecHits within the cell
0391   for (auto recHit = recHits.begin(); recHit != recHits.end(); recHit++) {
0392     float distTmp = recHitDistFromWire(*recHit, layer);
0393     if (fabs(distTmp - simHitDist) < res) {
0394       res = fabs(distTmp - simHitDist);
0395       theBestRecHit = &(*recHit);
0396     }
0397   }  // End of loop over RecHits within the cell
0398 
0399   return theBestRecHit;
0400 }
0401 
0402 // Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
0403 float DTRecHitQuality::recHitDistFromWire(const DTRecHit1DPair &hitPair, const DTLayer *layer) const {
0404   // Compute the rechit distance from wire
0405   return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
0406 }
0407 
0408 // Compute the distance from wire (cm) of a hits in a DTRecHit1D
0409 float DTRecHitQuality::recHitDistFromWire(const DTRecHit1D &recHit, const DTLayer *layer) const {
0410   return fabs(recHit.localPosition().x() - layer->specificTopology().wirePosition(recHit.wireId().wire()));
0411 }
0412 
0413 template <typename type>
0414 void DTRecHitQuality::compute(const DTGeometry &dtGeom,
0415                               const std::map<DTWireId, std::vector<PSimHit>> &simHitsPerWire,
0416                               const std::map<DTWireId, std::vector<type>> &recHitsPerWire,
0417                               Histograms const &histograms,
0418                               int step) const {
0419   // Loop over cells with a muon SimHit
0420   for (const auto &wireAndSHits : simHitsPerWire) {
0421     DTWireId wireId = wireAndSHits.first;
0422     int wheel = wireId.wheel();
0423     int sl = wireId.superLayer();
0424 
0425     vector<PSimHit> simHitsInCell = wireAndSHits.second;
0426 
0427     // Get the layer
0428     const DTLayer *layer = dtGeom.layer(wireId);
0429 
0430     // Look for a mu hit in the cell
0431     const PSimHit *muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
0432     if (muSimHit == nullptr) {
0433       if (debug_) {
0434         cout << "   No mu SimHit in channel: " << wireId << ", skipping! " << endl;
0435       }
0436       continue;  // Skip this cell
0437     }
0438 
0439     // Find the distance of the simhit from the wire
0440     float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
0441     // Skip simhits out of the cell
0442     if (simHitWireDist > 2.1) {
0443       if (debug_) {
0444         cout << "  [DTRecHitQuality]###Warning: The mu SimHit in out of the "
0445                 "cell, skipping!"
0446              << endl;
0447       }
0448       continue;  // Skip this cell
0449     }
0450     GlobalPoint simHitGlobalPos = layer->toGlobal(muSimHit->localPosition());
0451 
0452     // find SH impact angle
0453     float simHitTheta = simHitImpactAngle(layer, wireId, *muSimHit);
0454 
0455     // find SH distance from FE
0456     float simHitFEDist = simHitDistFromFE(layer, wireId, *muSimHit);
0457 
0458     bool recHitReconstructed = false;
0459 
0460     // Look for RecHits in the same cell
0461     if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
0462       // No RecHit found in this cell
0463       if (debug_) {
0464         cout << "   No RecHit found at Step: " << step << " in cell: " << wireId << endl;
0465       }
0466     } else {
0467       recHitReconstructed = true;
0468       // vector<type> recHits = (*wireAndRecHits).second;
0469       const vector<type> &recHits = recHitsPerWire.at(wireId);
0470       if (debug_) {
0471         cout << "   " << recHits.size() << " RecHits, Step " << step << " in channel: " << wireId << endl;
0472       }
0473 
0474       // Find the best RecHit
0475       const type *theBestRecHit = findBestRecHit(layer, wireId, recHits, simHitWireDist);
0476 
0477       float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
0478       if (debug_) {
0479         cout << "    SimHit distance from wire: " << simHitWireDist << endl
0480              << "    SimHit distance from FE:   " << simHitFEDist << endl
0481              << "    SimHit angle in layer RF:  " << simHitTheta << endl
0482              << "    RecHit distance from wire: " << recHitWireDist << endl;
0483       }
0484       float recHitErr = recHitPositionError(*theBestRecHit);
0485       HRes1DHit *hRes = nullptr;
0486       HRes1DHit *hResTot = nullptr;
0487 
0488       // Mirror angle in phi so that + and - wheels can be plotted together
0489       if (mirrorMinusWheels && wheel < 0 && sl != 2) {
0490         simHitTheta *= -1.;
0491         // Note: local X, if used, would have to be mirrored as well
0492       }
0493 
0494       // Fill residuals and pulls
0495       // Select the histo to be filled
0496       if (step == 1) {
0497         // Step 1
0498         if (sl != 2) {
0499           hResTot = histograms.hRes_S1RPhi.get();
0500           if (wheel == 0) {
0501             hRes = histograms.hRes_S1RPhi_W0.get();
0502           }
0503           if (abs(wheel) == 1) {
0504             hRes = histograms.hRes_S1RPhi_W1.get();
0505           }
0506           if (abs(wheel) == 2) {
0507             hRes = histograms.hRes_S1RPhi_W2.get();
0508           }
0509         } else {
0510           hResTot = histograms.hRes_S1RZ.get();
0511           if (wheel == 0) {
0512             hRes = histograms.hRes_S1RZ_W0.get();
0513           }
0514           if (abs(wheel) == 1) {
0515             hRes = histograms.hRes_S1RZ_W1.get();
0516           }
0517           if (abs(wheel) == 2) {
0518             hRes = histograms.hRes_S1RZ_W2.get();
0519           }
0520         }
0521 
0522       } else if (step == 2) {
0523         // Step 2
0524         if (sl != 2) {
0525           hRes = histograms.hRes_S2RPhi.get();
0526           if (wheel == 0) {
0527             hRes = histograms.hRes_S2RPhi_W0.get();
0528           }
0529           if (abs(wheel) == 1) {
0530             hRes = histograms.hRes_S2RPhi_W1.get();
0531           }
0532           if (abs(wheel) == 2) {
0533             hRes = histograms.hRes_S2RPhi_W2.get();
0534           }
0535         } else {
0536           hResTot = histograms.hRes_S2RZ.get();
0537           if (wheel == 0) {
0538             hRes = histograms.hRes_S2RZ_W0.get();
0539           }
0540           if (abs(wheel) == 1) {
0541             hRes = histograms.hRes_S2RZ_W1.get();
0542           }
0543           if (abs(wheel) == 2) {
0544             hRes = histograms.hRes_S2RZ_W2.get();
0545           }
0546         }
0547 
0548       } else if (step == 3) {
0549         // Step 3
0550         if (sl != 2) {
0551           hResTot = histograms.hRes_S3RPhi.get();
0552           if (wheel == 0) {
0553             hRes = histograms.hRes_S3RPhi_W0.get();
0554           }
0555           if (abs(wheel) == 1) {
0556             hRes = histograms.hRes_S3RPhi_W1.get();
0557           }
0558           if (abs(wheel) == 2) {
0559             hRes = histograms.hRes_S3RPhi_W2.get();
0560           }
0561           if (local_) {
0562             histograms.hRes_S3RPhiWS[abs(wheel)][wireId.station() - 1]->fill(simHitWireDist,
0563                                                                              simHitTheta,
0564                                                                              simHitFEDist,
0565                                                                              recHitWireDist,
0566                                                                              simHitGlobalPos.eta(),
0567                                                                              simHitGlobalPos.phi(),
0568                                                                              recHitErr,
0569                                                                              wireId.station());
0570           }
0571         } else {
0572           hResTot = histograms.hRes_S3RZ.get();
0573           if (wheel == 0) {
0574             hRes = histograms.hRes_S3RZ_W0.get();
0575           }
0576           if (abs(wheel) == 1) {
0577             hRes = histograms.hRes_S3RZ_W1.get();
0578           }
0579           if (abs(wheel) == 2) {
0580             hRes = histograms.hRes_S3RZ_W2.get();
0581           }
0582           if (local_) {
0583             histograms.hRes_S3RZWS[abs(wheel)][wireId.station() - 1]->fill(simHitWireDist,
0584                                                                            simHitTheta,
0585                                                                            simHitFEDist,
0586                                                                            recHitWireDist,
0587                                                                            simHitGlobalPos.eta(),
0588                                                                            simHitGlobalPos.phi(),
0589                                                                            recHitErr,
0590                                                                            wireId.station());
0591           }
0592         }
0593       }
0594 
0595       // Fill
0596       hRes->fill(simHitWireDist,
0597                  simHitTheta,
0598                  simHitFEDist,
0599                  recHitWireDist,
0600                  simHitGlobalPos.eta(),
0601                  simHitGlobalPos.phi(),
0602                  recHitErr,
0603                  wireId.station());
0604       if (hResTot != nullptr) {
0605         hResTot->fill(simHitWireDist,
0606                       simHitTheta,
0607                       simHitFEDist,
0608                       recHitWireDist,
0609                       simHitGlobalPos.eta(),
0610                       simHitGlobalPos.phi(),
0611                       recHitErr,
0612                       wireId.station());
0613       }
0614     }
0615 
0616     // Fill Efficiencies
0617     if (doall_) {
0618       HEff1DHit *hEff = nullptr;
0619       HEff1DHit *hEffTot = nullptr;
0620       if (step == 1) {
0621         // Step 1
0622         if (sl != 2) {
0623           hEff = histograms.hEff_S1RPhi.get();
0624           if (local_) {
0625             histograms.hEff_S1RPhiWS[abs(wheel)][wireId.station() - 1]->fill(
0626                 simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
0627           }
0628         } else {
0629           hEffTot = histograms.hEff_S1RZ.get();
0630           if (wheel == 0) {
0631             hEff = histograms.hEff_S1RZ_W0.get();
0632           }
0633           if (abs(wheel) == 1) {
0634             hEff = histograms.hEff_S1RZ_W1.get();
0635           }
0636           if (abs(wheel) == 2) {
0637             hEff = histograms.hEff_S1RZ_W2.get();
0638           }
0639           if (local_) {
0640             histograms.hEff_S1RZWS[abs(wheel)][wireId.station() - 1]->fill(
0641                 simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
0642           }
0643         }
0644 
0645       } else if (step == 2) {
0646         // Step 2
0647         if (sl != 2) {
0648           hEff = histograms.hEff_S2RPhi.get();
0649         } else {
0650           hEffTot = histograms.hEff_S2RZ.get();
0651           if (wheel == 0) {
0652             hEff = histograms.hEff_S2RZ_W0.get();
0653           }
0654           if (abs(wheel) == 1) {
0655             hEff = histograms.hEff_S2RZ_W1.get();
0656           }
0657           if (abs(wheel) == 2) {
0658             hEff = histograms.hEff_S2RZ_W2.get();
0659           }
0660         }
0661 
0662       } else if (step == 3) {
0663         // Step 3
0664         if (sl != 2) {
0665           hEff = histograms.hEff_S3RPhi.get();
0666           if (local_) {
0667             histograms.hEff_S3RPhiWS[abs(wheel)][wireId.station() - 1]->fill(
0668                 simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
0669           }
0670         } else {
0671           hEffTot = histograms.hEff_S3RZ.get();
0672           if (wheel == 0) {
0673             hEff = histograms.hEff_S3RZ_W0.get();
0674           }
0675           if (abs(wheel) == 1) {
0676             hEff = histograms.hEff_S3RZ_W1.get();
0677           }
0678           if (abs(wheel) == 2) {
0679             hEff = histograms.hEff_S3RZ_W2.get();
0680           }
0681           if (local_) {
0682             histograms.hEff_S3RZWS[abs(wheel)][wireId.station() - 1]->fill(
0683                 simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
0684           }
0685         }
0686       }
0687       // Fill
0688       hEff->fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
0689       if (hEffTot != nullptr) {
0690         hEffTot->fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
0691       }
0692     }
0693   }
0694 }
0695 
0696 // Return the error on the measured (cm) coordinate
0697 float DTRecHitQuality::recHitPositionError(const DTRecHit1DPair &recHit) const {
0698   return sqrt(recHit.localPositionError(DTEnums::Left).xx());
0699 }
0700 
0701 // Return the error on the measured (cm) coordinate
0702 float DTRecHitQuality::recHitPositionError(const DTRecHit1D &recHit) const {
0703   return sqrt(recHit.localPositionError().xx());
0704 }
0705 
0706 // declare this as a framework plugin
0707 #include "FWCore/Framework/interface/MakerMacros.h"
0708 DEFINE_FWK_MODULE(DTRecHitQuality);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!