Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-23 02:30:03

0001 // -*- C++ -*-
0002 //
0003 // Package:    MuonTimingValidator
0004 // Class:      MuonTimingValidator
0005 //
0006 /**\class MuonTimingValidator MuonTimingValidator.cc 
0007 
0008  Description: An example analyzer that fills muon timing information histograms 
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Piotr Traczyk
0015 //         Created:  Wed Sep 27 14:54:28 EDT 2006
0016 //
0017 //
0018 
0019 #include "RecoMuon/MuonIdentification/test/MuonTimingValidator.h"
0020 #include "FWCore/Utilities/interface/isFinite.h"
0021 
0022 // system include files
0023 #include <memory>
0024 #include <string>
0025 #include <iostream>
0026 #include <fstream>
0027 #include <iostream>
0028 #include <iomanip>
0029 
0030 // user include files
0031 #include "FWCore/Framework/interface/Frameworkfwd.h"
0032 
0033 #include "FWCore/Framework/interface/Event.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "FWCore/Utilities/interface/InputTag.h"
0038 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0039 
0040 #include "DataFormats/Common/interface/Ref.h"
0041 
0042 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0043 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0044 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0045 #include "Geometry/DTGeometry/interface/DTLayer.h"
0046 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0047 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0048 
0049 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
0050 #include <DataFormats/MuonDetId/interface/CSCDetId.h>
0051 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0052 #include <DataFormats/CSCRecHit/interface/CSCRangeMapAccessor.h>
0053 
0054 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0055 #include "DataFormats/MuonReco/interface/MuonCocktails.h"
0056 #include "DataFormats/MuonReco/interface/MuonTimeExtra.h"
0057 #include "DataFormats/MuonReco/interface/MuonTimeExtraMap.h"
0058 
0059 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0060 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0061 
0062 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0063 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0064 
0065 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0066 
0067 #include <TROOT.h>
0068 #include <TSystem.h>
0069 #include <TFile.h>
0070 #include <TH1.h>
0071 #include <TH2.h>
0072 #include <TLegend.h>
0073 #include <TStyle.h>
0074 #include <TCanvas.h>
0075 #include <TFrame.h>
0076 
0077 //
0078 // constructors and destructor
0079 //
0080 MuonTimingValidator::MuonTimingValidator(const edm::ParameterSet& iConfig)
0081     : TKtrackTags_(iConfig.getUntrackedParameter<edm::InputTag>("TKtracks")),
0082       TKtrackTokens_(consumes<reco::TrackCollection>(TKtrackTags_)),
0083       MuonTags_(iConfig.getUntrackedParameter<edm::InputTag>("Muons")),
0084       MuonTokens_(consumes<reco::MuonCollection>(MuonTags_)),
0085       CombinedTimeTags_(iConfig.getUntrackedParameter<edm::InputTag>("CombinedTiming")),
0086       CombinedTimeTokens_(consumes<reco::MuonTimeExtraMap>(CombinedTimeTags_)),
0087       DtTimeTags_(iConfig.getUntrackedParameter<edm::InputTag>("DtTiming")),
0088       DtTimeTokens_(consumes<reco::MuonTimeExtraMap>(DtTimeTags_)),
0089       CscTimeTags_(iConfig.getUntrackedParameter<edm::InputTag>("CscTiming")),
0090       CscTimeTokens_(consumes<reco::MuonTimeExtraMap>(CscTimeTags_)),
0091       trackingGeometryToken_(esConsumes()),
0092       out(iConfig.getParameter<std::string>("out")),
0093       open(iConfig.getParameter<std::string>("open")),
0094       theMinEta(iConfig.getParameter<double>("etaMin")),
0095       theMaxEta(iConfig.getParameter<double>("etaMax")),
0096       theMinPt(iConfig.getParameter<double>("simPtMin")),
0097       thePtCut(iConfig.getParameter<double>("PtCut")),
0098       theMinPtres(iConfig.getParameter<double>("PtresMin")),
0099       theMaxPtres(iConfig.getParameter<double>("PtresMax")),
0100       theScale(iConfig.getParameter<double>("PlotScale")),
0101       theDtCut(iConfig.getParameter<int>("DTcut")),
0102       theCscCut(iConfig.getParameter<int>("CSCcut")),
0103       theNBins(iConfig.getParameter<int>("nbins")) {
0104   //now do what ever initialization is needed
0105 }
0106 
0107 MuonTimingValidator::~MuonTimingValidator() {
0108   // do anything here that needs to be done at desctruction time
0109   // (e.g. close files, deallocate resources etc.)
0110   if (hFile != 0) {
0111     hFile->Close();
0112     delete hFile;
0113   }
0114 }
0115 
0116 //
0117 // member functions
0118 //
0119 
0120 // ------------ method called to for each event  ------------
0121 void MuonTimingValidator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0122   //  std::cout << "*** Begin Muon Timing Validatior " << std::endl;
0123   //  std::cout << " Event: " << iEvent.id() << "  Orbit: " << iEvent.orbitNumber() << "  BX: " << iEvent.bunchCrossing() << std::endl;
0124 
0125   iEvent.getByToken(TKtrackTokens_, TKTrackCollection);
0126   reco::TrackCollection tkTC;
0127   const reco::TrackCollection tkTC1 = *(TKTrackCollection.product());
0128 
0129   iEvent.getByToken(MuonTokens_, MuCollection);
0130   const reco::MuonCollection muonC = *(MuCollection.product());
0131   if (!muonC.size())
0132     return;
0133 
0134   iEvent.getByToken(CombinedTimeTokens_, timeMap1);
0135   const reco::MuonTimeExtraMap& timeMapCmb = *timeMap1;
0136   iEvent.getByToken(DtTimeTokens_, timeMap2);
0137   const reco::MuonTimeExtraMap& timeMapDT = *timeMap2;
0138   iEvent.getByToken(CscTimeTokens_, timeMap3);
0139   const reco::MuonTimeExtraMap& timeMapCSC = *timeMap3;
0140 
0141   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry = iSetup.getHandle(trackingGeometryToken_);
0142 
0143   reco::MuonCollection::const_iterator imuon;
0144 
0145   int imucount = 0;
0146   for (imuon = muonC.begin(); imuon != muonC.end(); ++imuon) {
0147     LogDebug("RecoMuonIdentification") << " Event: " << iEvent.id() << " Found muon. Pt: " << imuon->p();
0148 
0149     if ((fabs(imuon->eta()) < theMinEta) || (fabs(imuon->eta()) > theMaxEta))
0150       continue;
0151 
0152     reco::TrackRef trkTrack = imuon->track();
0153     if (trkTrack.isNonnull()) {
0154       hi_tk_pt->Fill(((*trkTrack).pt()));
0155       hi_tk_phi->Fill(((*trkTrack).phi()));
0156       hi_tk_eta->Fill(((*trkTrack).eta()));
0157     }
0158 
0159     reco::TrackRef staTrack = imuon->standAloneMuon();
0160     if (staTrack.isNonnull()) {
0161       hi_sta_pt->Fill((*staTrack).pt());
0162       hi_sta_phi->Fill((*staTrack).phi());
0163       hi_sta_eta->Fill(((*staTrack).eta()));
0164     }
0165 
0166     reco::TrackRef glbTrack = imuon->combinedMuon();
0167     if (glbTrack.isNonnull()) {
0168       hi_glb_pt->Fill((*glbTrack).pt());
0169       hi_glb_phi->Fill((*glbTrack).phi());
0170       hi_glb_eta->Fill((*glbTrack).eta());
0171 
0172       LogDebug("RecoMuonIdentification") << " Global Pt: " << (*glbTrack).pt();
0173     }
0174 
0175     // Analyze the short info stored directly in reco::Muon
0176 
0177     reco::MuonTime muonTime;
0178     if (imuon->isTimeValid()) {
0179       muonTime = imuon->time();
0180       //      std::cout << "    Time points: " << muonTime.nDof << "  time: " << muonTime.timeAtIpInOut << std::endl;
0181       if (muonTime.nDof) {
0182         hi_mutime_vtx->Fill(muonTime.timeAtIpInOut);
0183         hi_mutime_vtx_err->Fill(muonTime.timeAtIpInOutErr);
0184       }
0185     }
0186 
0187     reco::MuonEnergy muonE;
0188     //    if (imuon->isEnergyValid() && fabs(imuon->eta())<1.5) {
0189     if (imuon->isEnergyValid()) {
0190       muonE = imuon->calEnergy();
0191       if (muonE.emMax > 0.25) {
0192         hi_ecal_time->Fill(muonE.ecal_time);
0193         if (muonE.emMax > 1.)
0194           hi_ecal_time_ecut->Fill(muonE.ecal_time);
0195         hi_ecal_energy->Fill(muonE.emMax);
0196         if (muonE.hadMax > 0.25)
0197           hi_hcalecal_vtx->Fill(muonE.hcal_time - muonE.ecal_time);
0198 
0199         double emErr = 1.5 / muonE.emMax;
0200 
0201         if (emErr > 0.) {
0202           hi_ecal_time_err->Fill(emErr);
0203           hi_ecal_time_pull->Fill((muonE.ecal_time - 1.) / emErr);
0204           LogDebug("RecoMuonIdentification") << "     ECAL time: " << muonE.ecal_time << " +/- " << emErr;
0205         }
0206       }
0207 
0208       if (muonE.hadMax > 0.25) {
0209         hi_hcal_time->Fill(muonE.hcal_time);
0210         if (muonE.hadMax > 1.)
0211           hi_hcal_time_ecut->Fill(muonE.hcal_time);
0212         hi_hcal_energy->Fill(muonE.hadMax);
0213 
0214         double hadErr = 1.;  // DUMMY!!!
0215 
0216         hi_hcal_time_err->Fill(hadErr);
0217         hi_hcal_time_pull->Fill((muonE.hcal_time - 1.) / hadErr);
0218         LogDebug("RecoMuonIdentification") << "     HCAL time: " << muonE.hcal_time << " +/- " << hadErr;
0219       }
0220     }
0221 
0222     // Analyze the MuonTimeExtra information
0223     reco::MuonRef muonR(MuCollection, imucount);
0224     reco::MuonTimeExtra timec = timeMapCmb[muonR];
0225     reco::MuonTimeExtra timedt = timeMapDT[muonR];
0226     reco::MuonTimeExtra timecsc = timeMapCSC[muonR];
0227 
0228     hi_cmbtime_ndof->Fill(timec.nDof());
0229     hi_dttime_ndof->Fill(timedt.nDof());
0230     hi_csctime_ndof->Fill(timecsc.nDof());
0231 
0232     if (timedt.nDof() > theDtCut) {
0233       LogDebug("RecoMuonIdentification") << "          DT nDof: " << timedt.nDof();
0234       LogDebug("RecoMuonIdentification") << "          DT Time: " << timedt.timeAtIpInOut() << " +/- "
0235                                          << timedt.inverseBetaErr();
0236       LogDebug("RecoMuonIdentification") << "         DT FreeB: " << timedt.freeInverseBeta() << " +/- "
0237                                          << timedt.freeInverseBetaErr();
0238 
0239       hi_dttime_ibt->Fill(timedt.inverseBeta());
0240       hi_dttime_ibt_pt->Fill(imuon->pt(), timedt.inverseBeta());
0241       hi_dttime_ibt_err->Fill(timedt.inverseBetaErr());
0242       hi_dttime_fib->Fill(timedt.freeInverseBeta());
0243       hi_dttime_fib_err->Fill(timedt.freeInverseBetaErr());
0244       hi_dttime_vtx->Fill(timedt.timeAtIpInOut());
0245       hi_dttime_vtx_err->Fill(timedt.timeAtIpInOutErr());
0246       hi_dttime_vtxr->Fill(timedt.timeAtIpOutIn());
0247       hi_dttime_vtxr_err->Fill(timedt.timeAtIpOutInErr());
0248       hi_dttime_errdiff->Fill(timedt.timeAtIpInOutErr() - timedt.timeAtIpOutInErr());
0249 
0250       if (timedt.inverseBetaErr() > 0.)
0251         hi_dttime_ibt_pull->Fill((timedt.inverseBeta() - 1.) / timedt.inverseBetaErr());
0252       if (timedt.freeInverseBetaErr() > 0.)
0253         hi_dttime_fib_pull->Fill((timedt.freeInverseBeta() - 1.) / timedt.freeInverseBetaErr());
0254       if (timedt.timeAtIpInOutErr() > 0.)
0255         hi_dttime_vtx_pull->Fill(timedt.timeAtIpInOut() / timedt.timeAtIpInOutErr());
0256       if (timedt.timeAtIpOutInErr() > 0.)
0257         hi_dttime_vtxr_pull->Fill(timedt.timeAtIpOutIn() / timedt.timeAtIpOutInErr());
0258 
0259       if (timecsc.nDof() > theCscCut)
0260         hi_dtcsc_vtx->Fill(timedt.timeAtIpInOut() - timecsc.timeAtIpInOut());
0261       if (imuon->isEnergyValid()) {
0262         if (muonE.emMax > 0.25)
0263           hi_dtecal_vtx->Fill(timedt.timeAtIpInOut() - muonE.ecal_time);
0264         if (muonE.hadMax > 0.25)
0265           hi_dthcal_vtx->Fill(timedt.timeAtIpInOut() - muonE.hcal_time);
0266       }
0267     }
0268 
0269     if (timecsc.nDof() > theCscCut) {
0270       LogDebug("RecoMuonIdentification") << "         CSC nDof: " << timecsc.nDof();
0271       LogDebug("RecoMuonIdentification") << "         CSC Time: " << timecsc.timeAtIpInOut() << " +/- "
0272                                          << timecsc.inverseBetaErr();
0273       LogDebug("RecoMuonIdentification") << "        CSC FreeB: " << timecsc.freeInverseBeta() << " +/- "
0274                                          << timecsc.freeInverseBetaErr();
0275 
0276       hi_csctime_ibt->Fill(timecsc.inverseBeta());
0277       hi_csctime_ibt_pt->Fill(imuon->pt(), timecsc.inverseBeta());
0278       hi_csctime_ibt_err->Fill(timecsc.inverseBetaErr());
0279       hi_csctime_fib->Fill(timecsc.freeInverseBeta());
0280       hi_csctime_fib_err->Fill(timecsc.freeInverseBetaErr());
0281       hi_csctime_vtx->Fill(timecsc.timeAtIpInOut());
0282       hi_csctime_vtx_err->Fill(timecsc.timeAtIpInOutErr());
0283       hi_csctime_vtxr->Fill(timecsc.timeAtIpOutIn());
0284       hi_csctime_vtxr_err->Fill(timecsc.timeAtIpOutInErr());
0285 
0286       if (timec.inverseBetaErr() > 0.)
0287         hi_csctime_ibt_pull->Fill((timecsc.inverseBeta() - 1.) / timecsc.inverseBetaErr());
0288       if (timecsc.freeInverseBetaErr() > 0.)
0289         hi_csctime_fib_pull->Fill((timecsc.freeInverseBeta() - 1.) / timecsc.freeInverseBetaErr());
0290       if (timecsc.timeAtIpInOutErr() > 0.)
0291         hi_csctime_vtx_pull->Fill(timecsc.timeAtIpInOut() / timecsc.timeAtIpInOutErr());
0292       if (timecsc.timeAtIpOutInErr() > 0.)
0293         hi_csctime_vtxr_pull->Fill(timecsc.timeAtIpOutIn() / timecsc.timeAtIpOutInErr());
0294 
0295       if (imuon->isEnergyValid()) {
0296         if (muonE.emMax > 0.25)
0297           hi_ecalcsc_vtx->Fill(muonE.ecal_time - timecsc.timeAtIpInOut());
0298         if (muonE.hadMax > 0.25)
0299           hi_hcalcsc_vtx->Fill(muonE.hcal_time - timecsc.timeAtIpInOut());
0300       }
0301     }
0302 
0303     if (timec.nDof() > 0) {
0304       LogDebug("RecoMuonIdentification") << "    Combined nDof: " << timec.nDof();
0305       LogDebug("RecoMuonIdentification") << "    Combined Time: " << timec.timeAtIpInOut() << " +/- "
0306                                          << timec.inverseBetaErr();
0307       LogDebug("RecoMuonIdentification") << "   Combined FreeB: " << timec.freeInverseBeta() << " +/- "
0308                                          << timec.freeInverseBetaErr();
0309 
0310       hi_cmbtime_ibt->Fill(timec.inverseBeta());
0311       hi_cmbtime_ibt_pt->Fill(imuon->pt(), timec.inverseBeta());
0312       hi_cmbtime_ibt_err->Fill(timec.inverseBetaErr());
0313       hi_cmbtime_fib->Fill(timec.freeInverseBeta());
0314       hi_cmbtime_fib_err->Fill(timec.freeInverseBetaErr());
0315       hi_cmbtime_vtx->Fill(timec.timeAtIpInOut());
0316       hi_cmbtime_vtx_err->Fill(timec.timeAtIpInOutErr());
0317       hi_cmbtime_vtxr->Fill(timec.timeAtIpOutIn());
0318       hi_cmbtime_vtxr_err->Fill(timec.timeAtIpOutInErr());
0319 
0320       if (timec.inverseBetaErr() > 0.)
0321         hi_cmbtime_ibt_pull->Fill((timec.inverseBeta() - 1.) / timec.inverseBetaErr());
0322       if (timec.freeInverseBetaErr() > 0.)
0323         hi_cmbtime_fib_pull->Fill((timec.freeInverseBeta() - 1.) / timec.freeInverseBetaErr());
0324       if (timec.timeAtIpInOutErr() > 0.)
0325         hi_cmbtime_vtx_pull->Fill(timec.timeAtIpInOut() / timec.timeAtIpInOutErr());
0326       if (timec.timeAtIpOutInErr() > 0.)
0327         hi_cmbtime_vtxr_pull->Fill(timec.timeAtIpOutIn() / timec.timeAtIpOutInErr());
0328     }
0329 
0330     imucount++;
0331   }
0332 }
0333 
0334 // ------------ method called once each job just before starting event loop  ------------
0335 void MuonTimingValidator::beginJob() {
0336   hFile = new TFile(out.c_str(), open.c_str());
0337   hFile->cd();
0338 
0339   effStyle = new TStyle("effStyle", "Efficiency Study Style");
0340   effStyle->SetCanvasBorderMode(0);
0341   effStyle->SetPadBorderMode(1);
0342   effStyle->SetOptTitle(0);
0343   effStyle->SetStatFont(42);
0344   effStyle->SetTitleFont(22);
0345   effStyle->SetCanvasColor(10);
0346   effStyle->SetPadColor(0);
0347   effStyle->SetLabelFont(42, "x");
0348   effStyle->SetLabelFont(42, "y");
0349   effStyle->SetHistFillStyle(1001);
0350   effStyle->SetHistFillColor(0);
0351   effStyle->SetOptStat(0);
0352   effStyle->SetOptFit(0111);
0353   effStyle->SetStatH(0.05);
0354 
0355   hi_sta_pt = new TH1F("hi_sta_pt", "P_{T}^{STA}", theNBins, 0.0, theMaxPtres);
0356   hi_tk_pt = new TH1F("hi_tk_pt", "P_{T}^{TK}", theNBins, 0.0, theMaxPtres);
0357   hi_glb_pt = new TH1F("hi_glb_pt", "P_{T}^{GLB}", theNBins, 0.0, theMaxPtres);
0358 
0359   hi_sta_phi = new TH1F("hi_sta_phi", "#phi^{STA}", theNBins, -3.0, 3.);
0360   hi_tk_phi = new TH1F("hi_tk_phi", "#phi^{TK}", theNBins, -3.0, 3.);
0361   hi_glb_phi = new TH1F("hi_glb_phi", "#phi^{GLB}", theNBins, -3.0, 3.);
0362 
0363   hi_mutime_vtx = new TH1F("hi_mutime_vtx", "Time at Vertex (inout)", theNBins, -25. * theScale, 25. * theScale);
0364   hi_mutime_vtx_err = new TH1F("hi_mutime_vtx_err", "Time at Vertex Error (inout)", theNBins, 0., 25.0);
0365 
0366   hi_dtcsc_vtx = new TH1F("hi_dtcsc_vtx", "Time at Vertex (DT-CSC)", theNBins, -25. * theScale, 25. * theScale);
0367   hi_dtecal_vtx = new TH1F("hi_dtecal_vtx", "Time at Vertex (DT-ECAL)", theNBins, -25. * theScale, 25. * theScale);
0368   hi_ecalcsc_vtx = new TH1F("hi_ecalcsc_vtx", "Time at Vertex (ECAL-CSC)", theNBins, -25. * theScale, 25. * theScale);
0369   hi_dthcal_vtx = new TH1F("hi_dthcal_vtx", "Time at Vertex (DT-HCAL)", theNBins, -25. * theScale, 25. * theScale);
0370   hi_hcalcsc_vtx = new TH1F("hi_hcalcsc_vtx", "Time at Vertex (HCAL-CSC)", theNBins, -25. * theScale, 25. * theScale);
0371   hi_hcalecal_vtx =
0372       new TH1F("hi_hcalecal_vtx", "Time at Vertex (HCAL-ECAL)", theNBins, -25. * theScale, 25. * theScale);
0373 
0374   hi_cmbtime_ibt = new TH1F("hi_cmbtime_ibt", "Inverse Beta", theNBins, 0., 1.6);
0375   hi_cmbtime_ibt_pt =
0376       new TH2F("hi_cmbtime_ibt_pt", "P{T} vs Inverse Beta", theNBins, 0.0, theMaxPtres, theNBins, 0.7, 2.0);
0377   hi_cmbtime_ibt_err = new TH1F("hi_cmbtime_ibt_err", "Inverse Beta Error", theNBins, 0., 1.0);
0378   hi_cmbtime_fib = new TH1F("hi_cmbtime_fib", "Free Inverse Beta", theNBins, -5. * theScale, 7. * theScale);
0379   hi_cmbtime_fib_err = new TH1F("hi_cmbtime_fib_err", "Free Inverse Beta Error", theNBins, 0, 5.);
0380   hi_cmbtime_vtx = new TH1F("hi_cmbtime_vtx", "Time at Vertex (inout)", theNBins, -25. * theScale, 25. * theScale);
0381   hi_cmbtime_vtx_err = new TH1F("hi_cmbtime_vtx_err", "Time at Vertex Error (inout)", theNBins, 0., 25.0);
0382   hi_cmbtime_vtxr = new TH1F("hi_cmbtime_vtxR", "Time at Vertex (inout)", theNBins, 0., 75. * theScale);
0383   hi_cmbtime_vtxr_err = new TH1F("hi_cmbtime_vtxR_err", "Time at Vertex Error (inout)", theNBins, 0., 25.0);
0384   hi_cmbtime_ibt_pull = new TH1F("hi_cmbtime_ibt_pull", "Inverse Beta Pull", theNBins, -5., 5.0);
0385   hi_cmbtime_fib_pull = new TH1F("hi_cmbtime_fib_pull", "Free Inverse Beta Pull", theNBins, -5., 5.0);
0386   hi_cmbtime_vtx_pull = new TH1F("hi_cmbtime_vtx_pull", "Time at Vertex Pull (inout)", theNBins, -5., 5.0);
0387   hi_cmbtime_vtxr_pull = new TH1F("hi_cmbtime_vtxR_pull", "Time at Vertex Pull (inout)", theNBins, -5., 5.0);
0388 
0389   hi_cmbtime_ndof = new TH1F("hi_cmbtime_ndof", "Number of timing measurements", 48, 0., 48.0);
0390 
0391   hi_dttime_ibt = new TH1F("hi_dttime_ibt", "DT Inverse Beta", theNBins, 0., 1.6);
0392   hi_dttime_ibt_pt =
0393       new TH2F("hi_dttime_ibt_pt", "P{T} vs DT Inverse Beta", theNBins, 0.0, theMaxPtres, theNBins, 0.7, 2.0);
0394   hi_dttime_ibt_err = new TH1F("hi_dttime_ibt_err", "DT Inverse Beta Error", theNBins, 0., 0.3);
0395   hi_dttime_fib = new TH1F("hi_dttime_fib", "DT Free Inverse Beta", theNBins, -5. * theScale, 7. * theScale);
0396   hi_dttime_fib_err = new TH1F("hi_dttime_fib_err", "DT Free Inverse Beta Error", theNBins, 0, 5.);
0397   hi_dttime_vtx = new TH1F("hi_dttime_vtx", "DT Time at Vertex (inout)", theNBins, -25. * theScale, 25. * theScale);
0398   hi_dttime_vtx_err = new TH1F("hi_dttime_vtx_err", "DT Time at Vertex Error (inout)", theNBins, 0., 10.0);
0399   hi_dttime_vtxr = new TH1F("hi_dttime_vtxR", "DT Time at Vertex (inout)", theNBins, 0., 75. * theScale);
0400   hi_dttime_vtxr_err = new TH1F("hi_dttime_vtxR_err", "DT Time at Vertex Error (inout)", theNBins, 0., 10.0);
0401   hi_dttime_ibt_pull = new TH1F("hi_dttime_ibt_pull", "DT Inverse Beta Pull", theNBins, -5., 5.0);
0402   hi_dttime_fib_pull = new TH1F("hi_dttime_fib_pull", "DT Free Inverse Beta Pull", theNBins, -5., 5.0);
0403   hi_dttime_vtx_pull = new TH1F("hi_dttime_vtx_pull", "DT Time at Vertex Pull (inout)", theNBins, -5., 5.0);
0404   hi_dttime_vtxr_pull = new TH1F("hi_dttime_vtxR_pull", "DT Time at Vertex Pull (inout)", theNBins, -5., 5.0);
0405   hi_dttime_errdiff = new TH1F(
0406       "hi_dttime_errdiff", "DT Time at Vertex inout-outin error difference", theNBins, -2. * theScale, 2. * theScale);
0407 
0408   hi_dttime_ndof = new TH1F("hi_dttime_ndof", "Number of DT timing measurements", 48, 0., 48.0);
0409 
0410   hi_csctime_ibt = new TH1F("hi_csctime_ibt", "CSC Inverse Beta", theNBins, 0., 1.6);
0411   hi_csctime_ibt_pt =
0412       new TH2F("hi_csctime_ibt_pt", "P{T} vs CSC Inverse Beta", theNBins, 0.0, theMaxPtres, theNBins, 0.7, 2.0);
0413   hi_csctime_ibt_err = new TH1F("hi_csctime_ibt_err", "CSC Inverse Beta Error", theNBins, 0., 1.0);
0414   hi_csctime_fib = new TH1F("hi_csctime_fib", "CSC Free Inverse Beta", theNBins, -5. * theScale, 7. * theScale);
0415   hi_csctime_fib_err = new TH1F("hi_csctime_fib_err", "CSC Free Inverse Beta Error", theNBins, 0, 5.);
0416   hi_csctime_vtx = new TH1F("hi_csctime_vtx", "CSC Time at Vertex (inout)", theNBins, -25. * theScale, 25. * theScale);
0417   hi_csctime_vtx_err = new TH1F("hi_csctime_vtx_err", "CSC Time at Vertex Error (inout)", theNBins, 0., 25.0);
0418   hi_csctime_vtxr = new TH1F("hi_csctime_vtxR", "CSC Time at Vertex (inout)", theNBins, 0., 75. * theScale);
0419   hi_csctime_vtxr_err = new TH1F("hi_csctime_vtxR_err", "CSC Time at Vertex Error (inout)", theNBins, 0., 25.0);
0420   hi_csctime_ibt_pull = new TH1F("hi_csctime_ibt_pull", "CSC Inverse Beta Pull", theNBins, -5., 5.0);
0421   hi_csctime_fib_pull = new TH1F("hi_csctime_fib_pull", "CSC Free Inverse Beta Pull", theNBins, -5., 5.0);
0422   hi_csctime_vtx_pull = new TH1F("hi_csctime_vtx_pull", "CSC Time at Vertex Pull (inout)", theNBins, -5., 5.0);
0423   hi_csctime_vtxr_pull = new TH1F("hi_csctime_vtxR_pull", "CSC Time at Vertex Pull (inout)", theNBins, -5., 5.0);
0424 
0425   hi_csctime_ndof = new TH1F("hi_csctime_ndof", "Number of CSC timing measurements", 48, 0., 48.0);
0426 
0427   hi_ecal_time = new TH1F("hi_ecal_time", "ECAL Time at Vertex (inout)", theNBins, -40. * theScale, 40. * theScale);
0428   hi_ecal_time_err = new TH1F("hi_ecal_time_err", "ECAL Time at Vertex Error", theNBins, 0., 20.0);
0429   hi_ecal_time_pull = new TH1F("hi_ecal_time_pull", "ECAL Time at Vertex Pull", theNBins, -7.0, 7.0);
0430   hi_ecal_time_ecut = new TH1F(
0431       "hi_ecal_time_ecut", "ECAL Time at Vertex (inout) after energy cut", theNBins, -20. * theScale, 20. * theScale);
0432   hi_ecal_energy = new TH1F("hi_ecal_energy", "ECAL max energy in 5x5 crystals", theNBins, .0, 5.0);
0433 
0434   hi_hcal_time = new TH1F("hi_hcal_time", "HCAL Time at Vertex (inout)", theNBins, -40. * theScale, 40. * theScale);
0435   hi_hcal_time_err = new TH1F("hi_hcal_time_err", "HCAL Time at Vertex Error", theNBins, 0., 20.0);
0436   hi_hcal_time_pull = new TH1F("hi_hcal_time_pull", "HCAL Time at Vertex Pull", theNBins, -7.0, 7.0);
0437   hi_hcal_time_ecut = new TH1F(
0438       "hi_hcal_time_ecut", "HCAL Time at Vertex (inout) after energy cut", theNBins, -20. * theScale, 20. * theScale);
0439   hi_hcal_energy = new TH1F("hi_hcal_energy", "HCAL max energy in 5x5 crystals", theNBins, .0, 5.0);
0440 
0441   hi_sta_eta = new TH1F("hi_sta_eta", "#eta^{STA}", theNBins / 2, theMinEta, theMaxEta);
0442   hi_tk_eta = new TH1F("hi_tk_eta", "#eta^{TK}", theNBins / 2, theMinEta, theMaxEta);
0443   hi_glb_eta = new TH1F("hi_glb_eta", "#eta^{GLB}", theNBins / 2, theMinEta, theMaxEta);
0444 }
0445 
0446 // ------------ method called once each job just after ending the event loop  ------------
0447 void MuonTimingValidator::endJob() {
0448   hFile->cd();
0449 
0450   gROOT->SetStyle("effStyle");
0451 
0452   hi_sta_pt->Write();
0453   hi_tk_pt->Write();
0454   hi_glb_pt->Write();
0455 
0456   hi_sta_phi->Write();
0457   hi_tk_phi->Write();
0458   hi_glb_phi->Write();
0459 
0460   hi_sta_eta->Write();
0461   hi_tk_eta->Write();
0462   hi_glb_eta->Write();
0463 
0464   hi_mutime_vtx->Write();
0465   hi_mutime_vtx_err->Write();
0466 
0467   hFile->mkdir("differences");
0468   hFile->cd("differences");
0469 
0470   hi_dtcsc_vtx->Write();
0471   hi_dtecal_vtx->Write();
0472   hi_ecalcsc_vtx->Write();
0473   hi_dthcal_vtx->Write();
0474   hi_hcalcsc_vtx->Write();
0475   hi_hcalecal_vtx->Write();
0476 
0477   hFile->cd();
0478   hFile->mkdir("combined");
0479   hFile->cd("combined");
0480 
0481   hi_cmbtime_ibt->Write();
0482   hi_cmbtime_ibt_pt->Write();
0483   hi_cmbtime_ibt_err->Write();
0484   hi_cmbtime_fib->Write();
0485   hi_cmbtime_fib_err->Write();
0486   hi_cmbtime_vtx->Write();
0487   hi_cmbtime_vtx_err->Write();
0488   hi_cmbtime_vtxr->Write();
0489   hi_cmbtime_vtxr_err->Write();
0490   hi_cmbtime_ibt_pull->Write();
0491   hi_cmbtime_fib_pull->Write();
0492   hi_cmbtime_vtx_pull->Write();
0493   hi_cmbtime_vtxr_pull->Write();
0494   hi_cmbtime_ndof->Write();
0495 
0496   hFile->cd();
0497   hFile->mkdir("dt");
0498   hFile->cd("dt");
0499 
0500   hi_dttime_ibt->Write();
0501   hi_dttime_ibt_pt->Write();
0502   hi_dttime_ibt_err->Write();
0503   hi_dttime_fib->Write();
0504   hi_dttime_fib_err->Write();
0505   hi_dttime_vtx->Write();
0506   hi_dttime_vtx_err->Write();
0507   hi_dttime_vtxr->Write();
0508   hi_dttime_vtxr_err->Write();
0509   hi_dttime_ibt_pull->Write();
0510   hi_dttime_fib_pull->Write();
0511   hi_dttime_vtx_pull->Write();
0512   hi_dttime_vtxr_pull->Write();
0513   hi_dttime_errdiff->Write();
0514   hi_dttime_ndof->Write();
0515 
0516   hFile->cd();
0517   hFile->mkdir("csc");
0518   hFile->cd("csc");
0519 
0520   hi_csctime_ibt->Write();
0521   hi_csctime_ibt_pt->Write();
0522   hi_csctime_ibt_err->Write();
0523   hi_csctime_fib->Write();
0524   hi_csctime_fib_err->Write();
0525   hi_csctime_vtx->Write();
0526   hi_csctime_vtx_err->Write();
0527   hi_csctime_vtxr->Write();
0528   hi_csctime_vtxr_err->Write();
0529   hi_csctime_ibt_pull->Write();
0530   hi_csctime_fib_pull->Write();
0531   hi_csctime_vtx_pull->Write();
0532   hi_csctime_vtxr_pull->Write();
0533   hi_csctime_ndof->Write();
0534 
0535   hFile->cd();
0536   hFile->mkdir("ecal");
0537   hFile->cd("ecal");
0538 
0539   hi_ecal_time->Write();
0540   hi_ecal_time_ecut->Write();
0541   hi_ecal_time_err->Write();
0542   hi_ecal_time_pull->Write();
0543   hi_ecal_energy->Write();
0544 
0545   hFile->cd();
0546   hFile->mkdir("hcal");
0547   hFile->cd("hcal");
0548 
0549   hi_hcal_time->Write();
0550   hi_hcal_time_ecut->Write();
0551   hi_hcal_time_err->Write();
0552   hi_hcal_time_pull->Write();
0553   hi_hcal_energy->Write();
0554 
0555   hFile->Write();
0556 }
0557 
0558 float MuonTimingValidator::calculateDistance(const math::XYZVector& vect1, const math::XYZVector& vect2) {
0559   float dEta = vect1.eta() - vect2.eta();
0560   float dPhi = fabs(Geom::Phi<float>(vect1.phi()) - Geom::Phi<float>(vect2.phi()));
0561   float distance = sqrt(pow(dEta, 2) + pow(dPhi, 2));
0562 
0563   return distance;
0564 }
0565 
0566 //
0567 // return h1/h2 with recalculated errors
0568 //
0569 TH1F* MuonTimingValidator::divideErr(TH1F* h1, TH1F* h2, TH1F* hout) {
0570   hout->Reset();
0571   hout->Divide(h1, h2, 1., 1., "B");
0572 
0573   for (int i = 0; i <= hout->GetNbinsX() + 1; i++) {
0574     Float_t tot = h2->GetBinContent(i);
0575     Float_t tot_e = h2->GetBinError(i);
0576     Float_t eff = hout->GetBinContent(i);
0577     Float_t Err = 0.;
0578     if (tot > 0)
0579       Err = tot_e / tot * sqrt(eff * (1 - eff));
0580     if (eff == 1. || edm::isNotFinite(Err))
0581       Err = 1.e-3;
0582     hout->SetBinError(i, Err);
0583   }
0584   return hout;
0585 }
0586 
0587 //define this as a plug-in
0588 DEFINE_FWK_MODULE(MuonTimingValidator);