Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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 "DTSegment4DQuality.h"
0022 #include "Histograms.h"
0023 
0024 using namespace std;
0025 using namespace edm;
0026 
0027 namespace dtsegment4d {
0028   struct Histograms {
0029     std::unique_ptr<HRes4DHit> h4DHit;
0030     std::unique_ptr<HRes4DHit> h4DHit_W0;
0031     std::unique_ptr<HRes4DHit> h4DHit_W1;
0032     std::unique_ptr<HRes4DHit> h4DHit_W2;
0033     std::unique_ptr<HRes4DHit> h4DHitWS[3][4];
0034 
0035     std::unique_ptr<HEff4DHit> hEff_All;
0036     std::unique_ptr<HEff4DHit> hEff_W0;
0037     std::unique_ptr<HEff4DHit> hEff_W1;
0038     std::unique_ptr<HEff4DHit> hEff_W2;
0039     std::unique_ptr<HEff4DHit> hEffWS[3][4];
0040   };
0041 }  // namespace dtsegment4d
0042 using namespace dtsegment4d;
0043 
0044 // In phi SLs, The dependency on X and angle is specular in positive
0045 // and negative wheels. Since positive and negative wheels are filled
0046 // together into the same plots, it is useful to mirror negative wheels
0047 // so that the actual dependency can be observerd instead of an artificially
0048 // simmetrized one.
0049 // Set mirrorMinusWheels to avoid this.
0050 namespace {
0051   constexpr bool mirrorMinusWheels = true;
0052 }
0053 
0054 // Constructor
0055 DTSegment4DQuality::DTSegment4DQuality(const ParameterSet &pset) : muonGeomToken_(esConsumes()) {
0056   // Get the debug parameter for verbose output
0057   debug_ = pset.getUntrackedParameter<bool>("debug");
0058   DTHitQualityUtils::debug = debug_;
0059 
0060   // the name of the simhit collection
0061   simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
0062   simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
0063   // the name of the 2D rec hit collection
0064   segment4DLabel_ = pset.getUntrackedParameter<InputTag>("segment4DLabel");
0065   segment4DToken_ = consumes<DTRecSegment4DCollection>(pset.getUntrackedParameter<InputTag>("segment4DLabel"));
0066 
0067   // sigma resolution on position
0068   sigmaResX_ = pset.getParameter<double>("sigmaResX");
0069   sigmaResY_ = pset.getParameter<double>("sigmaResY");
0070   // sigma resolution on angle
0071   sigmaResAlpha_ = pset.getParameter<double>("sigmaResAlpha");
0072   sigmaResBeta_ = pset.getParameter<double>("sigmaResBeta");
0073   doall_ = pset.getUntrackedParameter<bool>("doall", false);
0074   local_ = pset.getUntrackedParameter<bool>("local", false);
0075 }
0076 
0077 void DTSegment4DQuality::bookHistograms(DQMStore::IBooker &booker,
0078                                         edm::Run const &run,
0079                                         edm::EventSetup const &setup,
0080                                         Histograms &histograms) const {
0081   histograms.h4DHit = std::make_unique<HRes4DHit>("All", booker, doall_, local_);
0082   histograms.h4DHit_W0 = std::make_unique<HRes4DHit>("W0", booker, doall_, local_);
0083   histograms.h4DHit_W1 = std::make_unique<HRes4DHit>("W1", booker, doall_, local_);
0084   histograms.h4DHit_W2 = std::make_unique<HRes4DHit>("W2", booker, doall_, local_);
0085 
0086   if (doall_) {
0087     histograms.hEff_All = std::make_unique<HEff4DHit>("All", booker);
0088     histograms.hEff_W0 = std::make_unique<HEff4DHit>("W0", booker);
0089     histograms.hEff_W1 = std::make_unique<HEff4DHit>("W1", booker);
0090     histograms.hEff_W2 = std::make_unique<HEff4DHit>("W2", booker);
0091   }
0092 
0093   if (local_) {
0094     // Plots with finer granularity, not to be included in DQM
0095     TString name = "W";
0096     for (long w = 0; w <= 2; ++w) {
0097       for (long s = 1; s <= 4; ++s) {
0098         // FIXME station 4 is not filled
0099         TString nameWS = (name + w + "_St" + s);
0100         histograms.h4DHitWS[w][s - 1] = std::make_unique<HRes4DHit>(nameWS.Data(), booker, doall_, local_);
0101         histograms.hEffWS[w][s - 1] = std::make_unique<HEff4DHit>(nameWS.Data(), booker);
0102       }
0103     }
0104   }
0105 };
0106 
0107 // The real analysis
0108 void DTSegment4DQuality::dqmAnalyze(edm::Event const &event,
0109                                     edm::EventSetup const &setup,
0110                                     Histograms const &histograms) const {
0111   const float epsilon = 5e-5;  // numerical accuracy on angles [rad}
0112 
0113   // Get the DT Geometry
0114   const DTGeometry &dtGeom = setup.getData(muonGeomToken_);
0115 
0116   // Get the SimHit collection from the event
0117   edm::Handle<PSimHitContainer> simHits;
0118   event.getByToken(simHitToken_, simHits);  // FIXME: second string to be removed
0119 
0120   // Map simHits by chamber
0121   map<DTChamberId, PSimHitContainer> simHitsPerCh;
0122   for (const auto &simHit : *simHits) {
0123     // Consider only muon simhits; the others are not considered elsewhere in
0124     // this class!
0125     if (abs(simHit.particleType()) == 13) {
0126       // Create the id of the chamber (the simHits in the DT known their wireId)
0127       DTChamberId chamberId = (((DTWireId(simHit.detUnitId())).layerId()).superlayerId()).chamberId();
0128       // Fill the map
0129       simHitsPerCh[chamberId].push_back(simHit);
0130     }
0131   }
0132 
0133   // Get the 4D rechits from the event
0134   Handle<DTRecSegment4DCollection> segment4Ds;
0135   event.getByToken(segment4DToken_, segment4Ds);
0136 
0137   if (!segment4Ds.isValid()) {
0138     if (debug_) {
0139       cout << "[DTSegment4DQuality]**Warning: no 4D Segments with label: " << segment4DLabel_
0140            << " in this event, skipping!" << endl;
0141     }
0142     return;
0143   }
0144 
0145   // Loop over all chambers containing a (muon) simhit
0146   for (auto &simHitsInChamber : simHitsPerCh) {
0147     DTChamberId chamberId = simHitsInChamber.first;
0148     int station = chamberId.station();
0149     if (station == 4 && !(local_)) {
0150       continue;  // use DTSegment2DSLPhiQuality to analyze MB4 performaces in DQM
0151     }
0152     int wheel = chamberId.wheel();
0153 
0154     //------------------------- simHits ---------------------------//
0155     // Get simHits of this chamber
0156     const PSimHitContainer &simHits = simHitsInChamber.second;
0157 
0158     // Map simhits per wire
0159     auto const &simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
0160     auto const &muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
0161     int nMuSimHit = muSimHitPerWire.size();
0162     if (nMuSimHit < 2) {  // Skip chamber with less than 2 cells with mu hits
0163       continue;
0164     }
0165     if (debug_) {
0166       cout << "=== Chamber " << chamberId << " has " << nMuSimHit << " SimHits" << endl;
0167     }
0168 
0169     // Find outer and inner mu SimHit to build a segment
0170     pair<const PSimHit *, const PSimHit *> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire);
0171 
0172     // Consider only sim segments crossing at least 2 SLs
0173     if ((DTWireId(inAndOutSimHit.first->detUnitId())).superlayer() ==
0174         (DTWireId(inAndOutSimHit.second->detUnitId())).superLayer()) {
0175       continue;
0176     }
0177 
0178     // Find direction and position of the sim Segment in Chamber RF
0179     pair<LocalVector, LocalPoint> dirAndPosSimSegm =
0180         DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit, chamberId, dtGeom);
0181 
0182     LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
0183     LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
0184     const DTChamber *chamber = dtGeom.chamber(chamberId);
0185     GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos);
0186     GlobalVector simSegmGlobalDir = chamber->toGlobal(simSegmLocalDir);
0187 
0188     // phi and theta angle of simulated segment in Chamber RF
0189     float alphaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
0190     float betaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).second;
0191     // x, y position of simulated segment in Chamber RF
0192     float xSimSeg = simSegmLocalPos.x();
0193     float ySimSeg = simSegmLocalPos.y();
0194     // Position (in eta, phi coordinates) in lobal RF
0195     float etaSimSeg = simSegmGlobalPos.eta();
0196     float phiSimSeg = simSegmGlobalPos.phi();
0197 
0198     double count_seg = 0;
0199 
0200     if (debug_) {
0201       cout << "  Simulated segment:  local direction " << simSegmLocalDir << endl
0202            << "                      local position  " << simSegmLocalPos << endl
0203            << "                      alpha           " << alphaSimSeg << endl
0204            << "                      beta            " << betaSimSeg << endl;
0205     }
0206 
0207     //---------------------------- recHits --------------------------//
0208     // Get the range of rechit for the corresponding chamberId
0209     bool recHitFound = false;
0210     DTRecSegment4DCollection::range range = segment4Ds->get(chamberId);
0211     int nsegm = distance(range.first, range.second);
0212     if (debug_) {
0213       cout << "   Chamber: " << chamberId << " has " << nsegm << " 4D segments" << endl;
0214     }
0215 
0216     if (nsegm != 0) {
0217       // Find the best RecHit: look for the 4D RecHit with the phi angle closest
0218       // to that of segment made of SimHits.
0219       // RecHits must have delta alpha and delta position within 5 sigma of
0220       // the residual distribution (we are looking for residuals of segments
0221       // usefull to the track fit) for efficency purpose
0222       const DTRecSegment4D *bestRecHit = nullptr;
0223       double deltaAlpha = 99999;
0224       double deltaBeta = 99999;
0225 
0226       // Loop over the recHits of this chamberId
0227       for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0228         // Consider only segments with both projections
0229         if (station != 4 && (*segment4D).dimension() != 4) {
0230           continue;
0231         }
0232         // Segment Local Direction and position (in Chamber RF)
0233         LocalVector recSegDirection = (*segment4D).localDirection();
0234         LocalPoint recSegPosition = (*segment4D).localPosition();
0235 
0236         pair<double, double> ab = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection);
0237         float recSegAlpha = ab.first;
0238         float recSegBeta = ab.second;
0239 
0240         if (debug_) {
0241           cout << &(*segment4D) << "  RecSegment direction: " << recSegDirection << endl
0242                << "             position : " << (*segment4D).localPosition() << endl
0243                << "             alpha    : " << recSegAlpha << endl
0244                << "             beta     : " << recSegBeta << endl
0245                << "             nhits    : " << (*segment4D).phiSegment()->recHits().size() << " "
0246                << (((*segment4D).zSegment() != nullptr) ? (*segment4D).zSegment()->recHits().size() : 0) << endl;
0247         }
0248 
0249         float dAlphaRecSim = fabs(recSegAlpha - alphaSimSeg);
0250         float dBetaRecSim = fabs(recSegBeta - betaSimSeg);
0251 
0252         if ((fabs(recSegPosition.x() - simSegmLocalPos.x()) <
0253              4)  // require rec and sim segments to be ~in the same cell in x
0254             && ((fabs(recSegPosition.y() - simSegmLocalPos.y()) < 4) ||
0255                 (*segment4D).dimension() < 4)) {  // ~in the same cell in y, if segment has the theta view
0256           ++count_seg;
0257 
0258           if (fabs(dAlphaRecSim - deltaAlpha) < epsilon) {  // Numerically equivalent alphas, choose based on beta
0259             if (dBetaRecSim < deltaBeta) {
0260               deltaAlpha = dAlphaRecSim;
0261               deltaBeta = dBetaRecSim;
0262               bestRecHit = &(*segment4D);
0263             }
0264 
0265           } else if (dAlphaRecSim < deltaAlpha) {
0266             deltaAlpha = dAlphaRecSim;
0267             deltaBeta = dBetaRecSim;
0268             bestRecHit = &(*segment4D);
0269           }
0270         }
0271 
0272       }  // End of Loop over all 4D RecHits
0273 
0274       if (bestRecHit) {
0275         if (debug_) {
0276           cout << endl << "Chosen: " << bestRecHit << endl;
0277         }
0278         // Best rechit direction and position in Chamber RF
0279         LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
0280         LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
0281         // Errors on x and y
0282         LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
0283         LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
0284 
0285         pair<double, double> ab = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir);
0286         float alphaBestRHit = ab.first;
0287         float betaBestRHit = ab.second;
0288         // Errors on alpha and beta
0289 
0290         // Get position and direction using the rx projection (so in SL
0291         // reference frame). Note that x (and y) are swapped wrt to Chamber
0292         // frame
0293         // if (bestRecHit->hasZed()) //
0294         const DTSLRecSegment2D *zedRecSeg = bestRecHit->zSegment();
0295         LocalPoint bestRecHitLocalPosRZ;
0296         LocalVector bestRecHitLocalDirRZ;
0297         LocalError bestRecHitLocalPosErrRZ;
0298         LocalError bestRecHitLocalDirErrRZ;
0299         LocalPoint simSegLocalPosRZ;  // FIXME: this is not set for segments with
0300                                       // only the phi view.
0301         float alphaBestRHitRZ = 0;    // angle measured in the RZ SL, in its own frame
0302         float alphaSimSegRZ = betaSimSeg;
0303         if (zedRecSeg) {
0304           bestRecHitLocalPosRZ = zedRecSeg->localPosition();
0305           bestRecHitLocalDirRZ = zedRecSeg->localDirection();
0306           // Errors on x and y
0307           bestRecHitLocalPosErrRZ = zedRecSeg->localPositionError();
0308           bestRecHitLocalDirErrRZ = zedRecSeg->localDirectionError();
0309 
0310           // angle measured in the RZ SL, in its own frame
0311           alphaBestRHitRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDirRZ).first;
0312 
0313           // Get SimSeg position and Direction in rZ SL frame
0314           const DTSuperLayer *sl = dtGeom.superLayer(zedRecSeg->superLayerId());
0315           LocalPoint simSegLocalPosRZTmp = sl->toLocal(simSegmGlobalPos);
0316           LocalVector simSegLocalDirRZ = sl->toLocal(simSegmGlobalDir);
0317           simSegLocalPosRZ =
0318               simSegLocalPosRZTmp + simSegLocalDirRZ * (-simSegLocalPosRZTmp.z() / (cos(simSegLocalDirRZ.theta())));
0319           alphaSimSegRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegLocalDirRZ).first;
0320 
0321           if (debug_) {
0322             cout << "RZ SL: recPos " << bestRecHitLocalPosRZ << "recDir " << bestRecHitLocalDirRZ << "recAlpha "
0323                  << alphaBestRHitRZ << endl
0324                  << "RZ SL: simPos " << simSegLocalPosRZ << "simDir " << simSegLocalDirRZ << "simAlpha "
0325                  << alphaSimSegRZ << endl;
0326           }
0327         }
0328 
0329         // get nhits and t0
0330         const DTChamberRecSegment2D *phiSeg = bestRecHit->phiSegment();
0331 
0332         float t0phi = -999;
0333         float t0theta = -999;
0334         int nHitPhi = 0;
0335         int nHitTheta = 0;
0336 
0337         if (phiSeg) {
0338           t0phi = phiSeg->t0();
0339           nHitPhi = phiSeg->recHits().size();
0340         }
0341 
0342         if (zedRecSeg) {
0343           t0theta = zedRecSeg->t0();
0344           nHitTheta = zedRecSeg->recHits().size();
0345         }
0346 
0347         recHitFound = true;
0348 
0349         // Mirror alpha in phi SLs so that + and - wheels can be plotted
0350         // together
0351         if (mirrorMinusWheels && wheel < 0) {
0352           alphaSimSeg *= -1.;
0353           alphaBestRHit *= -1.;
0354           // Note: local X (xSimSeg, bestRecHitLocalPos.x() would have to be
0355           // mirrored as well; but at the moment there is no plot of dependency
0356           // vs X, except for efficiency.
0357         }
0358 
0359         // Fill Residual histos
0360         HRes4DHit *histo = nullptr;
0361 
0362         if (wheel == 0) {
0363           histo = histograms.h4DHit_W0.get();
0364         } else if (abs(wheel) == 1) {
0365           histo = histograms.h4DHit_W1.get();
0366         } else if (abs(wheel) == 2) {
0367           histo = histograms.h4DHit_W2.get();
0368         }
0369 
0370         float sigmaAlphaBestRhit = sqrt(DTHitQualityUtils::sigmaAngle(alphaBestRHit, bestRecHitLocalDirErr.xx()));
0371         float sigmaBetaBestRhit =
0372             sqrt(DTHitQualityUtils::sigmaAngle(betaBestRHit,
0373                                                bestRecHitLocalDirErr.yy()));  // FIXME this misses the contribution
0374                                                                               // from uncertainty in extrapolation!
0375         float sigmaAlphaBestRhitRZ = sqrt(DTHitQualityUtils::sigmaAngle(alphaBestRHitRZ, bestRecHitLocalDirErrRZ.xx()));
0376 
0377         histo->fill(alphaSimSeg,
0378                     alphaBestRHit,
0379                     betaSimSeg,
0380                     betaBestRHit,
0381                     xSimSeg,
0382                     bestRecHitLocalPos.x(),
0383                     ySimSeg,
0384                     bestRecHitLocalPos.y(),
0385                     etaSimSeg,
0386                     phiSimSeg,
0387                     bestRecHitLocalPosRZ.x(),
0388                     simSegLocalPosRZ.x(),
0389                     alphaBestRHitRZ,
0390                     alphaSimSegRZ,
0391                     sigmaAlphaBestRhit,
0392                     sigmaBetaBestRhit,
0393                     sqrt(bestRecHitLocalPosErr.xx()),
0394                     sqrt(bestRecHitLocalPosErr.yy()),
0395                     sigmaAlphaBestRhitRZ,
0396                     sqrt(bestRecHitLocalPosErrRZ.xx()),
0397                     nHitPhi,
0398                     nHitTheta,
0399                     t0phi,
0400                     t0theta);
0401 
0402         histograms.h4DHit->fill(alphaSimSeg,
0403                                 alphaBestRHit,
0404                                 betaSimSeg,
0405                                 betaBestRHit,
0406                                 xSimSeg,
0407                                 bestRecHitLocalPos.x(),
0408                                 ySimSeg,
0409                                 bestRecHitLocalPos.y(),
0410                                 etaSimSeg,
0411                                 phiSimSeg,
0412                                 bestRecHitLocalPosRZ.x(),
0413                                 simSegLocalPosRZ.x(),
0414                                 alphaBestRHitRZ,
0415                                 alphaSimSegRZ,
0416                                 sigmaAlphaBestRhit,
0417                                 sigmaBetaBestRhit,
0418                                 sqrt(bestRecHitLocalPosErr.xx()),
0419                                 sqrt(bestRecHitLocalPosErr.yy()),
0420                                 sigmaAlphaBestRhitRZ,
0421                                 sqrt(bestRecHitLocalPosErrRZ.xx()),
0422                                 nHitPhi,
0423                                 nHitTheta,
0424                                 t0phi,
0425                                 t0theta);
0426 
0427         if (local_) {
0428           histograms.h4DHitWS[abs(wheel)][station - 1]->fill(alphaSimSeg,
0429                                                              alphaBestRHit,
0430                                                              betaSimSeg,
0431                                                              betaBestRHit,
0432                                                              xSimSeg,
0433                                                              bestRecHitLocalPos.x(),
0434                                                              ySimSeg,
0435                                                              bestRecHitLocalPos.y(),
0436                                                              etaSimSeg,
0437                                                              phiSimSeg,
0438                                                              bestRecHitLocalPosRZ.x(),
0439                                                              simSegLocalPosRZ.x(),
0440                                                              alphaBestRHitRZ,
0441                                                              alphaSimSegRZ,
0442                                                              sigmaAlphaBestRhit,
0443                                                              sigmaBetaBestRhit,
0444                                                              sqrt(bestRecHitLocalPosErr.xx()),
0445                                                              sqrt(bestRecHitLocalPosErr.yy()),
0446                                                              sigmaAlphaBestRhitRZ,
0447                                                              sqrt(bestRecHitLocalPosErrRZ.xx()),
0448                                                              nHitPhi,
0449                                                              nHitTheta,
0450                                                              t0phi,
0451                                                              t0theta);
0452         }
0453 
0454       }  // end of if (bestRecHit)
0455 
0456     }  // end of if (nsegm!= 0)
0457 
0458     // Fill Efficiency plot
0459     if (doall_) {
0460       HEff4DHit *heff = nullptr;
0461 
0462       if (wheel == 0) {
0463         heff = histograms.hEff_W0.get();
0464       } else if (abs(wheel) == 1) {
0465         heff = histograms.hEff_W1.get();
0466       } else if (abs(wheel) == 2) {
0467         heff = histograms.hEff_W2.get();
0468       }
0469       heff->fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
0470       histograms.hEff_All->fill(
0471           etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
0472       if (local_) {
0473         histograms.hEffWS[abs(wheel)][station - 1]->fill(
0474             etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
0475       }
0476     }
0477   }  // End of loop over chambers
0478 }
0479 
0480 // declare this as a framework plugin
0481 #include "FWCore/Framework/interface/MakerMacros.h"
0482 DEFINE_FWK_MODULE(DTSegment4DQuality);
[ 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!