Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:42:48 

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   float jet2_et = -9.0;
0539   unsigned int jetCollectionSize = jetCollection->size();
0540   int njets = 0;
0541   for (unsigned int i = 0; i < jetCollectionSize; i++) {
0542     const Jet& jet = jetCollection->at(i);
0543 
0544     float jet_current_et = jet.et();
0545     //  cout << "jet_current_et " << jet_current_et << endl;
0546     // if it overlaps with electron, it is not a jet
0547     if (electron_et > 0.0 && fabs(jet.eta() - electron_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron_phi) < 0.2)
0548       continue;
0549     if (electron2_et > 0.0 && fabs(jet.eta() - electron2_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron2_phi) < 0.2)
0550       continue;
0551 
0552     // if it has too low Et, throw away
0553     //  if (jet_current_et < eJetMin_) continue; //Keep if only want to plot
0554     // above jet cut
0555 
0556     if (jet.et() > eJetMin_) {
0557       njets++;
0558     }
0559     if (jet_current_et > jet_et) {
0560       jet2_et = jet_et;   // 2nd highest jet get's et from current highest
0561       jet_et = jet.et();  // current highest jet gets et from the new highest
0562       jet_eta = jet.eta();
0563     } else if (jet_current_et > jet2_et) {
0564       jet2_et = jet.et();
0565     }
0566   }
0567 
0568   // Fill After all electron cuts (or both before and after)
0569   if (jet_et > 10)  // don't want low energy "jets"
0570   {
0571     jet_et_before_->Fill(jet_et);
0572     //    jet2_et_before_  ->Fill(jet2_et);
0573     jet_eta_before_->Fill(jet_eta);
0574   }
0575 
0576   LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
0577   LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
0578   njets_before_->Fill(njets);
0579 
0580   // Start counting
0581   nall++;
0582 
0583   // Histograms per event should be done only once, so keep track of them
0584   //bool hlt_hist_done = false;
0585   // bool minv_hist_done = false;
0586   bool met_hist_done = false;
0587   //       bool nz1_hist_done = false;
0588   //       bool nz2_hist_done = false;
0589   bool njets_hist_done = false;
0590 
0591   // Central selection criteria
0592   // const int NFLAGS = 13; // number of individual selection criteria
0593   const int NFLAGS = 10;  // number of individual selection criteria
0594   // 0: pt cut           | rec
0595   // 1: eta cut          | rec
0596   // 2: sieie            | eid
0597   // 3: detain           | eid
0598   // 4: ecal iso         | iso
0599   // 5: hcal iso         | iso
0600   // 6: trk iso          | iso
0601   // 7: trigger fired    | hlt/all
0602   bool electron_sel[NFLAGS];
0603 
0604   // for invariant mass calculation
0605   // keep track of highest-pt electrons for initial (RECO) electrons
0606   // and "good" electrons (passing all cuts)
0607   // first array dimension is for first or second good electron
0608   // second array dimension is for relevant quantities of good electron
0609   //    [0]: 1 for electron found or 0 for not found (if 0, no other quantities
0610   // filled)
0611   //    [1]: mSqr
0612   //    [2]: E
0613   //    [3]: px
0614   //    [4]: py
0615   //    [5]: pz
0616   // inv mass = sqrt(m_1^2 + m_2^2 + 2*(E_1*E_2 - (px1*px2 + py1*py2 + pz1+pz2)
0617   // ) )
0618   double electron[2][6];
0619   double goodElectron[2][6];
0620   nGoodElectrons = 0;
0621   for (unsigned int i = 0; i < 2; i++) {
0622     for (unsigned int j = 0; j < 6; j++) {
0623       electron[i][j] = 0.;
0624       goodElectron[i][j] = 0.;
0625     }
0626   }
0627 
0628   for (unsigned int i = 0; i < electronCollectionSize; i++) {
0629     for (int j = 0; j < NFLAGS; ++j) {
0630       electron_sel[j] = false;
0631     }
0632 
0633     const GsfElectron& elec = electronCollection->at(i);
0634     // if (!mu.isGlobalMuon()) continue;
0635     // if (mu.globalTrack().isNull()) continue;
0636     // if (mu.innerTrack().isNull()) continue;
0637 
0638     LogTrace("") << "> elec: processing electron number " << i << "...";
0639     // reco::TrackRef gm = mu.globalTrack();
0640     // reco::TrackRef tk = mu.innerTrack();
0641     // should have stuff for electron track?
0642 
0643     if (i < 2) {
0644       electron[i][0] = 1.;
0645       electron[i][1] = elec.massSqr();
0646       electron[i][2] = elec.energy();
0647       electron[i][3] = elec.px();
0648       electron[i][4] = elec.py();
0649       electron[i][5] = elec.pz();
0650     }
0651 
0652     // Pt,eta cuts
0653     double pt = elec.pt();
0654     double px = elec.px();
0655     double py = elec.py();
0656     double eta = elec.eta();
0657     LogTrace("") << "\t... pt, eta: " << pt << " [GeV], " << eta;
0658     ;
0659     if (pt > ptCut_)
0660       electron_sel[0] = true;
0661     if (fabs(eta) < etaCut_)
0662       electron_sel[1] = true;
0663 
0664     bool isBarrel = false;
0665     bool isEndcap = false;
0666     if (eta < 1.4442 && eta > -1.4442) {
0667       isBarrel = true;
0668     } else if ((eta > 1.56 && eta < 2.4) || (eta < -1.56 && eta > -2.4)) {
0669       isEndcap = true;
0670     }
0671 
0672     //             // d0, chi2, nhits quality cuts
0673     //             double dxy = tk->dxy(beamSpotHandle->position());
0674     //             double normalizedChi2 = gm->normalizedChi2();
0675     //             double trackerHits = tk->numberOfValidHits();
0676     //             LogTrace("") << "\t... dxy, normalizedChi2, trackerHits,
0677     // isTrackerMuon?: " << dxy << " [cm], " << normalizedChi2 << ", " <<
0678     // trackerHits << ", " << mu.isTrackerMuon();
0679     //             if (fabs(dxy)<dxyCut_) muon_sel[2] = true;
0680     //             if (normalizedChi2<normalizedChi2Cut_) muon_sel[3] = true;
0681     //             if (trackerHits>=trackerHitsCut_) muon_sel[4] = true;
0682     //             if (mu.isTrackerMuon()) muon_sel[5] = true;
0683 
0684     pt_before_->Fill(pt);
0685     eta_before_->Fill(eta);
0686     //             dxy_before_->Fill(dxy);
0687     //             chi2_before_->Fill(normalizedChi2);
0688     //             nhits_before_->Fill(trackerHits);
0689     //             tkmu_before_->Fill(mu.isTrackerMuon());
0690 
0691     // Electron ID cuts
0692     double sieie = (double)elec.sigmaIetaIeta();
0693     double detain = (double)elec.deltaEtaSuperClusterTrackAtVtx();  // think this is detain
0694     if (sieie < sieieCutBarrel_ && isBarrel)
0695       electron_sel[2] = true;
0696     if (sieie < sieieCutEndcap_ && isEndcap)
0697       electron_sel[2] = true;
0698     if (detain < detainCutBarrel_ && isBarrel)
0699       electron_sel[3] = true;
0700     if (detain < detainCutEndcap_ && isEndcap)
0701       electron_sel[3] = true;
0702     if (isBarrel) {
0703       LogTrace("") << "\t... sieie value " << sieie << " (barrel), pass? " << electron_sel[2];
0704       LogTrace("") << "\t... detain value " << detain << " (barrel), pass? " << electron_sel[3];
0705     } else if (isEndcap) {
0706       LogTrace("") << "\t... sieie value " << sieie << " (endcap), pass? " << electron_sel[2];
0707       LogTrace("") << "\t... detain value " << detain << " (endcap), pass? " << electron_sel[2];
0708     }
0709 
0710     if (isBarrel) {
0711       sieiebarrel_before_->Fill(sieie);
0712       detainbarrel_before_->Fill(detain);
0713     } else if (isEndcap) {
0714       sieieendcap_before_->Fill(sieie);
0715       detainendcap_before_->Fill(detain);
0716     }
0717 
0718     // Isolation cuts
0719     // double isovar = mu.isolationR03().sumPt;
0720     double ecalisovar = elec.dr03EcalRecHitSumEt();  // picked one set!
0721     double hcalisovar = elec.dr03HcalTowerSumEt();   // try others if
0722     double trkisovar = elec.dr04TkSumPt();           // doesn't work
0723     // if (isCombinedIso_) {
0724     // isovar += mu.isolationR03().emEt;
0725     // isovar += mu.isolationR03().hadEt;
0726     //}
0727     // if (isRelativeIso_) isovar /= pt;
0728     if (ecalisovar < ecalIsoCutBarrel_ && isBarrel)
0729       electron_sel[4] = true;
0730     if (ecalisovar < ecalIsoCutEndcap_ && isEndcap)
0731       electron_sel[4] = true;
0732     if (hcalisovar < hcalIsoCutBarrel_ && isBarrel)
0733       electron_sel[5] = true;
0734     if (hcalisovar < hcalIsoCutEndcap_ && isEndcap)
0735       electron_sel[5] = true;
0736     if (trkisovar < trkIsoCutBarrel_ && isBarrel)
0737       electron_sel[6] = true;
0738     if (trkisovar < trkIsoCutEndcap_ && isEndcap)
0739       electron_sel[6] = true;
0740     if (isBarrel) {
0741       LogTrace("") << "\t... ecal isolation value " << ecalisovar << " (barrel), pass? " << electron_sel[4];
0742       LogTrace("") << "\t... hcal isolation value " << hcalisovar << " (barrel), pass? " << electron_sel[5];
0743       LogTrace("") << "\t... trk isolation value " << trkisovar << " (barrel), pass? " << electron_sel[6];
0744     } else if (isEndcap) {
0745       LogTrace("") << "\t... ecal isolation value " << ecalisovar << " (endcap), pass? " << electron_sel[4];
0746       LogTrace("") << "\t... hcal isolation value " << hcalisovar << " (endcap), pass? " << electron_sel[5];
0747       LogTrace("") << "\t... trk isolation value " << trkisovar << " (endcap), pass? " << electron_sel[6];
0748     }
0749 
0750     // iso_before_->Fill(isovar);
0751     if (isBarrel) {
0752       ecalisobarrel_before_->Fill(ecalisovar);
0753       hcalisobarrel_before_->Fill(hcalisovar);
0754       trkisobarrel_before_->Fill(trkisovar);
0755     } else if (isEndcap) {
0756       ecalisoendcap_before_->Fill(ecalisovar);
0757       hcalisoendcap_before_->Fill(hcalisovar);
0758       trkisoendcap_before_->Fill(trkisovar);
0759     }
0760 
0761     // HLT
0762     // if (trigger_fired) electron_sel[7] = true;
0763 
0764     //             // MET/MT cuts
0765     double w_et = met_et + pt;
0766     double w_px = met_px + px;
0767     double w_py = met_py + py;
0768 
0769     double massT = w_et * w_et - w_px * w_px - w_py * w_py;
0770     massT = (massT > 0) ? sqrt(massT) : 0;
0771 
0772     LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << w_px << ", " << w_py
0773                  << " [GeV]";
0774     if (massT > mtMin_ && massT < mtMax_)
0775       electron_sel[7] = true;
0776     mt_before_->Fill(massT);
0777     if (met_et > metMin_ && met_et < metMax_)
0778       electron_sel[8] = true;
0779 
0780     //             // Acoplanarity cuts
0781     //             Geom::Phi<double> deltaphi(mu.phi()-atan2(met_py,met_px));
0782     //             double acop = deltaphi.value();
0783     //             if (acop<0) acop = - acop;
0784     //             acop = M_PI - acop;
0785     //             LogTrace("") << "\t... acoplanarity: " << acop;
0786     //             if (acop<acopCut_) muon_sel[10] = true;
0787     //             acop_before_->Fill(acop);
0788 
0789     //             // Remaining flags (from global event information)
0790     //             if (nmuonsForZ1<1 || nmuonsForZ2<2) muon_sel[11] = true;
0791     if (njets <= nJetMax_)
0792       electron_sel[9] = true;
0793 
0794     // Collect necessary flags "per electron"
0795     int flags_passed = 0;
0796     bool rec_sel_this = true;
0797     bool eid_sel_this = true;
0798     bool iso_sel_this = true;
0799     bool all_sel_this = true;
0800     for (int j = 0; j < NFLAGS; ++j) {
0801       if (electron_sel[j])
0802         flags_passed += 1;
0803       if (j < 2 && !electron_sel[j])
0804         rec_sel_this = false;
0805       if (j < 4 && !electron_sel[j])
0806         eid_sel_this = false;
0807       if (j < 7 && !electron_sel[j])
0808         iso_sel_this = false;
0809       if (!electron_sel[j])
0810         all_sel_this = false;
0811     }
0812 
0813     if (all_sel_this) {
0814       if (nGoodElectrons < 2) {
0815         goodElectron[nGoodElectrons][0] = 1.;
0816         goodElectron[nGoodElectrons][1] = elec.massSqr();
0817         goodElectron[nGoodElectrons][2] = elec.energy();
0818         goodElectron[nGoodElectrons][3] = elec.px();
0819         goodElectron[nGoodElectrons][4] = elec.py();
0820         goodElectron[nGoodElectrons][5] = elec.pz();
0821       }
0822       nGoodElectrons++;
0823     }
0824 
0825     //             // "rec" => pt,eta and quality cuts are satisfied
0826     //             if (rec_sel_this) rec_sel = true;
0827     //             // "iso" => "rec" AND "muon is isolated"
0828     //             if (iso_sel_this) iso_sel = true;
0829     //             // "hlt" => "iso" AND "event is triggered"
0830     //             if (hlt_sel_this) hlt_sel = true;
0831     //           // "all" => "met" AND "Z/top rejection cuts"
0832     //             if (all_sel_this) all_sel = true;
0833 
0834     // "rec" => pt,eta cuts are satisfied
0835     if (rec_sel_this)
0836       rec_sel = true;
0837     // "eid" => "rec" AND "electron passes ID"
0838     if (eid_sel_this)
0839       iso_sel = true;
0840     // "iso" => "eid" AND "electron is isolated"
0841     if (iso_sel_this)
0842       iso_sel = true;
0843     // "met" => "iso" AND "MET/MT"
0844     // "all" => "met" AND "event is triggered"
0845     if (all_sel_this)
0846       all_sel = true;
0847 
0848     // Do N-1 histograms now (and only once for global event quantities)
0849     if (flags_passed >= (NFLAGS - 1)) {
0850       if (!electron_sel[0] || flags_passed == NFLAGS) {
0851         pt_after_->Fill(pt);
0852       }
0853       if (!electron_sel[1] || flags_passed == NFLAGS) {
0854         eta_after_->Fill(eta);
0855       }
0856       if (!electron_sel[2] || flags_passed == NFLAGS) {
0857         if (isBarrel) {
0858           sieiebarrel_after_->Fill(sieie);
0859         } else if (isEndcap) {
0860           sieieendcap_after_->Fill(sieie);
0861         }
0862       }
0863       if (!electron_sel[3] || flags_passed == NFLAGS) {
0864         if (isBarrel) {
0865           detainbarrel_after_->Fill(detain);
0866         } else if (isEndcap) {
0867           detainendcap_after_->Fill(detain);
0868         }
0869       }
0870       if (!electron_sel[4] || flags_passed == NFLAGS) {
0871         if (isBarrel) {
0872           ecalisobarrel_after_->Fill(ecalisovar);
0873         } else if (isEndcap) {
0874           ecalisoendcap_after_->Fill(ecalisovar);
0875         }
0876       }
0877       if (!electron_sel[5] || flags_passed == NFLAGS) {
0878         if (isBarrel) {
0879           hcalisobarrel_after_->Fill(hcalisovar);
0880         } else if (isEndcap) {
0881           hcalisoendcap_after_->Fill(hcalisovar);
0882         }
0883       }
0884       if (!electron_sel[6] || flags_passed == NFLAGS) {
0885         if (isBarrel) {
0886           trkisobarrel_after_->Fill(trkisovar);
0887         } else if (isEndcap) {
0888           trkisoendcap_after_->Fill(trkisovar);
0889         }
0890       }
0891       //        if (!electron_sel[3] || flags_passed==NFLAGS)
0892       //          {
0893       //            detain_after_->Fill(detain);
0894       //          }
0895       //        if (!electron_sel[4] || flags_passed==NFLAGS)
0896       //          {
0897       //            ecaliso_after_->Fill(trackerHits);
0898       //          }
0899       //        if (!electron_sel[5] || flags_passed==NFLAGS)
0900       //          {
0901       //            tkelectr_after_->Fill(electr.isTrackerElectron());
0902       //          }
0903       //        if (!electron_sel[6] || flags_passed==NFLAGS)
0904       //          {
0905       //            iso_after_->Fill(isovar);
0906       //          }
0907       /*  if (!electron_sel[7] || flags_passed == NFLAGS) {
0908         if (!hlt_hist_done) {
0909           trig_after_->Fill(trigger_fired);
0910         }
0911       }*/
0912       //hlt_hist_done = true;
0913       if (!electron_sel[7] || flags_passed == NFLAGS) {
0914         mt_after_->Fill(massT);
0915       }
0916       if (!electron_sel[8] || flags_passed == NFLAGS) {
0917         if (!met_hist_done) {
0918           met_after_->Fill(met_et);
0919         }
0920       }
0921       met_hist_done = true;
0922       //                   if (!muon_sel[10] || flags_passed==NFLAGS)
0923       //                         acop_after_->Fill(acop);
0924       //                   if (!muon_sel[11] || flags_passed==NFLAGS)
0925       //                         if (!nz1_hist_done)
0926       // nz1_after_->Fill(nmuonsForZ1);
0927       //                         nz1_hist_done = true;
0928       //                   if (!muon_sel[11] || flags_passed==NFLAGS)
0929       //                         if (!nz2_hist_done)
0930       // nz2_after_->Fill(nmuonsForZ2);
0931       //                         nz2_hist_done = true;
0932       if (!electron_sel[9] || flags_passed == NFLAGS) {
0933         if (!njets_hist_done) {
0934           njets_after_->Fill(njets);
0935           if (jet_et > 10)  // don't want low energy "jets"
0936           {
0937             jet_et_after_->Fill(jet_et);
0938             jet_eta_after_->Fill(jet_eta);
0939           }
0940         }
0941       }
0942       njets_hist_done = true;
0943 
0944     }  // end N-1 histos block
0945 
0946   }  // end loop through electrons
0947 
0948   // inv mass = sqrt(m_1^2 + m_2^2 + 2*(E_1*E_2 - (px1*px2 + py1*py2 + pz1+pz2)
0949   // ) )
0950   double invMass = 0;
0951 
0952   nelectrons_before_->Fill(electronCollectionSize);
0953   if (electronCollectionSize > 1) {
0954     invMass =
0955         sqrt(electron[0][1] + electron[1][1] +
0956              2 * (electron[0][2] * electron[1][2] - (electron[0][3] * electron[1][3] + electron[0][4] * electron[1][4] +
0957                                                      electron[0][5] * electron[1][5])));
0958     invmass_before_->Fill(invMass);
0959     invmassPU_before_->Fill(invMass, npvCount);
0960   }
0961 
0962   nelectrons_after_->Fill(nGoodElectrons);
0963   if (nGoodElectrons > 1) {
0964     invMass = sqrt(goodElectron[0][1] + goodElectron[1][1] +
0965                    2 * (goodElectron[0][2] * goodElectron[1][2] -
0966                         (goodElectron[0][3] * goodElectron[1][3] + goodElectron[0][4] * goodElectron[1][4] +
0967                          goodElectron[0][5] * goodElectron[1][5])));
0968     invmass_after_->Fill(invMass);
0969     invmassPU_afterZ_->Fill(invMass, npvCount);
0970     npvs_afterZ_->Fill(npvCount);
0971   }
0972 
0973   // Collect final flags
0974   if (rec_sel)
0975     nrec++;
0976   if (eid_sel)
0977     neid++;
0978   if (iso_sel)
0979     niso++;
0980   //      if (hlt_sel) nhlt++;
0981   //      if (met_sel) nmet++;
0982 
0983   if (all_sel) {
0984     nsel++;
0985     LogTrace("") << ">>>> Event ACCEPTED";
0986   } else {
0987     LogTrace("") << ">>>> Event REJECTED";
0988   }
0989 
0990   return;
0991 }
0992 
0993 // This always returns only a positive deltaPhi
0994 double EwkElecDQM::calcDeltaPhi(double phi1, double phi2) {
0995   double deltaPhi = phi1 - phi2;
0996 
0997   if (deltaPhi < 0)
0998     deltaPhi = -deltaPhi;
0999 
1000   if (deltaPhi > 3.1415926) {
1001     deltaPhi = 2 * 3.1415926 - deltaPhi;
1002   }
1003 
1004   return deltaPhi;
1005 }
1006 
1007 // Local Variables:
1008 // show-trailing-whitespace: t
1009 // truncate-lines: t
1010 // End:
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!