Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-09-20 02:50:18

0001 // system includes
0002 #include <fstream>
0003 #include <iostream>
0004 #include <map>
0005 #include <memory>
0006 #include <ostream>
0007 #include <set>
0008 #include <string>
0009 #include <vector>
0010 
0011 // user includes
0012 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0016 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0017 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0018 #include "DataFormats/TrackReco/interface/HitPattern.h"
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0034 
0035 // ROOT includes
0036 #include "TFile.h"
0037 #include "TLorentzVector.h"
0038 #include "TMath.h"
0039 
0040 class ZEEDetails : public DQMEDAnalyzer {
0041 public:
0042   ZEEDetails(const edm::ParameterSet&);
0043   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0044 
0045 protected:
0046   void analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) override;
0047   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0048 
0049 private:
0050   const std::string moduleName_;
0051   const std::string folderName_;
0052 
0053   const edm::InputTag electronTag_;
0054   const edm::InputTag bsTag_;
0055   const edm::InputTag puSummaryTag_;
0056   const edm::InputTag vertexTag_;
0057   const edm::EDGetTokenT<reco::GsfElectronCollection> electronToken_;
0058   const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0059   const edm::EDGetTokenT<std::vector<PileupSummaryInfo> > puSummaryToken_;
0060   const edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0061 
0062   const double maxEta_;
0063   const double minPt_;
0064   const double maxDeltaPhiInEB_;
0065   const double maxDeltaEtaInEB_;
0066   const double maxHOEEB_;
0067   const double maxSigmaiEiEEB_;
0068   const double maxDeltaPhiInEE_;
0069   const double maxDeltaEtaInEE_;
0070   const double maxHOEEE_;
0071   const double maxSigmaiEiEEE_;
0072   const double maxNormChi2_;
0073   const double maxD0_;
0074   const double maxDz_;
0075   const int minPixelHits_;
0076   const int minStripHits_;
0077   const double maxIso_;
0078   const double minPtHighest_;
0079   const double minInvMass_;
0080   const double maxInvMass_;
0081   const std::string trackQuality_;
0082   std::vector<float> vpu_;
0083   std::vector<float> vtrack_;
0084   const bool isMC_;
0085   const bool doPUCorrection_;
0086   const std::string puScaleFactorFile_;
0087 
0088   MonitorElement* Zpt_;
0089   MonitorElement* ZInvMass_;
0090   MonitorElement* EoverP_;
0091 };
0092 
0093 void ZEEDetails::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0094   edm::ParameterSetDescription desc;
0095   desc.addUntracked<std::string>("moduleName", "ZEEDetails");
0096   desc.addUntracked<std::string>("folderName", "ElectronTracks");
0097   desc.addUntracked<edm::InputTag>("electronInputTag", edm::InputTag("gedGsfElectrons"));
0098   desc.addUntracked<edm::InputTag>("offlineBeamSpot", edm::InputTag("offlineBeamSpot"));
0099   desc.addUntracked<edm::InputTag>("puTag", edm::InputTag("addPileupInfo"));
0100   desc.addUntracked<edm::InputTag>("vertexTag", edm::InputTag("offlinePrimaryVertices"));
0101   desc.addUntracked<double>("maxEta", 2.4);
0102   desc.addUntracked<double>("minPt", 5);
0103   desc.addUntracked<double>("maxDeltaPhiInEB", 0.15);
0104   desc.addUntracked<double>("maxDeltaEtaInEB", 0.007);
0105   desc.addUntracked<double>("maxHOEEB", 0.12);
0106   desc.addUntracked<double>("maxSigmaiEiEEB", 0.01);
0107   desc.addUntracked<double>("maxDeltaPhiInEE", 0.1);
0108   desc.addUntracked<double>("maxDeltaEtaInEE", 0.009);
0109   desc.addUntracked<double>("maxHOEEB_", .10);
0110   desc.addUntracked<double>("maxSigmaiEiEEE", 0.03);
0111   desc.addUntracked<double>("maxNormChi2", 10);
0112   desc.addUntracked<double>("maxD0", 0.02);
0113   desc.addUntracked<double>("maxDz", 20.);
0114   desc.addUntracked<uint32_t>("minPixelHits", 1);
0115   desc.addUntracked<uint32_t>("minStripHits", 8);
0116   desc.addUntracked<double>("maxIso", 0.3);
0117   desc.addUntracked<double>("minPtHighest", 24);
0118   desc.addUntracked<double>("minInvMass", 60);
0119   desc.addUntracked<double>("maxInvMass", 120);
0120   desc.addUntracked<std::string>("trackQuality", "highPurity");
0121   desc.addUntracked<bool>("isMC", false);
0122   desc.addUntracked<bool>("doPUCorrection", false);
0123   desc.addUntracked<std::string>("puScaleFactorFile", "PileupScaleFactor.root");
0124   descriptions.addWithDefaultLabel(desc);
0125 }
0126 
0127 ZEEDetails::ZEEDetails(const edm::ParameterSet& ps)
0128     : moduleName_(ps.getUntrackedParameter<std::string>("moduleName", "ZEEDetails")),
0129       folderName_(ps.getUntrackedParameter<std::string>("folderName", "ElectronTracks")),
0130       electronTag_(ps.getUntrackedParameter<edm::InputTag>("electronInputTag", edm::InputTag("gedGsfElectrons"))),
0131       bsTag_(ps.getUntrackedParameter<edm::InputTag>("offlineBeamSpot", edm::InputTag("offlineBeamSpot"))),
0132       puSummaryTag_(ps.getUntrackedParameter<edm::InputTag>("puTag", edm::InputTag("addPileupInfo"))),
0133       vertexTag_(ps.getUntrackedParameter<edm::InputTag>("vertexTag", edm::InputTag("offlinePrimaryVertices"))),
0134       electronToken_(consumes<reco::GsfElectronCollection>(electronTag_)),
0135       bsToken_(consumes<reco::BeamSpot>(bsTag_)),
0136       puSummaryToken_(consumes<std::vector<PileupSummaryInfo> >(puSummaryTag_)),
0137       vertexToken_(consumes<reco::VertexCollection>(vertexTag_)),
0138       maxEta_(ps.getUntrackedParameter<double>("maxEta", 2.4)),
0139       minPt_(ps.getUntrackedParameter<double>("minPt", 5)),
0140       maxDeltaPhiInEB_(ps.getUntrackedParameter<double>("maxDeltaPhiInEB", 0.15)),
0141       maxDeltaEtaInEB_(ps.getUntrackedParameter<double>("maxDeltaEtaInEB", 0.007)),
0142       maxHOEEB_(ps.getUntrackedParameter<double>("maxHOEEB", 0.12)),
0143       maxSigmaiEiEEB_(ps.getUntrackedParameter<double>("maxSigmaiEiEEB", 0.01)),
0144       maxDeltaPhiInEE_(ps.getUntrackedParameter<double>("maxDeltaPhiInEE", 0.1)),
0145       maxDeltaEtaInEE_(ps.getUntrackedParameter<double>("maxDeltaEtaInEE", 0.009)),
0146       maxHOEEE_(ps.getUntrackedParameter<double>("maxHOEEB_", .10)),
0147       maxSigmaiEiEEE_(ps.getUntrackedParameter<double>("maxSigmaiEiEEE", 0.03)),
0148       maxNormChi2_(ps.getUntrackedParameter<double>("maxNormChi2", 10)),
0149       maxD0_(ps.getUntrackedParameter<double>("maxD0", 0.02)),
0150       maxDz_(ps.getUntrackedParameter<double>("maxDz", 20.)),
0151       minPixelHits_(ps.getUntrackedParameter<uint32_t>("minPixelHits", 1)),
0152       minStripHits_(ps.getUntrackedParameter<uint32_t>("minStripHits", 8)),
0153       maxIso_(ps.getUntrackedParameter<double>("maxIso", 0.3)),
0154       minPtHighest_(ps.getUntrackedParameter<double>("minPtHighest", 24)),
0155       minInvMass_(ps.getUntrackedParameter<double>("minInvMass", 60)),
0156       maxInvMass_(ps.getUntrackedParameter<double>("maxInvMass", 120)),
0157       trackQuality_(ps.getUntrackedParameter<std::string>("trackQuality", "highPurity")),
0158       isMC_(ps.getUntrackedParameter<bool>("isMC", false)),
0159       doPUCorrection_(ps.getUntrackedParameter<bool>("doPUCorrection", false)),
0160       puScaleFactorFile_(ps.getUntrackedParameter<std::string>("puScaleFactorFile", "PileupScaleFactor.root")) {
0161   if (isMC_ && doPUCorrection_) {
0162     vpu_.clear();
0163     TFile* f1 = TFile::Open(puScaleFactorFile_.c_str());
0164     TH1F* h1 = dynamic_cast<TH1F*>(f1->Get("pileupweight"));
0165     for (int i = 1; i <= h1->GetNbinsX(); ++i)
0166       vpu_.push_back(h1->GetBinContent(i));
0167     f1->Close();
0168   }
0169 }
0170 
0171 void ZEEDetails::bookHistograms(DQMStore::IBooker& ibook, edm::Run const& iRun, edm::EventSetup const& iSetup) {
0172   std::string currentFolder = moduleName_ + "/" + folderName_;
0173   ibook.setCurrentFolder(currentFolder);
0174   Zpt_ = ibook.book1D("Zpt", "Z-Boson p_{T}", 100, 0.0, 100.0);
0175   ZInvMass_ = ibook.book1D("ZInvMass", "m_{ee}", 200, minInvMass_, maxInvMass_);
0176   EoverP_ = ibook.book3D("EoverP", "EoverP", 48, -2.4, 2.4, 36, -3.2, 3.2, 100, 0, 10);
0177 }
0178 
0179 void ZEEDetails::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0180   std::vector<TLorentzVector> list;
0181   std::vector<int> chrgeList;
0182   std::vector<reco::GsfElectron> finalelectrons;
0183 
0184   // Read Electron Collection
0185   edm::Handle<reco::GsfElectronCollection> electronColl;
0186   iEvent.getByToken(electronToken_, electronColl);
0187 
0188   edm::Handle<reco::BeamSpot> beamSpot;
0189   iEvent.getByToken(bsToken_, beamSpot);
0190 
0191   if (electronColl.isValid()) {
0192     for (auto const& ele : *electronColl) {
0193       if (!ele.ecalDriven())
0194         continue;
0195       if (ele.pt() < minPt_)
0196         continue;
0197       // set a max Eta cut
0198       if (!(ele.isEB() || ele.isEE()))
0199         continue;
0200 
0201       double hOverE = ele.hadronicOverEm();
0202       double sigmaee = ele.sigmaIetaIeta();
0203       double deltaPhiIn = ele.deltaPhiSuperClusterTrackAtVtx();
0204       double deltaEtaIn = ele.deltaEtaSuperClusterTrackAtVtx();
0205 
0206       // separate cut for barrel and endcap
0207       if (ele.isEB()) {
0208         if (fabs(deltaPhiIn) >= maxDeltaPhiInEB_ && fabs(deltaEtaIn) >= maxDeltaEtaInEB_ && hOverE >= maxHOEEB_ &&
0209             sigmaee >= maxSigmaiEiEEB_)
0210           continue;
0211       } else if (ele.isEE()) {
0212         if (fabs(deltaPhiIn) >= maxDeltaPhiInEE_ && fabs(deltaEtaIn) >= maxDeltaEtaInEE_ && hOverE >= maxHOEEE_ &&
0213             sigmaee >= maxSigmaiEiEEE_)
0214           continue;
0215       }
0216 
0217       reco::GsfTrackRef trk = ele.gsfTrack();
0218       reco::TrackRef tk = ele.closestCtfTrackRef();
0219       if (!trk.isNonnull())
0220         continue;  // only electrons with tracks
0221       if (!tk.isNonnull())
0222         continue;
0223       double chi2 = trk->chi2();
0224       double ndof = trk->ndof();
0225       double chbyndof = (ndof > 0) ? chi2 / ndof : 0;
0226       if (chbyndof >= maxNormChi2_)
0227         continue;
0228 
0229       double trkd0 = trk->d0();
0230       if (beamSpot.isValid()) {
0231         trkd0 = -(trk->dxy(beamSpot->position()));
0232       } else {
0233         edm::LogError("ElectronTrackProducer") << "Error >> Failed to get BeamSpot for label: " << bsTag_;
0234       }
0235       if (std::fabs(trkd0) >= maxD0_)
0236         continue;
0237 
0238       const reco::HitPattern& hitp = trk->hitPattern();
0239       int nPixelHits = hitp.numberOfValidPixelHits();
0240       if (nPixelHits < minPixelHits_)
0241         continue;
0242 
0243       int nStripHits = hitp.numberOfValidStripHits();
0244       if (nStripHits < minStripHits_)
0245         continue;
0246 
0247       // DB corrected PF Isolation
0248       reco::GsfElectron::PflowIsolationVariables pfIso = ele.pfIsolationVariables();
0249       const float eiso =
0250           pfIso.sumChargedHadronPt + std::max(0.0, pfIso.sumNeutralHadronEt + pfIso.sumPhotonEt - 0.5 * pfIso.sumPUPt);
0251       if (eiso > maxIso_ * ele.pt())
0252         continue;
0253 
0254       if (!tk->quality(reco::Track::qualityByName(trackQuality_)))
0255         continue;
0256 
0257       TLorentzVector lv;
0258       lv.SetPtEtaPhiE(ele.pt(), ele.eta(), ele.phi(), ele.energy());
0259       list.push_back(lv);
0260       chrgeList.push_back(ele.charge());
0261       finalelectrons.push_back(ele);
0262     }
0263   } else {
0264     edm::LogError("ElectronTrackProducer") << "Error >> Failed to get ElectronCollection for label: " << electronTag_;
0265   }
0266 
0267   edm::Handle<reco::VertexCollection> vertexColl;
0268   iEvent.getByToken(vertexToken_, vertexColl);
0269   if (!vertexColl.isValid()) {
0270     edm::LogError("DqmTrackStudy") << "Error! Failed to get reco::Vertex Collection, " << vertexTag_;
0271   }
0272   if (vertexColl->empty()) {
0273     edm::LogError("DqmTrackStudy") << "No good vertex in the event!!";
0274     return;
0275   }
0276 
0277   // Access PU information
0278   double wfac = 1.0;  // for data
0279   if (!iEvent.isRealData()) {
0280     edm::Handle<std::vector<PileupSummaryInfo> > PupInfo;
0281     iEvent.getByToken(puSummaryToken_, PupInfo);
0282 
0283     if (PupInfo.isValid()) {
0284       for (auto const& v : *PupInfo) {
0285         int bx = v.getBunchCrossing();
0286         if (bx == 0) {
0287           int nVertex = (vertexColl.isValid() ? vertexColl->size() : 0);
0288           if (doPUCorrection_) {
0289             if (nVertex > -1 && nVertex < int(vpu_.size()))
0290               wfac = vpu_.at(nVertex);
0291             else
0292               wfac = 0.0;
0293           }
0294         }
0295       }
0296     } else
0297       edm::LogError("DqmTrackStudy") << "PUSummary for input tag: " << puSummaryTag_ << " not found!!";
0298   }
0299 
0300   for (unsigned int I = 0; I != finalelectrons.size(); I++) {
0301     EoverP_->Fill(finalelectrons[I].superCluster()->eta(),
0302                   finalelectrons[I].superCluster()->phi(),
0303                   finalelectrons[I].eEleClusterOverPout(),
0304                   wfac);
0305   }
0306 
0307   if (list.size() >= 2) {
0308     if (chrgeList[0] + chrgeList[1] == 0) {
0309       if (list[0].Pt() >= minPtHighest_) {
0310         TLorentzVector zv = list[0] + list[1];
0311         if ((zv.M() >= minInvMass_) && (zv.M() <= maxInvMass_)) {
0312           Zpt_->Fill(zv.Pt(), wfac);
0313           ZInvMass_->Fill(zv.Mag(), wfac);
0314         }
0315       }
0316     }
0317   }
0318 }
0319 
0320 // Define this as a plug-in
0321 #include "FWCore/Framework/interface/MakerMacros.h"
0322 DEFINE_FWK_MODULE(ZEEDetails);