Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:07:07 

0001 
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author G. Mila - INFN Torino
0006  */
0007 
0008 #include "DQM/DTMonitorModule/interface/DTCalibValidationFromMuons.h"
0009 
0010 // Framework
0011 #include "DQMServices/Core/interface/DQMStore.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 
0018 // Geometry
0019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0020 
0021 // RecHit
0022 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0023 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0024 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0025 
0026 #include <iterator>
0027 
0028 using namespace edm;
0029 using namespace std;
0030 
0031 DTCalibValidationFromMuons::DTCalibValidationFromMuons(const ParameterSet &pset)
0032     : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
0033   parameters = pset;
0034 
0035   // the name of the 4D segments
0036   segment4DToken_ =
0037       consumes<DTRecSegment4DCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("segment4DLabel")));
0038   // muon collection for matching 4D segments to muons
0039   muonToken_ = consumes<reco::MuonCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("muonLabel")));
0040   // the counter of segments not used to compute residuals
0041   wrongSegment = 0;
0042   // the counter of segments used to compute residuals
0043   rightSegment = 0;
0044 
0045   nevent = 0;
0046 }
0047 
0048 DTCalibValidationFromMuons::~DTCalibValidationFromMuons() {
0049   // FR the following was previously in the endJob
0050 
0051   LogVerbatim("DTCalibValidationFromMuons") << "Segments used to compute residuals: " << rightSegment;
0052   LogVerbatim("DTCalibValidationFromMuons") << "Segments not used to compute residuals: " << wrongSegment;
0053 }
0054 
0055 void DTCalibValidationFromMuons::dqmBeginRun(const edm::Run &run, const edm::EventSetup &setup) {
0056   // get the geometry
0057   dtGeom = &setup.getData(muonGeomToken_);
0058 }
0059 
0060 void DTCalibValidationFromMuons::analyze(const edm::Event &event, const edm::EventSetup &setup) {
0061   ++nevent;
0062   LogTrace("DTCalibValidationFromMuons") << "[DTCalibValidationFromMuons] Analyze #Run: " << event.id().run()
0063                                          << " #Event: " << nevent;
0064 
0065   // RecHit mapping at Step 3 ---------------------------------
0066   LogTrace("DTCalibValidationFromMuons") << "  -- DTRecHit S3: begin analysis:";
0067   // Get the 4D rechits from the event
0068   Handle<reco::MuonCollection> muonH;
0069   event.getByToken(muonToken_, muonH);
0070   const vector<reco::Muon> *muons = muonH.product();
0071 
0072   // Get the 4D rechits from the event
0073   Handle<DTRecSegment4DCollection> segment4Ds;
0074   event.getByToken(segment4DToken_, segment4Ds);
0075 
0076   vector<const DTRecSegment4D *> selectedSegment4Ds;
0077 
0078   for (auto &imuon : *muons) {
0079     for (const auto &ch : imuon.matches()) {
0080       DetId chId(ch.id.rawId());
0081       if (chId.det() != DetId::Muon)
0082         continue;
0083       if (chId.subdetId() != MuonSubdetId::DT)
0084         continue;
0085       if (imuon.pt() < 15)
0086         continue;
0087       if (!imuon.isGlobalMuon())
0088         continue;
0089 
0090       int nsegs = ch.segmentMatches.size();
0091       if (!nsegs)
0092         continue;
0093 
0094       // get the DT segments that were used to construct the muon
0095       DTChamberId matchId = ch.id();
0096       DTRecSegment4DCollection::range segs = segment4Ds->get(matchId);
0097       for (DTRecSegment4DCollection::const_iterator segment = segs.first; segment != segs.second; ++segment) {
0098         LocalPoint posHit = segment->localPosition();
0099         float dx = (posHit.x() ? posHit.x() - ch.x : 0);
0100         float dy = (posHit.y() ? posHit.y() - ch.y : 0);
0101         float dr = sqrt(dx * dx + dy * dy);
0102         if (dr < 5)
0103           selectedSegment4Ds.push_back(&(*segment));
0104       }
0105     }
0106   }
0107 
0108   // Loop over all 4D segments
0109   for (auto segment : selectedSegment4Ds) {
0110     LogTrace("DTCalibValidationFromMuons") << "Anlysis on recHit at step 3";
0111     compute(dtGeom, *segment);
0112   }
0113 }
0114 
0115 // Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
0116 float DTCalibValidationFromMuons::recHitDistFromWire(const DTRecHit1DPair &hitPair, const DTLayer *layer) {
0117   return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
0118 }
0119 
0120 // Compute the distance from wire (cm) of a hits in a DTRecHit1D
0121 float DTCalibValidationFromMuons::recHitDistFromWire(const DTRecHit1D &recHit, const DTLayer *layer) {
0122   return fabs(recHit.localPosition().x() - layer->specificTopology().wirePosition(recHit.wireId().wire()));
0123 }
0124 
0125 // Compute the position (cm) of a hits in a DTRecHit1DPair
0126 float DTCalibValidationFromMuons::recHitPosition(
0127     const DTRecHit1DPair &hitPair, const DTLayer *layer, const DTChamber *chamber, float segmentPos, int sl) {
0128   // Get the layer and the wire position
0129   GlobalPoint hitPosGlob_right = layer->toGlobal(hitPair.localPosition(DTEnums::Right));
0130   LocalPoint hitPosInChamber_right = chamber->toLocal(hitPosGlob_right);
0131   GlobalPoint hitPosGlob_left = layer->toGlobal(hitPair.localPosition(DTEnums::Left));
0132   LocalPoint hitPosInChamber_left = chamber->toLocal(hitPosGlob_left);
0133 
0134   float recHitPos = -1;
0135   if (sl != 2) {
0136     if (fabs(hitPosInChamber_left.x() - segmentPos) < fabs(hitPosInChamber_right.x() - segmentPos))
0137       recHitPos = hitPosInChamber_left.x();
0138     else
0139       recHitPos = hitPosInChamber_right.x();
0140   } else {
0141     if (fabs(hitPosInChamber_left.y() - segmentPos) < fabs(hitPosInChamber_right.y() - segmentPos))
0142       recHitPos = hitPosInChamber_left.y();
0143     else
0144       recHitPos = hitPosInChamber_right.y();
0145   }
0146 
0147   return recHitPos;
0148 }
0149 
0150 // Compute the position (cm) of a hits in a  DTRecHit1D
0151 float DTCalibValidationFromMuons::recHitPosition(
0152     const DTRecHit1D &recHit, const DTLayer *layer, const DTChamber *chamber, float segmentPos, int sl) {
0153   // Get the layer and the wire position
0154   GlobalPoint recHitPosGlob = layer->toGlobal(recHit.localPosition());
0155   LocalPoint recHitPosInChamber = chamber->toLocal(recHitPosGlob);
0156 
0157   float recHitPos = -1;
0158   if (sl != 2)
0159     recHitPos = recHitPosInChamber.x();
0160   else
0161     recHitPos = recHitPosInChamber.y();
0162 
0163   return recHitPos;
0164 }
0165 
0166 // Compute the residuals
0167 void DTCalibValidationFromMuons::compute(const DTGeometry *dtGeom, const DTRecSegment4D &segment) {
0168   bool computeResidual = true;
0169 
0170   // Get all 1D RecHits at step 3 within the 4D segment
0171   vector<DTRecHit1D> recHits1D_S3;
0172 
0173   // Get 1D RecHits at Step 3 and select only events with
0174   // >=7 hits in phi and 4 hits in theta (if any)
0175   const DTChamberRecSegment2D *phiSeg = segment.phiSegment();
0176   if (phiSeg) {
0177     vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
0178     if (phiRecHits.size() < 7) {
0179       LogTrace("DTCalibValidationFromMuons") << "[DTCalibValidationFromMuons] Phi segments has: " << phiRecHits.size()
0180                                              << " hits, skipping";  // FIXME: info output
0181       computeResidual = false;
0182     }
0183     copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
0184   }
0185   if (!phiSeg) {
0186     LogTrace("DTCalibValidationFromMuons") << " [DTCalibValidationFromMuons] 4D segment has no phi segment! ";
0187     computeResidual = false;
0188   }
0189 
0190   if (segment.dimension() == 4) {
0191     const DTSLRecSegment2D *zSeg = segment.zSegment();
0192     if (zSeg) {
0193       vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
0194       if (zRecHits.size() != 4) {
0195         LogTrace("DTCalibValidationFromMuons") << "[DTCalibValidationFromMuons] Theta segments has: " << zRecHits.size()
0196                                                << " hits, skipping";  // FIXME: info output
0197         computeResidual = false;
0198       }
0199       copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
0200     }
0201     if (!zSeg) {
0202       LogTrace("DTCalibValidationFromMuons") << " [DTCalibValidationFromMuons] 4D segment has not the z segment! ";
0203       computeResidual = false;
0204     }
0205   }
0206 
0207   if (!computeResidual)
0208     ++wrongSegment;
0209 
0210   if (computeResidual) {
0211     ++rightSegment;
0212 
0213     // Loop over 1D RecHit inside 4D segment
0214     for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin(); recHit1D != recHits1D_S3.end();
0215          ++recHit1D) {
0216       const DTWireId wireId = (*recHit1D).wireId();
0217 
0218       // Get the layer and the wire position
0219       const DTLayer *layer = dtGeom->layer(wireId);
0220       float wireX = layer->specificTopology().wirePosition(wireId.wire());
0221 
0222       // Extrapolate the segment to the z of the wire
0223       // Get wire position in chamber RF
0224       // (y and z must be those of the hit to be coherent in the transf. of RF
0225       // in case of rotations of the layer alignment)
0226       LocalPoint wirePosInLay(wireX, (*recHit1D).localPosition().y(), (*recHit1D).localPosition().z());
0227       GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
0228       const DTChamber *chamber = dtGeom->chamber((*recHit1D).wireId().layerId().chamberId());
0229       LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
0230 
0231       // Segment position at Wire z in chamber local frame
0232       LocalPoint segPosAtZWire = segment.localPosition() + segment.localDirection() * wirePosInChamber.z() /
0233                                                                cos(segment.localDirection().theta());
0234 
0235       // Compute the distance of the segment from the wire
0236       int sl = wireId.superlayer();
0237       float SegmDistance = -1;
0238       if (sl == 1 || sl == 3) {
0239         // RPhi SL
0240         SegmDistance = fabs(wirePosInChamber.x() - segPosAtZWire.x());
0241         LogTrace("DTCalibValidationFromMuons") << "SegmDistance: " << SegmDistance;
0242       } else if (sl == 2) {
0243         // RZ SL
0244         SegmDistance = fabs(segPosAtZWire.y() - wirePosInChamber.y());
0245         LogTrace("DTCalibValidationFromMuons") << "SegmDistance: " << SegmDistance;
0246       }
0247 
0248       if (SegmDistance > 2.1)
0249         LogTrace("DTCalibValidationFromMuons") << "  Warning: dist segment-wire: " << SegmDistance;
0250 
0251       // Compute the distance of the recHit from the wire
0252       float recHitWireDist = recHitDistFromWire(*recHit1D, layer);
0253       LogTrace("DTCalibValidationFromMuons") << "recHitWireDist: " << recHitWireDist;
0254 
0255       // Compute the residuals
0256       float residualOnDistance = recHitWireDist - SegmDistance;
0257       LogTrace("DTCalibValidationFromMuons") << "WireId: " << wireId << "  ResidualOnDistance: " << residualOnDistance;
0258       float residualOnPosition = -1;
0259       float recHitPos = -1;
0260       if (sl == 1 || sl == 3) {
0261         recHitPos = recHitPosition(*recHit1D, layer, chamber, segPosAtZWire.x(), sl);
0262         residualOnPosition = recHitPos - segPosAtZWire.x();
0263       } else {
0264         recHitPos = recHitPosition(*recHit1D, layer, chamber, segPosAtZWire.y(), sl);
0265         residualOnPosition = recHitPos - segPosAtZWire.y();
0266       }
0267       LogTrace("DTCalibValidationFromMuons") << "WireId: " << wireId << "  ResidualOnPosition: " << residualOnPosition;
0268 
0269       // Fill the histos
0270       if (sl == 1 || sl == 3)
0271         fillHistos(wireId.superlayerId(),
0272                    SegmDistance,
0273                    residualOnDistance,
0274                    (wirePosInChamber.x() - segPosAtZWire.x()),
0275                    residualOnPosition,
0276                    3);
0277       else
0278         fillHistos(wireId.superlayerId(),
0279                    SegmDistance,
0280                    residualOnDistance,
0281                    (wirePosInChamber.y() - segPosAtZWire.y()),
0282                    residualOnPosition,
0283                    3);
0284     }
0285   }
0286 }
0287 
0288 void DTCalibValidationFromMuons::bookHistograms(DQMStore::IBooker &ibooker,
0289                                                 edm::Run const &iRun,
0290                                                 edm::EventSetup const &iSetup) {
0291   // FR substitute the DQMStore instance by ibooker
0292   ibooker.setCurrentFolder("DT/DTCalibValidationFromMuons");
0293 
0294   DTSuperLayerId slId;
0295 
0296   // Loop over all the chambers
0297   vector<const DTChamber *>::const_iterator ch_it = dtGeom->chambers().begin();
0298   vector<const DTChamber *>::const_iterator ch_end = dtGeom->chambers().end();
0299   for (; ch_it != ch_end; ++ch_it) {
0300     vector<const DTSuperLayer *>::const_iterator sl_it = (*ch_it)->superLayers().begin();
0301     vector<const DTSuperLayer *>::const_iterator sl_end = (*ch_it)->superLayers().end();
0302     // Loop over the SLs
0303     for (; sl_it != sl_end; ++sl_it) {
0304       slId = (*sl_it)->id();
0305 
0306       // TODO! fix this is a leftover
0307       int firstStep = 3;
0308       // Loop over the 3 steps
0309       for (int step = firstStep; step <= 3; ++step) {
0310         LogTrace("DTCalibValidationFromMuons") << "   Booking histos for SL: " << slId;
0311 
0312         // Compose the chamber name
0313         stringstream wheel;
0314         wheel << slId.wheel();
0315         stringstream station;
0316         station << slId.station();
0317         stringstream sector;
0318         sector << slId.sector();
0319         stringstream superLayer;
0320         superLayer << slId.superlayer();
0321         // Define the step
0322         stringstream Step;
0323         Step << step;
0324 
0325         string slHistoName = "_STEP" + Step.str() + "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str() +
0326                              "_SL" + superLayer.str();
0327 
0328         ibooker.setCurrentFolder("DT/DTCalibValidationFromMuons/Wheel" + wheel.str() + "/Station" + station.str() +
0329                                  "/Sector" + sector.str());
0330         // Create the monitor elements
0331         vector<MonitorElement *> histos;
0332         // Note the order matters
0333         histos.push_back(ibooker.book1D(
0334             "hResDist" + slHistoName, "Residuals on the distance from wire (rec_hit - segm_extr) (cm)", 200, -0.4, 0.4));
0335         histos.push_back(ibooker.book2D("hResDistVsDist" + slHistoName,
0336                                         "Residuals on the distance (cm) from wire (rec_hit "
0337                                         "- segm_extr) vs distance  (cm)",
0338                                         100,
0339                                         0,
0340                                         2.5,
0341                                         200,
0342                                         -0.4,
0343                                         0.4));
0344 
0345         histosPerSL[make_pair(slId, step)] = histos;
0346       }
0347     }
0348   }
0349 }
0350 
0351 // Fill a set of histograms for a given SL
0352 void DTCalibValidationFromMuons::fillHistos(
0353     DTSuperLayerId slId, float distance, float residualOnDistance, float position, float residualOnPosition, int step) {
0354   // FIXME: optimization of the number of searches
0355   vector<MonitorElement *> histos = histosPerSL[make_pair(slId, step)];
0356   histos[0]->Fill(residualOnDistance);
0357   histos[1]->Fill(distance, residualOnDistance);
0358 }
0359 
0360 // Local Variables:
0361 // show-trailing-whitespace: t
0362 // truncate-lines: t
0363 // End:
[ 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!