DTSegment/test/DTEffAnalyzer.cc

0001 /******* \class DTEffAnalyzer *******
0002  *
0003  * Description:
0004  *
0005  *  detailed description
0006  *
0007  * \author : Stefano Lacaprara - INFN LNL <stefano.lacaprara@pd.infn.it>
0008  *
0009  * Modification:
0010  *
0011  *********************************/
0012
0013 /* This Class Header */
0014 #include "RecoLocalMuon/DTSegment/test/DTEffAnalyzer.h"
0015
0016 /* Collaborating Class Header */
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023
0024 using namespace edm;
0025
0026 #include "TFile.h"
0027 #include "TH1F.h"
0028 #include "TH2F.h"
0029
0030 #include "DataFormats/LTCDigi/interface/LTCDigi.h"
0031
0032 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0033 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0034 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0035 #include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
0036
0037 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0038 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0039
0040 /* C++ Headers */
0041 #include <iostream>
0042 #include <cmath>
0043 using namespace std;
0044
0045 /* ====================================================================== */
0046
0047 /* Constructor */
0048 DTEffAnalyzer::DTEffAnalyzer(const ParameterSet& pset) : theDtGeomToken(esConsumes()) {
0049   // Get the debug parameter for verbose output
0050   debug = pset.getUntrackedParameter<bool>("debug");
0051   theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0052
0053   // the name of the 1D rec hits collection
0054   theRecHits1DLabel = pset.getParameter<string>("recHits1DLabel");
0055
0056   // the name of the 2D rec hits collection
0057   theRecHits2DLabel = pset.getParameter<string>("recHits2DLabel");
0058
0059   // the name of the 4D rec hits collection
0060   theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");
0061
0062   theMinHitsSegment = static_cast<unsigned int>(pset.getParameter<int>("minHitsSegment"));
0063   theMinChi2NormSegment = pset.getParameter<double>("minChi2NormSegment");
0064   theMinCloseDist = pset.getParameter<double>("minCloseDist");
0065
0066   consumes<DTRecSegment4DCollection>(theRecHits4DLabel);
0067 }
0068
0069 void DTEffAnalyzer::beginJob() {
0070   if (debug)
0071     cout << "beginOfJob" << endl;
0072
0073   // Create the root file
0074   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0075   bool dirStat = TH1::AddDirectoryStatus();
0076   TH1::AddDirectory(kTRUE);
0077
0078   // trigger Histos
0079   new TH1F("hTrigBits", "All trigger bits", 10, 0., 10.);
0080
0081   for (int w = -2; w <= 2; ++w) {
0082     stringstream nameWheel;
0083     nameWheel << "_Wh" << w;
0084     //cout << "Wheel " << nameWheel.str() << endl;
0085     for (int sec = 1; sec <= 14; ++sec) {  // section 1 to 14
0086       stringstream nameSector;
0087       nameSector << nameWheel.str() << "_Sec" << sec;
0088       //cout << "Sec " << nameSector.str() << endl;
0089       for (int st = 1; st <= 4; ++st) {  // station 1 to 4
0090
0091         stringstream nameChamber;
0092         nameChamber << nameSector.str() << "_St" << st;
0093
0094         //cout << nameChamber << endl;
0095         createTH1F("hDistSegFromExtrap", "Distance segments from extrap position ", nameChamber.str(), 200, 0., 200.);
0096         createTH1F("hNaiveEffSeg", "Naive eff ", nameChamber.str(), 10, 0., 10.);
0097         createTH2F(
0098             "hEffSegVsPosDen", "Eff vs local position (all) ", nameChamber.str(), 25, -250., 250., 25, -250., 250.);
0099         createTH2F(
0100             "hEffGoodSegVsPosDen", "Eff vs local position (good) ", nameChamber.str(), 25, -250., 250., 25, -250., 250.);
0101         createTH2F("hEffSegVsPosNum", "Eff vs local position ", nameChamber.str(), 25, -250., 250., 25, -250., 250.);
0102         createTH2F("hEffGoodSegVsPosNum",
0103                    "Eff vs local position (good segs) ",
0104                    nameChamber.str(),
0105                    25,
0106                    -250.,
0107                    250.,
0108                    25,
0109                    -250.,
0110                    250.);
0111         createTH2F("hEffGoodCloseSegVsPosNum",
0112                    "Eff vs local position (good aand close segs) ",
0113                    nameChamber.str(),
0114                    25,
0115                    -250.,
0116                    250.,
0117                    25,
0118                    -250.,
0119                    250.);
0120       }
0121     }
0122   }
0123   // cout << "List of created histograms " << endl;
0124   // theFile->ls();
0125   TH1::AddDirectory(dirStat);
0126 }
0127
0128 /* Destructor */
0129 DTEffAnalyzer::~DTEffAnalyzer() {
0130   theFile->cd();
0131   theFile->Write();
0132   theFile->Close();
0133 }
0134
0135 /* Operations */
0136 void DTEffAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0137   dtGeom = eventSetup.getHandle(theDtGeomToken);
0138
0139   if (debug)
0140     cout << endl
0141          << "--- [DTEffAnalyzer] Event analysed #Run: " << event.id().run() << " #Event: " << event.id().event()
0142          << endl;
0143
0144   effSegments(event, eventSetup);
0145 }
0146
0147 void DTEffAnalyzer::effSegments(const Event& event, const EventSetup& eventSetup) {
0148   // Get the 4D rechit collection from the event -------------------
0149   // Handle<DTRecSegment4DCollection> segs;
0150   event.getByLabel(theRecHits4DLabel, segs);
0151   if (debug) {
0152     cout << "4d " << segs->size() << endl;
0153     for (DTRecSegment4DCollection::const_iterator seg = segs->begin(); seg != segs->end(); ++seg)
0154       cout << *seg << endl;
0155   }
0156
0157   // Get events with 3 segments in different station and look what happen on
0158   // the other station. Note, must take into account geometrical acceptance
0159
0160   // trivial pattern recognition: get 3 segments in 3 different station of a
0161   // given wheel, sector
0162
0163   for (int wheel = -2; wheel <= 2; ++wheel) {
0164     for (int sector = 1; sector <= 12; ++sector) {
0165       evaluateEff(DTChamberId(wheel, 1, sector), 2, 3);  // get efficiency for MB1 using MB2 and MB3
0166       evaluateEff(DTChamberId(wheel, 2, sector), 1, 3);  // get efficiency for MB2 using MB1 and MB3
0167       evaluateEff(DTChamberId(wheel, 3, sector), 2, 4);  // get efficiency for MB3 using MB2 and MB4
0168       evaluateEff(DTChamberId(wheel, 4, sector), 2, 3);  // get efficiency for MB4 using MB2 and MB3
0169     }
0170   }
0171 }
0172
0173 void DTEffAnalyzer::evaluateEff(const DTChamberId& MidId, int bottom, int top) const {
0174   if (debug)
0175     cout << "evaluateEff " << MidId << " bott/top " << bottom << "/" << top << endl;
0176   // Select events with (good) segments in Bot and Top
0177   DTChamberId BotId(MidId.wheel(), bottom, MidId.sector());
0178   DTChamberId TopId(MidId.wheel(), top, MidId.sector());
0179
0180   // Get segments in the bottom chambers (if any)
0181   DTRecSegment4DCollection::range segsBot = segs->get(BotId);
0182   int nSegsBot = segsBot.second - segsBot.first;
0183   // check if any segments is there
0184   if (nSegsBot == 0)
0185     return;
0186
0187   // Get segments in the top chambers (if any)
0188   DTRecSegment4DCollection::range segsTop = segs->get(TopId);
0189   int nSegsTop = segsTop.second - segsTop.first;
0190
0191   // something more sophisticate check quality of segments
0192   const DTRecSegment4D& bestBotSeg = getBestSegment(segsBot);
0193   //cout << "BestBotSeg " << bestBotSeg << endl;
0194
0195   DTRecSegment4D* pBestTopSeg = 0;
0196   if (nSegsTop > 0)
0197     pBestTopSeg = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop));
0198   //if top chamber is MB4 sector 10, consider also sector 14
0199   if (TopId.station() == 4 && TopId.sector() == 10) {
0200     // cout << "nSegsTop " << nSegsTop << endl;
0201     // cout << "pBestTopSeg " << pBestTopSeg << endl;
0202     // cout << "Ch " << TopId << endl;
0203     DTChamberId TopId14(MidId.wheel(), top, 14);
0204     DTRecSegment4DCollection::range segsTop14 = segs->get(TopId14);
0205     int nSegsTop14 = segsTop14.second - segsTop14.first;
0206     //cout << "nSegsTop14 " << nSegsTop14 << endl;
0207     nSegsTop += nSegsTop;
0208     if (nSegsTop14) {
0209       DTRecSegment4D* pBestTopSeg14 = const_cast<DTRecSegment4D*>(&getBestSegment(segsTop14));
0210
0211       // get best between sector 10 and 14
0212       pBestTopSeg = const_cast<DTRecSegment4D*>(getBestSegment(pBestTopSeg, pBestTopSeg14));
0213       //cout << "pBestTopSeg " << pBestTopSeg << endl;
0214     }
0215   }
0216   if (!pBestTopSeg)
0217     return;
0218   const DTRecSegment4D& bestTopSeg = *pBestTopSeg;
0219   //cout << "BestTopSeg " << bestTopSeg << endl;
0220
0221   DTRecSegment4DCollection::range segsMid = segs->get(MidId);
0222   int nSegsMid = segsMid.second - segsMid.first;
0223   //cout << "nSegsMid " << nSegsMid << endl;
0224
0225   // very trivial efficiency, just count segments
0226   // cout << "MidId " << MidId << endl;
0227   // cout << "histo " << hName("hNaiveEffSeg",MidId) << endl;
0228   // cout << histo(hName("hNaiveEffSeg",MidId)) << endl;
0229   histo(hName("hNaiveEffSeg", MidId))->Fill(0);
0230   if (nSegsMid > 0)
0231     histo(hName("hNaiveEffSeg", MidId))->Fill(1);
0232
0233   // get position at Mid by interpolating the position (not direction) of best
0234   // segment in Bot and Top to Mid surface
0235   LocalPoint posAtMid = interpolate(bestBotSeg, bestTopSeg, MidId);
0236   // cout << "PosAtMid " << posAtMid << endl;
0237
0238   // is best segment good enough?
0239   //cout << "about to good " << endl;
0240   if (isGoodSegment(bestBotSeg) && isGoodSegment(bestTopSeg)) {
0241     histo2d(hName("hEffSegVsPosDen", MidId))->Fill(posAtMid.x(), posAtMid.y());
0242     //check if interpolation fall inside middle chamber
0243     if ((dtGeom->chamber(MidId))->surface().bounds().inside(posAtMid)) {
0244       // cout << "IsInside" << endl;
0245
0246       //cout << "good" << endl;
0247       histo2d(hName("hEffGoodSegVsPosDen", MidId))->Fill(posAtMid.x(), posAtMid.y());
0248
0249       if (nSegsMid > 0) {
0250         histo(hName("hNaiveEffSeg", MidId))->Fill(2);
0251         histo2d(hName("hEffSegVsPosNum", MidId))->Fill(posAtMid.x(), posAtMid.y());
0252         const DTRecSegment4D& bestMidSeg = getBestSegment(segsMid);
0253         // check if middle segments is good enough
0254         if (isGoodSegment(bestMidSeg)) {
0255           histo2d(hName("hEffGoodSegVsPosNum", MidId))->Fill(posAtMid.x(), posAtMid.y());
0256           LocalPoint midSegPos = bestMidSeg.localPosition();
0257           // check if middle segments is also close enough
0258           double dist;
0259           // cout << "bestBotSeg " << bestBotSeg.hasPhi() << " " <<
0260           //    bestBotSeg.hasZed() << " " << bestBotSeg << endl;
0261           // cout << "bestTopSeg " << bestTopSeg.hasPhi() << " " <<
0262           //   bestTopSeg.hasZed() << " " << bestTopSeg << endl;
0263           // cout << "midSegPos " << midSegPos << endl;
0264           // cout << "posAtMid " << posAtMid<< endl;
0265           // cout << "bestMidSeg " << bestMidSeg.hasPhi() << " " <<
0266           //   bestMidSeg.hasZed() << " " << bestMidSeg << endl;
0267           if (bestMidSeg.hasPhi()) {
0268             if (bestTopSeg.hasZed() && bestBotSeg.hasZed() && bestMidSeg.hasZed()) {
0269               dist = (midSegPos - posAtMid).mag();
0270             } else {
0271               dist = fabs((midSegPos - posAtMid).x());
0272             }
0273           } else {
0274             dist = fabs((midSegPos - posAtMid).y());
0275           }
0276           // cout << "dist " << dist << " theMinCloseDist " << theMinCloseDist<< endl;
0277           if (dist < theMinCloseDist) {
0278             histo2d(hName("hEffGoodCloseSegVsPosNum", MidId))->Fill(posAtMid.x(), posAtMid.y());
0279           }
0280           histo(hName("hDistSegFromExtrap", MidId))->Fill(dist);
0281         }
0282       }
0283     }
0284   }
0285   // else cout << "Outside " << endl;
0286 }
0287
0288 // as usual max number of hits and min chi2
0289 const DTRecSegment4D& DTEffAnalyzer::getBestSegment(const DTRecSegment4DCollection::range& segs) const {
0290   DTRecSegment4DCollection::const_iterator bestIter;
0291   unsigned int nHitBest = 0;
0292   double chi2Best = 99999.;
0293   for (DTRecSegment4DCollection::const_iterator seg = segs.first; seg != segs.second; ++seg) {
0294     unsigned int nHits = ((*seg).hasPhi() ? (*seg).phiSegment()->recHits().size() : 0);
0295     nHits += ((*seg).hasZed() ? (*seg).zSegment()->recHits().size() : 0);
0296
0297     if (nHits == nHitBest) {
0298       if ((*seg).chi2() / (*seg).degreesOfFreedom() < chi2Best) {
0299         chi2Best = (*seg).chi2() / (*seg).degreesOfFreedom();
0300         bestIter = seg;
0301       }
0302     } else if (nHits > nHitBest) {
0303       nHitBest = nHits;
0304       bestIter = seg;
0305     }
0306   }
0307   return *bestIter;
0308 }
0309
0310 const DTRecSegment4D* DTEffAnalyzer::getBestSegment(const DTRecSegment4D* s1, const DTRecSegment4D* s2) const {
0311   if (!s1)
0312     return s2;
0313   if (!s2)
0314     return s1;
0315   unsigned int nHits1 = (s1->hasPhi() ? s1->phiSegment()->recHits().size() : 0);
0316   nHits1 += (s1->hasZed() ? s1->zSegment()->recHits().size() : 0);
0317
0318   unsigned int nHits2 = (s2->hasPhi() ? s2->phiSegment()->recHits().size() : 0);
0319   nHits2 += (s2->hasZed() ? s2->zSegment()->recHits().size() : 0);
0320
0321   if (nHits1 == nHits2) {
0322     if (s1->chi2() / s1->degreesOfFreedom() < s2->chi2() / s2->degreesOfFreedom())
0323       return s1;
0324     else
0325       return s2;
0326   } else if (nHits1 > nHits2)
0327     return s1;
0328   return s2;
0329 }
0330
0331 bool DTEffAnalyzer::isGoodSegment(const DTRecSegment4D& seg) const {
0332   if (seg.chamberId().station() != 4 && !seg.hasZed())
0333     return false;
0334   unsigned int nHits = (seg.hasPhi() ? seg.phiSegment()->recHits().size() : 0);
0335   nHits += (seg.hasZed() ? seg.zSegment()->recHits().size() : 0);
0336   return (nHits >= theMinHitsSegment && seg.chi2() / seg.degreesOfFreedom() < theMinChi2NormSegment);
0337 }
0338
0339 LocalPoint DTEffAnalyzer::interpolate(const DTRecSegment4D& seg1,
0340                                       const DTRecSegment4D& seg3,
0341                                       const DTChamberId& id2) const {
0342   // Get GlobalPoition of Seg in MB1
0343   GlobalPoint gpos1 = (dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());
0344
0345   // Get GlobalPoition of Seg in MB3
0346   GlobalPoint gpos3 = (dtGeom->chamber(seg3.chamberId()))->toGlobal(seg3.localPosition());
0347
0348   // interpolate
0349   // get all in MB2 frame
0350   LocalPoint pos1 = (dtGeom->chamber(id2))->toLocal(gpos1);
0351   LocalPoint pos3 = (dtGeom->chamber(id2))->toLocal(gpos3);
0352   // cout << "pos1 " << pos1 << endl;
0353   // cout << "pos3 " << pos3 << endl;
0354
0355   // case 1: 1 and 3 has both projection. No problem
0356
0357   // case 2: one projection is missing for one of the segments. Keep the other's
0358   // segment position
0359   if (!seg1.hasZed())
0360     pos1 = LocalPoint(pos1.x(), pos3.y(), pos1.z());
0361   if (!seg3.hasZed())
0362     pos3 = LocalPoint(pos3.x(), pos1.y(), pos3.z());
0363
0364   if (!seg1.hasPhi())
0365     pos1 = LocalPoint(pos3.x(), pos1.y(), pos1.z());
0366   if (!seg3.hasPhi())
0367     pos3 = LocalPoint(pos1.x(), pos3.y(), pos3.z());
0368
0369   // cout << "pos1 " << pos1 << endl;
0370   // cout << "pos3 " << pos3 << endl;
0371   // direction
0372   LocalVector dir = (pos3 - pos1).unit();  // z points inward!
0373   // cout << "dir " << dir << endl;
0374   LocalPoint pos2 = pos1 + dir * pos1.z() / (-dir.z());
0375   // cout << "pos2 " << pos2 << endl;
0376
0377   return pos2;
0378 }
0379
0380 TH1F* DTEffAnalyzer::histo(const string& name) const {
0381   if (!theFile->Get(name.c_str()))
0382     throw cms::Exception("DTEffAnalyzer") << " TH1F not existing " << name;
0383   if (TH1F* h = dynamic_cast<TH1F*>(theFile->Get(name.c_str())))
0384     return h;
0385   else
0386     throw cms::Exception("DTEffAnalyzer") << " Not a TH1F " << name;
0387 }
0388
0389 TH2F* DTEffAnalyzer::histo2d(const string& name) const {
0390   if (!theFile->Get(name.c_str()))
0391     throw cms::Exception("DTEffAnalyzer") << " TH1F not existing " << name;
0392   if (TH2F* h = dynamic_cast<TH2F*>(theFile->Get(name.c_str())))
0393     return h;
0394   else
0395     throw cms::Exception("DTEffAnalyzer") << " Not a TH2F " << name;
0396 }
0397
0398 string DTEffAnalyzer::toString(const DTChamberId& id) const {
0399   stringstream result;
0400   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station();
0401   return result.str();
0402 }
0403
0404 template <class T>
0405 string DTEffAnalyzer::hName(const string& s, const T& id) const {
0406   string name(toString(id));
0407   stringstream hName;
0408   hName << s << name;
0409   return hName.str();
0410 }
0411
0412 void DTEffAnalyzer::createTH1F(const std::string& name,
0413                                const std::string& title,
0414                                const std::string& suffix,
0415                                int nbin,
0416                                const double& binMin,
0417                                const double& binMax) const {
0418   stringstream hName;
0419   stringstream hTitle;
0420   hName << name << suffix;
0421   hTitle << title << suffix;
0422   new TH1F(hName.str().c_str(), hTitle.str().c_str(), nbin, binMin, binMax);
0423 }
0424
0425 void DTEffAnalyzer::createTH2F(const std::string& name,
0426                                const std::string& title,
0427                                const std::string& suffix,
0428                                int nBinX,
0429                                const double& binXMin,
0430                                const double& binXMax,
0431                                int nBinY,
0432                                const double& binYMin,
0433                                const double& binYMax) const {
0434   stringstream hName;
0435   stringstream hTitle;
0436   hName << name << suffix;
0437   hTitle << title << suffix;
0438   new TH2F(hName.str().c_str(), hTitle.str().c_str(), nBinX, binXMin, binXMax, nBinY, binYMin, binYMax);
0439 }
0440
0441 DEFINE_FWK_MODULE(DTEffAnalyzer);