RecoEgamma/plugins/ElectronConversionRejectionValidator.cc

0001 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0002 #include "FWCore/ServiceRegistry/interface/Service.h"
0003 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0004 #include "Validation/RecoEgamma/plugins/ElectronConversionRejectionValidator.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0011 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0012 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0013 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0014 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0015 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0016 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0017 #include "DataFormats/Math/interface/deltaPhi.h"
0018 #include "CommonTools/Egamma/interface/ConversionTools.h"
0019 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0020
0021 #include "TFile.h"
0022 #include "TH1.h"
0023 #include "TH2.h"
0024 #include "TTree.h"
0025 #include "TVector3.h"
0026 #include "TProfile.h"
0027 #include "TMath.h"
0028
0029 #include <iostream>
0030
0031 /** \class ElectronConversionRejectionValidator
0032  **
0033  **
0034  **  $Id: ElectronConversionRejectionValidator
0035  **  \author J.Bendavid
0036  **
0037  ***/
0038
0039 using namespace std;
0040
0041 ElectronConversionRejectionValidator::ElectronConversionRejectionValidator(const edm::ParameterSet& pset) {
0042   fName_ = pset.getUntrackedParameter<std::string>("Name");
0043   verbosity_ = pset.getUntrackedParameter<int>("Verbosity");
0044   parameters_ = pset;
0045
0046   gsfElectronCollectionProducer_ = pset.getParameter<std::string>("gsfElectronProducer");
0047   gsfElectronCollection_ = pset.getParameter<std::string>("gsfElectronCollection");
0048
0049   conversionCollectionProducer_ = pset.getParameter<std::string>("convProducer");
0050   conversionCollection_ = pset.getParameter<std::string>("conversionCollection");
0051   // conversionTrackProducer_ = pset.getParameter<std::string>("trackProducer");
0052   gsfElecToken_ =
0053       consumes<reco::GsfElectronCollection>(edm::InputTag(gsfElectronCollectionProducer_, gsfElectronCollection_));
0054   convToken_ =
0055       consumes<reco::ConversionCollection>(edm::InputTag(conversionCollectionProducer_, conversionCollection_));
0056
0057   isRunCentrally_ = pset.getParameter<bool>("isRunCentrally");
0058
0059   elePtMin_ = pset.getParameter<double>("elePtMin");
0060   eleExpectedHitsInnerMax_ = pset.getParameter<int>("eleExpectedHitsInnerMax");
0061   eleD0Max_ = pset.getParameter<double>("eleD0Max");
0062   offline_pvToken_ = consumes<reco::VertexCollection>(
0063       pset.getUntrackedParameter<edm::InputTag>("offlinePV", edm::InputTag("offlinePrimaryVertices")));
0064   beamspotToken_ =
0065       consumes<reco::BeamSpot>(pset.getUntrackedParameter<edm::InputTag>("beamspot", edm::InputTag("offlineBeamSpot")));
0066
0067   nEvt_ = 0;
0068   nEntry_ = 0;
0069 }
0070
0071 ElectronConversionRejectionValidator::~ElectronConversionRejectionValidator() {}
0072
0073 void ElectronConversionRejectionValidator::bookHistograms(DQMStore::IBooker& ibooker,
0074                                                           edm::Run const&,
0075                                                           edm::EventSetup const&) {
0076   double ptMin = parameters_.getParameter<double>("ptMin");
0077   double ptMax = parameters_.getParameter<double>("ptMax");
0078   int ptBin = parameters_.getParameter<int>("ptBin");
0079
0080   double trackptMin = parameters_.getParameter<double>("trackptMin");
0081   double trackptMax = parameters_.getParameter<double>("trackptMax");
0082   int trackptBin = parameters_.getParameter<int>("trackptBin");
0083
0084   double etaMin = parameters_.getParameter<double>("etaMin");
0085   double etaMax = parameters_.getParameter<double>("etaMax");
0086   int etaBin = parameters_.getParameter<int>("etaBin");
0087
0088   double phiMin = -TMath::Pi();
0089   double phiMax = TMath::Pi();
0090   int phiBin = parameters_.getParameter<int>("phiBin");
0091
0092   double rhoMin = parameters_.getParameter<double>("rhoMin");
0093   double rhoMax = parameters_.getParameter<double>("rhoMax");
0094   int rhoBin = parameters_.getParameter<int>("rhoBin");
0095
0096   double zMin = parameters_.getParameter<double>("zMin");
0097   double zMax = parameters_.getParameter<double>("zMax");
0098   int zBin = parameters_.getParameter<int>("zBin");
0099
0100   //// All MC photons
0101   // SC from reco photons
0102
0103   //TString simfolder = TString(
0104   //std::string simpath = dqmpath_ + "SimulationInfo";
0105   ibooker.setCurrentFolder(dqmpath_);
0106   //
0107   // simulation information about conversions
0108   // Histograms for efficiencies
0109   h_elePtAll_ = ibooker.book1D("elePtAll", "# of Electrons", ptBin, ptMin, ptMax);
0110   h_eleEtaAll_ = ibooker.book1D("eleEtaAll", "# of Electrons", etaBin, etaMin, etaMax);
0111   h_elePhiAll_ = ibooker.book1D("elePhiAll", "# of Electrons", phiBin, phiMin, phiMax);
0112
0113   h_elePtPass_ = ibooker.book1D("elePtPass", "# of Electrons", ptBin, ptMin, ptMax);
0114   h_eleEtaPass_ = ibooker.book1D("eleEtaPass", "# of Electrons", etaBin, etaMin, etaMax);
0115   h_elePhiPass_ = ibooker.book1D("elePhiPass", "# of Electrons", phiBin, phiMin, phiMax);
0116
0117   h_elePtFail_ = ibooker.book1D("elePtFail", "# of Electrons", ptBin, ptMin, ptMax);
0118   h_eleEtaFail_ = ibooker.book1D("eleEtaFail", "# of Electrons", etaBin, etaMin, etaMax);
0119   h_elePhiFail_ = ibooker.book1D("elePhiFail", "# of Electrons", phiBin, phiMin, phiMax);
0120
0121   h_convPt_ = ibooker.book1D("convPt", "# of Electrons", ptBin, ptMin, ptMax);
0122   h_convEta_ = ibooker.book1D("convEta", "# of Electrons", etaBin, etaMin, etaMax);
0123   h_convPhi_ = ibooker.book1D("convPhi", "# of Electrons", phiBin, phiMin, phiMax);
0124   h_convRho_ = ibooker.book1D("convRho", "# of Electrons", rhoBin, rhoMin, rhoMax);
0125   h_convZ_ = ibooker.book1D("convZ", "# of Electrons", zBin, zMin, zMax);
0126   h_convProb_ = ibooker.book1D("convProb", "# of Electrons", 100, 0.0, 1.0);
0127
0128   h_convLeadTrackpt_ = ibooker.book1D("convLeadTrackpt", "# of Electrons", trackptBin, trackptMin, trackptMax);
0129
0130   h_convTrailTrackpt_ = ibooker.book1D("convTrailTrackpt", "# of Electrons", trackptBin, trackptMin, trackptMax);
0131
0132   h_convLog10TrailTrackpt_ = ibooker.book1D("convLog10TrailTrackpt", "# of Electrons", ptBin, -2.0, 3.0);
0133
0134   h_convLeadTrackAlgo_ = ibooker.book1D(
0135       "convLeadTrackAlgo", "# of Electrons", reco::TrackBase::algoSize, -0.5, reco::TrackBase::algoSize - 0.5);
0136   h_convTrailTrackAlgo_ = ibooker.book1D(
0137       "convLeadTrackAlgo", "# of Electrons", reco::TrackBase::algoSize, -0.5, reco::TrackBase::algoSize - 0.5);
0138 }
0139
0140 void ElectronConversionRejectionValidator::analyze(const edm::Event& e, const edm::EventSetup&) {
0141   using namespace edm;
0142
0143   nEvt_++;
0144   LogInfo("ElectronConversionRejectionValidator")
0145       << "ElectronConversionRejectionValidator Analyzing event number: " << e.id() << " Global Counter " << nEvt_
0146       << "\n";
0147
0148   ///// Get the recontructed  conversions
0149   auto convHandle = e.getHandle(convToken_);
0150   if (!convHandle.isValid()) {
0151     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the Conversion collection " << std::endl;
0152     return;
0153   }
0154
0155   ///// Get the recontructed  photons
0156   auto gsfElectronHandle = e.getHandle(gsfElecToken_);
0157   const reco::GsfElectronCollection& gsfElectronCollection = *(gsfElectronHandle.product());
0158   if (!gsfElectronHandle.isValid()) {
0159     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the Electron collection " << std::endl;
0160     return;
0161   }
0162
0163   // offline  Primary vertex
0164   auto vertexHandle = e.getHandle(offline_pvToken_);
0165   if (!vertexHandle.isValid()) {
0166     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the product primary Vertex Collection "
0167                                                           << "\n";
0168     return;
0169   }
0170   const reco::Vertex& thevtx = vertexHandle->at(0);
0171
0172   auto bsHandle = e.getHandle(beamspotToken_);
0173   if (!bsHandle.isValid()) {
0174     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the product beamspot Collection "
0175                                                           << "\n";
0176     return;
0177   }
0178   const reco::BeamSpot& thebs = *bsHandle.product();
0179
0180   //loop over electrons
0181   for (reco::GsfElectronCollection::const_iterator iele = gsfElectronCollection.begin();
0182        iele != gsfElectronCollection.end();
0183        ++iele) {
0184     //apply basic pre-selection cuts to remove the conversions with obviously displaced tracks which will anyways be
0185     //removed from the analysis by the hit pattern or impact parameter requirements
0186     if (iele->pt() < elePtMin_)
0187       continue;
0188     if (iele->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) >
0189         eleExpectedHitsInnerMax_)
0190       continue;
0191     if (std::abs(iele->gsfTrack()->dxy(thevtx.position())) > eleD0Max_)
0192       continue;
0193
0194     //fill information for all electrons
0195     h_elePtAll_->Fill(iele->pt());
0196     h_eleEtaAll_->Fill(iele->eta());
0197     h_elePhiAll_->Fill(iele->phi());
0198
0199     //find matching conversion if any
0200     reco::Conversion const* conv = ConversionTools::matchedConversion(*iele, *convHandle, thebs.position());
0201     //fill information on passing electrons only if there is no matching conversion (electron passed the conversion rejection cut!)
0202     if (conv == nullptr) {
0203       h_elePtPass_->Fill(iele->pt());
0204       h_eleEtaPass_->Fill(iele->eta());
0205       h_elePhiPass_->Fill(iele->phi());
0206     } else {
0207       //matching conversion, electron failed conversion rejection cut
0208       //fill information on electron and matching conversion
0209       //(Note that in case of multiple matching conversions passing the requirements, the conversion tools returns the one closest to the IP,
0210       //which is most likely to be the conversion of the primary photon in case there was one.)
0211
0212       //fill electron info
0213       h_elePtFail_->Fill(iele->pt());
0214       h_eleEtaFail_->Fill(iele->eta());
0215       h_elePhiFail_->Fill(iele->phi());
0216
0217       //fill conversion info
0218       math::XYZVectorF convmom = conv->refittedPairMomentum();
0219       h_convPt_->Fill(convmom.rho());
0220       h_convEta_->Fill(convmom.eta());
0221       h_convPhi_->Fill(convmom.phi());
0222       h_convRho_->Fill(conv->conversionVertex().position().rho());
0223       h_convZ_->Fill(conv->conversionVertex().position().z());
0224       h_convProb_->Fill(ChiSquaredProbability(conv->conversionVertex().chi2(), conv->conversionVertex().ndof()));
0225
0226       //fill information about conversion tracks
0227       if (conv->tracks().size() < 2)
0228         continue;
0229
0230       RefToBase<reco::Track> tk1 = conv->tracks().front();
0231       RefToBase<reco::Track> tk2 = conv->tracks().back();
0232
0233       RefToBase<reco::Track> tklead;
0234       RefToBase<reco::Track> tktrail;
0235       if (tk1->pt() >= tk2->pt()) {
0236         tklead = tk1;
0237         tktrail = tk2;
0238       } else {
0239         tklead = tk2;
0240         tktrail = tk1;
0241       }
0242       h_convLeadTrackpt_->Fill(tklead->pt());
0243       h_convTrailTrackpt_->Fill(tktrail->pt());
0244       h_convLog10TrailTrackpt_->Fill(log10(tktrail->pt()));
0245       h_convLeadTrackAlgo_->Fill(tklead->algo());
0246       h_convTrailTrackAlgo_->Fill(tktrail->algo());
0247     }
0248   }
0249 }