MomentumScaleCalibration/plugins/TestResolution.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestResolution
0004 // Class:      TestResolution
0005 //
0006 /**\class TestResolution TestResolution.cc MuonAnalysis/MomentumScaleCalibration/plugins/TestResolution.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Marco De Mattia
0015 //         Created:  Thu Sep 11 12:16:00 CEST 2008
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <string>
0022 #include <vector>
0023
0024 // user include files
0025 #include "DataFormats/Candidate/interface/Candidate.h"
0026 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0027 #include "DataFormats/MuonReco/interface/Muon.h"
0028 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0029 #include "DataFormats/TrackReco/interface/Track.h"
0030 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0036 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0037
0038 // For the momentum scale resolution
0039 #include "MuonAnalysis/MomentumScaleCalibration/interface/ResolutionFunction.h"
0040
0041 // ROOT includes
0042 #include "TCanvas.h"
0043 #include "TLegend.h"
0044 #include "TFile.h"
0045 #include "TProfile.h"
0046
0047 //
0048 // class decleration
0049 //
0050
0051 class TestResolution : public edm::one::EDAnalyzer<> {
0052 public:
0053   explicit TestResolution(const edm::ParameterSet&);
0054   ~TestResolution() override;
0055
0056 private:
0057   void analyze(const edm::Event&, const edm::EventSetup&) override;
0058   template <typename T>
0059   std::vector<reco::LeafCandidate> fillMuonCollection(const std::vector<T>& tracks) {
0060     std::vector<reco::LeafCandidate> muons;
0061     typename std::vector<T>::const_iterator track;
0062     for (track = tracks.begin(); track != tracks.end(); ++track) {
0063       // Where 0.011163612 is the squared muon mass.
0064       reco::Particle::LorentzVector mu(
0065           track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + 0.011163612));
0066       reco::LeafCandidate muon(track->charge(), mu);
0067       // Store muon
0068       // ----------
0069       muons.push_back(muon);
0070     }
0071     return muons;
0072   }
0073
0074   // ----------member data ---------------------------
0075
0076   // Collections labels
0077   // ------------------
0078   edm::InputTag theMuonLabel_;
0079   edm::EDGetTokenT<reco::MuonCollection> glbMuonsToken_;
0080   edm::EDGetTokenT<reco::TrackCollection> saMuonsToken_;
0081   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0082
0083   int theMuonType_;
0084   std::string theRootFileName_;
0085   TFile* outputFile_;
0086
0087   TProfile* sigmaPt_;
0088
0089   int eventCounter_;
0090
0091   std::unique_ptr<ResolutionFunction> resolutionFunction_;
0092 };
0093
0094 //
0095 // constructors and destructor
0096 //
0097 TestResolution::TestResolution(const edm::ParameterSet& iConfig)
0098     : theMuonLabel_(iConfig.getParameter<edm::InputTag>("MuonLabel")),
0099       glbMuonsToken_(mayConsume<reco::MuonCollection>(theMuonLabel_)),
0100       saMuonsToken_(mayConsume<reco::TrackCollection>(theMuonLabel_)),
0101       tracksToken_(mayConsume<reco::TrackCollection>(theMuonLabel_)),
0102       theMuonType_(iConfig.getParameter<int>("MuonType")),
0103       theRootFileName_(iConfig.getUntrackedParameter<std::string>("OutputFileName")) {
0104   //now do what ever initialization is needed
0105   outputFile_ = new TFile(theRootFileName_.c_str(), "RECREATE");
0106   outputFile_->cd();
0107   sigmaPt_ = new TProfile("sigmaPtOverPt", "sigmaPt/Pt vs muon Pt", 1000, 0, 100);
0108   eventCounter_ = 0;
0109   // Create the corrector and set the parameters
0110   resolutionFunction_ =
0111       std::make_unique<ResolutionFunction>(iConfig.getUntrackedParameter<std::string>("ResolutionsIdentifier"));
0112   std::cout << "resolutionFunction_ = " << &*resolutionFunction_ << std::endl;
0113 }
0114
0115 TestResolution::~TestResolution() {
0116   outputFile_->cd();
0117   TCanvas canvas("sigmaPtOverPt", "sigmaPt/Pt vs muon Pt", 1000, 800);
0118   canvas.cd();
0119   sigmaPt_->GetXaxis()->SetTitle("Pt(GeV)");
0120   //   TLegend * legend = new TLegend(0.7,0.71,0.98,1.);
0121   //   legend->SetTextSize(0.02);
0122   //   legend->SetFillColor(0); // Have a white background
0123   //   legend->AddEntry(uncorrectedPt_, "original pt");
0124   //   legend->AddEntry(correctedPt_, "corrected pt");
0125   sigmaPt_->Draw();
0126   //   legend->Draw("same");
0127
0128   canvas.Write();
0129   sigmaPt_->Write();
0130   outputFile_->Close();
0131
0132   std::cout << "Total analyzed events = " << eventCounter_ << std::endl;
0133 }
0134
0135 //
0136 // member functions
0137 //
0138
0139 // ------------ method called to for each event  ------------
0140 void TestResolution::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0141   using namespace edm;
0142
0143   ++eventCounter_;
0144   if (eventCounter_ % 100 == 0) {
0145     std::cout << "Event number " << eventCounter_ << std::endl;
0146   }
0147
0148   // Take the reco-muons, depending on the type selected in the cfg
0149   // --------------------------------------------------------------
0150
0151   std::vector<reco::LeafCandidate> muons;
0152
0153   if (theMuonType_ == 1) {  // GlobalMuons
0154     Handle<reco::MuonCollection> glbMuons;
0155     iEvent.getByToken(glbMuonsToken_, glbMuons);
0156     muons = fillMuonCollection(*glbMuons);
0157   } else if (theMuonType_ == 2) {  // StandaloneMuons
0158     Handle<reco::TrackCollection> saMuons;
0159     iEvent.getByToken(saMuonsToken_, saMuons);
0160     muons = fillMuonCollection(*saMuons);
0161   } else if (theMuonType_ == 3) {  // Tracker tracks
0162     Handle<reco::TrackCollection> tracks;
0163     iEvent.getByToken(tracksToken_, tracks);
0164     muons = fillMuonCollection(*tracks);
0165   }
0166
0167   // Loop on the recMuons
0168   std::vector<reco::LeafCandidate>::const_iterator recMuon = muons.begin();
0169   for (; recMuon != muons.end(); ++recMuon) {
0170     // Fill the histogram with uncorrected pt values
0171     sigmaPt_->Fill(resolutionFunction_->sigmaPt(*recMuon, 0), recMuon->pt());
0172   }
0173 }
0174
0175 //define this as a plug-in
0176 DEFINE_FWK_MODULE(TestResolution);