DTMonitorModule/src/DTCalibValidation.cc

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/DTCalibValidation.h"
0009
0010 // Framework
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "DQMServices/Core/interface/DQMStore.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 //Geometry
0019 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0020
0021 //RecHit
0022 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0023 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0024
0025 #include <iterator>
0026
0027 using namespace edm;
0028 using namespace std;
0029
0030 DTCalibValidation::DTCalibValidation(const ParameterSet& pset)
0031     : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
0032   parameters = pset;
0033
0034   //FR the following was previously in the beginJob
0035
0036   // the name of the rechits collection at step 1
0037   recHits1DToken_ =
0038       consumes<DTRecHitCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("recHits1DLabel")));
0039   // the name of the 2D segments
0040   segment2DToken_ =
0041       consumes<DTRecSegment2DCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("segment2DLabel")));
0042   // the name of the 4D segments
0043   segment4DToken_ =
0044       consumes<DTRecSegment4DCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("segment4DLabel")));
0045   // the counter of segments not used to compute residuals
0046   wrongSegment = 0;
0047   // the counter of segments used to compute residuals
0048   rightSegment = 0;
0049   // the analysis type
0050   detailedAnalysis = parameters.getUntrackedParameter<bool>("detailedAnalysis", false);
0051
0052   nevent = 0;
0053 }
0054
0055 DTCalibValidation::~DTCalibValidation() {
0056   //FR the following was previously in the endJob
0057
0058   LogVerbatim("DTCalibValidation") << "Segments used to compute residuals: " << rightSegment;
0059   LogVerbatim("DTCalibValidation") << "Segments not used to compute residuals: " << wrongSegment;
0060 }
0061
0062 void DTCalibValidation::dqmBeginRun(const edm::Run& run, const edm::EventSetup& setup) {
0063   // get the geometry
0064   dtGeom = &setup.getData(muonGeomToken_);
0065 }
0066
0067 void DTCalibValidation::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0068   ++nevent;
0069   LogTrace("DTCalibValidation") << "[DTCalibValidation] Analyze #Run: " << event.id().run() << " #Event: " << nevent;
0070
0071   // RecHit mapping at Step 1 -------------------------------
0072   map<DTWireId, vector<DTRecHit1DPair> > recHitsPerWire_1S;
0073
0074   // RecHit mapping at Step 2 ------------------------------
0075   map<DTWireId, vector<DTRecHit1D> > recHitsPerWire_2S;
0076
0077   if (detailedAnalysis) {
0078     LogTrace("DTCalibValidation") << "  -- DTRecHit S1: begin analysis:";
0079     // Get the rechit collection from the event
0080     Handle<DTRecHitCollection> dtRecHits;
0081     event.getByToken(recHits1DToken_, dtRecHits);
0082     recHitsPerWire_1S = map1DRecHitsPerWire(dtRecHits.product());
0083
0084     LogTrace("DTCalibValidation") << "  -- DTRecHit S2: begin analysis:";
0085     // Get the 2D rechits from the event
0086     Handle<DTRecSegment2DCollection> segment2Ds;
0087     event.getByToken(segment2DToken_, segment2Ds);
0088     recHitsPerWire_2S = map1DRecHitsPerWire(segment2Ds.product());
0089   }
0090
0091   // RecHit mapping at Step 3 ---------------------------------
0092   LogTrace("DTCalibValidation") << "  -- DTRecHit S3: begin analysis:";
0093   // Get the 4D rechits from the event
0094   Handle<DTRecSegment4DCollection> segment4Ds;
0095   event.getByToken(segment4DToken_, segment4Ds);
0096   map<DTWireId, vector<DTRecHit1D> > recHitsPerWire_3S = map1DRecHitsPerWire(segment4Ds.product());
0097
0098   // Loop over all 4D segments
0099   for (DTRecSegment4DCollection::const_iterator segment = segment4Ds->begin(); segment != segment4Ds->end();
0100        ++segment) {
0101     if (detailedAnalysis) {
0102       LogTrace("DTCalibValidation") << "Anlysis on recHit at step 1";
0103       compute(dtGeom, (*segment), recHitsPerWire_1S, 1);
0104
0105       LogTrace("DTCalibValidation") << "Anlysis on recHit at step 2";
0106       compute(dtGeom, (*segment), recHitsPerWire_2S, 2);
0107     }
0108
0109     LogTrace("DTCalibValidation") << "Anlysis on recHit at step 3";
0110     compute(dtGeom, (*segment), recHitsPerWire_3S, 3);
0111   }
0112 }
0113
0114 // Return a map between DTRecHit1DPair and wireId
0115 map<DTWireId, vector<DTRecHit1DPair> > DTCalibValidation::map1DRecHitsPerWire(
0116     const DTRecHitCollection* dt1DRecHitPairs) {
0117   map<DTWireId, vector<DTRecHit1DPair> > ret;
0118
0119   for (DTRecHitCollection::const_iterator rechit = dt1DRecHitPairs->begin(); rechit != dt1DRecHitPairs->end();
0120        ++rechit) {
0121     ret[(*rechit).wireId()].push_back(*rechit);
0122   }
0123
0124   return ret;
0125 }
0126
0127 // Return a map between DTRecHit1D at S2 and wireId
0128 map<DTWireId, vector<DTRecHit1D> > DTCalibValidation::map1DRecHitsPerWire(const DTRecSegment2DCollection* segment2Ds) {
0129   map<DTWireId, vector<DTRecHit1D> > ret;
0130
0131   // Loop over all 2D segments
0132   for (DTRecSegment2DCollection::const_iterator segment = segment2Ds->begin(); segment != segment2Ds->end();
0133        ++segment) {
0134     vector<DTRecHit1D> component1DHits = (*segment).specificRecHits();
0135     // Loop over all component 1D hits
0136     for (vector<DTRecHit1D>::const_iterator hit = component1DHits.begin(); hit != component1DHits.end(); ++hit) {
0137       ret[(*hit).wireId()].push_back(*hit);
0138     }
0139   }
0140   return ret;
0141 }
0142
0143 // Return a map between DTRecHit1D at S3 and wireId
0144 map<DTWireId, std::vector<DTRecHit1D> > DTCalibValidation::map1DRecHitsPerWire(
0145     const DTRecSegment4DCollection* segment4Ds) {
0146   map<DTWireId, vector<DTRecHit1D> > ret;
0147   // Loop over all 4D segments
0148   for (DTRecSegment4DCollection::const_iterator segment = segment4Ds->begin(); segment != segment4Ds->end();
0149        ++segment) {
0150     // Get component 2D segments
0151     vector<const TrackingRecHit*> segment2Ds = (*segment).recHits();
0152     // Loop over 2D segments:
0153     for (vector<const TrackingRecHit*>::const_iterator segment2D = segment2Ds.begin(); segment2D != segment2Ds.end();
0154          ++segment2D) {
0155       // Get 1D component rechits
0156       vector<const TrackingRecHit*> hits = (*segment2D)->recHits();
0157       // Loop over them
0158       for (vector<const TrackingRecHit*>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0159         const DTRecHit1D* hit1D = dynamic_cast<const DTRecHit1D*>(*hit);
0160         ret[hit1D->wireId()].push_back(*hit1D);
0161       }
0162     }
0163   }
0164
0165   return ret;
0166 }
0167
0168 // Find the RecHit closest to the segment4D
0169 template <typename type>
0170 const type* DTCalibValidation::findBestRecHit(const DTLayer* layer,
0171                                               DTWireId wireId,
0172                                               const vector<type>& recHits,
0173                                               const float segmDist) {
0174   float res = 99999;
0175   const type* theBestRecHit = nullptr;
0176   // Loop over RecHits within the cell
0177   for (typename vector<type>::const_iterator recHit = recHits.begin(); recHit != recHits.end(); ++recHit) {
0178     float distTmp = recHitDistFromWire(*recHit, layer);
0179     if (fabs(distTmp - segmDist) < res) {
0180       res = fabs(distTmp - segmDist);
0181       theBestRecHit = &(*recHit);
0182     }
0183   }  // End of loop over RecHits within the cell
0184
0185   return theBestRecHit;
0186 }
0187
0188 // Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
0189 float DTCalibValidation::recHitDistFromWire(const DTRecHit1DPair& hitPair, const DTLayer* layer) {
0190   return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
0191 }
0192
0193 // Compute the distance from wire (cm) of a hits in a DTRecHit1D
0194 float DTCalibValidation::recHitDistFromWire(const DTRecHit1D& recHit, const DTLayer* layer) {
0195   return fabs(recHit.localPosition().x() - layer->specificTopology().wirePosition(recHit.wireId().wire()));
0196 }
0197
0198 // Compute the position (cm) of a hits in a DTRecHit1DPair
0199 float DTCalibValidation::recHitPosition(
0200     const DTRecHit1DPair& hitPair, const DTLayer* layer, const DTChamber* chamber, float segmentPos, int sl) {
0201   // Get the layer and the wire position
0202   GlobalPoint hitPosGlob_right = layer->toGlobal(hitPair.localPosition(DTEnums::Right));
0203   LocalPoint hitPosInChamber_right = chamber->toLocal(hitPosGlob_right);
0204   GlobalPoint hitPosGlob_left = layer->toGlobal(hitPair.localPosition(DTEnums::Left));
0205   LocalPoint hitPosInChamber_left = chamber->toLocal(hitPosGlob_left);
0206
0207   float recHitPos = -1;
0208   if (sl != 2) {
0209     if (fabs(hitPosInChamber_left.x() - segmentPos) < fabs(hitPosInChamber_right.x() - segmentPos))
0210       recHitPos = hitPosInChamber_left.x();
0211     else
0212       recHitPos = hitPosInChamber_right.x();
0213   } else {
0214     if (fabs(hitPosInChamber_left.y() - segmentPos) < fabs(hitPosInChamber_right.y() - segmentPos))
0215       recHitPos = hitPosInChamber_left.y();
0216     else
0217       recHitPos = hitPosInChamber_right.y();
0218   }
0219
0220   return recHitPos;
0221 }
0222
0223 // Compute the position (cm) of a hits in a  DTRecHit1D
0224 float DTCalibValidation::recHitPosition(
0225     const DTRecHit1D& recHit, const DTLayer* layer, const DTChamber* chamber, float segmentPos, int sl) {
0226   // Get the layer and the wire position
0227   GlobalPoint recHitPosGlob = layer->toGlobal(recHit.localPosition());
0228   LocalPoint recHitPosInChamber = chamber->toLocal(recHitPosGlob);
0229
0230   float recHitPos = -1;
0231   if (sl != 2)
0232     recHitPos = recHitPosInChamber.x();
0233   else
0234     recHitPos = recHitPosInChamber.y();
0235
0236   return recHitPos;
0237 }
0238
0239 // Compute the residuals
0240 template <typename type>
0241 void DTCalibValidation::compute(const DTGeometry* dtGeom,
0242                                 const DTRecSegment4D& segment,
0243                                 const std::map<DTWireId, std::vector<type> >& recHitsPerWire,
0244                                 int step) {
0245   bool computeResidual = true;
0246
0247   // Get all 1D RecHits at step 3 within the 4D segment
0248   vector<DTRecHit1D> recHits1D_S3;
0249
0250   // Get 1D RecHits at Step 3 and select only events with
0251   // 8 hits in phi and 4 hits in theta (if any)
0252   const DTChamberRecSegment2D* phiSeg = segment.phiSegment();
0253   if (phiSeg) {
0254     vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
0255     if (phiRecHits.size() != 8) {
0256       LogTrace("DTCalibValidation") << "[DTCalibValidation] Phi segments has: " << phiRecHits.size()
0257                                     << " hits, skipping";  // FIXME: info output
0258       computeResidual = false;
0259     }
0260     copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
0261   }
0262   if (!phiSeg) {
0263     LogTrace("DTCalibValidation") << " [DTCalibValidation] 4D segment has not the phi segment! ";
0264     computeResidual = false;
0265   }
0266
0267   if (segment.dimension() == 4) {
0268     const DTSLRecSegment2D* zSeg = segment.zSegment();
0269     if (zSeg) {
0270       vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
0271       if (zRecHits.size() != 4) {
0272         LogTrace("DTCalibValidation") << "[DTCalibValidation] Theta segments has: " << zRecHits.size()
0273                                       << " hits, skipping";  // FIXME: info output
0274         computeResidual = false;
0275       }
0276       copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
0277     }
0278     if (!zSeg) {
0279       LogTrace("DTCalibValidation") << " [DTCalibValidation] 4D segment has not the z segment! ";
0280       computeResidual = false;
0281     }
0282   }
0283
0284   if (!computeResidual)
0285     ++wrongSegment;
0286   if (computeResidual) {
0287     ++rightSegment;
0288     // Loop over 1D RecHit inside 4D segment
0289     for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin(); recHit1D != recHits1D_S3.end();
0290          ++recHit1D) {
0291       const DTWireId wireId = (*recHit1D).wireId();
0292
0293       // Get the layer and the wire position
0294       const DTLayer* layer = dtGeom->layer(wireId);
0295       float wireX = layer->specificTopology().wirePosition(wireId.wire());
0296
0297       // Extrapolate the segment to the z of the wire
0298       // Get wire position in chamber RF
0299       // (y and z must be those of the hit to be coherent in the transf. of RF in case of rotations of the layer alignment)
0300       LocalPoint wirePosInLay(wireX, (*recHit1D).localPosition().y(), (*recHit1D).localPosition().z());
0301       GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
0302       const DTChamber* chamber = dtGeom->chamber((*recHit1D).wireId().layerId().chamberId());
0303       LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
0304
0305       // Segment position at Wire z in chamber local frame
0306       LocalPoint segPosAtZWire = segment.localPosition() + segment.localDirection() * wirePosInChamber.z() /
0307                                                                cos(segment.localDirection().theta());
0308
0309       // Compute the distance of the segment from the wire
0310       int sl = wireId.superlayer();
0311       float SegmDistance = -1;
0312       if (sl == 1 || sl == 3) {
0313         // RPhi SL
0314         SegmDistance = fabs(wirePosInChamber.x() - segPosAtZWire.x());
0315         LogTrace("DTCalibValidation") << "SegmDistance: " << SegmDistance;
0316       } else if (sl == 2) {
0317         // RZ SL
0318         SegmDistance = fabs(segPosAtZWire.y() - wirePosInChamber.y());
0319         LogTrace("DTCalibValidation") << "SegmDistance: " << SegmDistance;
0320       }
0321       if (SegmDistance > 2.1)
0322         LogTrace("DTCalibValidation") << "  Warning: dist segment-wire: " << SegmDistance;
0323
0324       // Look for RecHits in the same cell
0325       if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
0326         LogTrace("DTCalibValidation") << "   No RecHit found at Step: " << step << " in cell: " << wireId;
0327       } else {
0328         const vector<type>& recHits = recHitsPerWire.at(wireId);
0329         LogTrace("DTCalibValidation") << "   " << recHits.size() << " RecHits, Step " << step
0330                                       << " in channel: " << wireId;
0331
0332         // Get the layer
0333         const DTLayer* layer = dtGeom->layer(wireId);
0334         // Find the best RecHits
0335         const type* theBestRecHit = findBestRecHit(layer, wireId, recHits, SegmDistance);
0336         // Compute the distance of the recHit from the wire
0337         float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
0338         LogTrace("DTCalibValidation") << "recHitWireDist: " << recHitWireDist;
0339
0340         // Compute the residuals
0341         float residualOnDistance = recHitWireDist - SegmDistance;
0342         LogTrace("DTCalibValidation") << "WireId: " << wireId << "  ResidualOnDistance: " << residualOnDistance;
0343         float residualOnPosition = -1;
0344         float recHitPos = -1;
0345         if (sl == 1 || sl == 3) {
0346           recHitPos = recHitPosition(*theBestRecHit, layer, chamber, segPosAtZWire.x(), sl);
0347           residualOnPosition = recHitPos - segPosAtZWire.x();
0348         } else {
0349           recHitPos = recHitPosition(*theBestRecHit, layer, chamber, segPosAtZWire.y(), sl);
0350           residualOnPosition = recHitPos - segPosAtZWire.y();
0351         }
0352         LogTrace("DTCalibValidation") << "WireId: " << wireId << "  ResidualOnPosition: " << residualOnPosition;
0353
0354         // Fill the histos
0355         if (sl == 1 || sl == 3)
0356           fillHistos(wireId.superlayerId(),
0357                      SegmDistance,
0358                      residualOnDistance,
0359                      (wirePosInChamber.x() - segPosAtZWire.x()),
0360                      residualOnPosition,
0361                      step);
0362         else
0363           fillHistos(wireId.superlayerId(),
0364                      SegmDistance,
0365                      residualOnDistance,
0366                      (wirePosInChamber.y() - segPosAtZWire.y()),
0367                      residualOnPosition,
0368                      step);
0369       }
0370     }
0371   }
0372 }
0373
0374 void DTCalibValidation::bookHistograms(DQMStore::IBooker& ibooker,
0375                                        edm::Run const& iRun,
0376                                        edm::EventSetup const& iSetup) {
0377   //FR substitute the DQMStore instance by ibooker
0378   ibooker.setCurrentFolder("DT/DTCalibValidation");
0379
0380   DTSuperLayerId slId;
0381
0382   // Loop over all the chambers
0383   vector<const DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
0384   vector<const DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
0385   for (; ch_it != ch_end; ++ch_it) {
0386     vector<const DTSuperLayer*>::const_iterator sl_it = (*ch_it)->superLayers().begin();
0387     vector<const DTSuperLayer*>::const_iterator sl_end = (*ch_it)->superLayers().end();
0388     // Loop over the SLs
0389     for (; sl_it != sl_end; ++sl_it) {
0390       slId = (*sl_it)->id();
0391
0392       int firstStep = 1;
0393       if (!detailedAnalysis)
0394         firstStep = 3;
0395       // Loop over the 3 steps
0396       for (int step = firstStep; step <= 3; ++step) {
0397         LogTrace("DTCalibValidation") << "   Booking histos for SL: " << slId;
0398
0399         // Compose the chamber name
0400         stringstream wheel;
0401         wheel << slId.wheel();
0402         stringstream station;
0403         station << slId.station();
0404         stringstream sector;
0405         sector << slId.sector();
0406         stringstream superLayer;
0407         superLayer << slId.superlayer();
0408         // Define the step
0409         stringstream Step;
0410         Step << step;
0411
0412         string slHistoName = "_STEP" + Step.str() + "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str() +
0413                              "_SL" + superLayer.str();
0414
0415         ibooker.setCurrentFolder("DT/DTCalibValidation/Wheel" + wheel.str() + "/Station" + station.str() + "/Sector" +
0416                                  sector.str());
0417         // Create the monitor elements
0418         vector<MonitorElement*> histos;
0419         // Note the order matters
0420         histos.push_back(ibooker.book1D(
0421             "hResDist" + slHistoName, "Residuals on the distance from wire (rec_hit - segm_extr) (cm)", 200, -0.4, 0.4));
0422         histos.push_back(
0423             ibooker.book2D("hResDistVsDist" + slHistoName,
0424                            "Residuals on the distance (cm) from wire (rec_hit - segm_extr) vs distance  (cm)",
0425                            100,
0426                            0,
0427                            2.5,
0428                            200,
0429                            -0.4,
0430                            0.4));
0431         if (detailedAnalysis) {
0432           histos.push_back(ibooker.book1D("hResPos" + slHistoName,
0433                                           "Residuals on the position from wire (rec_hit - segm_extr) (cm)",
0434                                           200,
0435                                           -0.4,
0436                                           0.4));
0437           histos.push_back(
0438               ibooker.book2D("hResPosVsPos" + slHistoName,
0439                              "Residuals on the position (cm) from wire (rec_hit - segm_extr) vs distance  (cm)",
0440                              200,
0441                              -2.5,
0442                              2.5,
0443                              200,
0444                              -0.4,
0445                              0.4));
0446         }
0447
0448         histosPerSL[make_pair(slId, step)] = histos;
0449       }
0450     }
0451   }
0452 }
0453
0454 // Fill a set of histograms for a given SL
0455 void DTCalibValidation::fillHistos(
0456     DTSuperLayerId slId, float distance, float residualOnDistance, float position, float residualOnPosition, int step) {
0457   // FIXME: optimization of the number of searches
0458   vector<MonitorElement*> histos = histosPerSL[make_pair(slId, step)];
0459   histos[0]->Fill(residualOnDistance);
0460   histos[1]->Fill(distance, residualOnDistance);
0461   if (detailedAnalysis) {
0462     histos[2]->Fill(residualOnPosition);
0463     histos[3]->Fill(position, residualOnPosition);
0464   }
0465 }
0466
0467 // Local Variables:
0468 // show-trailing-whitespace: t
0469 // truncate-lines: t
0470 // End: