L1TTrackMatch/test/L1TkMETAnalyser.cc

0001 #include <cmath>
0002 #include <iostream>
0003 #include <map>
0004 #include <string>
0005 #include <vector>
0006
0007 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0008 #include "DataFormats/Candidate/interface/Candidate.h"
0009 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0010 #include "DataFormats/L1TCorrelator/interface/TkEtMiss.h"
0011 #include "DataFormats/L1TCorrelator/interface/TkEtMissFwd.h"
0012 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0013 #include "DataFormats/L1Trigger/interface/EtSum.h"
0014 #include "DataFormats/L1Trigger/interface/Vertex.h"
0015 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
0016 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0017 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0025 #include "L1Trigger/L1TTrackMatch/interface/L1TkEtMissEmuAlgo.h"
0026 #include "TEfficiency.h"
0027 #include "TGraphAsymmErrors.h"
0028 #include "TGraphErrors.h"
0029 #include "TH1F.h"
0030 #include "TH2F.h"
0031 #include "TPad.h"
0032 #include "TProfile.h"
0033 #include "TTree.h"
0034
0035 using namespace std;
0036
0037 class L1TkMETAnalyser : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0038 public:
0039   explicit L1TkMETAnalyser(const edm::ParameterSet& iConfig);
0040   ~L1TkMETAnalyser() override = default;
0041
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043
0044 private:
0045   void beginJob() override;
0046   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0047   void endJob() override;
0048
0049   edm::ParameterSet config;
0050
0051   edm::InputTag TrackMETSimInputTag;
0052   edm::InputTag TrackMETEmuInputTag;
0053   edm::InputTag TrackMETHWInputTag;
0054
0055   edm::EDGetTokenT<std::vector<l1t::TkEtMiss>> TrackMETSimToken_;
0056   edm::EDGetTokenT<std::vector<l1t::EtSum>> TrackMETEmuToken_;
0057   edm::EDGetTokenT<std::vector<l1t::EtSum>> TrackMETHWToken_;
0058
0059   float EtScale;
0060   float EtphiScale;
0061
0062   bool available_;
0063   TTree* eventTree;
0064
0065   std::vector<float>* m_SimMET;
0066   std::vector<float>* m_EmuMET;
0067   std::vector<float>* m_HwMET;
0068
0069   std::vector<float>* m_SimMETphi;
0070   std::vector<float>* m_EmuMETphi;
0071   std::vector<float>* m_HwMETphi;
0072
0073   std::vector<float>* m_SimNtrk;
0074   std::vector<float>* m_EmuNtrk;
0075   std::vector<float>* m_HwNtrk;
0076
0077   bool HW_analysis_;
0078
0079   TH1F* hisTkSimMET_;
0080   TH1F* hisTkEmuMET_;
0081   TH1F* hisTkHWMET_;
0082
0083   TH1F* hisTkSimPhi_;
0084   TH1F* hisTkEmuPhi_;
0085   TH1F* hisTkHWPhi_;
0086
0087   TH1F* hisTkSimNtrk_;
0088   TH1F* hisTkEmuNtrk_;
0089   TH1F* hisTkHWNtrk_;
0090
0091   TH1F* hisMETResidual_;
0092   TH1F* hisPhiResidual_;
0093   TH1F* hisNtrkResidual_;
0094 };
0095
0096 L1TkMETAnalyser::L1TkMETAnalyser(edm::ParameterSet const& iConfig) : config(iConfig) {
0097   HW_analysis_ = iConfig.getParameter<bool>("HW_Analysis");
0098
0099   TrackMETSimInputTag = iConfig.getParameter<edm::InputTag>("TrackMETInputTag");
0100   TrackMETEmuInputTag = iConfig.getParameter<edm::InputTag>("TrackMETEmuInputTag");
0101   if (HW_analysis_) {
0102     TrackMETHWInputTag = iConfig.getParameter<edm::InputTag>("TrackMETHWInputTag");
0103   }
0104   TrackMETSimToken_ = consumes<std::vector<l1t::TkEtMiss>>(TrackMETSimInputTag);
0105   TrackMETEmuToken_ = consumes<std::vector<l1t::EtSum>>(TrackMETEmuInputTag);
0106   if (HW_analysis_) {
0107     TrackMETHWToken_ = consumes<std::vector<l1t::EtSum>>(TrackMETHWInputTag);
0108   }
0109   usesResource(TFileService::kSharedResource);
0110 }
0111
0112 void L1TkMETAnalyser::beginJob() {
0113   edm::Service<TFileService> fs;
0114   TFileDirectory inputDir = fs->mkdir("TkMETAnalysis");
0115   available_ = fs.isAvailable();
0116   if (not available_)
0117     return;  // No ROOT file open.
0118
0119   m_SimMET = new std::vector<float>;
0120   m_EmuMET = new std::vector<float>;
0121   m_HwMET = new std::vector<float>;
0122
0123   m_SimMETphi = new std::vector<float>;
0124   m_EmuMETphi = new std::vector<float>;
0125   m_HwMETphi = new std::vector<float>;
0126
0127   m_SimNtrk = new std::vector<float>;
0128   m_EmuNtrk = new std::vector<float>;
0129   m_HwNtrk = new std::vector<float>;
0130
0131   eventTree = fs->make<TTree>("eventTree", "Event tree");
0132
0133   eventTree->Branch("SimMET", &m_SimMET);
0134   eventTree->Branch("EmuMET", &m_EmuMET);
0135   eventTree->Branch("HwMET", &m_HwMET);
0136
0137   eventTree->Branch("SimMETphi", &m_SimMETphi);
0138   eventTree->Branch("EmuMETphi", &m_EmuMETphi);
0139   eventTree->Branch("HwMETphi", &m_HwMETphi);
0140
0141   eventTree->Branch("SimNtrk", &m_SimNtrk);
0142   eventTree->Branch("EmuNtrk", &m_EmuNtrk);
0143   eventTree->Branch("HwNtrk", &m_HwNtrk);
0144
0145   hisTkSimMET_ = inputDir.make<TH1F>("hisTkSimMET_", "sim TkMET [GeV]", 101, 0, 500);
0146   hisTkEmuMET_ = inputDir.make<TH1F>("hisTkEmuMET_", "emu TkMET [GeV]", 101, 0, 500);
0147
0148   hisTkSimPhi_ = inputDir.make<TH1F>("hisTkSimPhi_", "sim phi [rad]", 101, -M_PI, M_PI);
0149   hisTkEmuPhi_ = inputDir.make<TH1F>("hisTkEmuPhi_", "emu phi [rad]", 101, -M_PI, M_PI);
0150
0151   hisTkSimNtrk_ = inputDir.make<TH1F>("hisTkSimNtrk_", "sim ntrks", 101, 0, 256);
0152   hisTkEmuNtrk_ = inputDir.make<TH1F>("hisTkEmuNtrk_", "emu ntrks", 101, 0, 256);
0153
0154   if (!HW_analysis_) {
0155     hisMETResidual_ = inputDir.make<TH1F>("hisMETResidual_", "sim - emu TkMET [GeV]", 101, -100, 100);
0156     hisPhiResidual_ = inputDir.make<TH1F>("hisPhiResidual_", "sim - emu phi [rad]", 101, -1, 1);
0157     hisNtrkResidual_ = inputDir.make<TH1F>("hisNtrkResidual_", "sim - emu ntrks", 101, -10, 10);
0158   }
0159
0160   if (HW_analysis_) {
0161     hisTkHWMET_ = inputDir.make<TH1F>("hisTkHWMET_", "hw TkMET [GeV]", 101, 0, 500);
0162     hisTkHWPhi_ = inputDir.make<TH1F>("hisTkHWPhi_", "hw phi [rad]", 101, -M_PI, M_PI);
0163     hisTkHWNtrk_ = inputDir.make<TH1F>("hisTkEmuNtrk_", "hw ntrks", 101, 0, 256);
0164
0165     hisMETResidual_ = inputDir.make<TH1F>("hisMETResidual_", "emu - hw TkMET [GeV]", 101, -100, 100);
0166     hisPhiResidual_ = inputDir.make<TH1F>("hisPhiResidual_", "emu - hw phi [rad]", 101, -1, 1);
0167     hisNtrkResidual_ = inputDir.make<TH1F>("hisNtrkResidual_", "emu - hw ntrks", 101, -10, 10);
0168   }
0169 }
0170
0171 void L1TkMETAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0172   if (not available_)
0173     return;  // No ROOT file open.
0174
0175   edm::Handle<std::vector<l1t::TkEtMiss>> L1TkMETSimHandle;
0176   edm::Handle<std::vector<l1t::EtSum>> L1TkMETEmuHandle;
0177
0178   iEvent.getByToken(TrackMETSimToken_, L1TkMETSimHandle);
0179   iEvent.getByToken(TrackMETEmuToken_, L1TkMETEmuHandle);
0180
0181   m_SimMET->clear();
0182   m_EmuMET->clear();
0183   m_HwMET->clear();
0184   m_SimMETphi->clear();
0185   m_EmuMETphi->clear();
0186   m_HwMETphi->clear();
0187   m_SimNtrk->clear();
0188   m_EmuNtrk->clear();
0189   m_HwNtrk->clear();
0190
0191   float SimEtmiss = L1TkMETSimHandle->begin()->etMiss();
0192   float EmuEtmiss = L1TkMETEmuHandle->begin()->hwPt() * l1tmetemu::kStepMETwordEt;
0193
0194   float SimEtPhi = L1TkMETSimHandle->begin()->etPhi();
0195   float EmuEtPhi = L1TkMETEmuHandle->begin()->hwPhi() * l1tmetemu::kStepMETwordPhi;
0196
0197   int SimEtNtrk = L1TkMETSimHandle->begin()->etQual();
0198   int EmuEtNtrk = L1TkMETEmuHandle->begin()->hwQual();
0199
0200   if (!HW_analysis_) {
0201     hisMETResidual_->Fill(EmuEtmiss - SimEtmiss);
0202     hisPhiResidual_->Fill(EmuEtPhi - SimEtPhi);
0203     hisNtrkResidual_->Fill(EmuEtNtrk - SimEtNtrk);
0204   }
0205
0206   m_SimMET->push_back(SimEtmiss);
0207   m_EmuMET->push_back(EmuEtmiss);
0208   m_SimMETphi->push_back(SimEtPhi);
0209   m_EmuMETphi->push_back(EmuEtPhi);
0210   m_SimNtrk->push_back(SimEtNtrk);
0211   m_EmuNtrk->push_back(EmuEtNtrk);
0212
0213   hisTkSimMET_->Fill(SimEtmiss);
0214   hisTkEmuMET_->Fill(EmuEtmiss);
0215
0216   hisTkSimPhi_->Fill(SimEtPhi);
0217   hisTkEmuPhi_->Fill(EmuEtPhi);
0218
0219   hisTkSimNtrk_->Fill(SimEtNtrk);
0220   hisTkEmuNtrk_->Fill(EmuEtNtrk);
0221
0222   if (HW_analysis_) {
0223     edm::Handle<std::vector<l1t::EtSum>> L1TkMETHWHandle;
0224     iEvent.getByToken(TrackMETHWToken_, L1TkMETHWHandle);
0225     float HWEtmiss = L1TkMETHWHandle->begin()->hwPt();
0226     float HWEtPhi = L1TkMETHWHandle->begin()->hwPhi();
0227     int HWEtNtrk = L1TkMETHWHandle->begin()->hwQual();
0228
0229     hisTkHWMET_->Fill(HWEtmiss);
0230     hisTkHWPhi_->Fill(HWEtPhi);
0231     hisTkHWNtrk_->Fill(HWEtNtrk);
0232
0233     hisMETResidual_->Fill(EmuEtmiss - HWEtmiss);
0234     hisPhiResidual_->Fill(EmuEtPhi - HWEtPhi);
0235     hisNtrkResidual_->Fill(EmuEtNtrk - HWEtNtrk);
0236
0237     m_HwMET->push_back(HWEtmiss);
0238     m_HwMETphi->push_back(HWEtPhi);
0239     m_HwNtrk->push_back(HWEtNtrk);
0240   }
0241   eventTree->Fill();
0242 }
0243
0244 //////////
0245 // END JOB
0246 void L1TkMETAnalyser::endJob() {
0247   // things to be done at the exit of the event Loop
0248   edm::LogInfo("L1TkMETAnalyser") << "==================== TkMET RECONSTRUCTION ======================\n"
0249                                   << "MET Residual Bias: " << hisMETResidual_->GetMean() << " GeV\n"
0250                                   << "MET Resolution: " << hisMETResidual_->GetStdDev() << " GeV\n"
0251                                   << "Phi Residual Bias: " << hisPhiResidual_->GetMean() << " rad\n"
0252                                   << "Phi Resolution: " << hisPhiResidual_->GetStdDev() << " rad\n"
0253                                   << "NTrk Residual Bias: " << hisNtrkResidual_->GetMean() << " Tracks\n"
0254                                   << "Ntrk Resolution: " << hisNtrkResidual_->GetStdDev() << " Tracks\n";
0255 }
0256
0257 void L1TkMETAnalyser::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0258   //The following says we do not know what parameters are allowed so do no validation
0259   // Please change this to state exactly what you do use, even if it is no parameters
0260   edm::ParameterSetDescription desc;
0261   desc.setUnknown();
0262   descriptions.addDefault(desc);
0263 }
0264
0265 ///////////////////////////
0266 // DEFINE THIS AS A PLUG-IN
0267 DEFINE_FWK_MODULE(L1TkMETAnalyser);