Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:10:54

0001 #include "DQM/Physics/src/EwkElecDQM.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006 
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0009 
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011 
0012 #include "DataFormats/TrackReco/interface/Track.h"
0013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0014 #include "DataFormats/VertexReco/interface/Vertex.h"
0015 
0016 //#include "DataFormats/MuonReco/interface/Muon.h"
0017 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"  // I guess this is the right one??
0018 // also need Fwd.h file ???
0019 #include "DataFormats/METReco/interface/MET.h"
0020 #include "DataFormats/JetReco/interface/Jet.h"
0021 
0022 #include "DataFormats/GeometryVector/interface/Phi.h"
0023 
0024 #include "FWCore/Common/interface/TriggerNames.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "DataFormats/Common/interface/TriggerResults.h"
0027 
0028 #include "DataFormats/Common/interface/View.h"
0029 
0030 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0031 
0032 using namespace edm;
0033 using namespace std;
0034 using namespace reco;
0035 
0036 EwkElecDQM::EwkElecDQM(const ParameterSet& cfg)
0037     :  // Input collections
0038       metTag_(cfg.getUntrackedParameter<edm::InputTag>("METTag", edm::InputTag("met"))),
0039       jetTag_(cfg.getUntrackedParameter<edm::InputTag>("JetTag", edm::InputTag("sisCone5CaloJets"))),
0040       //      trigTag_(consumes<edm::TriggerResults>(
0041       //        cfg.getUntrackedParameter<edm::InputTag>(
0042       //         "TrigTag", edm::InputTag("TriggerResults::HLT")))),
0043       elecTag_(consumes<edm::View<reco::GsfElectron> >(
0044           cfg.getUntrackedParameter<edm::InputTag>("ElecTag", edm::InputTag("gsfElectrons")))),
0045       metToken_(
0046           consumes<edm::View<reco::MET> >(cfg.getUntrackedParameter<edm::InputTag>("METTag", edm::InputTag("met")))),
0047       jetToken_(consumes<edm::View<reco::Jet> >(
0048           cfg.getUntrackedParameter<edm::InputTag>("JetTag", edm::InputTag("sisCone5CaloJets")))),
0049       vertexTag_(consumes<edm::View<reco::Vertex> >(
0050           cfg.getUntrackedParameter<edm::InputTag>("VertexTag", edm::InputTag("offlinePrimaryVertices")))),
0051       beamSpotTag_(
0052           consumes<reco::BeamSpot>(cfg.getUntrackedParameter<edm::InputTag>("BeamSpotTag", edm::InputTag("BeamSpot")))),
0053 
0054       // Main cuts
0055       //      muonTrig_(cfg.getUntrackedParameter<std::string> ("MuonTrig",
0056       // "HLT_Mu9")),
0057       // elecTrig_(cfg.getUntrackedParameter<std::vector< std::string >
0058       // >("ElecTrig", "HLT_Ele10_SW_L1R")),
0059       elecTrig_(cfg.getUntrackedParameter<std::vector<std::string> >("ElecTrig")),
0060       //      ptCut_(cfg.getUntrackedParameter<double>("PtCut", 25.)),
0061       ptCut_(cfg.getUntrackedParameter<double>("PtCut", 10.)),
0062       //      etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.1)),
0063       etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.4)),
0064       sieieCutBarrel_(cfg.getUntrackedParameter<double>("SieieBarrel", 0.01)),
0065       sieieCutEndcap_(cfg.getUntrackedParameter<double>("SieieEndcap", 0.028)),
0066       detainCutBarrel_(cfg.getUntrackedParameter<double>("DetainBarrel", 0.0071)),
0067       detainCutEndcap_(cfg.getUntrackedParameter<double>("DetainEndcap", 0.0066)),
0068       //      isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso",
0069       // true)),
0070       //      isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso",
0071       // false)),
0072       //      isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
0073       ecalIsoCutBarrel_(cfg.getUntrackedParameter<double>("EcalIsoCutBarrel", 5.7)),
0074       ecalIsoCutEndcap_(cfg.getUntrackedParameter<double>("EcalIsoCutEndcap", 5.0)),
0075       hcalIsoCutBarrel_(cfg.getUntrackedParameter<double>("HcalIsoCutBarrel", 8.1)),
0076       hcalIsoCutEndcap_(cfg.getUntrackedParameter<double>("HcalIsoCutEndcap", 3.4)),
0077       trkIsoCutBarrel_(cfg.getUntrackedParameter<double>("TrkIsoCutBarrel", 7.2)),
0078       trkIsoCutEndcap_(cfg.getUntrackedParameter<double>("TrkIsoCutEndcap", 5.1)),
0079       mtMin_(cfg.getUntrackedParameter<double>("MtMin", -999999)),
0080       mtMax_(cfg.getUntrackedParameter<double>("MtMax", 999999.)),
0081       metMin_(cfg.getUntrackedParameter<double>("MetMin", -999999.)),
0082       metMax_(cfg.getUntrackedParameter<double>("MetMax", 999999.)),
0083       //      acopCut_(cfg.getUntrackedParameter<double>("AcopCut", 2.)),
0084 
0085       // Muon quality cuts
0086       //      dxyCut_(cfg.getUntrackedParameter<double>("DxyCut", 0.2)),
0087       //      normalizedChi2Cut_(cfg.getUntrackedParameter<double>("NormalizedChi2Cut",
0088       // 10.)),
0089       //      trackerHitsCut_(cfg.getUntrackedParameter<int>("TrackerHitsCut",
0090       // 11)),
0091       //      isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon",
0092       // true)),
0093 
0094       // Z rejection
0095       //      ptThrForZ1_(cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
0096       //      ptThrForZ2_(cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
0097 
0098       // Top rejection
0099       eJetMin_(cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
0100       nJetMax_(cfg.getUntrackedParameter<int>("NJetMax", 999999)),
0101       PUMax_(cfg.getUntrackedParameter<unsigned int>("PUMax", 60)),
0102       PUBinCount_(cfg.getUntrackedParameter<unsigned int>("PUBinCount", 12)),
0103       hltPrescaleProvider_(cfg, consumesCollector(), *this)
0104 //       caloJetCollection_(cfg.getUntrackedParameter<edm:InputTag>("CaloJetCollection","sisCone5CaloJets"))
0105 {
0106   isValidHltConfig_ = false;
0107 }
0108 
0109 void EwkElecDQM::dqmBeginRun(const Run& iRun, const EventSetup& iSet) {
0110   nall = 0;
0111   nsel = 0;
0112 
0113   nrec = 0;
0114   neid = 0;
0115   niso = 0;
0116   //       nhlt = 0;
0117   //       nmet = 0;
0118 
0119   // passed as parameter to HLTConfigProvider::init(), not yet used
0120   bool isConfigChanged = false;
0121   // isValidHltConfig_ could be used to short-circuit analyze() in case of
0122   // problems
0123   isValidHltConfig_ = hltPrescaleProvider_.init(iRun, iSet, "HLT", isConfigChanged);
0124 
0125   LogTrace("") << "isValidHltConfig_=" << isValidHltConfig_ << "\n";
0126 }
0127 
0128 void EwkElecDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0129   ibooker.setCurrentFolder("Physics/EwkElecDQM");
0130 
0131   char chtitle[256] = "";
0132 
0133   pt_before_ = ibooker.book1D("PT_BEFORECUTS", "Electron transverse momentum [GeV]", 100, 0., 100.);
0134   pt_after_ = ibooker.book1D("PT_LASTCUT", "Electron transverse momentum [GeV]", 100, 0., 100.);
0135 
0136   eta_before_ = ibooker.book1D("ETA_BEFORECUTS", "Electron pseudo-rapidity", 50, -2.5, 2.5);
0137   eta_after_ = ibooker.book1D("ETA_LASTCUT", "Electron pseudo-rapidity", 50, -2.5, 2.5);
0138 
0139   sieiebarrel_before_ = ibooker.book1D("SIEIEBARREL_BEFORECUTS", "Electron #sigma_{i#etai#eta} (barrel)", 70, 0., 0.07);
0140   sieiebarrel_after_ = ibooker.book1D("SIEIEBARREL_LASTCUT", "Electron #sigma_{i#etai#eta} (barrel)", 70, 0., 0.07);
0141 
0142   sieieendcap_before_ = ibooker.book1D("SIEIEENDCAP_BEFORECUTS", "Electron #sigma_{i#etai#eta} (endcap)", 70, 0., 0.07);
0143   sieieendcap_after_ = ibooker.book1D("SIEIEENDCAP_LASTCUT", "Electron #sigma_{i#etai#eta} (endcap)", 70, 0., 0.07);
0144 
0145   detainbarrel_before_ =
0146       ibooker.book1D("DETAINBARREL_BEFORECUTS", "Electron #Delta#eta_{in} (barrel)", 40, -0.02, 0.02);
0147   detainbarrel_after_ = ibooker.book1D("DETAINBARREL_LASTCUT", "Electron #Delta#eta_{in} (barrel)", 40, -0.02, 0.02);
0148 
0149   detainendcap_before_ =
0150       ibooker.book1D("DETAINENDCAP_BEFORECUTS", "Electron #Delta#eta_{in} (endcap)", 40, -0.02, 0.02);
0151   detainendcap_after_ = ibooker.book1D("DETAINENDCAP_LASTCUT", "Electron #Delta#eta_{in} (endcap)", 40, -0.02, 0.02);
0152 
0153   ecalisobarrel_before_ = ibooker.book1D(
0154       "ECALISOBARREL_BEFORECUTS", "Absolute electron ECAL isolation variable (barrel) [GeV]", 50, 0., 50.);
0155   ecalisobarrel_after_ =
0156       ibooker.book1D("ECALISOBARREL_LASTCUT", "Absolute electron ECAL isolation variable (barrel) [GeV]", 50, 0., 50.);
0157 
0158   ecalisoendcap_before_ = ibooker.book1D(
0159       "ECALISOENDCAP_BEFORECUTS", "Absolute electron ECAL isolation variable (endcap) [GeV]", 50, 0., 50.);
0160   ecalisoendcap_after_ =
0161       ibooker.book1D("ECALISOENDCAP_LASTCUT", "Absolute electron ECAL isolation variable (endcap) [GeV]", 50, 0., 50.);
0162 
0163   hcalisobarrel_before_ = ibooker.book1D(
0164       "HCALISOBARREL_BEFORECUTS", "Absolute electron HCAL isolation variable (barrel) [GeV]", 50, 0., 50.);
0165   hcalisobarrel_after_ =
0166       ibooker.book1D("HCALISOBARREL_LASTCUT", "Absolute electron HCAL isolation variable (barrel) [GeV]", 50, 0., 50.);
0167 
0168   hcalisoendcap_before_ = ibooker.book1D(
0169       "HCALISOENDCAP_BEFORECUTS", "Absolute electron HCAL isolation variable (endcap) [GeV]", 50, 0., 50.);
0170   hcalisoendcap_after_ =
0171       ibooker.book1D("HCALISOENDCAP_LASTCUT", "Absolute electron HCAL isolation variable (endcap) [GeV]", 50, 0., 50.);
0172 
0173   trkisobarrel_before_ = ibooker.book1D(
0174       "TRKISOBARREL_BEFORECUTS", "Absolute electron track isolation variable (barrel) [GeV]", 50, 0., 50.);
0175   trkisobarrel_after_ =
0176       ibooker.book1D("TRKISOBARREL_LASTCUT", "Absolute electron track isolation variable (barrel) [GeV]", 50, 0., 50.);
0177 
0178   trkisoendcap_before_ = ibooker.book1D(
0179       "TRKISOENDCAP_BEFORECUTS", "Absolute electron track isolation variable (endcap) [GeV]", 50, 0., 50.);
0180   trkisoendcap_after_ =
0181       ibooker.book1D("TRKISOENDCAP_LASTCUT", "Absolute electron track isolation variable (endcap) [GeV]", 50, 0., 50.);
0182 
0183   // trig_before_ = ibooker.book1D("TRIG_BEFORECUTS", "Trigger response", 2, -0.5,
0184   /// 1.5);  // elecTrig_ is now a vector of strings!
0185   //  trig_after_ = ibooker.book1D("TRIG_LASTCUT", "Trigger response", 2, -0.5, 1.5);
0186 
0187   invmass_before_ = ibooker.book1D("INVMASS_BEFORECUTS", "Di-electron invariant mass [GeV]", 100, 0., 200.);
0188   invmass_after_ = ibooker.book1D("INVMASS_AFTERCUTS", "Di-electron invariant mass [GeV]", 100, 0., 200.);
0189 
0190   invmassPU_before_ = ibooker.book2D("INVMASS_PU_BEFORECUTS",
0191                                      "Di-electron invariant mass [GeV] vs PU; mass [GeV]; PU count",
0192                                      100,
0193                                      0.,
0194                                      200.,
0195                                      PUBinCount_,
0196                                      -0.5,
0197                                      PUMax_ + 0.5);
0198   invmassPU_afterZ_ = ibooker.book2D("INVMASS_PU_AFTERZCUTS",
0199                                      "Di-electron invariant mass [GeV] vs PU; mass [GeV]; PU count",
0200                                      100,
0201                                      0.,
0202                                      200.,
0203                                      PUBinCount_,
0204                                      -0.5,
0205                                      PUMax_ + 0.5);
0206 
0207   npvs_before_ =
0208       ibooker.book1D("NPVs_BEFORECUTS", "Number of Valid Primary Vertices; nGoodPVs", PUMax_ + 1, -0.5, PUMax_ + 0.5);
0209 
0210   npvs_afterZ_ =
0211       ibooker.book1D("NPVs_AFTERZCUTS", "Number of Valid Primary Vertices; nGoodPVs", PUMax_ + 1, -0.5, PUMax_ + 0.5);
0212 
0213   nelectrons_before_ = ibooker.book1D("NELECTRONS_BEFORECUTS", "Number of electrons in event", 10, -0.5, 9.5);
0214   nelectrons_after_ = ibooker.book1D("NELECTRONS_AFTERCUTS", "Number of electrons in event", 10, -0.5, 9.5);
0215 
0216   snprintf(chtitle, 255, "Transverse mass (%s) [GeV]", metTag_.label().data());
0217   mt_before_ = ibooker.book1D("MT_BEFORECUTS", chtitle, 150, 0., 300.);
0218   mt_after_ = ibooker.book1D("MT_LASTCUT", chtitle, 150, 0., 300.);
0219 
0220   snprintf(chtitle, 255, "Missing transverse energy (%s) [GeV]", metTag_.label().data());
0221   met_before_ = ibooker.book1D("MET_BEFORECUTS", chtitle, 100, 0., 200.);
0222   met_after_ = ibooker.book1D("MET_LASTCUT", chtitle, 100, 0., 200.);
0223 
0224   snprintf(chtitle, 255, "Number of jets (%s) above %.2f GeV", jetTag_.label().data(), eJetMin_);
0225   njets_before_ = ibooker.book1D("NJETS_BEFORECUTS", chtitle, 10, -0.5, 9.5);
0226   njets_after_ = ibooker.book1D("NJETS_LASTCUT", chtitle, 10, -0.5, 9.5);
0227 
0228   snprintf(chtitle, 255, "Jet with highest E_{T} (%s)", jetTag_.label().data());
0229   jet_et_before_ = ibooker.book1D("JETET1_BEFORECUTS", chtitle, 20, 0., 200.0);
0230   jet_et_after_ = ibooker.book1D("JETET1_AFTERCUTS", chtitle, 20, 0., 200.0);
0231 
0232   snprintf(chtitle, 255, "Eta of Jet with highest E_{T} (%s)", jetTag_.label().data());
0233   jet_eta_before_ = ibooker.book1D("JETETA1_BEFORECUTS", chtitle, 20, -5, 5);
0234   jet_eta_after_ = ibooker.book1D("JETETA1_AFTERCUTS", chtitle, 20, -5, 5);
0235 }
0236 
0237 void EwkElecDQM::dqmEndRun(const Run& r, const EventSetup&) {
0238   // overall
0239   double all = nall;
0240   double esel = nsel / all;
0241   LogVerbatim("") << "\n>>>>>> SELECTION SUMMARY BEGIN >>>>>>>>>>>>>>>";
0242   LogVerbatim("") << "Total number of events analyzed: " << nall << " [events]";
0243   LogVerbatim("") << "Total number of events selected: " << nsel << " [events]";
0244   LogVerbatim("") << "Overall efficiency:             "
0245                   << "(" << setprecision(4) << esel * 100. << " +/- " << setprecision(2)
0246                   << sqrt(esel * (1 - esel) / all) * 100. << ")%";
0247 
0248   double erec = nrec / all;
0249   double eeid = neid / all;
0250   double eiso = niso / all;
0251   //   double ehlt = nhlt/all;
0252   //   double emet = nmet/all;
0253 
0254   // general reconstruction step??
0255   double num = nrec;
0256   double eff = erec;
0257   double err = sqrt(eff * (1 - eff) / all);
0258   LogVerbatim("") << "Passing Pt/Eta/Quality cuts:    " << num << " [events], (" << setprecision(4) << eff * 100.
0259                   << " +/- " << setprecision(2) << err * 100. << ")%";
0260 
0261   // electron ID step
0262   num = neid;
0263   eff = eeid;
0264   err = sqrt(eff * (1 - eff) / all);
0265   double effstep = 0.;
0266   double errstep = 0.;
0267   if (nrec > 0)
0268     effstep = eeid / erec;
0269   if (nrec > 0)
0270     errstep = sqrt(effstep * (1 - effstep) / nrec);
0271   LogVerbatim("") << "Passing eID cuts:         " << num << " [events], (" << setprecision(4) << eff * 100. << " +/- "
0272                   << setprecision(2) << err * 100. << ")%, to previous step: (" << setprecision(4) << effstep * 100.
0273                   << " +/- " << setprecision(2) << errstep * 100. << ")%";
0274 
0275   // isolation step
0276   num = niso;
0277   eff = eiso;
0278   err = sqrt(eff * (1 - eff) / all);
0279   effstep = 0.;
0280   errstep = 0.;
0281   if (neid > 0)
0282     effstep = eiso / eeid;
0283   if (neid > 0)
0284     errstep = sqrt(effstep * (1 - effstep) / neid);
0285   LogVerbatim("") << "Passing isolation cuts:         " << num << " [events], (" << setprecision(4) << eff * 100.
0286                   << " +/- " << setprecision(2) << err * 100. << ")%, to previous step: (" << setprecision(4)
0287                   << effstep * 100. << " +/- " << setprecision(2) << errstep * 100. << ")%";
0288 
0289   //   // trigger step
0290   //   num = nhlt;
0291   //   eff = ehlt;
0292   //   err = sqrt(eff*(1-eff)/all);
0293   //   effstep = 0.;
0294   //   errstep = 0.;
0295   //   if (niso>0) effstep = ehlt/eiso;
0296   //   if (niso>0) errstep = sqrt(effstep*(1-effstep)/niso);
0297   //   LogVerbatim("") << "Passing HLT criteria:           " << num << "
0298   // [events], (" << setprecision(4) << eff*100. <<" +/- "<< setprecision(2) <<
0299   // err*100. << ")%, to previous step: (" <<  setprecision(4) << effstep*100.
0300   // << " +/- "<< setprecision(2) << errstep*100. <<")%";
0301 
0302   // trigger step
0303   num = nsel;
0304   eff = esel;
0305   err = sqrt(eff * (1 - eff) / all);
0306   effstep = 0.;
0307   errstep = 0.;
0308   if (niso > 0)
0309     effstep = esel / eiso;
0310   if (niso > 0)
0311     errstep = sqrt(effstep * (1 - effstep) / niso);
0312   LogVerbatim("") << "Passing HLT criteria:           " << num << " [events], (" << setprecision(4) << eff * 100.
0313                   << " +/- " << setprecision(2) << err * 100. << ")%, to previous step: (" << setprecision(4)
0314                   << effstep * 100. << " +/- " << setprecision(2) << errstep * 100. << ")%";
0315 
0316   //   // met/acoplanarity cuts
0317   //   num = nmet;
0318   //   eff = emet;
0319   //   err = sqrt(eff*(1-eff)/all);
0320   //   effstep = 0.;
0321   //   errstep = 0.;
0322   //   if (nhlt>0) effstep = emet/ehlt;
0323   //   if (nhlt>0) errstep = sqrt(effstep*(1-effstep)/nhlt);
0324   //   LogVerbatim("") << "Passing MET/acoplanarity cuts:  " << num << "
0325   // [events], (" << setprecision(4) << eff*100. <<" +/- "<< setprecision(2) <<
0326   // err*100. << ")%, to previous step: (" <<  setprecision(4) << effstep*100.
0327   // << " +/- "<< setprecision(2) << errstep*100. <<")%";
0328 
0329   //   // Z/top selection cuts ALSO LAST STEP so "sel" for "selection"
0330   //   num = nsel;
0331   //   eff = esel;
0332   //   err = sqrt(eff*(1-eff)/all);
0333   //   effstep = 0.;
0334   //   errstep = 0.;
0335   //   if (nmet>0) effstep = esel/emet;
0336   //   if (nmet>0) errstep = sqrt(effstep*(1-effstep)/nmet);
0337   //   LogVerbatim("") << "Passing Z/top rejection cuts:   " << num << "
0338   // [events], (" << setprecision(4) << eff*100. <<" +/- "<< setprecision(2) <<
0339   // err*100. << ")%, to previous step: (" <<  setprecision(4) << effstep*100.
0340   // << " +/- "<< setprecision(2) << errstep*100. <<")%";
0341 
0342   LogVerbatim("") << ">>>>>> SELECTION SUMMARY END   >>>>>>>>>>>>>>>\n";
0343 }
0344 
0345 inline void HERE(const char* msg) { std::cout << msg << "\n"; }
0346 
0347 void EwkElecDQM::analyze(const Event& ev, const EventSetup& iSet) {
0348   // Reset global event selection flags
0349   bool rec_sel = false;
0350   bool eid_sel = false;
0351   bool iso_sel = false;
0352   bool all_sel = false;
0353 
0354   // Electron collection
0355   Handle<View<GsfElectron> > electronCollection;
0356   if (!ev.getByToken(elecTag_, electronCollection)) {
0357     // LogWarning("") << ">>> Electron collection does not exist !!!";
0358     return;
0359   }
0360   unsigned int electronCollectionSize = electronCollection->size();
0361 
0362   // Beam spot
0363   Handle<reco::BeamSpot> beamSpotHandle;
0364   if (!ev.getByToken(beamSpotTag_, beamSpotHandle)) {
0365     // LogWarning("") << ">>> No beam spot found !!!";
0366     return;
0367   }
0368 
0369   // MET
0370   double met_px = 0.;
0371   double met_py = 0.;
0372   Handle<View<MET> > metCollection;
0373   if (!ev.getByToken(metToken_, metCollection)) {
0374     // LogWarning("") << ">>> MET collection does not exist !!!";
0375     return;
0376   }
0377 
0378   const MET& met = metCollection->at(0);
0379   met_px = met.px();
0380   met_py = met.py();
0381   //       if (!metIncludesMuons_) {
0382   //             for (unsigned int i=0; i<muonCollectionSize; i++) {
0383   //                   const Muon& mu = muonCollection->at(i);
0384   //                   if (!mu.isGlobalMuon()) continue;
0385   //                   met_px -= mu.px();
0386   //                   met_py -= mu.py();
0387   //             }
0388   //       }
0389   double met_et = sqrt(met_px * met_px + met_py * met_py);
0390   LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met_px << ", " << met_py << " [GeV]";
0391   met_before_->Fill(met_et);
0392 
0393   // Vertices in the event
0394   int npvCount = 0;
0395   Handle<View<reco::Vertex> > vertexCollection;
0396   if (!ev.getByToken(vertexTag_, vertexCollection)) {
0397     LogError("") << ">>> Vertex collection does not exist !!!";
0398     return;
0399   }
0400 
0401   for (unsigned int i = 0; i < vertexCollection->size(); i++) {
0402     const Vertex& vertex = vertexCollection->at(i);
0403     if (vertex.isValid())
0404       npvCount++;
0405   }
0406   npvs_before_->Fill(npvCount);
0407 
0408   // Trigger
0409   // Handle<TriggerResults> triggerResults;
0410   // if (!ev.getByToken(trigTag_, triggerResults)) {
0411   // LogWarning("") << ">>> TRIGGER collection does not exist !!!";
0412   return;
0413   // }
0414   // const edm::TriggerNames& trigNames = ev.triggerNames(*triggerResults);
0415   // bool trigger_fired = false;
0416 
0417   // HLTConfigProvider const&  hltConfigProvider = hltPrescaleProvider_.hltConfigProvider();
0418 
0419   /* very old code
0420   for (unsigned int i=0; i<triggerResults->size(); i++) {
0421         if (triggerResults->accept(i)) {
0422               LogTrace("") << "Accept by: " << i << ", Trigger: " <<
0423   trigNames.triggerName(i);
0424         }
0425   }
0426   */
0427 
0428   // the following gives error on CRAFT08 data where itrig1=19 (vector index out
0429   // of range)
0430   /*
0431   int itrig1 = trigNames.triggerIndex(muonTrig_);
0432   if (triggerResults->accept(itrig1)) trigger_fired = true;
0433   */
0434   // suggested replacement: lm250909
0435   /*  Fix buggy trigger logic
0436   for (unsigned int i=0; i<triggerResults->size(); i++)
0437     {
0438       std::string trigName = trigNames.triggerName(i);
0439       bool found=false;
0440 
0441 //    for (unsigned int j = 0; j < elecTrig_.size(); j++)
0442 //      {
0443 //        if ( trigName == elecTrig_.at(j) && triggerResults->accept(i))
0444 //      {
0445 //        trigger_fired = true;
0446 //      }
0447 //      }
0448 
0449 
0450      for(unsigned int index=0; index<elecTrig_.size() && found==false; index++)
0451 {
0452        size_t trigPath = trigName.find(elecTrig_.at(index)); // 0 if found, pos
0453 if not
0454                if (trigPath==0) found=true;
0455           }
0456           if(!found) continue;
0457 
0458           bool prescaled=false;
0459           for (unsigned int ps= 0; ps<  hltConfigProvider.prescaleSize();
0460 ps++){
0461               const unsigned int prescaleValue =
0462 hltConfigProvider.prescaleValue(ps, trigName) ;
0463               if (prescaleValue != 1) prescaled =true;
0464           }
0465 
0466           if(triggerResults->accept(i) && !prescaled) trigger_fired=true;
0467 
0468     }
0469   */
0470 
0471   // get the prescale set for this event
0472   const int prescaleSet = hltPrescaleProvider_.prescaleSet(ev, iSet);
0473   if (prescaleSet == -1) {
0474     LogTrace("") << "Failed to determine prescaleSet\n";
0475     // std::cout << "Failed to determine prescaleSet. Check cmsRun GlobalTag\n";
0476     return;
0477   }
0478 
0479   // for (unsigned int i = 0;
0480   //    (i < triggerResults->size()) && (trigger_fired == false); i++) {
0481   // skip trigger, if it did not fire
0482   //  if (!triggerResults->accept(i)) continue;
0483 
0484   // skip trigger, if it is not on our list
0485   // bool found = false;
0486   // const std::string trigName = trigNames.triggerName(i);
0487   // for (unsigned int index = 0; index < elecTrig_.size() && found == false;
0488   //    index++) {
0489   // if (trigName.find(elecTrig_.at(index)) == 0) found = true;
0490   // }
0491   // if (!found) continue;
0492 
0493   // skip trigger, if it is prescaled
0494   // if (hltConfigProvider.prescaleValue(prescaleSet, trigName) != 1) continue;
0495 
0496   // std::cout << "found unprescaled trigger that fired: " << trigName <<
0497   // "\n";
0498   // trigger_fired = true;
0499   // }
0500 
0501   /*  LogTrace("") << ">>> Trigger bit: " << trigger_fired << " for one of ( ";
0502   for (unsigned int k = 0; k < elecTrig_.size(); k++) {
0503     LogTrace("") << elecTrig_.at(k) << " ";
0504   }
0505   LogTrace("") << ")";
0506   trig_before_->Fill(trigger_fired);
0507 */
0508   // Jet collection
0509   Handle<View<Jet> > jetCollection;
0510   if (!ev.getByToken(jetToken_, jetCollection)) {
0511     // LogError("") << ">>> JET collection does not exist !!!";
0512     return;
0513   }
0514   float electron_et = -8.0;
0515   float electron_eta = -8.0;
0516   float electron_phi = -8.0;
0517   float electron2_et = -9.0;
0518   float electron2_eta = -9.0;
0519   float electron2_phi = -9.0;
0520   // need to get some electron info so jets can be cleaned of them
0521   for (unsigned int i = 0; i < electronCollectionSize; i++) {
0522     const GsfElectron& elec = electronCollection->at(i);
0523 
0524     if (i < 1) {
0525       electron_et = elec.pt();
0526       electron_eta = elec.eta();
0527       electron_phi = elec.phi();
0528     }
0529     if (i == 2) {
0530       electron2_et = elec.pt();
0531       electron2_eta = elec.eta();
0532       electron2_phi = elec.phi();
0533     }
0534   }
0535 
0536   float jet_et = -8.0;
0537   float jet_eta = -8.0;
0538   int jet_count = 0;
0539   float jet2_et = -9.0;
0540   unsigned int jetCollectionSize = jetCollection->size();
0541   int njets = 0;
0542   for (unsigned int i = 0; i < jetCollectionSize; i++) {
0543     const Jet& jet = jetCollection->at(i);
0544 
0545     float jet_current_et = jet.et();
0546     //  cout << "jet_current_et " << jet_current_et << endl;
0547     // if it overlaps with electron, it is not a jet
0548     if (electron_et > 0.0 && fabs(jet.eta() - electron_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron_phi) < 0.2)
0549       continue;
0550     if (electron2_et > 0.0 && fabs(jet.eta() - electron2_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron2_phi) < 0.2)
0551       continue;
0552 
0553     // if it has too low Et, throw away
0554     //  if (jet_current_et < eJetMin_) continue; //Keep if only want to plot
0555     // above jet cut
0556 
0557     if (jet.et() > eJetMin_) {
0558       njets++;
0559       jet_count++;
0560     }
0561     if (jet_current_et > jet_et) {
0562       jet2_et = jet_et;   // 2nd highest jet get's et from current highest
0563       jet_et = jet.et();  // current highest jet gets et from the new highest
0564       jet_eta = jet.eta();
0565     } else if (jet_current_et > jet2_et) {
0566       jet2_et = jet.et();
0567     }
0568   }
0569 
0570   // Fill After all electron cuts (or both before and after)
0571   if (jet_et > 10)  // don't want low energy "jets"
0572   {
0573     jet_et_before_->Fill(jet_et);
0574     //    jet2_et_before_  ->Fill(jet2_et);
0575     jet_eta_before_->Fill(jet_eta);
0576   }
0577 
0578   LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
0579   LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
0580   njets_before_->Fill(njets);
0581 
0582   // Start counting
0583   nall++;
0584 
0585   // Histograms per event should be done only once, so keep track of them
0586   //bool hlt_hist_done = false;
0587   // bool minv_hist_done = false;
0588   bool met_hist_done = false;
0589   //       bool nz1_hist_done = false;
0590   //       bool nz2_hist_done = false;
0591   bool njets_hist_done = false;
0592 
0593   // Central selection criteria
0594   // const int NFLAGS = 13; // number of individual selection criteria
0595   const int NFLAGS = 10;  // number of individual selection criteria
0596   // 0: pt cut           | rec
0597   // 1: eta cut          | rec
0598   // 2: sieie            | eid
0599   // 3: detain           | eid
0600   // 4: ecal iso         | iso
0601   // 5: hcal iso         | iso
0602   // 6: trk iso          | iso
0603   // 7: trigger fired    | hlt/all
0604   bool electron_sel[NFLAGS];
0605 
0606   // for invariant mass calculation
0607   // keep track of highest-pt electrons for initial (RECO) electrons
0608   // and "good" electrons (passing all cuts)
0609   // first array dimension is for first or second good electron
0610   // second array dimension is for relevant quantities of good electron
0611   //    [0]: 1 for electron found or 0 for not found (if 0, no other quantities
0612   // filled)
0613   //    [1]: mSqr
0614   //    [2]: E
0615   //    [3]: px
0616   //    [4]: py
0617   //    [5]: pz
0618   // inv mass = sqrt(m_1^2 + m_2^2 + 2*(E_1*E_2 - (px1*px2 + py1*py2 + pz1+pz2)
0619   // ) )
0620   double electron[2][6];
0621   double goodElectron[2][6];
0622   nGoodElectrons = 0;
0623   for (unsigned int i = 0; i < 2; i++) {
0624     for (unsigned int j = 0; j < 6; j++) {
0625       electron[i][j] = 0.;
0626       goodElectron[i][j] = 0.;
0627     }
0628   }
0629 
0630   for (unsigned int i = 0; i < electronCollectionSize; i++) {
0631     for (int j = 0; j < NFLAGS; ++j) {
0632       electron_sel[j] = false;
0633     }
0634 
0635     const GsfElectron& elec = electronCollection->at(i);
0636     // if (!mu.isGlobalMuon()) continue;
0637     // if (mu.globalTrack().isNull()) continue;
0638     // if (mu.innerTrack().isNull()) continue;
0639 
0640     LogTrace("") << "> elec: processing electron number " << i << "...";
0641     // reco::TrackRef gm = mu.globalTrack();
0642     // reco::TrackRef tk = mu.innerTrack();
0643     // should have stuff for electron track?
0644 
0645     if (i < 2) {
0646       electron[i][0] = 1.;
0647       electron[i][1] = elec.massSqr();
0648       electron[i][2] = elec.energy();
0649       electron[i][3] = elec.px();
0650       electron[i][4] = elec.py();
0651       electron[i][5] = elec.pz();
0652     }
0653 
0654     // Pt,eta cuts
0655     double pt = elec.pt();
0656     double px = elec.px();
0657     double py = elec.py();
0658     double eta = elec.eta();
0659     LogTrace("") << "\t... pt, eta: " << pt << " [GeV], " << eta;
0660     ;
0661     if (pt > ptCut_)
0662       electron_sel[0] = true;
0663     if (fabs(eta) < etaCut_)
0664       electron_sel[1] = true;
0665 
0666     bool isBarrel = false;
0667     bool isEndcap = false;
0668     if (eta < 1.4442 && eta > -1.4442) {
0669       isBarrel = true;
0670     } else if ((eta > 1.56 && eta < 2.4) || (eta < -1.56 && eta > -2.4)) {
0671       isEndcap = true;
0672     }
0673 
0674     //             // d0, chi2, nhits quality cuts
0675     //             double dxy = tk->dxy(beamSpotHandle->position());
0676     //             double normalizedChi2 = gm->normalizedChi2();
0677     //             double trackerHits = tk->numberOfValidHits();
0678     //             LogTrace("") << "\t... dxy, normalizedChi2, trackerHits,
0679     // isTrackerMuon?: " << dxy << " [cm], " << normalizedChi2 << ", " <<
0680     // trackerHits << ", " << mu.isTrackerMuon();
0681     //             if (fabs(dxy)<dxyCut_) muon_sel[2] = true;
0682     //             if (normalizedChi2<normalizedChi2Cut_) muon_sel[3] = true;
0683     //             if (trackerHits>=trackerHitsCut_) muon_sel[4] = true;
0684     //             if (mu.isTrackerMuon()) muon_sel[5] = true;
0685 
0686     pt_before_->Fill(pt);
0687     eta_before_->Fill(eta);
0688     //             dxy_before_->Fill(dxy);
0689     //             chi2_before_->Fill(normalizedChi2);
0690     //             nhits_before_->Fill(trackerHits);
0691     //             tkmu_before_->Fill(mu.isTrackerMuon());
0692 
0693     // Electron ID cuts
0694     double sieie = (double)elec.sigmaIetaIeta();
0695     double detain = (double)elec.deltaEtaSuperClusterTrackAtVtx();  // think this is detain
0696     if (sieie < sieieCutBarrel_ && isBarrel)
0697       electron_sel[2] = true;
0698     if (sieie < sieieCutEndcap_ && isEndcap)
0699       electron_sel[2] = true;
0700     if (detain < detainCutBarrel_ && isBarrel)
0701       electron_sel[3] = true;
0702     if (detain < detainCutEndcap_ && isEndcap)
0703       electron_sel[3] = true;
0704     if (isBarrel) {
0705       LogTrace("") << "\t... sieie value " << sieie << " (barrel), pass? " << electron_sel[2];
0706       LogTrace("") << "\t... detain value " << detain << " (barrel), pass? " << electron_sel[3];
0707     } else if (isEndcap) {
0708       LogTrace("") << "\t... sieie value " << sieie << " (endcap), pass? " << electron_sel[2];
0709       LogTrace("") << "\t... detain value " << detain << " (endcap), pass? " << electron_sel[2];
0710     }
0711 
0712     if (isBarrel) {
0713       sieiebarrel_before_->Fill(sieie);
0714       detainbarrel_before_->Fill(detain);
0715     } else if (isEndcap) {
0716       sieieendcap_before_->Fill(sieie);
0717       detainendcap_before_->Fill(detain);
0718     }
0719 
0720     // Isolation cuts
0721     // double isovar = mu.isolationR03().sumPt;
0722     double ecalisovar = elec.dr03EcalRecHitSumEt();  // picked one set!
0723     double hcalisovar = elec.dr03HcalTowerSumEt();   // try others if
0724     double trkisovar = elec.dr04TkSumPt();           // doesn't work
0725     // if (isCombinedIso_) {
0726     // isovar += mu.isolationR03().emEt;
0727     // isovar += mu.isolationR03().hadEt;
0728     //}
0729     // if (isRelativeIso_) isovar /= pt;
0730     if (ecalisovar < ecalIsoCutBarrel_ && isBarrel)
0731       electron_sel[4] = true;
0732     if (ecalisovar < ecalIsoCutEndcap_ && isEndcap)
0733       electron_sel[4] = true;
0734     if (hcalisovar < hcalIsoCutBarrel_ && isBarrel)
0735       electron_sel[5] = true;
0736     if (hcalisovar < hcalIsoCutEndcap_ && isEndcap)
0737       electron_sel[5] = true;
0738     if (trkisovar < trkIsoCutBarrel_ && isBarrel)
0739       electron_sel[6] = true;
0740     if (trkisovar < trkIsoCutEndcap_ && isEndcap)
0741       electron_sel[6] = true;
0742     if (isBarrel) {
0743       LogTrace("") << "\t... ecal isolation value " << ecalisovar << " (barrel), pass? " << electron_sel[4];
0744       LogTrace("") << "\t... hcal isolation value " << hcalisovar << " (barrel), pass? " << electron_sel[5];
0745       LogTrace("") << "\t... trk isolation value " << trkisovar << " (barrel), pass? " << electron_sel[6];
0746     } else if (isEndcap) {
0747       LogTrace("") << "\t... ecal isolation value " << ecalisovar << " (endcap), pass? " << electron_sel[4];
0748       LogTrace("") << "\t... hcal isolation value " << hcalisovar << " (endcap), pass? " << electron_sel[5];
0749       LogTrace("") << "\t... trk isolation value " << trkisovar << " (endcap), pass? " << electron_sel[6];
0750     }
0751 
0752     // iso_before_->Fill(isovar);
0753     if (isBarrel) {
0754       ecalisobarrel_before_->Fill(ecalisovar);
0755       hcalisobarrel_before_->Fill(hcalisovar);
0756       trkisobarrel_before_->Fill(trkisovar);
0757     } else if (isEndcap) {
0758       ecalisoendcap_before_->Fill(ecalisovar);
0759       hcalisoendcap_before_->Fill(hcalisovar);
0760       trkisoendcap_before_->Fill(trkisovar);
0761     }
0762 
0763     // HLT
0764     // if (trigger_fired) electron_sel[7] = true;
0765 
0766     //             // MET/MT cuts
0767     double w_et = met_et + pt;
0768     double w_px = met_px + px;
0769     double w_py = met_py + py;
0770 
0771     double massT = w_et * w_et - w_px * w_px - w_py * w_py;
0772     massT = (massT > 0) ? sqrt(massT) : 0;
0773 
0774     LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << w_px << ", " << w_py
0775                  << " [GeV]";
0776     if (massT > mtMin_ && massT < mtMax_)
0777       electron_sel[7] = true;
0778     mt_before_->Fill(massT);
0779     if (met_et > metMin_ && met_et < metMax_)
0780       electron_sel[8] = true;
0781 
0782     //             // Acoplanarity cuts
0783     //             Geom::Phi<double> deltaphi(mu.phi()-atan2(met_py,met_px));
0784     //             double acop = deltaphi.value();
0785     //             if (acop<0) acop = - acop;
0786     //             acop = M_PI - acop;
0787     //             LogTrace("") << "\t... acoplanarity: " << acop;
0788     //             if (acop<acopCut_) muon_sel[10] = true;
0789     //             acop_before_->Fill(acop);
0790 
0791     //             // Remaining flags (from global event information)
0792     //             if (nmuonsForZ1<1 || nmuonsForZ2<2) muon_sel[11] = true;
0793     if (njets <= nJetMax_)
0794       electron_sel[9] = true;
0795 
0796     // Collect necessary flags "per electron"
0797     int flags_passed = 0;
0798     bool rec_sel_this = true;
0799     bool eid_sel_this = true;
0800     bool iso_sel_this = true;
0801     bool all_sel_this = true;
0802     for (int j = 0; j < NFLAGS; ++j) {
0803       if (electron_sel[j])
0804         flags_passed += 1;
0805       if (j < 2 && !electron_sel[j])
0806         rec_sel_this = false;
0807       if (j < 4 && !electron_sel[j])
0808         eid_sel_this = false;
0809       if (j < 7 && !electron_sel[j])
0810         iso_sel_this = false;
0811       if (!electron_sel[j])
0812         all_sel_this = false;
0813     }
0814 
0815     if (all_sel_this) {
0816       if (nGoodElectrons < 2) {
0817         goodElectron[nGoodElectrons][0] = 1.;
0818         goodElectron[nGoodElectrons][1] = elec.massSqr();
0819         goodElectron[nGoodElectrons][2] = elec.energy();
0820         goodElectron[nGoodElectrons][3] = elec.px();
0821         goodElectron[nGoodElectrons][4] = elec.py();
0822         goodElectron[nGoodElectrons][5] = elec.pz();
0823       }
0824       nGoodElectrons++;
0825     }
0826 
0827     //             // "rec" => pt,eta and quality cuts are satisfied
0828     //             if (rec_sel_this) rec_sel = true;
0829     //             // "iso" => "rec" AND "muon is isolated"
0830     //             if (iso_sel_this) iso_sel = true;
0831     //             // "hlt" => "iso" AND "event is triggered"
0832     //             if (hlt_sel_this) hlt_sel = true;
0833     //           // "all" => "met" AND "Z/top rejection cuts"
0834     //             if (all_sel_this) all_sel = true;
0835 
0836     // "rec" => pt,eta cuts are satisfied
0837     if (rec_sel_this)
0838       rec_sel = true;
0839     // "eid" => "rec" AND "electron passes ID"
0840     if (eid_sel_this)
0841       iso_sel = true;
0842     // "iso" => "eid" AND "electron is isolated"
0843     if (iso_sel_this)
0844       iso_sel = true;
0845     // "met" => "iso" AND "MET/MT"
0846     // "all" => "met" AND "event is triggered"
0847     if (all_sel_this)
0848       all_sel = true;
0849 
0850     // Do N-1 histograms now (and only once for global event quantities)
0851     if (flags_passed >= (NFLAGS - 1)) {
0852       if (!electron_sel[0] || flags_passed == NFLAGS) {
0853         pt_after_->Fill(pt);
0854       }
0855       if (!electron_sel[1] || flags_passed == NFLAGS) {
0856         eta_after_->Fill(eta);
0857       }
0858       if (!electron_sel[2] || flags_passed == NFLAGS) {
0859         if (isBarrel) {
0860           sieiebarrel_after_->Fill(sieie);
0861         } else if (isEndcap) {
0862           sieieendcap_after_->Fill(sieie);
0863         }
0864       }
0865       if (!electron_sel[3] || flags_passed == NFLAGS) {
0866         if (isBarrel) {
0867           detainbarrel_after_->Fill(detain);
0868         } else if (isEndcap) {
0869           detainendcap_after_->Fill(detain);
0870         }
0871       }
0872       if (!electron_sel[4] || flags_passed == NFLAGS) {
0873         if (isBarrel) {
0874           ecalisobarrel_after_->Fill(ecalisovar);
0875         } else if (isEndcap) {
0876           ecalisoendcap_after_->Fill(ecalisovar);
0877         }
0878       }
0879       if (!electron_sel[5] || flags_passed == NFLAGS) {
0880         if (isBarrel) {
0881           hcalisobarrel_after_->Fill(hcalisovar);
0882         } else if (isEndcap) {
0883           hcalisoendcap_after_->Fill(hcalisovar);
0884         }
0885       }
0886       if (!electron_sel[6] || flags_passed == NFLAGS) {
0887         if (isBarrel) {
0888           trkisobarrel_after_->Fill(trkisovar);
0889         } else if (isEndcap) {
0890           trkisoendcap_after_->Fill(trkisovar);
0891         }
0892       }
0893       //        if (!electron_sel[3] || flags_passed==NFLAGS)
0894       //          {
0895       //            detain_after_->Fill(detain);
0896       //          }
0897       //        if (!electron_sel[4] || flags_passed==NFLAGS)
0898       //          {
0899       //            ecaliso_after_->Fill(trackerHits);
0900       //          }
0901       //        if (!electron_sel[5] || flags_passed==NFLAGS)
0902       //          {
0903       //            tkelectr_after_->Fill(electr.isTrackerElectron());
0904       //          }
0905       //        if (!electron_sel[6] || flags_passed==NFLAGS)
0906       //          {
0907       //            iso_after_->Fill(isovar);
0908       //          }
0909       /*  if (!electron_sel[7] || flags_passed == NFLAGS) {
0910         if (!hlt_hist_done) {
0911           trig_after_->Fill(trigger_fired);
0912         }
0913       }*/
0914       //hlt_hist_done = true;
0915       if (!electron_sel[7] || flags_passed == NFLAGS) {
0916         mt_after_->Fill(massT);
0917       }
0918       if (!electron_sel[8] || flags_passed == NFLAGS) {
0919         if (!met_hist_done) {
0920           met_after_->Fill(met_et);
0921         }
0922       }
0923       met_hist_done = true;
0924       //                   if (!muon_sel[10] || flags_passed==NFLAGS)
0925       //                         acop_after_->Fill(acop);
0926       //                   if (!muon_sel[11] || flags_passed==NFLAGS)
0927       //                         if (!nz1_hist_done)
0928       // nz1_after_->Fill(nmuonsForZ1);
0929       //                         nz1_hist_done = true;
0930       //                   if (!muon_sel[11] || flags_passed==NFLAGS)
0931       //                         if (!nz2_hist_done)
0932       // nz2_after_->Fill(nmuonsForZ2);
0933       //                         nz2_hist_done = true;
0934       if (!electron_sel[9] || flags_passed == NFLAGS) {
0935         if (!njets_hist_done) {
0936           njets_after_->Fill(njets);
0937           if (jet_et > 10)  // don't want low energy "jets"
0938           {
0939             jet_et_after_->Fill(jet_et);
0940             jet_eta_after_->Fill(jet_eta);
0941           }
0942         }
0943       }
0944       njets_hist_done = true;
0945 
0946     }  // end N-1 histos block
0947 
0948   }  // end loop through electrons
0949 
0950   // inv mass = sqrt(m_1^2 + m_2^2 + 2*(E_1*E_2 - (px1*px2 + py1*py2 + pz1+pz2)
0951   // ) )
0952   double invMass = 0;
0953 
0954   nelectrons_before_->Fill(electronCollectionSize);
0955   if (electronCollectionSize > 1) {
0956     invMass =
0957         sqrt(electron[0][1] + electron[1][1] +
0958              2 * (electron[0][2] * electron[1][2] - (electron[0][3] * electron[1][3] + electron[0][4] * electron[1][4] +
0959                                                      electron[0][5] * electron[1][5])));
0960     invmass_before_->Fill(invMass);
0961     invmassPU_before_->Fill(invMass, npvCount);
0962   }
0963 
0964   nelectrons_after_->Fill(nGoodElectrons);
0965   if (nGoodElectrons > 1) {
0966     invMass = sqrt(goodElectron[0][1] + goodElectron[1][1] +
0967                    2 * (goodElectron[0][2] * goodElectron[1][2] -
0968                         (goodElectron[0][3] * goodElectron[1][3] + goodElectron[0][4] * goodElectron[1][4] +
0969                          goodElectron[0][5] * goodElectron[1][5])));
0970     invmass_after_->Fill(invMass);
0971     invmassPU_afterZ_->Fill(invMass, npvCount);
0972     npvs_afterZ_->Fill(npvCount);
0973   }
0974 
0975   // Collect final flags
0976   if (rec_sel)
0977     nrec++;
0978   if (eid_sel)
0979     neid++;
0980   if (iso_sel)
0981     niso++;
0982   //      if (hlt_sel) nhlt++;
0983   //      if (met_sel) nmet++;
0984 
0985   if (all_sel) {
0986     nsel++;
0987     LogTrace("") << ">>>> Event ACCEPTED";
0988   } else {
0989     LogTrace("") << ">>>> Event REJECTED";
0990   }
0991 
0992   return;
0993 }
0994 
0995 // This always returns only a positive deltaPhi
0996 double EwkElecDQM::calcDeltaPhi(double phi1, double phi2) {
0997   double deltaPhi = phi1 - phi2;
0998 
0999   if (deltaPhi < 0)
1000     deltaPhi = -deltaPhi;
1001 
1002   if (deltaPhi > 3.1415926) {
1003     deltaPhi = 2 * 3.1415926 - deltaPhi;
1004   }
1005 
1006   return deltaPhi;
1007 }
1008 
1009 // Local Variables:
1010 // show-trailing-whitespace: t
1011 // truncate-lines: t
1012 // End: