MomentumScaleCalibration/plugins/TestCorrection.cc

0001 #ifndef TESTCORRECTION_HH
0002 #define TESTCORRECTION_HH
0003
0004 // -*- C++ -*-
0005 //
0006 // Package:    TestCorrection
0007 // Class:      TestCorrection
0008 //
0009 /**\class TestCorrection TestCorrection.cc MuonAnalysis/MomentumScaleCalibration/plugins/TestCorrection.cc
0010
0011  Description: <one line class summary>
0012
0013  Implementation:
0014      <Notes on implementation>
0015 */
0016 //
0017 // Original Author:  Marco De Mattia
0018 //         Created:  Thu Sep 11 12:16:00 CEST 2008
0019 //
0020 //
0021
0022 // system include files
0023 #include <memory>
0024 #include <string>
0025 #include <vector>
0026
0027 // user include files
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032
0033 #include "DataFormats/TrackReco/interface/Track.h"
0034 #include "DataFormats/MuonReco/interface/Muon.h"
0035 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0036 #include "DataFormats/Candidate/interface/Candidate.h"
0037 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0038 #include "DataFormats/PatCandidates/interface/Muon.h"
0039 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0040 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0041
0042 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0043
0044 // For the momentum scale correction
0045 #include "MuScleFitBase.h"
0046 #include "MuonAnalysis/MomentumScaleCalibration/interface/MomentumScaleCorrector.h"
0047 #include "MuonAnalysis/MomentumScaleCalibration/interface/ResolutionFunction.h"
0048 #include "MuonAnalysis/MomentumScaleCalibration/interface/BackgroundFunction.h"
0049
0050 #include "TFile.h"
0051 #include "TProfile.h"
0052 #include "TH1F.h"
0053 #include "TCanvas.h"
0054 #include "TLegend.h"
0055
0056 //
0057 // class decleration
0058 //
0059
0060 class TestCorrection : public edm::one::EDAnalyzer<>, MuScleFitBase {
0061 public:
0062   explicit TestCorrection(const edm::ParameterSet&);
0063   ~TestCorrection() override;
0064
0065 private:
0066   virtual void initialize(const edm::EventSetup&);
0067   void analyze(const edm::Event&, const edm::EventSetup&) override;
0068   void endJob() override{};
0069   template <typename T>
0070   std::vector<MuScleFitMuon> fillMuonCollection(const std::vector<T>& tracks) {
0071     std::vector<MuScleFitMuon> muons;
0072     typename std::vector<T>::const_iterator track;
0073     for (track = tracks.begin(); track != tracks.end(); ++track) {
0074       reco::Particle::LorentzVector mu;
0075       mu = reco::Particle::LorentzVector(
0076           track->px(), track->py(), track->pz(), sqrt(track->p() * track->p() + +0.011163612));
0077
0078       Double_t hitsTk(0), hitsMuon(0), ptError(0);
0079       if (const reco::Muon* myMu = dynamic_cast<const reco::Muon*>(&(*track))) {
0080         hitsTk = myMu->innerTrack()->hitPattern().numberOfValidTrackerHits();
0081         hitsMuon = myMu->innerTrack()->hitPattern().numberOfValidMuonHits();
0082         ptError = myMu->innerTrack()->ptError();
0083       } else if (const pat::Muon* myMu = dynamic_cast<const pat::Muon*>(&(*track))) {
0084         hitsTk = myMu->innerTrack()->hitPattern().numberOfValidTrackerHits();
0085         hitsMuon = myMu->innerTrack()->hitPattern().numberOfValidMuonHits();
0086         ptError = myMu->innerTrack()->ptError();
0087       } else if (const reco::Track* myMu = dynamic_cast<const reco::Track*>(&(*track))) {
0088         hitsTk = myMu->hitPattern().numberOfValidTrackerHits();
0089         hitsMuon = myMu->hitPattern().numberOfValidMuonHits();
0090         ptError = myMu->ptError();
0091       }
0092
0093       MuScleFitMuon muon(mu, track->charge(), ptError, hitsTk, hitsMuon);
0094
0095       if (debug_ > 0) {
0096         edm::LogPrint("TestCorrection") << "[TestCorrection::fillMuonCollection] after MuScleFitMuon initialization"
0097                                         << std::endl;
0098         edm::LogPrint("TestCorrection") << "  muon = " << muon << std::endl;
0099       }
0100
0101       muons.push_back(muon);
0102     }
0103     return muons;
0104   }
0105
0106   lorentzVector correctMuon(const lorentzVector& muon);
0107
0108   // ----------member data ---------------------------
0109
0110   // Collections labels
0111   // ------------------
0112   TH1F* uncorrectedPt_;
0113   TProfile* uncorrectedPtVsEta_;
0114   TH1F* correctedPt_;
0115   TProfile* correctedPtVsEta_;
0116
0117   int eventCounter_;
0118
0119   std::unique_ptr<MomentumScaleCorrector> corrector_;
0120   std::unique_ptr<ResolutionFunction> resolution_;
0121   std::unique_ptr<BackgroundFunction> background_;
0122
0123   edm::EDGetTokenT<reco::MuonCollection> glbMuonsToken_;
0124   edm::EDGetTokenT<reco::TrackCollection> saMuonsToken_;
0125   edm::EDGetTokenT<reco::TrackCollection> tracksToken_;
0126 };
0127
0128 #endif  // TESTCORRECTION_HH
0129
0130 //
0131 // constructors and destructor
0132 //
0133 TestCorrection::TestCorrection(const edm::ParameterSet& iConfig)
0134     : MuScleFitBase(iConfig),
0135       glbMuonsToken_(mayConsume<reco::MuonCollection>(theMuonLabel_)),
0136       saMuonsToken_(mayConsume<reco::TrackCollection>(theMuonLabel_)),
0137       tracksToken_(mayConsume<reco::TrackCollection>(theMuonLabel_)) {
0138   //now do what ever initialization is needed
0139   TFile* outputFile = new TFile(theRootFileName_.c_str(), "RECREATE");
0140   theFiles_.push_back(outputFile);
0141   // outputFile_ = new TFile(theRootFileName_.c_str(), "RECREATE");
0142   // outputFile_->cd();
0143   outputFile->cd();
0144   MuScleFitUtils::resfind = iConfig.getParameter<std::vector<int> >("resfind");
0145   fillHistoMap(outputFile, 0);
0146   uncorrectedPt_ = new TH1F("uncorrectedPt", "uncorrected pt", 1000, 0, 100);
0147   uncorrectedPtVsEta_ = new TProfile("uncorrectedPtVsEta", "uncorrected pt vs eta", 1000, 0, 100, -3., 3.);
0148   correctedPt_ = new TH1F("correctedPt", "corrected pt", 1000, 0, 100);
0149   correctedPtVsEta_ = new TProfile("correctedPtVsEta", "corrected pt vs eta", 1000, 0, 100, -3., 3.);
0150   eventCounter_ = 0;
0151   // Create the corrector and set the parameters
0152   corrector_ =
0153       std::make_unique<MomentumScaleCorrector>(iConfig.getUntrackedParameter<std::string>("CorrectionsIdentifier"));
0154   edm::LogPrint("TestCorrection") << "corrector_ = " << &*corrector_ << std::endl;
0155   resolution_ =
0156       std::make_unique<ResolutionFunction>(iConfig.getUntrackedParameter<std::string>("ResolutionsIdentifier"));
0157   edm::LogPrint("TestCorrection") << "resolution_ = " << &*resolution_ << std::endl;
0158   background_ =
0159       std::make_unique<BackgroundFunction>(iConfig.getUntrackedParameter<std::string>("BackgroundIdentifier"));
0160
0161   // Initialize the parameters of MuScleFitUtils from those saved in the functions.
0162   // MuScleFitUtils::parScale = corrector_.getFunction(0)->parameters();
0163   MuScleFitUtils::resolutionFunction = resolution_->function(0);
0164   MuScleFitUtils::resolutionFunctionForVec = resolutionFunctionVecService(resolution_->identifiers()[0]);
0165
0166   MuScleFitUtils::parResol = resolution_->parameters();
0167 }
0168
0169 TestCorrection::~TestCorrection() {
0170   theFiles_[0]->cd();
0171   TCanvas canvas("ptComparison", "pt comparison", 1000, 800);
0172   canvas.cd();
0173   uncorrectedPt_->GetXaxis()->SetTitle("Pt(GeV)");
0174   correctedPt_->SetLineColor(kRed);
0175   TLegend* legend = new TLegend(0.7, 0.71, 0.98, 1.);
0176   legend->SetTextSize(0.02);
0177   legend->SetFillColor(0);  // Have a white background
0178   legend->AddEntry(uncorrectedPt_, "original pt");
0179   legend->AddEntry(correctedPt_, "corrected pt");
0180   uncorrectedPt_->Draw();
0181   correctedPt_->Draw("same");
0182   legend->Draw("same");
0183
0184   canvas.Write();
0185   uncorrectedPt_->Write();
0186   uncorrectedPtVsEta_->Write();
0187   correctedPt_->Write();
0188   correctedPtVsEta_->Write();
0189
0190   writeHistoMap(0);
0191   theFiles_[0]->Close();
0192
0193   edm::LogPrint("TestCorrection") << "Total analyzed events = " << eventCounter_ << std::endl;
0194 }
0195
0196 //
0197 // member functions
0198 //
0199
0200 // ------------ method called to for each event  ------------
0201 void TestCorrection::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0202   using namespace edm;
0203
0204   initialize(iSetup);
0205
0206   ++eventCounter_;
0207   if (eventCounter_ % 100 == 0) {
0208     edm::LogPrint("TestCorrection") << "Event number " << eventCounter_ << std::endl;
0209   }
0210
0211   // Take the reco-muons, depending on the type selected in the cfg
0212   // --------------------------------------------------------------
0213
0214   std::vector<MuScleFitMuon> muons;
0215
0216   if (theMuonType_ == 1) {  // GlobalMuons
0217     Handle<reco::MuonCollection> glbMuons;
0218     iEvent.getByToken(glbMuonsToken_, glbMuons);
0219     muons = fillMuonCollection(*glbMuons);
0220   } else if (theMuonType_ == 2) {  // StandaloneMuons
0221     Handle<reco::TrackCollection> saMuons;
0222     iEvent.getByToken(saMuonsToken_, saMuons);
0223     muons = fillMuonCollection(*saMuons);
0224   } else if (theMuonType_ == 3) {  // Tracker tracks
0225     Handle<reco::TrackCollection> tracks;
0226     iEvent.getByToken(tracksToken_, tracks);
0227     muons = fillMuonCollection(*tracks);
0228   }
0229
0230   // Find the two muons from the resonance, and set ResFound bool
0231   // ------------------------------------------------------------
0232   std::pair<MuScleFitMuon, MuScleFitMuon> recMuFromBestRes = MuScleFitUtils::findBestRecoRes(muons);
0233   if (MuScleFitUtils::ResFound) {
0234     MuScleFitUtils::SavedPair.push_back(std::make_pair(recMuFromBestRes.first.p4(), recMuFromBestRes.second.p4()));
0235   } else {
0236     MuScleFitUtils::SavedPair.push_back(std::make_pair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.)));
0237   }
0238
0239   // If resonance found, do the hard work
0240   // ------------------------------------
0241   if (MuScleFitUtils::ResFound) {
0242     // Find weight and reference mass for this muon pair
0243     // -------------------------------------------------
0244     // double weight = MuScleFitUtils::computeWeight ((recMu1+recMu2).mass());
0245
0246     // Use the correction function to correct the pt scale of the muons. Note that this takes into
0247     // account the corrections from all iterations.
0248     lorentzVector recMu1;
0249     recMu1 = correctMuon(recMu1);
0250     lorentzVector recMu2;
0251     recMu2 = correctMuon(recMu2);
0252
0253     reco::Particle::LorentzVector bestRecRes(recMu1 + recMu2);
0254
0255     //Fill histograms
0256     //------------------
0257     mapHisto_["hRecBestMu"]->Fill(recMu1);
0258     if ((std::abs(recMu1.eta()) < 2.5) && (recMu1.pt() > 2.5)) {
0259       mapHisto_["hRecBestMu_Acc"]->Fill(recMu1);
0260     }
0261     mapHisto_["hRecBestMu"]->Fill(recMu2);
0262     if ((std::abs(recMu2.eta()) < 2.5) && (recMu2.pt() > 2.5)) {
0263       mapHisto_["hRecBestMu_Acc"]->Fill(recMu2);
0264     }
0265     mapHisto_["hDeltaRecBestMu"]->Fill(recMu1, recMu2);
0266
0267     mapHisto_["hRecBestRes"]->Fill(bestRecRes);
0268     if ((std::abs(recMu1.eta()) < 2.5) && (recMu1.pt() > 2.5) && (std::abs(recMu2.eta()) < 2.5) &&
0269         (recMu2.pt() > 2.5)) {
0270       mapHisto_["hRecBestRes_Acc"]->Fill(bestRecRes);
0271       // Fill histogram of Res mass vs muon variable
0272       mapHisto_["hRecBestResVSMu"]->Fill(recMu1, bestRecRes, -1);
0273       mapHisto_["hRecBestResVSMu"]->Fill(recMu2, bestRecRes, +1);
0274     }
0275   }
0276
0277   // Loop on the recMuons
0278   std::vector<MuScleFitMuon>::const_iterator recMuon = muons.begin();
0279   int muonCount = 0;
0280   for (; recMuon != muons.end(); ++recMuon, ++muonCount) {
0281     // Fill the histogram with uncorrected pt values
0282     uncorrectedPt_->Fill(recMuon->pt());
0283     uncorrectedPtVsEta_->Fill(recMuon->pt(), recMuon->eta());
0284
0285     // Fill the histogram with corrected pt values
0286     edm::LogPrint("TestCorrection") << "correcting muon[" << muonCount << "] with pt = " << recMuon->pt() << std::endl;
0287     double corrPt = (*corrector_)(*recMuon);
0288     edm::LogPrint("TestCorrection") << "to pt = " << corrPt << std::endl;
0289     correctedPt_->Fill(corrPt);
0290     correctedPtVsEta_->Fill(corrPt, recMuon->eta());
0291     // correctedPt_->Fill(recMuon->pt());
0292   }
0293 }
0294
0295 lorentzVector TestCorrection::correctMuon(const lorentzVector& muon) {
0296   double corrPt = corrector_->correct(muon);
0297   double ptEtaPhiE[4] = {corrPt, muon.Eta(), muon.Phi(), muon.E()};
0298   return MuScleFitUtils::fromPtEtaPhiToPxPyPz(ptEtaPhiE);
0299 }
0300
0301 // ------------ method called once each job just before starting event loop  ------------
0302 void TestCorrection::initialize(const edm::EventSetup&) {
0303   // Read the pdf from root file. They are used by massProb when finding the muon pair, needed
0304   // for the mass histograms.
0305   readProbabilityDistributionsFromFile();
0306 }
0307
0308 #include "FWCore/Framework/interface/MakerMacros.h"
0309 DEFINE_FWK_MODULE(TestCorrection);