File indexing completed on 2024-09-11 04:32:28
0001
0002
0003
0004
0005
0006
0007
0008 #include "DTResolutionAnalysisTask.h"
0009
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "DQMServices/Core/interface/DQMStore.h"
0020
0021
0022 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0023
0024
0025 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0026
0027 #include <iterator>
0028
0029 using namespace edm;
0030 using namespace std;
0031
0032 DTResolutionAnalysisTask::DTResolutionAnalysisTask(const ParameterSet& pset)
0033 : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
0034 edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0035 << "[DTResolutionAnalysisTask] Constructor called!" << endl;
0036
0037
0038 recHits4DToken_ =
0039 consumes<DTRecSegment4DCollection>(edm::InputTag(pset.getUntrackedParameter<string>("recHits4DLabel")));
0040
0041 prescaleFactor = pset.getUntrackedParameter<int>("diagnosticPrescale", 1);
0042 resetCycle = pset.getUntrackedParameter<int>("ResetCycle", -1);
0043
0044 topHistoFolder = pset.getUntrackedParameter<string>("topHistoFolder", "DT/02-Segments");
0045
0046 thePhiHitsCut = pset.getUntrackedParameter<u_int32_t>("phiHitsCut", 8);
0047 theZHitsCut = pset.getUntrackedParameter<u_int32_t>("zHitsCut", 4);
0048 }
0049
0050 DTResolutionAnalysisTask::~DTResolutionAnalysisTask() {
0051 edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0052 << "[DTResolutionAnalysisTask] Destructor called!" << endl;
0053 }
0054
0055 void DTResolutionAnalysisTask::dqmBeginRun(const Run& run, const EventSetup& setup) {
0056
0057 dtGeom = &setup.getData(muonGeomToken_);
0058 }
0059
0060 void DTResolutionAnalysisTask::bookHistograms(DQMStore::IBooker& ibooker,
0061 edm::Run const& iRun,
0062 edm::EventSetup const& ) {
0063
0064 vector<const DTChamber*> chambers = dtGeom->chambers();
0065 for (vector<const DTChamber*>::const_iterator chamber = chambers.begin(); chamber != chambers.end();
0066 ++chamber) {
0067 DTChamberId dtChId = (*chamber)->id();
0068 for (int sl = 1; sl <= 3; ++sl) {
0069 if (dtChId.station() == 4 && sl == 2)
0070 continue;
0071 const DTSuperLayerId dtSLId(dtChId, sl);
0072 bookHistos(ibooker, dtSLId);
0073 }
0074 }
0075 }
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093 void DTResolutionAnalysisTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0094 edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0095 << "[DTResolutionAnalysisTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event() << endl;
0096
0097
0098 edm::Handle<DTRecSegment4DCollection> all4DSegments;
0099 event.getByToken(recHits4DToken_, all4DSegments);
0100
0101
0102 if (!all4DSegments.isValid())
0103 return;
0104
0105
0106 DTRecSegment4DCollection::id_iterator chamberId;
0107 for (chamberId = all4DSegments->id_begin(); chamberId != all4DSegments->id_end(); ++chamberId) {
0108
0109 DTRecSegment4DCollection::range range = all4DSegments->get(*chamberId);
0110
0111
0112 const DTChamber* chamber = dtGeom->chamber(*chamberId);
0113
0114
0115 for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0116
0117
0118 if ((*chamberId).station() != 4 && (*segment4D).dimension() != 4) {
0119 edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0120 << "[DTResolutionAnalysisTask]***Warning: RecSegment dimension is not 4 but " << (*segment4D).dimension()
0121 << "!" << endl;
0122 continue;
0123 } else if ((*chamberId).station() == 4 && (*segment4D).dimension() != 2) {
0124 edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0125 << "[DTResolutionAnalysisTask]***Warning: RecSegment dimension is not 2 but " << (*segment4D).dimension()
0126 << "!" << endl;
0127 continue;
0128 }
0129
0130
0131 vector<DTRecHit1D> recHits1D_S3;
0132
0133
0134
0135
0136 if ((*segment4D).hasPhi()) {
0137 const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
0138 vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
0139
0140 if (phiRecHits.size() < thePhiHitsCut) {
0141 continue;
0142 }
0143 copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
0144 } else {
0145 }
0146
0147 if ((*segment4D).hasZed()) {
0148 const DTSLRecSegment2D* zSeg = (*segment4D).zSegment();
0149 vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
0150 if (zRecHits.size() < theZHitsCut) {
0151 continue;
0152 }
0153 copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
0154 }
0155
0156
0157 for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin(); recHit1D != recHits1D_S3.end();
0158 recHit1D++) {
0159 const DTWireId wireId = (*recHit1D).wireId();
0160
0161
0162 const DTLayer* layer = chamber->superLayer(wireId.superlayerId())->layer(wireId.layerId());
0163 float wireX = layer->specificTopology().wirePosition(wireId.wire());
0164
0165
0166 float distRecHitToWire = fabs(wireX - (*recHit1D).localPosition().x());
0167
0168
0169
0170
0171 LocalPoint wirePosInLay(wireX, (*recHit1D).localPosition().y(), (*recHit1D).localPosition().z());
0172 GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
0173 LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
0174
0175
0176 LocalPoint segPosAtZWire = (*segment4D).localPosition() + (*segment4D).localDirection() * wirePosInChamber.z() /
0177 cos((*segment4D).localDirection().theta());
0178
0179
0180 int sl = wireId.superlayer();
0181
0182 double distSegmToWire = -1;
0183 if (sl == 1 || sl == 3) {
0184
0185 distSegmToWire = fabs(wirePosInChamber.x() - segPosAtZWire.x());
0186 } else if (sl == 2) {
0187
0188 distSegmToWire = fabs(wirePosInChamber.y() - segPosAtZWire.y());
0189 }
0190
0191 if (distSegmToWire > 2.1)
0192 edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0193 << " Warning: dist segment-wire: " << distSegmToWire << endl;
0194
0195 double residual = distRecHitToWire - distSegmToWire;
0196
0197 fillHistos(wireId.superlayerId(), distSegmToWire, residual);
0198
0199 }
0200 }
0201 }
0202
0203 }
0204
0205
0206 void DTResolutionAnalysisTask::bookHistos(DQMStore::IBooker& ibooker, DTSuperLayerId slId) {
0207 edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << " Booking histos for SL: " << slId << endl;
0208
0209
0210 stringstream wheel;
0211 wheel << slId.wheel();
0212 stringstream station;
0213 station << slId.station();
0214 stringstream sector;
0215 sector << slId.sector();
0216 stringstream superLayer;
0217 superLayer << slId.superlayer();
0218
0219 string slHistoName = "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str() + "_SL" + superLayer.str();
0220
0221 ibooker.setCurrentFolder(topHistoFolder + "/Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" +
0222 station.str());
0223
0224 vector<MonitorElement*> histos;
0225
0226 histos.push_back(ibooker.book1D(
0227 "hResDist" + slHistoName, "Residuals on the distance from wire (rec_hit - segm_extr) (cm)", 200, -0.4, 0.4));
0228 histosPerSL[slId] = histos;
0229 }
0230
0231
0232 void DTResolutionAnalysisTask::fillHistos(DTSuperLayerId slId, float distExtr, float residual) {
0233 vector<MonitorElement*> histos = histosPerSL[slId];
0234 histos[0]->Fill(residual);
0235 }
0236
0237
0238
0239
0240