Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-21 22:39:44

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_ = consumes<DTRecSegment4DCollection>(edm::InputTag(pset.getParameter<string>("recHits4DLabel")));
0039 
0040   prescaleFactor = pset.getUntrackedParameter<int>("diagnosticPrescale", 1);
0041   resetCycle = pset.getUntrackedParameter<int>("ResetCycle", -1);
0042   // top folder for the histograms in DQMStore
0043   topHistoFolder = pset.getUntrackedParameter<string>("topHistoFolder", "DT/02-Segments");
0044 
0045   thePhiHitsCut = pset.getUntrackedParameter<u_int32_t>("phiHitsCut", 8);
0046   theZHitsCut = pset.getUntrackedParameter<u_int32_t>("zHitsCut", 4);
0047 }
0048 
0049 DTResolutionAnalysisTask::~DTResolutionAnalysisTask() {
0050   edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0051       << "[DTResolutionAnalysisTask] Destructor called!" << endl;
0052 }
0053 
0054 void DTResolutionAnalysisTask::dqmBeginRun(const Run& run, const EventSetup& setup) {
0055   // Get the DT Geometry
0056   dtGeom = &setup.getData(muonGeomToken_);
0057 }
0058 
0059 void DTResolutionAnalysisTask::bookHistograms(DQMStore::IBooker& ibooker,
0060                                               edm::Run const& iRun,
0061                                               edm::EventSetup const& /* iSetup */) {
0062   // Book the histograms
0063   vector<const DTChamber*> chambers = dtGeom->chambers();
0064   for (vector<const DTChamber*>::const_iterator chamber = chambers.begin(); chamber != chambers.end();
0065        ++chamber) {  // Loop over all chambers
0066     DTChamberId dtChId = (*chamber)->id();
0067     for (int sl = 1; sl <= 3; ++sl) {  // Loop over SLs
0068       if (dtChId.station() == 4 && sl == 2)
0069         continue;
0070       const DTSuperLayerId dtSLId(dtChId, sl);
0071       bookHistos(ibooker, dtSLId);
0072     }
0073   }
0074 }
0075 /*
0076 void DTResolutionAnalysisTask::beginLuminosityBlock(const LuminosityBlock& lumiSeg, const EventSetup& context) {
0077   edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0078       << "[DTResolutionTask]: Begin of LS transition" << endl;
0079 
0080   if (resetCycle != -1 && lumiSeg.id().luminosityBlock() % resetCycle == 0) {
0081     for (map<DTSuperLayerId, vector<MonitorElement*> >::const_iterator histo = histosPerSL.begin();
0082          histo != histosPerSL.end();
0083          histo++) {
0084       int size = (*histo).second.size();
0085       for (int i = 0; i < size; i++) {
0086         (*histo).second[i]->Reset();
0087       }
0088     }
0089   }
0090 }
0091 */
0092 void DTResolutionAnalysisTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0093   edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0094       << "[DTResolutionAnalysisTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event() << endl;
0095 
0096   // Get the 4D segment collection from the event
0097   edm::Handle<DTRecSegment4DCollection> all4DSegments;
0098   event.getByToken(recHits4DToken_, all4DSegments);
0099 
0100   // check the validity of the collection
0101   if (!all4DSegments.isValid())
0102     return;
0103 
0104   // Loop over all chambers containing a segment
0105   DTRecSegment4DCollection::id_iterator chamberId;
0106   for (chamberId = all4DSegments->id_begin(); chamberId != all4DSegments->id_end(); ++chamberId) {
0107     // Get the range for the corresponding ChamerId
0108     DTRecSegment4DCollection::range range = all4DSegments->get(*chamberId);
0109 
0110     // Get the chamber
0111     const DTChamber* chamber = dtGeom->chamber(*chamberId);
0112 
0113     // Loop over the rechits of this ChamerId
0114     for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0115       // If Statio != 4 skip RecHits with dimension != 4
0116       // For the Station 4 consider 2D RecHits
0117       if ((*chamberId).station() != 4 && (*segment4D).dimension() != 4) {
0118         edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0119             << "[DTResolutionAnalysisTask]***Warning: RecSegment dimension is not 4 but " << (*segment4D).dimension()
0120             << "!" << endl;
0121         continue;
0122       } else if ((*chamberId).station() == 4 && (*segment4D).dimension() != 2) {
0123         edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0124             << "[DTResolutionAnalysisTask]***Warning: RecSegment dimension is not 2 but " << (*segment4D).dimension()
0125             << "!" << endl;
0126         continue;
0127       }
0128 
0129       // Get all 1D RecHits at step 3 within the 4D segment
0130       vector<DTRecHit1D> recHits1D_S3;
0131 
0132       // Get 1D RecHits at Step 3 and select only events with
0133       // 8 hits in phi and 4 hits in theta (if any)
0134 
0135       if ((*segment4D).hasPhi()) {  // has phi component
0136         const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
0137         vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
0138 
0139         if (phiRecHits.size() < thePhiHitsCut) {
0140           continue;
0141         }
0142         copy(phiRecHits.begin(), phiRecHits.end(), back_inserter(recHits1D_S3));
0143       } else {
0144       }
0145 
0146       if ((*segment4D).hasZed()) {
0147         const DTSLRecSegment2D* zSeg = (*segment4D).zSegment();
0148         vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
0149         if (zRecHits.size() < theZHitsCut) {
0150           continue;
0151         }
0152         copy(zRecHits.begin(), zRecHits.end(), back_inserter(recHits1D_S3));
0153       }
0154 
0155       // Loop over 1D RecHit inside 4D segment
0156       for (vector<DTRecHit1D>::const_iterator recHit1D = recHits1D_S3.begin(); recHit1D != recHits1D_S3.end();
0157            recHit1D++) {
0158         const DTWireId wireId = (*recHit1D).wireId();
0159 
0160         // Get the layer and the wire position
0161         const DTLayer* layer = chamber->superLayer(wireId.superlayerId())->layer(wireId.layerId());
0162         float wireX = layer->specificTopology().wirePosition(wireId.wire());
0163 
0164         // Distance of the 1D rechit from the wire
0165         float distRecHitToWire = fabs(wireX - (*recHit1D).localPosition().x());
0166 
0167         // Extrapolate the segment to the z of the wire
0168 
0169         // Get wire position in chamber RF
0170         LocalPoint wirePosInLay(wireX, (*recHit1D).localPosition().y(), (*recHit1D).localPosition().z());
0171         GlobalPoint wirePosGlob = layer->toGlobal(wirePosInLay);
0172         LocalPoint wirePosInChamber = chamber->toLocal(wirePosGlob);
0173 
0174         // Segment position at Wire z in chamber local frame
0175         LocalPoint segPosAtZWire = (*segment4D).localPosition() + (*segment4D).localDirection() * wirePosInChamber.z() /
0176                                                                       cos((*segment4D).localDirection().theta());
0177 
0178         // Compute the distance of the segment from the wire
0179         int sl = wireId.superlayer();
0180 
0181         double distSegmToWire = -1;
0182         if (sl == 1 || sl == 3) {
0183           // RPhi SL
0184           distSegmToWire = fabs(wirePosInChamber.x() - segPosAtZWire.x());
0185         } else if (sl == 2) {
0186           // RZ SL
0187           distSegmToWire = fabs(wirePosInChamber.y() - segPosAtZWire.y());
0188         }
0189 
0190         if (distSegmToWire > 2.1)
0191           edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask")
0192               << "  Warning: dist segment-wire: " << distSegmToWire << endl;
0193 
0194         double residual = distRecHitToWire - distSegmToWire;
0195         // FIXME: Fill the histos
0196         fillHistos(wireId.superlayerId(), distSegmToWire, residual);
0197 
0198       }  // End of loop over 1D RecHit inside 4D segment
0199     }    // End of loop over the rechits of this ChamerId
0200   }
0201   // -----------------------------------------------------------------------------
0202 }
0203 
0204 // Book a set of histograms for a given SL
0205 void DTResolutionAnalysisTask::bookHistos(DQMStore::IBooker& ibooker, DTSuperLayerId slId) {
0206   edm::LogVerbatim("DTDQM|DTMonitorModule|DTResolutionAnalysisTask") << "   Booking histos for SL: " << slId << endl;
0207 
0208   // Compose the chamber name
0209   stringstream wheel;
0210   wheel << slId.wheel();
0211   stringstream station;
0212   station << slId.station();
0213   stringstream sector;
0214   sector << slId.sector();
0215   stringstream superLayer;
0216   superLayer << slId.superlayer();
0217 
0218   string slHistoName = "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str() + "_SL" + superLayer.str();
0219 
0220   ibooker.setCurrentFolder(topHistoFolder + "/Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" +
0221                            station.str());
0222   // Create the monitor elements
0223   vector<MonitorElement*> histos;
0224   // Note the order matters
0225   histos.push_back(ibooker.book1D(
0226       "hResDist" + slHistoName, "Residuals on the distance from wire (rec_hit - segm_extr) (cm)", 200, -0.4, 0.4));
0227   histosPerSL[slId] = histos;
0228 }
0229 
0230 // Fill a set of histograms for a given SL
0231 void DTResolutionAnalysisTask::fillHistos(DTSuperLayerId slId, float distExtr, float residual) {
0232   vector<MonitorElement*> histos = histosPerSL[slId];
0233   histos[0]->Fill(residual);
0234 }
0235 
0236 // Local Variables:
0237 // show-trailing-whitespace: t
0238 // truncate-lines: t
0239 // End: