Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-19 10:28:03

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