DTMonitorModule/src/DTResolutionAnalysisTask.cc

0001
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author G. Cerminara - INFN Torino
0006  */
0007
0008 #include "DTResolutionAnalysisTask.h"
0009
0010 // Framework
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 //#include "FWCore/Framework/interface/LuminosityBlock.h"
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 //Geometry
0022 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0023
0024 //RecHit
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   // the name of the 4D rec hits collection
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   // top folder for the histograms in DQMStore
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   // Get the DT Geometry
0057   dtGeom = &setup.getData(muonGeomToken_);
0058 }
0059
0060 void DTResolutionAnalysisTask::bookHistograms(DQMStore::IBooker& ibooker,
0061                                               edm::Run const& iRun,
0062                                               edm::EventSetup const& /* iSetup */) {
0063   // Book the histograms
0064   vector<const DTChamber*> chambers = dtGeom->chambers();
0065   for (vector<const DTChamber*>::const_iterator chamber = chambers.begin(); chamber != chambers.end();
0066        ++chamber) {  // Loop over all chambers
0067     DTChamberId dtChId = (*chamber)->id();
0068     for (int sl = 1; sl <= 3; ++sl) {  // Loop over SLs
0069       if (dtChId.station() == 4 && sl == 2)
0070         continue;
0071       const DTSuperLayerId dtSLId(dtChId, sl);
0072       bookHistos(ibooker, dtSLId);
0073     }
0074   }
0075 }
0076 /*
0077 void DTResolutionAnalysisTask::beginLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& context) {
0078   edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0079       << "[DTResolutionTask]: Begin of LS transition" << endl;
0080
0081   if (resetCycle != -1 && lumiSeg.id().luminosityBlock() % resetCycle == 0) {
0082     for (map<DTSuperLayerId, vector<MonitorElement*> >::const_iterator histo = histosPerSL.begin();
0083          histo != histosPerSL.end();
0084          histo++) {
0085       int size = (*histo).second.size();
0086       for (int i = 0; i < size; i++) {
0087         (*histo).second[i]->Reset();
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   // Get the 4D segment collection from the event
0098   edm::Handle<DTRecSegment4DCollection> all4DSegments;
0099   event.getByToken(recHits4DToken_, all4DSegments);
0100
0101   // check the validity of the collection
0102   if (!all4DSegments.isValid())
0103     return;
0104
0105   // Loop over all chambers containing a segment
0106   DTRecSegment4DCollection::id_iterator chamberId;
0107   for (chamberId = all4DSegments->id_begin(); chamberId != all4DSegments->id_end(); ++chamberId) {
0108     // Get the range for the corresponding ChamerId
0109     DTRecSegment4DCollection::range range = all4DSegments->get(*chamberId);
0110
0111     // Get the chamber
0112     const DTChamber* chamber = dtGeom->chamber(*chamberId);
0113
0114     // Loop over the rechits of this ChamerId
0115     for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0116       // If Statio != 4 skip RecHits with dimension != 4
0117       // For the Station 4 consider 2D RecHits
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       // Get all 1D RecHits at step 3 within the 4D segment
0131       vector<DTRecHit1D> recHits1D_S3;
0132
0133       // Get 1D RecHits at Step 3 and select only events with
0134       // 8 hits in phi and 4 hits in theta (if any)
0135
0136       if ((*segment4D).hasPhi()) {  // has phi component
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       // Loop over 1D RecHit inside 4D segment
0157       for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin(); recHit1D != recHits1D_S3.end();
0158            recHit1D++) {
0159         const DTWireId wireId = (*recHit1D).wireId();
0160
0161         // Get the layer and the wire position
0162         const DTLayer* layer = chamber->superLayer(wireId.superlayerId())->layer(wireId.layerId());
0163         float wireX = layer->specificTopology().wirePosition(wireId.wire());
0164
0165         // Distance of the 1D rechit from the wire
0166         float distRecHitToWire = fabs(wireX - (*recHit1D).localPosition().x());
0167
0168         // Extrapolate the segment to the z of the wire
0169
0170         // Get wire position in chamber RF
0171         LocalPoint wirePosInLay(wireX, (*recHit1D).localPosition().y(), (*recHit1D).localPosition().z());
0172         GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
0173         LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
0174
0175         // Segment position at Wire z in chamber local frame
0176         LocalPoint segPosAtZWire = (*segment4D).localPosition() + (*segment4D).localDirection() * wirePosInChamber.z() /
0177                                                                       cos((*segment4D).localDirection().theta());
0178
0179         // Compute the distance of the segment from the wire
0180         int sl = wireId.superlayer();
0181
0182         double distSegmToWire = -1;
0183         if (sl == 1 || sl == 3) {
0184           // RPhi SL
0185           distSegmToWire = fabs(wirePosInChamber.x() - segPosAtZWire.x());
0186         } else if (sl == 2) {
0187           // RZ SL
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         // FIXME: Fill the histos
0197         fillHistos(wireId.superlayerId(), distSegmToWire, residual);
0198
0199       }  // End of loop over 1D RecHit inside 4D segment
0200     }    // End of loop over the rechits of this ChamerId
0201   }
0202   // -----------------------------------------------------------------------------
0203 }
0204
0205 // Book a set of histograms for a given SL
0206 void DTResolutionAnalysisTask::bookHistos(DQMStore::IBooker& ibooker, DTSuperLayerId slId) {
0207   edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "   Booking histos for SL: " << slId << endl;
0208
0209   // Compose the chamber name
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   // Create the monitor elements
0224   vector<MonitorElement*> histos;
0225   // Note the order matters
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 // Fill a set of histograms for a given SL
0232 void DTResolutionAnalysisTask::fillHistos(DTSuperLayerId slId, float distExtr, float residual) {
0233   vector<MonitorElement*> histos = histosPerSL[slId];
0234   histos[0]->Fill(residual);
0235 }
0236
0237 // Local Variables:
0238 // show-trailing-whitespace: t
0239 // truncate-lines: t
0240 // End: