Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:07:10

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Cerminara - INFN Torino
0005  *  revised by G. Mila - INFN Torino
0006  */
0007 
0008 #include "DTSegmentAnalysisTask.h"
0009 
0010 // Framework
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 
0017 #include "DQMServices/Core/interface/DQMStore.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 
0020 //Geometry
0021 #include "DataFormats/GeometryVector/interface/Pi.h"
0022 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0023 #include "Geometry/DTGeometry/interface/DTTopology.h"
0024 
0025 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0026 
0027 #include <iterator>
0028 #include "TMath.h"
0029 
0030 using namespace edm;
0031 using namespace std;
0032 
0033 DTSegmentAnalysisTask::DTSegmentAnalysisTask(const edm::ParameterSet& pset)
0034     : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()), statusMapToken_(esConsumes()), nevents(0) {
0035   edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "[DTSegmentAnalysisTask] Constructor called!";
0036 
0037   // switch for detailed analysis
0038   detailedAnalysis = pset.getUntrackedParameter<bool>("detailedAnalysis", false);
0039   // the name of the 4D rec hits collection
0040   recHits4DToken_ =
0041       consumes<DTRecSegment4DCollection>(edm::InputTag(pset.getUntrackedParameter<string>("recHits4DLabel")));
0042   // Get the map of noisy channels
0043   checkNoisyChannels = pset.getUntrackedParameter<bool>("checkNoisyChannels", false);
0044   phiSegmCut = pset.getUntrackedParameter<double>("phiSegmCut", 30.);
0045   nhitsCut = pset.getUntrackedParameter<int>("nhitsCut", 12);
0046 
0047   // top folder for the histograms in DQMStore
0048   topHistoFolder = pset.getUntrackedParameter<string>("topHistoFolder", "DT/02-Segments");
0049   // hlt DQM mode
0050   hltDQMMode = pset.getUntrackedParameter<bool>("hltDQMMode", false);
0051 }
0052 
0053 DTSegmentAnalysisTask::~DTSegmentAnalysisTask() {
0054   //FR moved fron endjob
0055   edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "[DTSegmentAnalysisTask] Destructor called!";
0056 }
0057 
0058 void DTSegmentAnalysisTask::dqmBeginRun(const Run& run, const edm::EventSetup& context) {
0059   // Get the DT Geometry
0060   dtGeom = &context.getData(muonGeomToken_);
0061 }
0062 
0063 void DTSegmentAnalysisTask::bookHistograms(DQMStore::IBooker& ibooker,
0064                                            edm::Run const& iRun,
0065                                            edm::EventSetup const& context) {
0066   if (!hltDQMMode) {
0067     ibooker.setCurrentFolder("DT/EventInfo/Counters");
0068     nEventMonitor = ibooker.bookFloat("nProcessedEventsSegment");
0069   }
0070 
0071   for (int wh = -2; wh <= 2; wh++) {
0072     stringstream wheel;
0073     wheel << wh;
0074     ibooker.setCurrentFolder(topHistoFolder + "/Wheel" + wheel.str());
0075     string histoName = "numberOfSegments_W" + wheel.str();
0076 
0077     summaryHistos[wh] = ibooker.book2D(histoName.c_str(), histoName.c_str(), 12, 1, 13, 4, 1, 5);
0078     summaryHistos[wh]->setAxisTitle("Sector", 1);
0079     summaryHistos[wh]->setBinLabel(1, "1", 1);
0080     summaryHistos[wh]->setBinLabel(2, "2", 1);
0081     summaryHistos[wh]->setBinLabel(3, "3", 1);
0082     summaryHistos[wh]->setBinLabel(4, "4", 1);
0083     summaryHistos[wh]->setBinLabel(5, "5", 1);
0084     summaryHistos[wh]->setBinLabel(6, "6", 1);
0085     summaryHistos[wh]->setBinLabel(7, "7", 1);
0086     summaryHistos[wh]->setBinLabel(8, "8", 1);
0087     summaryHistos[wh]->setBinLabel(9, "9", 1);
0088     summaryHistos[wh]->setBinLabel(10, "10", 1);
0089     summaryHistos[wh]->setBinLabel(11, "11", 1);
0090     summaryHistos[wh]->setBinLabel(12, "12", 1);
0091     summaryHistos[wh]->setBinLabel(1, "MB1", 2);
0092     summaryHistos[wh]->setBinLabel(2, "MB2", 2);
0093     summaryHistos[wh]->setBinLabel(3, "MB3", 2);
0094     summaryHistos[wh]->setBinLabel(4, "MB4", 2);
0095   }
0096 
0097   // loop over all the DT chambers & book the histos
0098   const vector<const DTChamber*>& chambers = dtGeom->chambers();
0099   vector<const DTChamber*>::const_iterator ch_it = chambers.begin();
0100   vector<const DTChamber*>::const_iterator ch_end = chambers.end();
0101   for (; ch_it != ch_end; ++ch_it) {
0102     bookHistos(ibooker, (*ch_it)->id());
0103   }
0104 }
0105 
0106 void DTSegmentAnalysisTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0107   nevents++;
0108   nEventMonitor->Fill(nevents);
0109 
0110   edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask")
0111       << "[DTSegmentAnalysisTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event();
0112   if (!(event.id().event() % 1000))
0113     edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask")
0114         << "[DTSegmentAnalysisTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event();
0115 
0116   if (checkNoisyChannels) {
0117     statusMap = &setup.getData(statusMapToken_);
0118   }
0119 
0120   // -- 4D segment analysis  -----------------------------------------------------
0121 
0122   // Get the 4D segment collection from the event
0123   edm::Handle<DTRecSegment4DCollection> all4DSegments;
0124   event.getByToken(recHits4DToken_, all4DSegments);
0125 
0126   if (!all4DSegments.isValid())
0127     return;
0128 
0129   // Loop over all chambers containing a segment
0130   DTRecSegment4DCollection::id_iterator chamberId;
0131   for (chamberId = all4DSegments->id_begin(); chamberId != all4DSegments->id_end(); ++chamberId) {
0132     // Get the range for the corresponding ChamerId
0133     DTRecSegment4DCollection::range range = all4DSegments->get(*chamberId);
0134 
0135     edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask")
0136         << "   Chamber: " << *chamberId << " has " << distance(range.first, range.second) << " 4D segments";
0137 
0138     // Loop over the rechits of this ChamerId
0139     for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D != range.second; ++segment4D) {
0140       //FOR NOISY CHANNELS////////////////////////////////
0141       bool segmNoisy = false;
0142       if (checkNoisyChannels) {
0143         if ((*segment4D).hasPhi()) {
0144           const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment();
0145           vector<DTRecHit1D> phiHits = phiSeg->specificRecHits();
0146           map<DTSuperLayerId, vector<DTRecHit1D> > hitsBySLMap;
0147           for (vector<DTRecHit1D>::const_iterator hit = phiHits.begin(); hit != phiHits.end(); ++hit) {
0148             DTWireId wireId = (*hit).wireId();
0149 
0150             // Check for noisy channels to skip them
0151             bool isNoisy = false;
0152             bool isFEMasked = false;
0153             bool isTDCMasked = false;
0154             bool isTrigMask = false;
0155             bool isDead = false;
0156             bool isNohv = false;
0157             statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
0158             if (isNoisy) {
0159               edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask")
0160                   << "Wire: " << wireId << " is noisy, skipping!";
0161               segmNoisy = true;
0162             }
0163           }
0164         }
0165 
0166         if ((*segment4D).hasZed()) {
0167           const DTSLRecSegment2D* zSeg = (*segment4D).zSegment();  // zSeg lives in the SL RF
0168           // Check for noisy channels to skip them
0169           vector<DTRecHit1D> zHits = zSeg->specificRecHits();
0170           for (vector<DTRecHit1D>::const_iterator hit = zHits.begin(); hit != zHits.end(); ++hit) {
0171             DTWireId wireId = (*hit).wireId();
0172             bool isNoisy = false;
0173             bool isFEMasked = false;
0174             bool isTDCMasked = false;
0175             bool isTrigMask = false;
0176             bool isDead = false;
0177             bool isNohv = false;
0178             statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
0179             if (isNoisy) {
0180               edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask")
0181                   << "Wire: " << wireId << " is noisy, skipping!";
0182               segmNoisy = true;
0183             }
0184           }
0185         }
0186 
0187       }  // end of switch on noisy channels
0188       if (segmNoisy) {
0189         edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask")
0190             << "skipping the segment: it contains noisy cells";
0191         continue;
0192       }
0193       //END FOR NOISY CHANNELS////////////////////////////////
0194 
0195       int nHits = 0;
0196       if ((*segment4D).hasPhi())
0197         nHits = (((*segment4D).phiSegment())->specificRecHits()).size();
0198       if ((*segment4D).hasZed())
0199         nHits = nHits + ((((*segment4D).zSegment())->specificRecHits()).size());
0200 
0201       double anglePhiSegm(0.);
0202       if ((*segment4D).hasPhi()) {
0203         double xdir = (*segment4D).phiSegment()->localDirection().x();
0204         double zdir = (*segment4D).phiSegment()->localDirection().z();
0205 
0206         anglePhiSegm = atan(xdir / zdir) * 180. / TMath::Pi();
0207       }
0208       if (fabs(anglePhiSegm) > phiSegmCut)
0209         continue;
0210       // If the segment is in Wh+-2/SecX/MB1, get the DT chambers just above and check if there is a segment
0211       // to validate the segment present in MB1
0212       if (fabs((*chamberId).wheel()) == 2 && (*chamberId).station() == 1) {
0213         bool segmOk = false;
0214         int mb(2);
0215         while (mb < 4) {
0216           DTChamberId checkMB((*chamberId).wheel(), mb, (*chamberId).sector());
0217           DTRecSegment4DCollection::range ckrange = all4DSegments->get(checkMB);
0218 
0219           for (DTRecSegment4DCollection::const_iterator cksegment4D = ckrange.first; cksegment4D != ckrange.second;
0220                ++cksegment4D) {
0221             int nHits = 0;
0222             if ((*cksegment4D).hasPhi())
0223               nHits = (((*cksegment4D).phiSegment())->specificRecHits()).size();
0224             if ((*cksegment4D).hasZed())
0225               nHits = nHits + ((((*cksegment4D).zSegment())->specificRecHits()).size());
0226 
0227             if (nHits >= nhitsCut)
0228               segmOk = true;
0229           }
0230           mb++;
0231         }
0232 
0233         if (!segmOk)
0234           continue;
0235       }
0236       fillHistos(*chamberId, nHits, (*segment4D).chi2() / (*segment4D).degreesOfFreedom());
0237     }
0238   }
0239 
0240   // -----------------------------------------------------------------------------
0241 }
0242 
0243 // Book a set of histograms for a give chamber
0244 void DTSegmentAnalysisTask::bookHistos(DQMStore::IBooker& ibooker, DTChamberId chamberId) {
0245   edm::LogVerbatim("DTDQM|DTMonitorModule|DTSegmentAnalysisTask") << "   Booking histos for chamber: " << chamberId;
0246 
0247   // Compose the chamber name
0248   stringstream wheel;
0249   wheel << chamberId.wheel();
0250   stringstream station;
0251   station << chamberId.station();
0252   stringstream sector;
0253   sector << chamberId.sector();
0254 
0255   string chamberHistoName = "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str();
0256 
0257   ibooker.setCurrentFolder(topHistoFolder + "/Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" +
0258                            station.str());
0259 
0260   // Create the monitor elements
0261   vector<MonitorElement*> histos;
0262   histos.push_back(ibooker.book1D("h4DSegmNHits" + chamberHistoName, "# of hits per segment", 16, 0.5, 16.5));
0263   if (detailedAnalysis) {
0264     histos.push_back(ibooker.book1D("h4DChi2" + chamberHistoName, "4D Segment reduced Chi2", 20, 0, 20));
0265   }
0266   histosPerCh[chamberId] = histos;
0267 }
0268 
0269 // Fill a set of histograms for a give chamber
0270 void DTSegmentAnalysisTask::fillHistos(DTChamberId chamberId, int nHits, float chi2) {
0271   int sector = chamberId.sector();
0272   if (chamberId.sector() == 13) {
0273     sector = 4;
0274   } else if (chamberId.sector() == 14) {
0275     sector = 10;
0276   }
0277 
0278   summaryHistos[chamberId.wheel()]->Fill(sector, chamberId.station());
0279 
0280   vector<MonitorElement*> histos = histosPerCh[chamberId];
0281   histos[0]->Fill(nHits);
0282   if (detailedAnalysis) {
0283     histos[1]->Fill(chi2);
0284   }
0285 }
0286 
0287 // Local Variables:
0288 // show-trailing-whitespace: t
0289 // truncate-lines: t
0290 // End: