DTSegment/test/DTAnalyzerDetailed.cc

0001 /******* \class DTAnalyzerDetailed *******
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/DTAnalyzerDetailed.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 DTAnalyzerDetailed::DTAnalyzerDetailed(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
0094   for (int w = -2; w <= 2; ++w) {
0095     stringstream nameWheel;
0096     nameWheel << "_Wh" << w;
0097     //cout << "Wheel " << nameWheel.str() << endl;
0098     for (int sec = 1; sec <= 14; ++sec) {  // section 1 to 14
0099       stringstream nameSector;
0100       nameSector << nameWheel.str() << "_Sec" << sec;
0101       //cout << "Sec " << nameSector.str() << endl;
0102       for (int st = 1; st <= 4; ++st) {  // station 1 to 4
0103         stringstream nameChamber;
0104         nameChamber << nameSector.str() << "_St" << st;
0105         //cout << "Ch " << nameChamber.str() << endl;
0106         for (int sl = 1; sl <= 3; ++sl) {  // SL 1 to 3 (2 missing for st==4)
0107           if (st == 4 && sl == 2)
0108             continue;
0109           stringstream nameSL;
0110           nameSL << nameChamber.str() << "_Sl" << sl;
0111           //cout << "SL " << nameSL.str() << endl;
0112           for (int l = 1; l <= 4; ++l) {  // layer 1 to 4
0113             stringstream nameLayer;
0114             nameLayer << nameSL.str() << "_Lay" << l;
0115             //cout << "Lay " << nameLayer.str() << endl;
0116
0117             // Create hist for each Layer
0118             createTH1F("hDigiTime", "Digi Time (ns) ", nameLayer.str(), 100, -100., 600.);
0119
0120             createTH1F("hPosLeft", "Pos of Left hit (cm) in local frame", nameLayer.str(), 100, -220., 220.);
0121             createTH1F("hPosRight", "Pos of Right hit (cm) in local frame", nameLayer.str(), 100, -220., 220.);
0122           }
0123           // Create hist for each SuperLayer
0124           createTH1F("hDigiTime", "Digi Time (ns) ", nameSL.str(), 100, -100., 600.);
0125           createTH1F("hPosLeft", "Pos of Left hit (cm) in local frame", nameSL.str(), 100, -220., 220.);
0126           createTH1F("hPosRight", "Pos of Right hit (cm) in local frame", nameSL.str(), 100, -220., 220.);
0127           createTH1F("hMeanTimer", "Tmax ", nameSL.str(), 100, 200., 600.);
0128           createTH1F("hMeanTimerSeg", "Tmax from segments hits ", nameSL.str(), 100, 200., 600.);
0129           createTH2F("hMeanTimerSegAlongWire",
0130                      "Tmax from segments hits vs pos along wire ",
0131                      nameSL.str(),
0132                      40,
0133                      -150.,
0134                      150,
0135                      100,
0136                      200.,
0137                      600.);
0138           createTH2F(
0139               "hMeanTimerSegVsAngle", "Tmax from segments hits vs angle ", nameSL.str(), 40, -0., 1.0, 100, 200., 600.);
0140           createTH2F("hMeanTimerSegVsNHits",
0141                      "Tmax from segments hits vs n hits segment ",
0142                      nameSL.str(),
0143                      10,
0144                      0.5,
0145                      10.5,
0146                      100,
0147                      200.,
0148                      600.);
0149           createTH1F("hHitResidualSeg", "Hits residual wrt segments ", nameSL.str(), 100, -.2, +.2);
0150           createTH1F("hHitResidualSegCellDX", "Hits residual wrt segments semicell DX ", nameSL.str(), 100, -.2, +.2);
0151           createTH1F("hHitResidualSegCellSX", "Hits residual wrt segments semicell SX", nameSL.str(), 100, -.2, +.2);
0152           createTH2F("hHitResidualSegAlongWire",
0153                      "Hits residual wrt segments vs pos along wire ",
0154                      nameSL.str(),
0155                      40,
0156                      -150.,
0157                      150,
0158                      100,
0159                      -.2,
0160                      +.2);
0161           createTH2F("hHitResidualSegVsWireDis",
0162                      "Hits residual wrt segments vs wire distance ",
0163                      nameSL.str(),
0164                      40,
0165                      0.,
0166                      2.1,
0167                      100,
0168                      -.2,
0169                      +.2);
0170           createTH2F("hHitResidualSegVsAngle",
0171                      "Hits residual wrt segments vs impact angle ",
0172                      nameSL.str(),
0173                      40,
0174                      .0,
0175                      1.0,
0176                      100,
0177                      -.2,
0178                      +.2);
0179           createTH2F("hHitResidualSegVsNHits",
0180                      "Hits residual wrt segments vs num hits ",
0181                      nameSL.str(),
0182                      10,
0183                      -0.5,
0184                      9.5,
0185                      100,
0186                      -.2,
0187                      +.2);
0188           createTH2F(
0189               "hHitResidualSegVsChi2", "Hits residual wrt segments vs chi2 ", nameSL.str(), 25, .0, 25.0, 100, -.2, +.2);
0190
0191           createTH2F("hNsegs2dVsNhits", "N segs 2d vs n Hits ", nameSL.str(), 20, -0.5, 19.5, 5, -0.5, 4.5);
0192         }
0193         // Create hist for each Chamber
0194         createTH1F("hDigiTime", "Digi Time (ns) ", nameChamber.str(), 100, -100., 600.);
0195         createTH1F("hPosLeft", "Pos of Left hit (cm) in local frame", nameChamber.str(), 100, -220., 220.);
0196         createTH1F("hPosRight", "Pos of Right hit (cm) in local frame", nameChamber.str(), 100, -220., 220.);
0197         //segments
0198         createTH1F("hNSegs", "N segments ", nameChamber.str(), 20, 0., 20.);
0199         createTH1F("hChi2Seg", "Chi2 segments ", nameChamber.str(), 25, 0., 25.);
0200         createTH1F("hChi2SegPhi", "Chi2 segments phi ", nameChamber.str(), 25, 0., 25.);
0201         createTH1F("hChi2SegZed", "Chi2 segments zed ", nameChamber.str(), 25, 0., 25.);
0202         createTH1F("hNHitsSeg", "N hits segment ", nameChamber.str(), 15, -0.5, 14.5);
0203         createTH1F("hNHitsSegPhi", "N hits segment phi ", nameChamber.str(), 12, -0.5, 11.5);
0204         createTH2F(
0205             "hNHitsSegPhiVsAngle", "N hits segment phi vs angle ", nameChamber.str(), 40, .0, 1.0, 12, -0.5, 11.5);
0206         createTH2F("hNHitsSegPhiVsOtherHits",
0207                    "N hits segment phi vs hits in SLs ",
0208                    nameChamber.str(),
0209                    20,
0210                    -0.5,
0211                    19.5,
0212                    12,
0213                    -0.5,
0214                    11.5);
0215         createTH2F("hNHitsSegPhiVsNumSegs",
0216                    "N hits segment zed vs num segs in ch ",
0217                    nameChamber.str(),
0218                    6,
0219                    -0.5,
0220                    5.5,
0221                    12,
0222                    -0.5,
0223                    11.5);
0224         createTH2F(
0225             "hNHitsSegPhiVsChi2", "N hits segment zed vs chi2/NDoF ", nameChamber.str(), 25, 0., 25., 12, -0.5, 11.5);
0226         createTH1F("hNHitsSegZed", "N hits segment zed ", nameChamber.str(), 10, -0.5, 9.5);
0227         createTH2F(
0228             "hNHitsSegZedVsAngle", "N hits segment zed vs angle ", nameChamber.str(), 40, .0, 1.0, 12, -0.5, 11.5);
0229         createTH2F("hNHitsSegZedVsOtherHits",
0230                    "N hits segment zed vs hits in SLs ",
0231                    nameChamber.str(),
0232                    20,
0233                    -0.5,
0234                    19.5,
0235                    12,
0236                    -0.5,
0237                    11.5);
0238         createTH2F("hNHitsSegZedVsNumSegs",
0239                    "N hits segment zed vs num segs in ch ",
0240                    nameChamber.str(),
0241                    6,
0242                    -0.5,
0243                    5.5,
0244                    12,
0245                    -0.5,
0246                    11.5);
0247         createTH2F(
0248             "hNHitsSegZedVsChi2", "N hits segment zed vs chi2/NDoF ", nameChamber.str(), 25, 0., 25., 12, -0.5, 11.5);
0249
0250         createTH1F("hHitResidualSeg", "Hits residual wrt segments ", nameChamber.str(), 100, -.2, +.2);
0251         createTH1F(
0252             "hHitResidualSegCellDX", "Hits residual wrt segments semicell DX ", nameChamber.str(), 100, -.2, +.2);
0253         createTH1F("hHitResidualSegCellSX", "Hits residual wrt segments semicell SX", nameChamber.str(), 100, -.2, +.2);
0254         createTH2F("hHitResidualSegAlongWire",
0255                    "Hits residual wrt segments vs pos along wire ",
0256                    nameChamber.str(),
0257                    40,
0258                    -150.,
0259                    150,
0260                    100,
0261                    -.2,
0262                    +.2);
0263         createTH2F("hHitResidualSegVsWireDis",
0264                    "Hits residual wrt segments vs wire distance ",
0265                    nameChamber.str(),
0266                    40,
0267                    0.,
0268                    2.1,
0269                    100,
0270                    -.2,
0271                    +.2);
0272         createTH2F("hHitResidualSegVsAngle",
0273                    "Hits residual wrt segments vs impact angle ",
0274                    nameChamber.str(),
0275                    40,
0276                    .0,
0277                    1.0,
0278                    100,
0279                    -.2,
0280                    +.2);
0281         createTH2F("hHitResidualSegVsNHits",
0282                    "Hits residual wrt segments vs num hits ",
0283                    nameChamber.str(),
0284                    10,
0285                    -0.5,
0286                    9.5,
0287                    100,
0288                    -.2,
0289                    +.2);
0290         createTH2F("hHitResidualSegVsChi2",
0291                    "Hits residual wrt segments vs chi2 ",
0292                    nameChamber.str(),
0293                    25,
0294                    .0,
0295                    25.0,
0296                    100,
0297                    -.2,
0298                    +.2);
0299
0300         // eff
0301         createTH2F("hNsegs4dVsNhits", "N segs 4d vs n Hits ", nameChamber.str(), 20, 0.5, 20.5, 5, -0.5, 4.5);
0302         createTH2F("hNsegs4dVsNhitsPhi", "N segs 4d vs n HitsPhi ", nameChamber.str(), 20, 0.5, 20.5, 5, -0.5, 4.5);
0303         createTH2F("hNsegs4dVsNhitsZed", "N segs 4d vs n HitsZed ", nameChamber.str(), 20, 0.5, 20.5, 5, -0.5, 4.5);
0304
0305         createTH2F("hNsegs4dVsNsegs2d", "N segs 4d vs n segs2d ", nameChamber.str(), 4, 0.5, 4.5, 5, -0.5, 4.5);
0306         createTH2F("hNsegs4dVsNsegs2dPhi", "N segs 4d vs n segs2d Phi ", nameChamber.str(), 4, 0.5, 4.5, 5, -0.5, 4.5);
0307         createTH2F("hNsegs4dVsNsegs2dZed", "N segs 4d vs n segs2d Zed ", nameChamber.str(), 4, 0.5, 4.5, 5, -0.5, 4.5);
0308
0309         createTH2F("hNsegs2dSL1VsNsegs2dSL3", "N segs 2d SL1 vs SL3 ", nameChamber.str(), 5, -0.5, 4.5, 5, -0.5, 4.5);
0310         createTH2F("hNsegs2dSL1VsNsegs2dSL2", "N segs 2d SL1 vs SL2 ", nameChamber.str(), 5, -0.5, 4.5, 5, -0.5, 4.5);
0311         createTH2F("hNsegs2dSL2VsNsegs2dSL3", "N segs 2d SL2 vs SL3 ", nameChamber.str(), 5, -0.5, 4.5, 5, -0.5, 4.5);
0312
0313         // trigger eff
0314         createTH1F("hSegEff", "Eff LocaTrig if Seg", nameChamber.str(), 5, 0., 5.);
0315       }
0316     }
0317   }
0318
0319   // segs
0320   new TH1F("hnSegDT", "Num seg DT", 50, 0., 50.);
0321   TH1::AddDirectory(dirStat);
0322 }
0323
0324 /* Destructor */
0325 DTAnalyzerDetailed::~DTAnalyzerDetailed() {
0326   theFile->cd();
0327   theFile->Write();
0328   theFile->Close();
0329 }
0330
0331 /* Operations */
0332 void DTAnalyzerDetailed::analyze(const Event& event, const EventSetup& eventSetup) {
0333   theSync->setES(eventSetup);
0334   _ev++;
0335
0336   if (debug)
0337     cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() << " Num " << _ev << endl;
0338
0339   static int j = 1;
0340   if ((_ev % j) == 0) {
0341     if ((_ev / j) == 9)
0342       j *= 10;
0343     cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() << " Num " << _ev << endl;
0344   }
0345
0346   if (doHits)
0347     analyzeDTHits(event, eventSetup);
0348   if (doSegs)
0349     analyzeDTSegments(event, eventSetup);
0350 }
0351
0352 void DTAnalyzerDetailed::analyzeDTHits(const Event& event, const EventSetup& eventSetup) {
0353   // Get the DT Geometry
0354   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(theDTGeomToken);
0355
0356   // Get the 1D rechits from the event --------------
0357   Handle<DTRecHitCollection> dtRecHits;
0358   event.getByLabel(theRecHits1DLabel, dtRecHits);
0359
0360   int nHitDT = dtRecHits->size();
0361   histo("hnHitDT")->Fill(nHitDT);
0362
0363   //float ttrigg = 1895.; // should get this from CondDB...
0364   for (DTRecHitCollection::const_iterator hit = dtRecHits->begin(); hit != dtRecHits->end(); ++hit) {
0365     // Get the wireId of the rechit
0366     DTWireId wireId = (*hit).wireId();
0367
0368     float ttrig = theSync->offset(wireId);
0369     //cout << "TTrig " << ttrig << endl;
0370
0371     float time = (*hit).digiTime() - ttrig;
0372     double xLeft = (*hit).localPosition(DTEnums::Left).x();
0373     double xRight = (*hit).localPosition(DTEnums::Right).x();
0374
0375     histo("hDigiTime")->Fill(time);
0376     {  // per layer
0377       const DTLayerId& id((*hit).wireId().layerId());
0378       histo(hName("hDigiTime", id))->Fill(time);
0379       histo(hName("hPosLeft", id))->Fill(xLeft);
0380       histo(hName("hPosRight", id))->Fill(xRight);
0381     }
0382
0383     {  // per SL
0384       const DTSuperLayerId& id((*hit).wireId().superlayerId());
0385       histo(hName("hDigiTime", id))->Fill(time);
0386       histo(hName("hPosLeft", id))->Fill(xLeft);
0387       histo(hName("hPosRight", id))->Fill(xRight);
0388     }
0389
0390     {  // per Chamber
0391       const DTChamberId& id((*hit).wireId().chamberId());
0392       histo(hName("hDigiTime", id))->Fill(time);
0393       histo(hName("hPosLeft", id))->Fill(xLeft);
0394       histo(hName("hPosRight", id))->Fill(xRight);
0395     }
0396   }
0397
0398   // MeanTimer analysis
0399   // loop on SL
0400   //cout << "MeanTimer analysis" << endl;
0401   const std::vector<const DTSuperLayer*>& sls = dtGeom->superLayers();
0402   for (auto sl = sls.begin(); sl != sls.end(); ++sl) {
0403     DTSuperLayerId slid = (*sl)->id();
0404
0405     DTMeanTimer meanTimer(dtGeom->superLayer(slid), dtRecHits, theSync.get());
0406     vector<double> tMaxs = meanTimer.run();
0407     for (vector<double>::const_iterator tMax = tMaxs.begin(); tMax != tMaxs.end(); ++tMax) {
0408       //cout << "Filling " << hName("hMeanTimer", slid) << " with " << *tMax << endl;
0409       histo(hName("hMeanTimer", slid))->Fill(*tMax);
0410     }
0411   }
0412 }
0413
0414 void DTAnalyzerDetailed::analyzeDTSegments(const Event& event, const EventSetup& eventSetup) {
0415   if (debug)
0416     cout << "analyzeDTSegments" << endl;
0417   // Get the DT Geometry
0418   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(theDTGeomToken);
0419
0420   // Get the 4D rechit collection from the event -------------------
0421   edm::Handle<DTRecSegment4DCollection> segs;
0422   event.getByLabel(theRecHits4DLabel, segs);
0423   if (debug)
0424     cout << "4d " << segs->size() << endl;
0425
0426   // Get the 2D rechit collection from the event -------------------
0427   edm::Handle<DTRecSegment2DCollection> segs2d;
0428   event.getByLabel(theRecHits2DLabel, segs2d);
0429   if (debug)
0430     cout << "2d " << segs2d->size() << endl;
0431
0432   // Get the 1D rechits from the event --------------
0433   Handle<DTRecHitCollection> dtRecHits;
0434   event.getByLabel(theRecHits1DLabel, dtRecHits);
0435   if (debug)
0436     cout << "1d " << dtRecHits->size() << endl;
0437
0438   int nsegs = segs->size();
0439   histo("hnSegDT")->Fill(nsegs);
0440   const std::vector<const DTChamber*>& chs = dtGeom->chambers();
0441
0442   for (auto ch = chs.begin(); ch != chs.end(); ++ch) {
0443     DTChamberId chid((*ch)->id());
0444     //cout << "chid " << chid << endl;
0445     //int w= chid.wheel();
0446     //int se= chid.sector();
0447     //int st= chid.station();
0448     DTRecSegment4DCollection::range segsch = segs->get(chid);
0449     int nSegsCh = segsch.second - segsch.first;
0450     histo(hName("hNSegs", chid))->Fill(nSegsCh);
0451
0452     // some quality on segments
0453     //bool hasGoodSegment=false;
0454
0455     for (DTRecSegment4DCollection::const_iterator seg = segsch.first; seg != segsch.second; ++seg) {
0456       // first the the two projection separately
0457       // phi
0458       if (debug)
0459         cout << *seg << endl;
0460       const DTChamberRecSegment2D* phiSeg = (*seg).phiSegment();
0461
0462       vector<DTRecHit1D> phiHits;
0463       if (phiSeg) {
0464         if (debug)
0465           cout << "Phi " << *phiSeg << endl;
0466         phiHits = phiSeg->specificRecHits();
0467         //if (phiHits.size()>=6) hasGoodSegment = true;
0468
0469         DTSuperLayerId slid1(phiSeg->chamberId(), 1);
0470         /// Mean timer analysis
0471         DTMeanTimer meanTimer1(dtGeom->superLayer(slid1), phiHits, theSync.get());
0472         vector<double> tMaxs1 = meanTimer1.run();
0473         for (vector<double>::const_iterator tMax = tMaxs1.begin(); tMax != tMaxs1.end(); ++tMax) {
0474           histo(hName("hMeanTimerSeg", slid1))->Fill(*tMax);
0475           histo2d(hName("hMeanTimerSegVsNHits", slid1))->Fill(phiHits.size(), *tMax);
0476           if ((*seg).hasZed()) {
0477             histo2d(hName("hMeanTimerSegAlongWire", slid1))->Fill((*seg).localPosition().y(), *tMax);
0478             histo2d(hName("hMeanTimerSegVsAngle", slid1))->Fill(M_PI - (*seg).localDirection().theta(), *tMax);
0479           }
0480         }
0481
0482         DTSuperLayerId slid3(phiSeg->chamberId(), 3);
0483         /// Mean timer analysis
0484         DTMeanTimer meanTimer3(dtGeom->superLayer(slid3), phiHits, theSync.get());
0485         vector<double> tMaxs3 = meanTimer3.run();
0486         for (vector<double>::const_iterator tMax = tMaxs3.begin(); tMax != tMaxs3.end(); ++tMax) {
0487           histo(hName("hMeanTimerSeg", slid3))->Fill(*tMax);
0488           histo2d(hName("hMeanTimerSegVsNHits", slid1))->Fill(phiHits.size(), *tMax);
0489           if ((*seg).hasZed()) {
0490             histo2d(hName("hMeanTimerSegAlongWire", slid1))->Fill((*seg).localPosition().y(), *tMax);
0491             histo2d(hName("hMeanTimerSegVsAngle", slid1))->Fill(M_PI - (*seg).localDirection().theta(), *tMax);
0492           }
0493         }
0494       }
0495       // zed
0496       const DTSLRecSegment2D* zedSeg = (*seg).zSegment();
0497       //cout << "zedSeg " << zedSeg << endl;
0498       vector<DTRecHit1D> zedHits;
0499       if (zedSeg) {
0500         if (debug)
0501           cout << "Zed " << *zedSeg << endl;
0502         zedHits = zedSeg->specificRecHits();
0503         DTSuperLayerId slid(zedSeg->superLayerId());
0504         /// Mean timer analysis
0505         DTMeanTimer meanTimer(dtGeom->superLayer(slid), zedHits, theSync.get());
0506         vector<double> tMaxs = meanTimer.run();
0507         for (vector<double>::const_iterator tMax = tMaxs.begin(); tMax != tMaxs.end(); ++tMax) {
0508           histo(hName("hMeanTimerSeg", slid))->Fill(*tMax);
0509           histo2d(hName("hMeanTimerSegVsNHits", slid))->Fill(zedHits.size(), *tMax);
0510           if ((*seg).hasPhi()) {
0511             histo2d(hName("hMeanTimerSegAlongWire", slid))->Fill((*seg).localPosition().x(), *tMax);
0512             histo2d(hName("hMeanTimerSegVsAngle", slid))->Fill(M_PI - (*seg).localDirection().theta(), *tMax);
0513           }
0514         }
0515       }
0516
0517       histo(hName("hChi2Seg", chid))->Fill((*seg).chi2() / (*seg).degreesOfFreedom());
0518       if (phiSeg) {
0519         histo(hName("hNHitsSegPhi", chid))->Fill(phiSeg->recHits().size());
0520         histo2d(hName("hNHitsSegPhiVsAngle", chid))
0521             ->Fill(M_PI - phiSeg->localDirection().theta(), phiSeg->recHits().size());
0522         // get the total numer of hits in this SL(s)
0523         DTSuperLayerId slid1(chid, 1);
0524         DTRecHitCollection::range rangeH1 = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid1));
0525         DTSuperLayerId slid3(chid, 3);
0526         DTRecHitCollection::range rangeH3 = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid3));
0527         histo2d(hName("hNHitsSegPhiVsOtherHits", chid))
0528             ->Fill((rangeH1.second - rangeH1.first) + (rangeH3.second - rangeH3.first), phiSeg->recHits().size());
0529         histo2d(hName("hNHitsSegPhiVsNumSegs", chid))->Fill(nSegsCh, phiSeg->recHits().size());
0530         histo2d(hName("hNHitsSegPhiVsChi2", chid))
0531             ->Fill(phiSeg->chi2() / phiSeg->degreesOfFreedom(), phiSeg->recHits().size());
0532         histo(hName("hChi2SegPhi", chid))->Fill(phiSeg->chi2() / phiSeg->degreesOfFreedom());
0533       }
0534
0535       if (zedSeg) {
0536         histo(hName("hNHitsSegZed", chid))->Fill(zedSeg->recHits().size());
0537         histo2d(hName("hNHitsSegZedVsAngle", chid))
0538             ->Fill(M_PI - zedSeg->localDirection().theta(), zedSeg->recHits().size());
0539         DTSuperLayerId slid(chid, 2);
0540         DTRecHitCollection::range rangeH = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid));
0541         histo2d(hName("hNHitsSegZedVsOtherHits", chid))->Fill(rangeH.second - rangeH.first, zedSeg->recHits().size());
0542         histo2d(hName("hNHitsSegZedVsNumSegs", chid))->Fill(nSegsCh, zedSeg->recHits().size());
0543         histo2d(hName("hNHitsSegZedVsChi2", chid))
0544             ->Fill(zedSeg->chi2() / zedSeg->degreesOfFreedom(), zedSeg->recHits().size());
0545         histo(hName("hChi2SegZed", chid))->Fill(zedSeg->chi2() / zedSeg->degreesOfFreedom());
0546       }
0547       if (phiSeg && zedSeg)
0548         histo(hName("hNHitsSeg", chid))->Fill(phiSeg->recHits().size() + zedSeg->recHits().size());
0549
0550       // residual analysis
0551       if (phiSeg) {
0552         DTSegmentResidual res(phiSeg, *ch);
0553         res.run();
0554         vector<DTSegmentResidual::DTResidual> deltas = res.residuals();
0555         for (vector<DTSegmentResidual::DTResidual>::const_iterator delta = deltas.begin(); delta != deltas.end();
0556              ++delta) {
0557           histo(hName("hHitResidualSeg", (*ch)->id()))->Fill((*delta).value);
0558           if ((*delta).side == DTEnums::Right)
0559             histo(hName("hHitResidualSegCellDX", (*ch)->id()))->Fill((*delta).value);
0560           else if ((*delta).side == DTEnums::Left)
0561             histo(hName("hHitResidualSegCellSX", (*ch)->id()))->Fill((*delta).value);
0562
0563           histo2d(hName("hHitResidualSegVsWireDis", (*ch)->id()))->Fill((*delta).wireDistance, (*delta).value);
0564
0565           histo2d(hName("hHitResidualSegVsAngle", (*ch)->id()))->Fill((*delta).angle, (*delta).value);
0566           histo2d(hName("hHitResidualSegVsNHits", (*ch)->id()))->Fill(phiSeg->recHits().size(), (*delta).value);
0567           histo2d(hName("hHitResidualSegVsChi2", (*ch)->id()))->Fill(phiSeg->chi2(), (*delta).value);
0568
0569           if ((*seg).hasPhi())
0570             histo2d(hName("hHitResidualSegAlongWire", (*ch)->id()))->Fill((*seg).localPosition().x(), (*delta).value);
0571         }
0572       }
0573
0574       if (zedSeg) {
0575         const DTSuperLayer* sl = (*ch)->superLayer(2);
0576         DTSegmentResidual res(zedSeg, sl);
0577         res.run();
0578         vector<DTSegmentResidual::DTResidual> deltas = res.residuals();
0579         for (vector<DTSegmentResidual::DTResidual>::const_iterator delta = deltas.begin(); delta != deltas.end();
0580              ++delta) {
0581           histo(hName("hHitResidualSeg", sl->id()))->Fill((*delta).value);
0582           if ((*delta).side == DTEnums::Right)
0583             histo(hName("hHitResidualSegCellDX", sl->id()))->Fill((*delta).value);
0584           else if ((*delta).side == DTEnums::Left)
0585             histo(hName("hHitResidualSegCellSX", sl->id()))->Fill((*delta).value);
0586
0587           histo2d(hName("hHitResidualSegVsWireDis", sl->id()))->Fill((*delta).wireDistance, (*delta).value);
0588
0589           histo2d(hName("hHitResidualSegVsAngle", sl->id()))->Fill((*delta).angle, (*delta).value);
0590           histo2d(hName("hHitResidualSegVsNHits", sl->id()))->Fill(zedSeg->recHits().size(), (*delta).value);
0591           histo2d(hName("hHitResidualSegVsChi2", sl->id()))->Fill(zedSeg->chi2(), (*delta).value);
0592
0593           if ((*seg).hasPhi())
0594             histo2d(hName("hHitResidualSegAlongWire", sl->id()))->Fill((*seg).localPosition().x(), (*delta).value);
0595         }
0596       }  // if ZedSeg
0597
0598     }  // end loop segment in chamber
0599
0600     // Efficiency analysis
0601     // first get number of segments vs number of hits
0602
0603     // get No hits in this chamber
0604     // and no of 2d segments in chamber
0605     int nHitSl[3] = {0, 0, 0};
0606     int nSeg2dSl[3] = {0, 0, 0};
0607     for (int sl = 1; sl <= 3; ++sl) {
0608       if (sl == 2 && chid.station() == 4)
0609         continue;
0610       DTSuperLayerId slid(chid, sl);
0611       DTRecHitCollection::range rangeH = dtRecHits->get(DTRangeMapAccessor::layersBySuperLayer(slid));
0612       DTRecSegment2DCollection::range rangeS = segs2d->get(slid);
0613       nHitSl[sl - 1] = rangeH.second - rangeH.first;
0614       nSeg2dSl[sl - 1] = rangeS.second - rangeS.first;
0615       histo2d(hName("hNsegs2dVsNhits", slid))->Fill(nHitSl[sl - 1], nSeg2dSl[sl - 1]);
0616     }
0617     histo2d(hName("hNsegs4dVsNhits", chid))->Fill(nHitSl[0] + nHitSl[1] + nHitSl[2], nSegsCh);
0618     histo2d(hName("hNsegs4dVsNhitsPhi", chid))->Fill(nHitSl[0] + nHitSl[2], nSegsCh);
0619     histo2d(hName("hNsegs4dVsNhitsZed", chid))->Fill(nHitSl[1], nSegsCh);
0620
0621     histo2d(hName("hNsegs4dVsNsegs2d", chid))->Fill(nSeg2dSl[0] + nSeg2dSl[1] + nSeg2dSl[2], nSegsCh);
0622     histo2d(hName("hNsegs4dVsNsegs2dPhi", chid))->Fill(nSeg2dSl[0] + nSeg2dSl[2], nSegsCh);
0623     histo2d(hName("hNsegs4dVsNsegs2dZed", chid))->Fill(nSeg2dSl[1], nSegsCh);
0624
0625     histo2d(hName("hNsegs2dSL1VsNsegs2dSL3", chid))->Fill(nSeg2dSl[2], nSeg2dSl[0]);
0626     histo2d(hName("hNsegs2dSL1VsNsegs2dSL2", chid))->Fill(nSeg2dSl[1], nSeg2dSl[0]);
0627     histo2d(hName("hNsegs2dSL2VsNsegs2dSL3", chid))->Fill(nSeg2dSl[2], nSeg2dSl[1]);
0628
0629   }  // loop on all chamber
0630 }
0631
0632 TH1F* DTAnalyzerDetailed::histo(const string& name) const {
0633   if (TH1F* h = dynamic_cast<TH1F*>(theFile->Get(name.c_str())))
0634     return h;
0635   else
0636     throw cms::Exception("DTAnalyzerDetailed") << " Not a TH1F " << name;
0637 }
0638
0639 TH2F* DTAnalyzerDetailed::histo2d(const string& name) const {
0640   if (TH2F* h = dynamic_cast<TH2F*>(theFile->Get(name.c_str())))
0641     return h;
0642   else
0643     throw cms::Exception("DTAnalyzerDetailed") << " Not a TH2F " << name;
0644 }
0645
0646 string DTAnalyzerDetailed::toString(const DTLayerId& id) const {
0647   stringstream result;
0648   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer() << "_Lay"
0649          << id.layer();
0650   return result.str();
0651 }
0652
0653 string DTAnalyzerDetailed::toString(const DTSuperLayerId& id) const {
0654   stringstream result;
0655   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer();
0656   return result.str();
0657 }
0658
0659 string DTAnalyzerDetailed::toString(const DTChamberId& id) const {
0660   stringstream result;
0661   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station();
0662   return result.str();
0663 }
0664
0665 template <class T>
0666 string DTAnalyzerDetailed::hName(const string& s, const T& id) const {
0667   string name(toString(id));
0668   stringstream hName;
0669   hName << s << name;
0670   return hName.str();
0671 }
0672
0673 void DTAnalyzerDetailed::createTH1F(const std::string& name,
0674                                     const std::string& title,
0675                                     const std::string& suffix,
0676                                     int nbin,
0677                                     const double& binMin,
0678                                     const double& binMax) const {
0679   stringstream hName;
0680   stringstream hTitle;
0681   hName << name << suffix;
0682   hTitle << title << suffix;
0683   new TH1F(hName.str().c_str(), hTitle.str().c_str(), nbin, binMin, binMax);
0684 }
0685
0686 void DTAnalyzerDetailed::createTH2F(const std::string& name,
0687                                     const std::string& title,
0688                                     const std::string& suffix,
0689                                     int nBinX,
0690                                     const double& binXMin,
0691                                     const double& binXMax,
0692                                     int nBinY,
0693                                     const double& binYMin,
0694                                     const double& binYMax) const {
0695   stringstream hName;
0696   stringstream hTitle;
0697   hName << name << suffix;
0698   hTitle << title << suffix;
0699   new TH2F(hName.str().c_str(), hTitle.str().c_str(), nBinX, binXMin, binXMax, nBinY, binYMin, binYMax);
0700 }
0701
0702 DEFINE_FWK_MODULE(DTAnalyzerDetailed);