DTCalibration/plugins/DTResidualCalibration.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  */
0006
0007 #include "DTResidualCalibration.h"
0008
0009 // Framework
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 //Geometry
0017 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0018 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0019
0020 //RecHit
0021 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0022 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0023
0024 #include "CommonTools/Utils/interface/TH1AddDirectorySentry.h"
0025 #include "CalibMuon/DTCalibration/interface/DTSegmentSelector.h"
0026 #include "CalibMuon/DTCalibration/interface/DTRecHitSegmentResidual.h"
0027
0028 #include "TFile.h"
0029 #include "TH1F.h"
0030 #include "TH2F.h"
0031
0032 #include <algorithm>
0033
0034 DTResidualCalibration::DTResidualCalibration(const edm::ParameterSet& pset)
0035     : histRange_(pset.getParameter<double>("histogramRange")),
0036       segment4DToken_(consumes<DTRecSegment4DCollection>(pset.getParameter<edm::InputTag>("segment4DLabel"))),
0037       rootBaseDir_(pset.getUntrackedParameter<std::string>("rootBaseDir", "DT/Residuals")),
0038       detailedAnalysis_(pset.getUntrackedParameter<bool>("detailedAnalysis", false)),
0039       dtGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
0040   edm::ConsumesCollector collector(consumesCollector());
0041   select_ = new DTSegmentSelector(pset, collector);
0042
0043   LogDebug("Calibration") << "[DTResidualCalibration] Constructor called.";
0044   std::string rootFileName = pset.getUntrackedParameter<std::string>("rootFileName", "residuals.root");
0045   rootFile_ = new TFile(rootFileName.c_str(), "RECREATE");
0046   rootFile_->cd();
0047
0048   segmok = 0;
0049   segmbad = 0;
0050   nevent = 0;
0051 }
0052
0053 DTResidualCalibration::~DTResidualCalibration() {
0054   delete select_;
0055   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Destructor called.";
0056   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Analyzed events: " << nevent;
0057   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Good segments: " << segmok;
0058   edm::LogVerbatim("Calibration") << "[DTResidualCalibration] Bad segments: " << segmbad;
0059 }
0060
0061 void DTResidualCalibration::beginJob() { TH1::SetDefaultSumw2(true); }
0062
0063 void DTResidualCalibration::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
0064   // get the geometry
0065   edm::ESHandle<DTGeometry> dtGeomH;
0066   dtGeomH = setup.getHandle(dtGeomToken_);
0067   dtGeom_ = dtGeomH.product();
0068
0069   // Loop over all the chambers
0070   if (histoMapTH1F_.empty()) {
0071     for (auto ch_it : dtGeom_->chambers()) {
0072       // Loop over the SLs
0073       for (auto sl_it : ch_it->superLayers()) {
0074         DTSuperLayerId slId = (sl_it)->id();
0075         bookHistos(slId);
0076         if (detailedAnalysis_) {
0077           for (auto layer_it : (sl_it)->layers()) {
0078             DTLayerId layerId = (layer_it)->id();
0079             bookHistos(layerId);
0080           }
0081         }
0082       }
0083     }
0084   }
0085 }
0086
0087 void DTResidualCalibration::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0088   rootFile_->cd();
0089   ++nevent;
0090
0091   // Get the 4D rechits from the event
0092   const edm::Handle<DTRecSegment4DCollection>& segments4D = event.getHandle(segment4DToken_);
0093
0094   // Loop over segments by chamber
0095   DTRecSegment4DCollection::id_iterator chamberIdIt;
0096   for (chamberIdIt = segments4D->id_begin(); chamberIdIt != segments4D->id_end(); ++chamberIdIt) {
0097     const DTChamber* chamber = dtGeom_->chamber(*chamberIdIt);
0098
0099     // Get the range for the corresponding ChamberId
0100     DTRecSegment4DCollection::range range = segments4D->get((*chamberIdIt));
0101     // Loop over the rechits of this DetUnit
0102     for (DTRecSegment4DCollection::const_iterator segment = range.first; segment != range.second; ++segment) {
0103       LogTrace("Calibration") << "Segment local pos (in chamber RF): " << (*segment).localPosition()
0104                               << "\nSegment global pos: " << chamber->toGlobal((*segment).localPosition());
0105
0106       if (!(*select_)(*segment, event, setup)) {
0107         segmbad++;
0108         continue;
0109       }
0110       segmok++;
0111
0112       // Get all 1D RecHits at step 3 within the 4D segment
0113       std::vector<DTRecHit1D> recHits1D_S3;
0114
0115       if ((*segment).hasPhi()) {
0116         const DTChamberRecSegment2D* phiSeg = (*segment).phiSegment();
0117         const std::vector<DTRecHit1D>& phiRecHits = phiSeg->specificRecHits();
0118         std::copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
0119       }
0120
0121       if ((*segment).hasZed()) {
0122         const DTSLRecSegment2D* zSeg = (*segment).zSegment();
0123         const std::vector<DTRecHit1D>& zRecHits = zSeg->specificRecHits();
0124         std::copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
0125       }
0126
0127       // Loop over 1D RecHit inside 4D segment
0128       for (std::vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin(); recHit1D != recHits1D_S3.end();
0129            ++recHit1D) {
0130         const DTWireId wireId = recHit1D->wireId();
0131
0132         float segmDistance = segmentToWireDistance(*recHit1D, *segment);
0133         if (segmDistance > 2.1)
0134           LogTrace("Calibration") << "WARNING: segment-wire distance: " << segmDistance;
0135         else
0136           LogTrace("Calibration") << "segment-wire distance: " << segmDistance;
0137
0138         float residualOnDistance = DTRecHitSegmentResidual().compute(dtGeom_, *recHit1D, *segment);
0139         LogTrace("Calibration") << "Wire Id " << wireId << " residual on distance: " << residualOnDistance;
0140
0141         fillHistos(wireId.superlayerId(), segmDistance, residualOnDistance);
0142         if (detailedAnalysis_)
0143           fillHistos(wireId.layerId(), segmDistance, residualOnDistance);
0144       }
0145     }
0146   }
0147 }
0148
0149 float DTResidualCalibration::segmentToWireDistance(const DTRecHit1D& recHit1D, const DTRecSegment4D& segment) {
0150   // Get the layer and the wire position
0151   const DTWireId wireId = recHit1D.wireId();
0152   const DTLayer* layer = dtGeom_->layer(wireId);
0153   float wireX = layer->specificTopology().wirePosition(wireId.wire());
0154
0155   // Extrapolate the segment to the z of the wire
0156   // Get wire position in chamber RF
0157   // (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)
0158   LocalPoint wirePosInLay(wireX, recHit1D.localPosition().y(), recHit1D.localPosition().z());
0159   GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
0160   const DTChamber* chamber = dtGeom_->chamber(wireId.layerId().chamberId());
0161   LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
0162
0163   // Segment position at Wire z in chamber local frame
0164   LocalPoint segPosAtZWire =
0165       segment.localPosition() + segment.localDirection() * wirePosInChamber.z() / cos(segment.localDirection().theta());
0166
0167   // Compute the distance of the segment from the wire
0168   int sl = wireId.superlayer();
0169   float segmDistance = -1;
0170   if (sl == 1 || sl == 3)
0171     segmDistance = fabs(wirePosInChamber.x() - segPosAtZWire.x());
0172   else if (sl == 2)
0173     segmDistance = fabs(segPosAtZWire.y() - wirePosInChamber.y());
0174
0175   return segmDistance;
0176 }
0177
0178 void DTResidualCalibration::endJob() {
0179   LogDebug("Calibration") << "[DTResidualCalibration] Writing histos to file.";
0180   rootFile_->cd();
0181   rootFile_->Write();
0182   rootFile_->Close();
0183
0184   /*std::map<DTSuperLayerId, TH1F* >::const_iterator itSlHistos = histoMapTH1F_.begin();
0185   std::map<DTSuperLayerId, TH1F* >::const_iterator itSlHistos_end = histoMapTH1F_.end();
0186   for(; itSlHistos != itSlHistos_end; ++itSlHistos){
0187      std::vector<TH1F*>::const_iterator itHistTH1F = (*itSlHistos).second.begin();
0188      std::vector<TH1F*>::const_iterator itHistTH1F_end = (*itSlHistos).second.end();
0189      for(; itHistTH1F != itHistTH1F_end; ++itHistTH1F) (*itHistTH1F)->Write();
0190
0191      std::vector<TH2F*>::const_iterator itHistTH2F = histoMapTH2F_[(*itSlHistos).first].begin();
0192      std::vector<TH2F*>::const_iterator itHistTH2F_end = histoMapTH2F_[(*itSlHistos).first].end();
0193      for(; itHistTH2F != itHistTH2F_end; ++itHistTH2F) (*itHistTH2F)->Write();
0194   }*/
0195 }
0196
0197 void DTResidualCalibration::bookHistos(DTSuperLayerId slId) {
0198   TH1AddDirectorySentry addDir;
0199   rootFile_->cd();
0200
0201   LogDebug("Calibration") << "[DTResidualCalibration] Booking histos for SL: " << slId;
0202
0203   // Compose the chamber name
0204   // Define the step
0205   int step = 3;
0206
0207   std::string wheelStr = std::to_string(slId.wheel());
0208   std::string stationStr = std::to_string(slId.station());
0209   std::string sectorStr = std::to_string(slId.sector());
0210
0211   std::string slHistoName = "_STEP" + std::to_string(step) + "_W" + wheelStr + "_St" + stationStr + "_Sec" + sectorStr +
0212                             "_SL" + std::to_string(slId.superlayer());
0213
0214   LogDebug("Calibration") << "Accessing " << rootBaseDir_;
0215   TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
0216   if (!baseDir)
0217     baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
0218   LogDebug("Calibration") << "Accessing " << ("Wheel" + wheelStr);
0219   TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr).c_str());
0220   if (!wheelDir)
0221     wheelDir = baseDir->mkdir(("Wheel" + wheelStr).c_str());
0222   LogDebug("Calibration") << "Accessing " << ("Station" + stationStr);
0223   TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr).c_str());
0224   if (!stationDir)
0225     stationDir = wheelDir->mkdir(("Station" + stationStr).c_str());
0226   LogDebug("Calibration") << "Accessing " << ("Sector" + sectorStr);
0227   TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr).c_str());
0228   if (!sectorDir)
0229     sectorDir = stationDir->mkdir(("Sector" + sectorStr).c_str());
0230
0231   sectorDir->cd();
0232
0233   // Create the monitor elements
0234   TH1F* histosTH1F = new TH1F(("hResDist" + slHistoName).c_str(),
0235                               "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
0236                               200,
0237                               -histRange_,
0238                               histRange_);
0239   TH2F* histosTH2F = new TH2F(("hResDistVsDist" + slHistoName).c_str(),
0240                               "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
0241                               100,
0242                               0,
0243                               2.5,
0244                               200,
0245                               -histRange_,
0246                               histRange_);
0247   histoMapTH1F_[slId] = histosTH1F;
0248   histoMapTH2F_[slId] = histosTH2F;
0249 }
0250
0251 void DTResidualCalibration::bookHistos(DTLayerId layerId) {
0252   TH1AddDirectorySentry addDir;
0253   rootFile_->cd();
0254
0255   LogDebug("Calibration") << "[DTResidualCalibration] Booking histos for layer: " << layerId;
0256
0257   // Compose the chamber name
0258   std::string wheelStr = std::to_string(layerId.wheel());
0259   std::string stationStr = std::to_string(layerId.station());
0260   std::string sectorStr = std::to_string(layerId.sector());
0261   std::string superLayerStr = std::to_string(layerId.superlayer());
0262   std::string layerStr = std::to_string(layerId.layer());
0263   // Define the step
0264   int step = 3;
0265
0266   std::string layerHistoName = "_STEP" + std::to_string(step) + "_W" + wheelStr + "_St" + stationStr + "_Sec" +
0267                                sectorStr + "_SL" + superLayerStr + "_Layer" + layerStr;
0268
0269   LogDebug("Calibration") << "Accessing " << rootBaseDir_;
0270   TDirectory* baseDir = rootFile_->GetDirectory(rootBaseDir_.c_str());
0271   if (!baseDir)
0272     baseDir = rootFile_->mkdir(rootBaseDir_.c_str());
0273   LogDebug("Calibration") << "Accessing " << ("Wheel" + wheelStr);
0274   TDirectory* wheelDir = baseDir->GetDirectory(("Wheel" + wheelStr).c_str());
0275   if (!wheelDir)
0276     wheelDir = baseDir->mkdir(("Wheel" + wheelStr).c_str());
0277   LogDebug("Calibration") << "Accessing " << ("Station" + stationStr);
0278   TDirectory* stationDir = wheelDir->GetDirectory(("Station" + stationStr).c_str());
0279   if (!stationDir)
0280     stationDir = wheelDir->mkdir(("Station" + stationStr).c_str());
0281   LogDebug("Calibration") << "Accessing " << ("Sector" + sectorStr);
0282   TDirectory* sectorDir = stationDir->GetDirectory(("Sector" + sectorStr).c_str());
0283   if (!sectorDir)
0284     sectorDir = stationDir->mkdir(("Sector" + sectorStr).c_str());
0285   LogDebug("Calibration") << "Accessing " << ("SL" + superLayerStr);
0286   TDirectory* superLayerDir = sectorDir->GetDirectory(("SL" + superLayerStr).c_str());
0287   if (!superLayerDir)
0288     superLayerDir = sectorDir->mkdir(("SL" + superLayerStr).c_str());
0289
0290   superLayerDir->cd();
0291   // Create histograms
0292   TH1F* histosTH1F = new TH1F(("hResDist" + layerHistoName).c_str(),
0293                               "Residuals on the distance from wire (rec_hit - segm_extr) (cm)",
0294                               200,
0295                               -histRange_,
0296                               histRange_);
0297   TH2F* histosTH2F = new TH2F(("hResDistVsDist" + layerHistoName).c_str(),
0298                               "Residuals on the dist. (cm) from wire (rec_hit - segm_extr) vs dist. (cm)",
0299                               100,
0300                               0,
0301                               2.5,
0302                               200,
0303                               -histRange_,
0304                               histRange_);
0305   histoMapPerLayerTH1F_[layerId] = histosTH1F;
0306   histoMapPerLayerTH2F_[layerId] = histosTH2F;
0307 }
0308
0309 // Fill a set of histograms for a given SL
0310 void DTResidualCalibration::fillHistos(DTSuperLayerId slId, float distance, float residualOnDistance) {
0311   histoMapTH1F_[slId]->Fill(residualOnDistance);
0312   histoMapTH2F_[slId]->Fill(distance, residualOnDistance);
0313 }
0314
0315 // Fill a set of histograms for a given layer
0316 void DTResidualCalibration::fillHistos(DTLayerId layerId, float distance, float residualOnDistance) {
0317   histoMapPerLayerTH1F_[layerId]->Fill(residualOnDistance);
0318   histoMapPerLayerTH2F_[layerId]->Fill(distance, residualOnDistance);
0319 }