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("convLeadTrackAlgo",
0135                                         "# of Electrons",
0136                                         reco::TrackBase::algoSize,
0137                                         -0.5,
0138                                         static_cast<double>(reco::TrackBase::algoSize) - 0.5);
0139   h_convTrailTrackAlgo_ = ibooker.book1D("convLeadTrackAlgo",
0140                                          "# of Electrons",
0141                                          reco::TrackBase::algoSize,
0142                                          -0.5,
0143                                          static_cast<double>(reco::TrackBase::algoSize) - 0.5);
0144 }
0145
0146 void ElectronConversionRejectionValidator::analyze(const edm::Event& e, const edm::EventSetup&) {
0147   using namespace edm;
0148
0149   nEvt_++;
0150   LogInfo("ElectronConversionRejectionValidator")
0151       << "ElectronConversionRejectionValidator Analyzing event number: " << e.id() << " Global Counter " << nEvt_
0152       << "\n";
0153
0154   ///// Get the recontructed  conversions
0155   auto convHandle = e.getHandle(convToken_);
0156   if (!convHandle.isValid()) {
0157     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the Conversion collection " << std::endl;
0158     return;
0159   }
0160
0161   ///// Get the recontructed  photons
0162   auto gsfElectronHandle = e.getHandle(gsfElecToken_);
0163   const reco::GsfElectronCollection& gsfElectronCollection = *(gsfElectronHandle.product());
0164   if (!gsfElectronHandle.isValid()) {
0165     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the Electron collection " << std::endl;
0166     return;
0167   }
0168
0169   // offline  Primary vertex
0170   auto vertexHandle = e.getHandle(offline_pvToken_);
0171   if (!vertexHandle.isValid()) {
0172     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the product primary Vertex Collection "
0173                                                           << "\n";
0174     return;
0175   }
0176   const reco::Vertex& thevtx = vertexHandle->at(0);
0177
0178   auto bsHandle = e.getHandle(beamspotToken_);
0179   if (!bsHandle.isValid()) {
0180     edm::LogError("ElectronConversionRejectionValidator") << "Error! Can't get the product beamspot Collection "
0181                                                           << "\n";
0182     return;
0183   }
0184   const reco::BeamSpot& thebs = *bsHandle.product();
0185
0186   //loop over electrons
0187   for (reco::GsfElectronCollection::const_iterator iele = gsfElectronCollection.begin();
0188        iele != gsfElectronCollection.end();
0189        ++iele) {
0190     //apply basic pre-selection cuts to remove the conversions with obviously displaced tracks which will anyways be
0191     //removed from the analysis by the hit pattern or impact parameter requirements
0192     if (iele->pt() < elePtMin_)
0193       continue;
0194     if (iele->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) >
0195         eleExpectedHitsInnerMax_)
0196       continue;
0197     if (std::abs(iele->gsfTrack()->dxy(thevtx.position())) > eleD0Max_)
0198       continue;
0199
0200     //fill information for all electrons
0201     h_elePtAll_->Fill(iele->pt());
0202     h_eleEtaAll_->Fill(iele->eta());
0203     h_elePhiAll_->Fill(iele->phi());
0204
0205     //find matching conversion if any
0206     reco::Conversion const* conv = ConversionTools::matchedConversion(*iele, *convHandle, thebs.position());
0207     //fill information on passing electrons only if there is no matching conversion (electron passed the conversion rejection cut!)
0208     if (conv == nullptr) {
0209       h_elePtPass_->Fill(iele->pt());
0210       h_eleEtaPass_->Fill(iele->eta());
0211       h_elePhiPass_->Fill(iele->phi());
0212     } else {
0213       //matching conversion, electron failed conversion rejection cut
0214       //fill information on electron and matching conversion
0215       //(Note that in case of multiple matching conversions passing the requirements, the conversion tools returns the one closest to the IP,
0216       //which is most likely to be the conversion of the primary photon in case there was one.)
0217
0218       //fill electron info
0219       h_elePtFail_->Fill(iele->pt());
0220       h_eleEtaFail_->Fill(iele->eta());
0221       h_elePhiFail_->Fill(iele->phi());
0222
0223       //fill conversion info
0224       math::XYZVectorF convmom = conv->refittedPairMomentum();
0225       h_convPt_->Fill(convmom.rho());
0226       h_convEta_->Fill(convmom.eta());
0227       h_convPhi_->Fill(convmom.phi());
0228       h_convRho_->Fill(conv->conversionVertex().position().rho());
0229       h_convZ_->Fill(conv->conversionVertex().position().z());
0230       h_convProb_->Fill(ChiSquaredProbability(conv->conversionVertex().chi2(), conv->conversionVertex().ndof()));
0231
0232       //fill information about conversion tracks
0233       if (conv->tracks().size() < 2)
0234         continue;
0235
0236       RefToBase<reco::Track> tk1 = conv->tracks().front();
0237       RefToBase<reco::Track> tk2 = conv->tracks().back();
0238
0239       RefToBase<reco::Track> tklead;
0240       RefToBase<reco::Track> tktrail;
0241       if (tk1->pt() >= tk2->pt()) {
0242         tklead = tk1;
0243         tktrail = tk2;
0244       } else {
0245         tklead = tk2;
0246         tktrail = tk1;
0247       }
0248       h_convLeadTrackpt_->Fill(tklead->pt());
0249       h_convTrailTrackpt_->Fill(tktrail->pt());
0250       h_convLog10TrailTrackpt_->Fill(log10(tktrail->pt()));
0251       h_convLeadTrackAlgo_->Fill(tklead->algo());
0252       h_convTrailTrackAlgo_->Fill(tktrail->algo());
0253     }
0254   }
0255 }