DTSegment/test/DTSegAnalyzer.cc

0001 /******* \class DTSegAnalyzer *******
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/DTSegAnalyzer.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/Framework/interface/ConsumesCollector.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024
0025 using namespace edm;
0026
0027 #include "TFile.h"
0028 #include "TH1F.h"
0029 #include "TH2F.h"
0030
0031 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0032 #include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
0033 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0034 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0035 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0036 #include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
0037 #include "RecoLocalMuon/DTSegment/test/DTMeanTimer.h"
0038 #include "RecoLocalMuon/DTSegment/test/DTSegmentResidual.h"
0039
0040 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0041 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0045 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0046
0047 #include "DataFormats/DTDigi/interface/DTLocalTriggerCollection.h"
0048 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0049 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0050
0051 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0052 #include "DataFormats/TrackReco/interface/Track.h"
0053
0054 /* C++ Headers */
0055 #include <iostream>
0056 #include <cmath>
0057 using namespace std;
0058
0059 /* ====================================================================== */
0060
0061 /* Constructor */
0062 DTSegAnalyzer::DTSegAnalyzer(const ParameterSet& pset)
0063     : _ev(0),
0064       theSync{DTTTrigSyncFactory::get()->create(pset.getUntrackedParameter<string>("tTrigMode"),
0065                                                 pset.getUntrackedParameter<ParameterSet>("tTrigModeConfig"),
0066                                                 consumesCollector())},
0067       theDTGeomToken{esConsumes()} {
0068   // Get the debug parameter for verbose output
0069   debug = pset.getUntrackedParameter<bool>("debug");
0070   theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0071
0072   // the name of the 1D rec hits collection
0073   theRecHits1DLabel = pset.getParameter<string>("recHits1DLabel");
0074
0075   // the name of the 2D rec hits collection
0076   theRecHits2DLabel = pset.getParameter<string>("recHits2DLabel");
0077
0078   // the name of the 4D rec hits collection
0079   theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");
0080
0081   doHits = pset.getParameter<bool>("doHits");
0082   doSegs = pset.getParameter<bool>("doSegs");
0083
0084   // Create the root file
0085   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0086   bool dirStat = TH1::AddDirectoryStatus();
0087   TH1::AddDirectory(kTRUE);
0088
0089   /// DT histos
0090   // 1d hits
0091   new TH1F("hnHitDT", "Num 1d hits DT", 200, 0., 200.);
0092   new TH1F("hDigiTime", "Digi time (ns)", 700, -100., 600.);
0093   // new TH1F("hMeanTimerSL","Mean timer per SL (ns)",100,200.,600.);
0094   // new TH1F("hDigiTimeSL", "Digi Time (ns) ",100,-100.,600.);
0095
0096   new TH1F("hPosLeft", "Pos of Left hit (cm) in local frame", 100, -220., 220.);
0097   new TH1F("hPosRight", "Pos of Right hit (cm) in local frame", 100, -220., 220.);
0098   new TH1F("hMeanTimer", "Tmax ", 100, 200., 600.);
0099   new TH1F("hMeanTimerSeg", "Tmax from segments hits ", 100, 200., 600.);
0100   new TH2F("hMeanTimerSegAlongWire", "Tmax from segments hits vs pos along wire ", 40, -150., 150, 100, 200., 600.);
0101   new TH2F("hMeanTimerSegVsAngle", "Tmax from segments hits vs angle ", 40, -0., 1.0, 100, 200., 600.);
0102   new TH2F("hMeanTimerSegVsNHits", "Tmax from segments hits vs n hits segment ", 10, 0.5, 10.5, 100, 200., 600.);
0103
0104   // segs
0105   new TH1F("hnSegDT", "Num seg DT", 50, 0., 50.);
0106   new TH1F("hNSegs", "N segments ", 20, 0., 20.);
0107   new TH1F("hChi2Seg", "Chi2 segments ", 25, 0., 25.);
0108   new TH1F("hChi2SegPhi", "Chi2 segments phi ", 25, 0., 25.);
0109   new TH1F("hChi2SegZed", "Chi2 segments zed ", 25, 0., 25.);
0110   new TH1F("hNHitsSeg", "N hits segment ", 15, -0.5, 14.5);
0111   new TH1F("hNHitsSegPhi", "N hits segment phi ", 12, -0.5, 11.5);
0112   new TH2F("hNHitsSegPhiVsAngle", "N hits segment phi vs angle ", 40, .0, 1.0, 12, -0.5, 11.5);
0113   new TH2F("hNHitsSegPhiVsOtherHits", "N hits segment phi vs hits in SLs ", 20, -0.5, 19.5, 12, -0.5, 11.5);
0114   new TH2F("hNHitsSegPhiVsNumSegs", "N hits segment zed vs num segs in ch ", 6, -0.5, 5.5, 12, -0.5, 11.5);
0115   new TH2F("hNHitsSegPhiVsChi2", "N hits segment zed vs chi2/NDoF ", 25, 0., 25., 12, -0.5, 11.5);
0116   new TH1F("hNHitsSegZed", "N hits segment zed ", 10, -0.5, 9.5);
0117   new TH2F("hNHitsSegZedVsAngle", "N hits segment zed vs angle ", 40, .0, 1.0, 12, -0.5, 11.5);
0118   new TH2F("hNHitsSegZedVsOtherHits", "N hits segment zed vs hits in SLs ", 20, -0.5, 19.5, 12, -0.5, 11.5);
0119   new TH2F("hNHitsSegZedVsNumSegs", "N hits segment zed vs num segs in ch ", 6, -0.5, 5.5, 12, -0.5, 11.5);
0120   new TH2F("hNHitsSegZedVsChi2", "N hits segment zed vs chi2/NDoF ", 25, 0., 25., 12, -0.5, 11.5);
0121
0122   new TH1F("hHitResidualSeg", "Hits residual wrt segments ", 100, -.2, +.2);
0123   new TH1F("hHitResidualSegCellDX", "Hits residual wrt segments semicell DX ", 100, -.2, +.2);
0124   new TH1F("hHitResidualSegCellSX", "Hits residual wrt segments semicell SX", 100, -.2, +.2);
0125   new TH2F("hHitResidualSegAlongWire", "Hits residual wrt segments vs pos along wire ", 40, -150., 150, 100, -.2, +.2);
0126   new TH2F("hHitResidualSegVsWireDis", "Hits residual wrt segments vs wire distance ", 40, 0., 2.1, 100, -.2, +.2);
0127   new TH2F("hHitResidualSegVsAngle", "Hits residual wrt segments vs impact angle ", 40, .0, 1.0, 100, -.2, +.2);
0128   new TH2F("hHitResidualSegVsNHits", "Hits residual wrt segments vs num hits ", 10, -0.5, 9.5, 100, -.2, +.2);
0129   new TH2F("hHitResidualSegVsChi2", "Hits residual wrt segments vs chi2 ", 25, .0, 25.0, 100, -.2, +.2);
0130   new TH2F("hNsegs2dVsNhits", "N segs 2d vs n Hits ", 20, -0.5, 19.5, 5, -0.5, 4.5);
0131
0132   // eff
0133   new TH2F("hNsegs4dVsNhits", "N segs 4d vs n Hits ", 20, 0.5, 20.5, 5, -0.5, 4.5);
0134   new TH2F("hNsegs4dVsNhitsPhi", "N segs 4d vs n HitsPhi ", 20, 0.5, 20.5, 5, -0.5, 4.5);
0135   new TH2F("hNsegs4dVsNhitsZed", "N segs 4d vs n HitsZed ", 20, 0.5, 20.5, 5, -0.5, 4.5);
0136
0137   new TH2F("hNsegs4dVsNsegs2d", "N segs 4d vs n segs2d ", 4, 0.5, 4.5, 5, -0.5, 4.5);
0138   new TH2F("hNsegs4dVsNsegs2dPhi", "N segs 4d vs n segs2d Phi ", 4, 0.5, 4.5, 5, -0.5, 4.5);
0139   new TH2F("hNsegs4dVsNsegs2dZed", "N segs 4d vs n segs2d Zed ", 4, 0.5, 4.5, 5, -0.5, 4.5);
0140
0141   new TH2F("hNsegs2dSL1VsNsegs2dSL3", "N segs 2d SL1 vs SL3 ", 5, -0.5, 4.5, 5, -0.5, 4.5);
0142   new TH2F("hNsegs2dSL1VsNsegs2dSL2", "N segs 2d SL1 vs SL2 ", 5, -0.5, 4.5, 5, -0.5, 4.5);
0143   new TH2F("hNsegs2dSL2VsNsegs2dSL3", "N segs 2d SL2 vs SL3 ", 5, -0.5, 4.5, 5, -0.5, 4.5);
0144
0145   TH1::AddDirectory(dirStat);
0146 }
0147
0148 /* Destructor */
0149 DTSegAnalyzer::~DTSegAnalyzer() {
0150   theFile->cd();
0151   theFile->Write();
0152   theFile->Close();
0153 }
0154
0155 /* Operations */
0156 void DTSegAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0157   theSync->setES(eventSetup);
0158   _ev++;
0159
0160   if (debug)
0161     cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() << " Num " << _ev << endl;
0162
0163   static int j = 1;
0164   if ((_ev % j) == 0) {
0165     if ((_ev / j) == 9)
0166       j *= 10;
0167     cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() << " Num " << _ev << endl;
0168   }
0169
0170   if (doHits)
0171     analyzeDTHits(event, eventSetup);
0172   if (doSegs)
0173     analyzeDTSegments(event, eventSetup);
0174 }
0175
0176 void DTSegAnalyzer::analyzeDTHits(const Event& event, const EventSetup& eventSetup) {
0177   // Get the DT Geometry
0178   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(theDTGeomToken);
0179
0180   // Get the 1D rechits from the event --------------
0181   Handle<DTRecHitCollection> dtRecHits;
0182   event.getByLabel(theRecHits1DLabel, dtRecHits);
0183
0184   int nHitDT = dtRecHits->size();
0185   histo("hnHitDT")->Fill(nHitDT);
0186
0187   //float ttrigg = 1895.; // should get this from CondDB...
0188   for (DTRecHitCollection::const_iterator hit = dtRecHits->begin(); hit != dtRecHits->end(); ++hit) {
0189     // Get the wireId of the rechit
0190     DTWireId wireId = (*hit).wireId();
0191
0192     float ttrig = theSync->offset(wireId);
0193     //cout << "TTrig " << ttrig << endl;
0194
0195     float time = (*hit).digiTime() - ttrig;
0196     double xLeft = (*hit).localPosition(DTEnums::Left).x();
0197     double xRight = (*hit).localPosition(DTEnums::Right).x();
0198
0199     histo("hDigiTime")->Fill(time);
0200     histo("hPosLeft")->Fill(xLeft);
0201     histo("hPosRight")->Fill(xRight);
0202   }
0203
0204   // MeanTimer analysis
0205   // loop on SL
0206   //cout << "MeanTimer analysis" << endl;
0207   const std::vector<const DTSuperLayer*>& sls = dtGeom->superLayers();
0208   for (auto sl = sls.begin(); sl != sls.end(); ++sl) {
0209     DTSuperLayerId slid = (*sl)->id();
0210
0211     DTMeanTimer meanTimer(dtGeom->superLayer(slid), dtRecHits, theSync.get());
0212     vector<double> tMaxs = meanTimer.run();
0213     for (vector<double>::const_iterator tMax = tMaxs.begin(); tMax != tMaxs.end(); ++tMax) {
0214       //cout << "Filling " << hName("hMeanTimer", slid) << " with " << *tMax << endl;
0215       histo("hMeanTimer")->Fill(*tMax);
0216     }
0217   }
0218 }
0219
0220 void DTSegAnalyzer::analyzeDTSegments(const Event& event, const EventSetup& eventSetup) {
0221   if (debug)
0222     cout << "analyzeDTSegments" << endl;
0223   // Get the DT Geometry
0224   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(theDTGeomToken);
0225
0226   // Get the 4D rechit collection from the event -------------------
0227   edm::Handle<DTRecSegment4DCollection> segs;
0228   event.getByLabel(theRecHits4DLabel, segs);
0229   if (debug)
0230     cout << "4d " << segs->size() << endl;
0231
0232   // Get the 2D rechit collection from the event -------------------
0233   edm::Handle<DTRecSegment2DCollection> segs2d;
0234   event.getByLabel(theRecHits2DLabel, segs2d);
0235   if (debug)
0236     cout << "2d " << segs2d->size() << endl;
0237
0238   // Get the 1D rechits from the event --------------
0239   Handle<DTRecHitCollection> dtRecHits;
0240   event.getByLabel(theRecHits1DLabel, dtRecHits);
0241   if (debug)
0242     cout << "1d " << dtRecHits->size() << endl;
0243
0244   int nsegs = segs->size();
0245   histo("hnSegDT")->Fill(nsegs);
0246   const std::vector<const DTChamber*>& chs = dtGeom->chambers();
0247
0248   for (auto ch = chs.begin(); ch != chs.end(); ++ch) {
0249     DTChamberId chid((*ch)->id());
0250     // Segment 4d in this chamber
0251     DTRecSegment4DCollection::range segsch = segs->get(chid);
0252     int nSegsCh = segsch.second - segsch.first;
0253     histo("hNSegs")->Fill(nSegsCh);
0254
0255     // some quality on segments
0256     //bool hasGoodSegment=false;
0257
0258     for (DTRecSegment4DCollection::const_iterator seg = segsch.first; seg != segsch.second; ++seg) {
0259       // first the the two projection separately
0260       // phi
0261       if (debug)
0262         cout << *seg << endl;
0263       const DTChamberRecSegment2D* phiSeg = (*seg).phiSegment();
0264
0265       vector<DTRecHit1D> phiHits;
0266       if (phiSeg) {
0267         if (debug)
0268           cout << "Phi " << *phiSeg << endl;
0269         phiHits = phiSeg->specificRecHits();
0270         //if (phiHits.size()>=6) hasGoodSegment = true;
0271
0272         DTSuperLayerId slid1(phiSeg->chamberId(), 1);
0273         /// Mean timer analysis
0274         DTMeanTimer meanTimer1(dtGeom->superLayer(slid1), phiHits, theSync.get());
0275         vector<double> tMaxs1 = meanTimer1.run();
0276         for (vector<double>::const_iterator tMax = tMaxs1.begin(); tMax != tMaxs1.end(); ++tMax) {
0277           histo("hMeanTimerSeg")->Fill(*tMax);
0278           histo2d("hMeanTimerSegVsNHits")->Fill(phiHits.size(), *tMax);
0279           if ((*seg).hasZed()) {
0280             histo2d("hMeanTimerSegAlongWire")->Fill((*seg).localPosition().y(), *tMax);
0281             histo2d("hMeanTimerSegVsAngle")->Fill(M_PI - (*seg).localDirection().theta(), *tMax);
0282           }
0283         }
0284
0285         DTSuperLayerId slid3(phiSeg->chamberId(), 3);
0286         /// Mean timer analysis
0287         DTMeanTimer meanTimer3(dtGeom->superLayer(slid3), phiHits, theSync.get());
0288         vector<double> tMaxs3 = meanTimer3.run();
0289         for (vector<double>::const_iterator tMax = tMaxs3.begin(); tMax != tMaxs3.end(); ++tMax) {
0290           histo("hMeanTimerSeg")->Fill(*tMax);
0291           histo2d("hMeanTimerSegVsNHits")->Fill(phiHits.size(), *tMax);
0292           if ((*seg).hasZed()) {
0293             histo2d("hMeanTimerSegAlongWire")->Fill((*seg).localPosition().y(), *tMax);
0294             histo2d("hMeanTimerSegVsAngle")->Fill(M_PI - (*seg).localDirection().theta(), *tMax);
0295           }
0296         }
0297       }
0298       // zed
0299       const DTSLRecSegment2D* zedSeg = (*seg).zSegment();
0300       //cout << "zedSeg " << zedSeg << endl;
0301       vector<DTRecHit1D> zedHits;
0302       if (zedSeg) {
0303         if (debug)
0304           cout << "Zed " << *zedSeg << endl;
0305         zedHits = zedSeg->specificRecHits();
0306         DTSuperLayerId slid(zedSeg->superLayerId());
0307         /// Mean timer analysis
0308         DTMeanTimer meanTimer(dtGeom->superLayer(slid), zedHits, theSync.get());
0309         vector<double> tMaxs = meanTimer.run();
0310         for (vector<double>::const_iterator tMax = tMaxs.begin(); tMax != tMaxs.end(); ++tMax) {
0311           histo("hMeanTimerSeg")->Fill(*tMax);
0312           histo2d("hMeanTimerSegVsNHits")->Fill(zedHits.size(), *tMax);
0313           if ((*seg).hasPhi()) {
0314             histo2d("hMeanTimerSegAlongWire")->Fill((*seg).localPosition().x(), *tMax);
0315             histo2d("hMeanTimerSegVsAngle")->Fill(M_PI - (*seg).localDirection().theta(), *tMax);
0316           }
0317         }
0318       }
0319
0320       histo("hChi2Seg")->Fill((*seg).chi2() / (*seg).degreesOfFreedom());
0321       if (phiSeg) {
0322         histo("hNHitsSegPhi")->Fill(phiSeg->recHits().size());
0323         histo2d("hNHitsSegPhiVsAngle")->Fill(M_PI - phiSeg->localDirection().theta(), phiSeg->recHits().size());
0324         // get the total numer of hits in this SL(s)
0325         DTSuperLayerId slid1(chid, 1);
0326         DTRecHitCollection::range rangeH1 = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid1));
0327         DTSuperLayerId slid3(chid, 3);
0328         DTRecHitCollection::range rangeH3 = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid3));
0329         histo2d("hNHitsSegPhiVsOtherHits")
0330             ->Fill((rangeH1.second - rangeH1.first) + (rangeH3.second - rangeH3.first), phiSeg->recHits().size());
0331         histo2d("hNHitsSegPhiVsNumSegs")->Fill(nSegsCh, phiSeg->recHits().size());
0332         histo2d("hNHitsSegPhiVsChi2")->Fill(phiSeg->chi2() / phiSeg->degreesOfFreedom(), phiSeg->recHits().size());
0333         histo("hChi2SegPhi")->Fill(phiSeg->chi2() / phiSeg->degreesOfFreedom());
0334       }
0335
0336       if (zedSeg) {
0337         histo("hNHitsSegZed")->Fill(zedSeg->recHits().size());
0338         histo2d("hNHitsSegZedVsAngle")->Fill(M_PI - zedSeg->localDirection().theta(), zedSeg->recHits().size());
0339         DTSuperLayerId slid(chid, 2);
0340         DTRecHitCollection::range rangeH = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid));
0341         histo2d("hNHitsSegZedVsOtherHits")->Fill(rangeH.second - rangeH.first, zedSeg->recHits().size());
0342         histo2d("hNHitsSegZedVsNumSegs")->Fill(nSegsCh, zedSeg->recHits().size());
0343         histo2d("hNHitsSegZedVsChi2")->Fill(zedSeg->chi2() / zedSeg->degreesOfFreedom(), zedSeg->recHits().size());
0344         histo("hChi2SegZed")->Fill(zedSeg->chi2() / zedSeg->degreesOfFreedom());
0345       }
0346       if (phiSeg && zedSeg)
0347         histo("hNHitsSeg")->Fill(phiSeg->recHits().size() + zedSeg->recHits().size());
0348
0349       // residual analysis
0350       if (phiSeg) {
0351         DTSegmentResidual res(phiSeg, *ch);
0352         res.run();
0353         vector<DTSegmentResidual::DTResidual> deltas = res.residuals();
0354         for (vector<DTSegmentResidual::DTResidual>::const_iterator delta = deltas.begin(); delta != deltas.end();
0355              ++delta) {
0356           histo("hHitResidualSeg")->Fill((*delta).value);
0357           if ((*delta).side == DTEnums::Right)
0358             histo("hHitResidualSegCellDX")->Fill((*delta).value);
0359           else if ((*delta).side == DTEnums::Left)
0360             histo("hHitResidualSegCellSX")->Fill((*delta).value);
0361
0362           histo2d("hHitResidualSegVsWireDis")->Fill((*delta).wireDistance, (*delta).value);
0363
0364           histo2d("hHitResidualSegVsAngle")->Fill((*delta).angle, (*delta).value);
0365           histo2d("hHitResidualSegVsNHits")->Fill(phiSeg->recHits().size(), (*delta).value);
0366           histo2d("hHitResidualSegVsChi2")->Fill(phiSeg->chi2(), (*delta).value);
0367
0368           if ((*seg).hasPhi())
0369             histo2d("hHitResidualSegAlongWire")->Fill((*seg).localPosition().x(), (*delta).value);
0370         }
0371       }
0372
0373       if (zedSeg) {
0374         const DTSuperLayer* sl = (*ch)->superLayer(2);
0375         DTSegmentResidual res(zedSeg, sl);
0376         res.run();
0377         vector<DTSegmentResidual::DTResidual> deltas = res.residuals();
0378         for (vector<DTSegmentResidual::DTResidual>::const_iterator delta = deltas.begin(); delta != deltas.end();
0379              ++delta) {
0380           histo("hHitResidualSeg")->Fill((*delta).value);
0381           if ((*delta).side == DTEnums::Right)
0382             histo("hHitResidualSegCellDX")->Fill((*delta).value);
0383           else if ((*delta).side == DTEnums::Left)
0384             histo("hHitResidualSegCellSX")->Fill((*delta).value);
0385
0386           histo2d("hHitResidualSegVsWireDis")->Fill((*delta).wireDistance, (*delta).value);
0387
0388           histo2d("hHitResidualSegVsAngle")->Fill((*delta).angle, (*delta).value);
0389           histo2d("hHitResidualSegVsNHits")->Fill(zedSeg->recHits().size(), (*delta).value);
0390           histo2d("hHitResidualSegVsChi2")->Fill(zedSeg->chi2(), (*delta).value);
0391
0392           if ((*seg).hasPhi())
0393             histo2d("hHitResidualSegAlongWire")->Fill((*seg).localPosition().x(), (*delta).value);
0394         }
0395       }  // if ZedSeg
0396
0397     }  // end loop segment in chamber
0398
0399     // Efficiency analysis
0400     // first get number of segments vs number of hits
0401
0402     // get No hits in this chamber
0403     // and no of 2d segments in chamber
0404     int nHitSl[3] = {0, 0, 0};
0405     int nSeg2dSl[3] = {0, 0, 0};
0406     for (int sl = 1; sl <= 3; ++sl) {
0407       if (sl == 2 && chid.station() == 4)
0408         continue;
0409       DTSuperLayerId slid(chid, sl);
0410       DTRecHitCollection::range rangeH = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid));
0411       DTRecSegment2DCollection::range rangeS = segs2d->get(slid);
0412       nHitSl[sl - 1] = rangeH.second - rangeH.first;
0413       nSeg2dSl[sl - 1] = rangeS.second - rangeS.first;
0414       histo2d("hNsegs2dVsNhits")->Fill(nHitSl[sl - 1], nSeg2dSl[sl - 1]);
0415     }
0416     histo2d("hNsegs4dVsNhits")->Fill(nHitSl[0] + nHitSl[1] + nHitSl[2], nSegsCh);
0417     histo2d("hNsegs4dVsNhitsPhi")->Fill(nHitSl[0] + nHitSl[2], nSegsCh);
0418     histo2d("hNsegs4dVsNhitsZed")->Fill(nHitSl[1], nSegsCh);
0419
0420     histo2d("hNsegs4dVsNsegs2d")->Fill(nSeg2dSl[0] + nSeg2dSl[1] + nSeg2dSl[2], nSegsCh);
0421     histo2d("hNsegs4dVsNsegs2dPhi")->Fill(nSeg2dSl[0] + nSeg2dSl[2], nSegsCh);
0422     histo2d("hNsegs4dVsNsegs2dZed")->Fill(nSeg2dSl[1], nSegsCh);
0423
0424     histo2d("hNsegs2dSL1VsNsegs2dSL3")->Fill(nSeg2dSl[2], nSeg2dSl[0]);
0425     histo2d("hNsegs2dSL1VsNsegs2dSL2")->Fill(nSeg2dSl[1], nSeg2dSl[0]);
0426     histo2d("hNsegs2dSL2VsNsegs2dSL3")->Fill(nSeg2dSl[2], nSeg2dSl[1]);
0427
0428   }  // loop on all chamber
0429 }
0430
0431 TH1F* DTSegAnalyzer::histo(const string& name) const {
0432   if (TH1F* h = dynamic_cast<TH1F*>(theFile->Get(name.c_str())))
0433     return h;
0434   else
0435     throw cms::Exception("DTSegAnalyzer") << " Not a TH1F " << name;
0436 }
0437
0438 TH2F* DTSegAnalyzer::histo2d(const string& name) const {
0439   if (TH2F* h = dynamic_cast<TH2F*>(theFile->Get(name.c_str())))
0440     return h;
0441   else
0442     throw cms::Exception("DTSegAnalyzer") << " Not a TH2F " << name;
0443 }
0444
0445 string DTSegAnalyzer::toString(const DTLayerId& id) const {
0446   stringstream result;
0447   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer() << "_Lay"
0448          << id.layer();
0449   return result.str();
0450 }
0451
0452 string DTSegAnalyzer::toString(const DTSuperLayerId& id) const {
0453   stringstream result;
0454   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer();
0455   return result.str();
0456 }
0457
0458 string DTSegAnalyzer::toString(const DTChamberId& id) const {
0459   stringstream result;
0460   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station();
0461   return result.str();
0462 }
0463
0464 template <class T>
0465 string DTSegAnalyzer::hName(const string& s, const T& id) const {
0466   string name(toString(id));
0467   stringstream hName;
0468   hName << s << name;
0469   return hName.str();
0470 }
0471
0472 void DTSegAnalyzer::createTH1F(const std::string& name,
0473                                const std::string& title,
0474                                const std::string& suffix,
0475                                int nbin,
0476                                const double& binMin,
0477                                const double& binMax) const {
0478   stringstream hName;
0479   stringstream hTitle;
0480   hName << name << suffix;
0481   hTitle << title << suffix;
0482   new TH1F(hName.str().c_str(), hTitle.str().c_str(), nbin, binMin, binMax);
0483 }
0484
0485 void DTSegAnalyzer::createTH2F(const std::string& name,
0486                                const std::string& title,
0487                                const std::string& suffix,
0488                                int nBinX,
0489                                const double& binXMin,
0490                                const double& binXMax,
0491                                int nBinY,
0492                                const double& binYMin,
0493                                const double& binYMax) const {
0494   stringstream hName;
0495   stringstream hTitle;
0496   hName << name << suffix;
0497   hTitle << title << suffix;
0498   new TH2F(hName.str().c_str(), hTitle.str().c_str(), nBinX, binXMin, binXMax, nBinY, binYMin, binYMax);
0499 }
0500
0501 DEFINE_FWK_MODULE(DTSegAnalyzer);