Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-13 02:31:35 

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Mila - INFN Torino
0005  */
0006 
0007 #include "DTChamberEfficiencyTask.h"
0008 
0009 //Framework
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 
0014 #include "DQMServices/Core/interface/DQMStore.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include <iterator>
0018 #include <iostream>
0019 #include <cmath>
0020 using namespace edm;
0021 using namespace std;
0022 
0023 DTChamberEfficiencyTask::DTChamberEfficiencyTask(const ParameterSet& pset)
0024     : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
0025   debug = pset.getUntrackedParameter<bool>("debug", false);
0026 
0027   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "[DTChamberEfficiencyTask] Constructor called!";
0028 
0029   parameters = pset;
0030 
0031   // the name of the 4D rec hits collection
0032   recHits4DToken_ =
0033       consumes<DTRecSegment4DCollection>(edm::InputTag(parameters.getUntrackedParameter<string>("recHits4DLabel")));
0034 
0035   // parameters to use for the segment quality check
0036   theMinHitsSegment = static_cast<unsigned int>(parameters.getParameter<int>("minHitsSegment"));
0037   theMinChi2NormSegment = parameters.getParameter<double>("minChi2NormSegment");
0038   // parameter to use for the exstrapolated segment check
0039   theMinCloseDist = parameters.getParameter<double>("minCloseDist");
0040 
0041   // the running modality
0042   onlineMonitor = parameters.getUntrackedParameter<bool>("onlineMonitor");
0043 
0044   // the analysis mode
0045   detailedAnalysis = parameters.getUntrackedParameter<bool>("detailedAnalysis");
0046 }
0047 
0048 DTChamberEfficiencyTask::~DTChamberEfficiencyTask() {
0049   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "[DTChamberEfficiencyTask] Destructor called!";
0050 }
0051 
0052 void DTChamberEfficiencyTask::beginLuminosityBlock(LuminosityBlock const& lumiSeg, EventSetup const& context) {
0053   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask")
0054       << "[DTChamberEfficiencyTask]: Begin of LS transition";
0055 
0056   if (lumiSeg.id().luminosityBlock() % parameters.getUntrackedParameter<int>("ResetCycle", 3) == 0 && onlineMonitor) {
0057     for (map<DTChamberId, vector<MonitorElement*> >::const_iterator histo = histosPerCh.begin();
0058          histo != histosPerCh.end();
0059          histo++) {
0060       int size = (*histo).second.size();
0061       for (int i = 0; i < size; i++) {
0062         (*histo).second[i]->Reset();
0063       }
0064     }
0065   }
0066 }
0067 
0068 void DTChamberEfficiencyTask::dqmBeginRun(const edm::Run& run, const edm::EventSetup& setup) {
0069   // Get the DT Geometry
0070   dtGeom = &setup.getData(muonGeomToken_);
0071 }
0072 
0073 void DTChamberEfficiencyTask::bookHistograms(DQMStore::IBooker& ibooker,
0074                                              edm::Run const& iRun,
0075                                              edm::EventSetup const& context) {
0076   ibooker.setCurrentFolder("DT/DTChamberEfficiencyTask");
0077 
0078   // Loop over all the chambers
0079   vector<const DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
0080   vector<const DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
0081   for (; ch_it != ch_end; ++ch_it) {
0082     // histo booking
0083     bookHistos(ibooker, (*ch_it)->id());
0084   }
0085 }
0086 
0087 // Book a set of histograms for a given Layer
0088 void DTChamberEfficiencyTask::bookHistos(DQMStore::IBooker& ibooker, DTChamberId chId) {
0089   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask") << "   Booking histos for CH : " << chId;
0090 
0091   // Compose the chamber name
0092   stringstream wheel;
0093   wheel << chId.wheel();
0094   stringstream station;
0095   station << chId.station();
0096   stringstream sector;
0097   sector << chId.sector();
0098 
0099   string HistoName = "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str();
0100 
0101   ibooker.setCurrentFolder("DT/01-DTChamberEfficiency/Task/Wheel" + wheel.str() + "/Sector" + sector.str() +
0102                            "/Station" + station.str());
0103 
0104   // Create the monitor elements
0105   vector<MonitorElement*> histos;
0106 
0107   //efficiency selection cuts
0108   // a- number of segments of the top chamber > 0 && number of segments of the bottom chamber > 0
0109   // b- number of segments of the middle chamber > 0
0110   // c- check of the top and bottom segment quality
0111   // d- check if interpolation falls inside the middle chamber
0112   // e- check of the middle segment quality
0113   // f- check if the distance between the reconstructed and the exstrapolated segments is ok
0114 
0115   // histo for efficiency with cuts a-/c-/d-
0116   histos.push_back(ibooker.book2D(
0117       "hEffGoodSegVsPosDen" + HistoName, "Eff vs local position (good) ", 25, -250., 250., 25, -250., 250.));
0118   // histo for efficiency with cuts a-/b-/c-/d-/e-/f-
0119   histos.push_back(ibooker.book2D("hEffGoodCloseSegVsPosNum" + HistoName,
0120                                   "Eff vs local position (good and close segs) ",
0121                                   25,
0122                                   -250.,
0123                                   250.,
0124                                   25,
0125                                   -250.,
0126                                   250.));
0127   if (detailedAnalysis) {
0128     histos.push_back(
0129         ibooker.book1D("hDistSegFromExtrap" + HistoName, "Distance segments from extrap position ", 200, 0., 200.));
0130     // histo for efficiency from segment counting
0131     histos.push_back(ibooker.book1D("hNaiveEffSeg" + HistoName, "Naive eff ", 10, 0., 10.));
0132     // histo for efficiency with cuts a-/c-
0133     histos.push_back(ibooker.book2D(
0134         "hEffSegVsPosDen" + HistoName, "Eff vs local position (all) ", 25, -250., 250., 25, -250., 250.));
0135     // histo for efficiency with cuts a-/b-/c-/d-
0136     histos.push_back(
0137         ibooker.book2D("hEffSegVsPosNum" + HistoName, "Eff vs local position ", 25, -250., 250., 25, -250., 250.));
0138     // histo for efficiency with cuts a-/b-/c-/d-/e-
0139     histos.push_back(ibooker.book2D(
0140         "hEffGoodSegVsPosNum" + HistoName, "Eff vs local position (good segs) ", 25, -250., 250., 25, -250., 250.));
0141   }
0142   histosPerCh[chId] = histos;
0143 }
0144 
0145 void DTChamberEfficiencyTask::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0146   edm::LogVerbatim("DTDQM|DTMonitorModule|DTChamberEfficiencyTask")
0147       << "[DTChamberEfficiencyTask] Analyze #Run: " << event.id().run() << " #Event: " << event.id().event();
0148 
0149   // Get the 4D rechit collection from the event
0150   event.getByToken(recHits4DToken_, segs);
0151 
0152   int bottom = 0, top = 0;
0153 
0154   // Loop over all the chambers
0155   vector<const DTChamber*>::const_iterator ch_it = dtGeom->chambers().begin();
0156   vector<const DTChamber*>::const_iterator ch_end = dtGeom->chambers().end();
0157   for (; ch_it != ch_end; ++ch_it) {
0158     DTChamberId ch = (*ch_it)->id();
0159     int wheel = ch.wheel();
0160     int sector = ch.sector();
0161     int station = ch.station();
0162 
0163     DTChamberId MidId(wheel, station, sector);
0164 
0165     // get efficiency for MB1 using MB2 and MB3
0166     if (station == 1) {
0167       bottom = 2;
0168       top = 3;
0169     }
0170 
0171     // get efficiency for MB2 using MB1 and MB3
0172     if (station == 2) {
0173       bottom = 1;
0174       top = 3;
0175     }
0176 
0177     // get efficiency for MB2 using MB2 and MB4
0178     if (station == 3) {
0179       bottom = 2;
0180       top = 4;
0181     }
0182 
0183     // get efficiency for MB4 using MB2 and MB3
0184     if (station == 4) {
0185       bottom = 2;
0186       top = 3;
0187     }
0188 
0189     // Select events with (good) segments in Bot and Top
0190     DTChamberId BotId(wheel, bottom, sector);
0191     DTChamberId TopId(wheel, top, sector);
0192 
0193     // Get segments in the bottom chambers (if any)
0194     DTRecSegment4DCollection::range segsBot = segs->get(BotId);
0195     int nSegsBot = segsBot.second - segsBot.first;
0196     // check if any segments is there
0197     if (nSegsBot == 0)
0198       continue;
0199 
0200     vector<MonitorElement*> histos = histosPerCh[MidId];
0201 
0202     // Get segments in the top chambers (if any)
0203     DTRecSegment4DCollection::range segsTop = segs->get(TopId);
0204     int nSegsTop = segsTop.second - segsTop.first;
0205 
0206     // Select one segment for the bottom chamber
0207     const DTRecSegment4D& bestBotSeg = getBestSegment(segsBot);
0208 
0209     // Select one segment for the top chamber
0210     DTRecSegment4D* pBestTopSeg = nullptr;
0211     if (nSegsTop > 0)
0212       pBestTopSeg = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop));
0213     //if top chamber is MB4 sector 10, consider also sector 14
0214     if (TopId.station() == 4 && TopId.sector() == 10) {
0215       DTChamberId TopId14(wheel, top, 14);
0216       DTRecSegment4DCollection::range segsTop14 = segs->get(TopId14);
0217       int nSegsTop14 = segsTop14.second - segsTop14.first;
0218       if (nSegsTop14) {
0219         DTRecSegment4D* pBestTopSeg14 = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop14));
0220         // get best between sector 10 and 14
0221         pBestTopSeg = const_cast<DTRecSegment4D*>(getBestSegment(pBestTopSeg, pBestTopSeg14));
0222       }
0223     }
0224     if (!pBestTopSeg)
0225       continue;
0226     const DTRecSegment4D& bestTopSeg = *pBestTopSeg;
0227 
0228     DTRecSegment4DCollection::range segsMid = segs->get(MidId);
0229     int nSegsMid = segsMid.second - segsMid.first;
0230 
0231     if (detailedAnalysis) {
0232       // very trivial efficiency, just count segments
0233       histos[3]->Fill(0);
0234       if (nSegsMid > 0)
0235         histos[3]->Fill(1);
0236     }
0237 
0238     // get position at Mid by interpolating the position (not direction) of best
0239     // segment in Bot and Top to Mid surface
0240     LocalPoint posAtMid = interpolate(bestBotSeg, bestTopSeg, MidId);
0241 
0242     // is best segment good enough?
0243     if (isGoodSegment(bestBotSeg) && isGoodSegment(bestTopSeg)) {
0244       if (detailedAnalysis)
0245         histos[4]->Fill(posAtMid.x(), posAtMid.y());
0246       //check if interpolation fall inside middle chamber
0247       if ((dtGeom->chamber(MidId))->surface().bounds().inside(posAtMid)) {
0248         histos[0]->Fill(posAtMid.x(), posAtMid.y());
0249         if (nSegsMid > 0) {
0250           if (detailedAnalysis) {
0251             histos[3]->Fill(2);
0252             histos[5]->Fill(posAtMid.x(), posAtMid.y());
0253           }
0254 
0255           const DTRecSegment4D& bestMidSeg = getBestSegment(segsMid);
0256           // check if middle segments is good enough
0257           if (isGoodSegment(bestMidSeg)) {
0258             if (detailedAnalysis)
0259               histos[6]->Fill(posAtMid.x(), posAtMid.y());
0260             LocalPoint midSegPos = bestMidSeg.localPosition();
0261 
0262             // check if middle segments is also close enough
0263             double dist;
0264             if (bestMidSeg.hasPhi()) {
0265               if (bestTopSeg.hasZed() && bestBotSeg.hasZed() && bestMidSeg.hasZed()) {
0266                 dist = (midSegPos - posAtMid).mag();
0267               } else {
0268                 dist = fabs((midSegPos - posAtMid).x());
0269               }
0270             } else {
0271               dist = fabs((midSegPos - posAtMid).y());
0272             }
0273             if (dist < theMinCloseDist) {
0274               histos[1]->Fill(posAtMid.x(), posAtMid.y());
0275             }
0276             if (detailedAnalysis)
0277               histos[2]->Fill(dist);
0278           }
0279         }
0280       }
0281     }
0282   }  // loop over stations
0283 }
0284 
0285 // requirements : max number of hits and min chi2
0286 const DTRecSegment4D& DTChamberEfficiencyTask::getBestSegment(const DTRecSegment4DCollection::range& segs) const {
0287   DTRecSegment4DCollection::const_iterator bestIter;
0288   unsigned int nHitBest = 0;
0289   double chi2Best = 99999.;
0290   for (DTRecSegment4DCollection::const_iterator seg = segs.first; seg != segs.second; ++seg) {
0291     unsigned int nHits = ((*seg).hasPhi() ? (*seg).phiSegment()->recHits().size() : 0);
0292     nHits += ((*seg).hasZed() ? (*seg).zSegment()->recHits().size() : 0);
0293 
0294     if (nHits == nHitBest) {
0295       if ((*seg).chi2() / (*seg).degreesOfFreedom() < chi2Best) {
0296         chi2Best = (*seg).chi2() / (*seg).degreesOfFreedom();
0297         bestIter = seg;
0298       }
0299     } else if (nHits > nHitBest) {
0300       nHitBest = nHits;
0301       bestIter = seg;
0302     }
0303   }
0304   return *bestIter;
0305 }
0306 
0307 const DTRecSegment4D* DTChamberEfficiencyTask::getBestSegment(const DTRecSegment4D* s1,
0308                                                               const DTRecSegment4D* s2) const {
0309   if (!s1)
0310     return s2;
0311   if (!s2)
0312     return s1;
0313   unsigned int nHits1 = (s1->hasPhi() ? s1->phiSegment()->recHits().size() : 0);
0314   nHits1 += (s1->hasZed() ? s1->zSegment()->recHits().size() : 0);
0315 
0316   unsigned int nHits2 = (s2->hasPhi() ? s2->phiSegment()->recHits().size() : 0);
0317   nHits2 += (s2->hasZed() ? s2->zSegment()->recHits().size() : 0);
0318 
0319   if (nHits1 == nHits2) {
0320     if (s1->chi2() / s1->degreesOfFreedom() < s2->chi2() / s2->degreesOfFreedom())
0321       return s1;
0322     else
0323       return s2;
0324   } else if (nHits1 > nHits2)
0325     return s1;
0326   return s2;
0327 }
0328 
0329 LocalPoint DTChamberEfficiencyTask::interpolate(const DTRecSegment4D& seg1,
0330                                                 const DTRecSegment4D& seg3,
0331                                                 const DTChamberId& id2) const {
0332   // Get GlobalPoition of Seg in MB1
0333   GlobalPoint gpos1 = (dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());
0334 
0335   // Get GlobalPoition of Seg in MB3
0336   GlobalPoint gpos3 = (dtGeom->chamber(seg3.chamberId()))->toGlobal(seg3.localPosition());
0337 
0338   // interpolate
0339   // get all in MB2 frame
0340   LocalPoint pos1 = (dtGeom->chamber(id2))->toLocal(gpos1);
0341   LocalPoint pos3 = (dtGeom->chamber(id2))->toLocal(gpos3);
0342 
0343   // case 1: 1 and 3 has both projection. No problem
0344 
0345   // case 2: one projection is missing for one of the segments. Keep the other's segment position
0346   if (!seg1.hasZed())
0347     pos1 = LocalPoint(pos1.x(), pos3.y(), pos1.z());
0348   if (!seg3.hasZed())
0349     pos3 = LocalPoint(pos3.x(), pos1.y(), pos3.z());
0350 
0351   if (!seg1.hasPhi())
0352     pos1 = LocalPoint(pos3.x(), pos1.y(), pos1.z());
0353   if (!seg3.hasPhi())
0354     pos3 = LocalPoint(pos1.x(), pos3.y(), pos3.z());
0355 
0356   // direction
0357   LocalVector dir = (pos3 - pos1).unit();  // z points inward!
0358   LocalPoint pos2 = pos1 + dir * pos1.z() / (-dir.z());
0359 
0360   return pos2;
0361 }
0362 
0363 bool DTChamberEfficiencyTask::isGoodSegment(const DTRecSegment4D& seg) const {
0364   if (seg.chamberId().station() != 4 && !seg.hasZed())
0365     return false;
0366   unsigned int nHits = (seg.hasPhi() ? seg.phiSegment()->recHits().size() : 0);
0367   nHits += (seg.hasZed() ? seg.zSegment()->recHits().size() : 0);
0368   return (nHits >= theMinHitsSegment && seg.chi2() / seg.degreesOfFreedom() < theMinChi2NormSegment);
0369 }
0370 
0371 // Local Variables:
0372 // show-trailing-whitespace: t
0373 // truncate-lines: t
0374 // End:
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!