text: #include "DQM/Physics/src/EwkMuDQM.h"
  
  #include "FWCore/MessageLogger/interface/MessageLogger.h"
  #include "FWCore/ParameterSet/interface/ParameterSet.h"
  #include "DataFormats/Common/interface/Handle.h"
  
  #include "FWCore/ServiceRegistry/interface/Service.h"
  #include "CommonTools/UtilAlgos/interface/TFileService.h"
  
  #include "DQMServices/Core/interface/DQMStore.h"
  #include "DQMServices/Core/interface/MonitorElement.h"
  
  #include "DataFormats/TrackReco/interface/Track.h"
  #include "DataFormats/BeamSpot/interface/BeamSpot.h"
  #include "DataFormats/VertexReco/interface/Vertex.h"
  
  #include "DataFormats/MuonReco/interface/Muon.h"
  #include "DataFormats/MuonReco/interface/MuonSelectors.h"
  #include "DataFormats/METReco/interface/MET.h"
  #include "DataFormats/JetReco/interface/Jet.h"
  #include "DataFormats/EgammaCandidates/interface/Photon.h"
  
  #include "DataFormats/GeometryVector/interface/Phi.h"
  
  #include "FWCore/Common/interface/TriggerNames.h"
  #include "FWCore/Framework/interface/Event.h"
  #include "DataFormats/Common/interface/TriggerResults.h"
  
  #include "DataFormats/Common/interface/View.h"
  
  #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
  
  using namespace edm;
  using namespace std;
  using namespace reco;
  
  EwkMuDQM::EwkMuDQM(const ParameterSet& cfg)
      :  // Input collections
        metTag_(cfg.getUntrackedParameter<edm::InputTag>("METTag", edm::InputTag("pfmet"))),
        jetTag_(cfg.getUntrackedParameter<edm::InputTag>("JetTag", edm::InputTag("ak4PFJets"))),
        //trigTag_(consumes<edm::TriggerResults>(
        // cfg.getUntrackedParameter<edm::InputTag>(
        // "TrigTag", edm::InputTag("TriggerResults::HLT")))),
        muonTag_(consumes<edm::View<reco::Muon> >(
            cfg.getUntrackedParameter<edm::InputTag>("MuonTag", edm::InputTag("muons")))),
        metToken_(
            consumes<edm::View<reco::MET> >(cfg.getUntrackedParameter<edm::InputTag>("METTag", edm::InputTag("pfmet")))),
        jetToken_(consumes<edm::View<reco::Jet> >(
            cfg.getUntrackedParameter<edm::InputTag>("JetTag", edm::InputTag("ak4PFJets")))),
        phoTag_(consumes<edm::View<reco::Photon> >(
            cfg.getUntrackedParameter<edm::InputTag>("phoTag", edm::InputTag("photons")))),
        vertexTag_(consumes<edm::View<reco::Vertex> >(
            cfg.getUntrackedParameter<edm::InputTag>("VertexTag", edm::InputTag("offlinePrimaryVertices")))),
        beamSpotTag_(consumes<reco::BeamSpot>(
            cfg.getUntrackedParameter<edm::InputTag>("beamSpotTag", edm::InputTag("offlineBeamSpot")))),
        // trigPathNames_(cfg.getUntrackedParameter<std::vector<std::string> >(
        //   "TrigPathNames")),
  
        // Muon quality cuts
        isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon", true)),  // Glb muon also tracker muon
        dxyCut_(cfg.getUntrackedParameter<double>("DxyCut", 0.2)),                       // dxy < 0.2 cm
        normalizedChi2Cut_(
            cfg.getUntrackedParameter<double>("NormalizedChi2Cut", 10.)),  // chi2/ndof (of global fit) <10.0
        trackerHitsCut_(cfg.getUntrackedParameter<int>("TrackerHitsCut",
                                                       11)),               // Tracker Hits >10
        pixelHitsCut_(cfg.getUntrackedParameter<int>("PixelHitsCut", 1)),  // Pixel Hits >0
        muonHitsCut_(cfg.getUntrackedParameter<int>("MuonHitsCut",
                                                    1)),                 // Valid Muon Hits >0
        nMatchesCut_(cfg.getUntrackedParameter<int>("NMatchesCut", 2)),  // At least 2 Chambers with matches
  
        // W-boson cuts
        isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso", true)),
        isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso", false)),
        isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
        acopCut_(cfg.getUntrackedParameter<double>("AcopCut", 999.)),
        metMin_(cfg.getUntrackedParameter<double>("MetMin", -999999.)),
        metMax_(cfg.getUntrackedParameter<double>("MetMax", 999999.)),
        mtMin_(cfg.getUntrackedParameter<double>("MtMin", 50.)),
        mtMax_(cfg.getUntrackedParameter<double>("MtMax", 200.)),
        ptCut_(cfg.getUntrackedParameter<double>("PtCut", 20.)),
        etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.4)),
  
        // Z rejection
        ptThrForZ1_(cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
        ptThrForZ2_(cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
  
        // Z selection
        dimuonMassMin_(cfg.getUntrackedParameter<double>("dimuonMassMin", 80.)),
        dimuonMassMax_(cfg.getUntrackedParameter<double>("dimuonMassMax", 120.)),
  
        // Top rejection
        eJetMin_(cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
        nJetMax_(cfg.getUntrackedParameter<int>("NJetMax", 999999)),
  
        // Photon cuts
        ptThrForPhoton_(cfg.getUntrackedParameter<double>("ptThrForPhoton", 5.)),
        nPhoMax_(cfg.getUntrackedParameter<int>("nPhoMax", 999999)),
        hltPrescaleProvider_(cfg, consumesCollector(), *this) {
    isValidHltConfig_ = false;
  }
  
  void EwkMuDQM::dqmBeginRun(const Run& iRun, const EventSetup& iSet) {
    nall = 0;
    nsel = 0;
    nz = 0;
  
    nrec = 0;
    niso = 0;
    nhlt = 0;
    nmet = 0;
  
    // passed as parameter to HLTConfigProvider::init(), not yet used
    bool isConfigChanged = false;
    // isValidHltConfig_ used to short-circuit analyze() in case of problems
    isValidHltConfig_ = hltPrescaleProvider_.init(iRun, iSet, "HLT", isConfigChanged);
  }
  
  void EwkMuDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
    ibooker.setCurrentFolder("Physics/EwkMuDQM");
  
    char chtitle[256] = "";
  
    pt_before_ = ibooker.book1D("PT_BEFORECUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
    pt_after_ = ibooker.book1D("PT_AFTERWCUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
  
    eta_before_ = ibooker.book1D("ETA_BEFORECUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
    eta_after_ = ibooker.book1D("ETA_AFTERWCUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
  
    dxy_before_ = ibooker.book1D("DXY_BEFORECUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
    dxy_after_ = ibooker.book1D("DXY_AFTERWCUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
  
    goodewkmuon_before_ = ibooker.book1D("GOODEWKMUON_BEFORECUTS", "Quality-muon flag", 2, -0.5, 1.5);
    goodewkmuon_after_ = ibooker.book1D("GOODEWKMUON_AFTERWCUTS", "Quality-muon flag", 2, -0.5, 1.5);
  
    if (isRelativeIso_) {
      if (isCombinedIso_) {
        iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Relative (combined) isolation variable", 100, 0., 1.);
        iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Relative (combined) isolation variable", 100, 0., 1.);
      } else {
        iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
        iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
      }
    } else {
      if (isCombinedIso_) {
        iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
        iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
      } else {
        iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
        iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
      }
    }
  
    /*  trig_before_ = ibooker.book1D("TRIG_BEFORECUTS",
        "Trigger response (boolean of muon triggers)", 2, -0.5, 1.5);
    trig_after_ = ibooker.book1D("TRIG_AFTERWCUTS",
        "Trigger response (boolean of muon triggers)", 2, -0.5, 1.5);
  */
    snprintf(chtitle, 255, "Transverse mass (%s) [GeV]", metTag_.label().data());
    mt_before_ = ibooker.book1D("MT_BEFORECUTS", chtitle, 150, 0., 300.);
    mt_after_ = ibooker.book1D("MT_AFTERWCUTS", chtitle, 150, 0., 300.);
  
    snprintf(chtitle, 255, "Missing transverse energy (%s) [GeV]", metTag_.label().data());
    met_before_ = ibooker.book1D("MET_BEFORECUTS", chtitle, 100, 0., 200.);
    met_after_ = ibooker.book1D("MET_AFTERWCUTS", chtitle, 100, 0., 200.);
    met_afterZ_ = ibooker.book1D("MET_AFTERZCUTS", chtitle, 100, 0., 200.);
  
    snprintf(chtitle, 255, "MU-MET (%s) acoplanarity", metTag_.label().data());
    acop_before_ = ibooker.book1D("ACOP_BEFORECUTS", chtitle, 50, 0., M_PI);
    acop_after_ = ibooker.book1D("ACOP_AFTERWCUTS", chtitle, 50, 0., M_PI);
  
    snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ1_);
    n_zselPt1thr_ = ibooker.book1D("NZSELPT1THR", chtitle, 10, -0.5, 9.5);
    snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ2_);
    n_zselPt2thr_ = ibooker.book1D("NZSELPT2THR", chtitle, 10, -0.5, 9.5);
  
    snprintf(chtitle, 255, "Number of jets (%s) above %.2f GeV", jetTag_.label().data(), eJetMin_);
    njets_before_ = ibooker.book1D("NJETS_BEFORECUTS", chtitle, 16, -0.5, 15.5);
    njets_after_ = ibooker.book1D("NJETS_AFTERWCUTS", chtitle, 16, -0.5, 15.5);
    njets_afterZ_ = ibooker.book1D("NJETS_AFTERZCUTS", chtitle, 16, -0.5, 15.5);
  
    leadingjet_pt_before_ = ibooker.book1D("LEADINGJET_PT_BEFORECUTS", "Leading Jet transverse momentum", 300, 0., 300.);
    leadingjet_pt_after_ = ibooker.book1D("LEADINGJET_PT_AFTERWCUTS", "Leading Jet transverse momentum", 300, 0., 300.);
    leadingjet_pt_afterZ_ = ibooker.book1D("LEADINGJET_PT_AFTERZCUTS", "Leading Jet transverse momentum", 300, 0., 300.);
  
    leadingjet_eta_before_ = ibooker.book1D("LEADINGJET_ETA_BEFORECUTS", "Leading Jet pseudo-rapidity", 50, -2.5, 2.5);
    leadingjet_eta_after_ = ibooker.book1D("LEADINGJET_ETA_AFTERWCUTS", "Leading Jet pseudo-rapidity", 50, -2.5, 2.5);
    leadingjet_eta_afterZ_ = ibooker.book1D("LEADINGJET_ETA_AFTERZCUTS", "Leading Jet pseudo-rapidity", 50, -2.5, 2.5);
  
    ptDiffPM_before_ = ibooker.book1D("PTDIFFPM_BEFORE_CUTS", "pt(Muon+)-pt(Muon-) after Z cuts [GeV]", 200, -100., 100.);
    ptDiffPM_afterZ_ = ibooker.book1D("PTDIFFPM_AFTERZ_CUTS", "pt(Muon+)-pt(Muon-) after Z cuts [GeV]", 200, -100., 100.);
  
    /**\ For Z-boson events  */
  
    pt1_afterZ_ = ibooker.book1D("PT1_AFTERZCUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
    eta1_afterZ_ = ibooker.book1D("ETA1_AFTERZCUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
    dxy1_afterZ_ = ibooker.book1D("DXY1_AFTERZCUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
    goodewkmuon1_afterZ_ = ibooker.book1D("GOODEWKMUON1_AFTERZCUTS", "Quality-muon flag", 2, -0.5, 1.5);
  
    if (isRelativeIso_) {
      if (isCombinedIso_) {
        iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Relative (combined) isolation variable", 100, 0., 1.);
        iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Relative (combined) isolation variable", 100, 0., 1.);
      } else {
        iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
        iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
      }
    } else {
      if (isCombinedIso_) {
        iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
        iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
      } else {
        iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
        iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
      }
    }
  
    pt2_afterZ_ = ibooker.book1D("PT2_AFTERZCUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
    eta2_afterZ_ = ibooker.book1D("ETA2_AFTERZCUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
    dxy2_afterZ_ = ibooker.book1D("DXY2_AFTERZCUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
    goodewkmuon2_afterZ_ = ibooker.book1D("GOODEWKMUON2_AFTERZCUTS", "Quality-muon flag", 2, -0.5, 1.5);
    /*  ztrig_afterZ_ = ibooker.book1D("ZTRIG_AFTERZCUTS",
        "Trigger response (boolean of muon triggers)", 2, -0.5, 1.5);
   */
    dimuonmass_before_ = ibooker.book1D("DIMUONMASS_BEFORECUTS", "DiMuonMass (2 globals)", 100, 0, 200);
    dimuonmass_afterZ_ = ibooker.book1D("DIMUONMASS_AFTERZCUTS", "DiMuonMass (2 globals)", 100, 0, 200);
    npvs_before_ = ibooker.book1D("NPVs_BEFORECUTS", "Number of Valid Primary Vertices", 51, -0.5, 50.5);
    npvs_after_ = ibooker.book1D("NPVs_AFTERWCUTS", "Number of Valid Primary Vertices", 51, -0.5, 50.5);
    npvs_afterZ_ = ibooker.book1D("NPVs_AFTERZCUTS", "Number of Valid Primary Vertices", 51, -0.5, 50.5);
    muoncharge_before_ = ibooker.book1D("MUONCHARGE_BEFORECUTS", "Muon Charge", 3, -1.5, 1.5);
    muoncharge_after_ = ibooker.book1D("MUONCHARGE_AFTERWCUTS", "Muon Charge", 3, -1.5, 1.5);
    muoncharge_afterZ_ = ibooker.book1D("MUONCHARGE_AFTERZCUTS", "Muon Charge", 3, -1.5, 1.5);
  
    // Adding these to replace the NZ ones (more useful, since they are more
    // general?)
    nmuons_ = ibooker.book1D("NMuons", "Number of muons in the event", 10, -0.5, 9.5);
    ngoodmuons_ = ibooker.book1D("NGoodMuons", "Number of muons passing the quality criteria", 10, -0.5, 9.5);
  
    nph_ = ibooker.book1D("nph", "Number of photons in the event", 20, 0., 20.);
    phPt_ = ibooker.book1D("phPt", "Photon transverse momentum [GeV]", 100, 0., 1000.);
    snprintf(chtitle, 255, "Photon pseudorapidity (pT>%4.1f)", ptThrForPhoton_);
    phEta_ = ibooker.book1D("phEta", chtitle, 100, -2.5, 2.5);
  }
  
  void EwkMuDQM::endRun(const Run& r, const EventSetup& iSet) {}
  
  void EwkMuDQM::analyze(const Event& ev, const EventSetup& iSet) {
    // Muon collection
    Handle<View<Muon> > muonCollection;
    if (!ev.getByToken(muonTag_, muonCollection)) {
      // LogWarning("") << ">>> Muon collection does not exist !!!";
      return;
    }
    unsigned int muonCollectionSize = muonCollection->size();
  
    // Beam spot
    Handle<reco::BeamSpot> beamSpotHandle;
    if (!ev.getByToken(beamSpotTag_, beamSpotHandle)) {
      // LogWarning("") << ">>> No beam spot found !!!";
      return;
    }
  
    // Loop to reject/control Z->mumu is done separately
    unsigned int nmuonsForZ1 = 0;
    unsigned int nmuonsForZ2 = 0;
    bool cosmic = false;
    for (unsigned int i = 0; i < muonCollectionSize; i++) {
      const Muon& mu = muonCollection->at(i);
      if (!mu.isGlobalMuon())
        continue;
      double pt = mu.pt();
      double dxy = mu.innerTrack()->dxy(beamSpotHandle->position());
  
      if (fabs(dxy) > 1) {
        cosmic = true;
        break;
      }
  
      if (pt > ptThrForZ1_)
        nmuonsForZ1++;
      if (pt > ptThrForZ2_)
        nmuonsForZ2++;
  
      for (unsigned int j = i + 1; j < muonCollectionSize; j++) {
        const Muon& mu2 = muonCollection->at(j);
        if (mu2.isGlobalMuon() && (mu.charge() * mu2.charge() == -1)) {
          const math::XYZTLorentzVector ZRecoGlb(
              mu.px() + mu2.px(), mu.py() + mu2.py(), mu.pz() + mu2.pz(), mu.p() + mu2.p());
          dimuonmass_before_->Fill(ZRecoGlb.mass());
          if (mu.charge() > 0) {
            ptDiffPM_before_->Fill(mu.pt() - mu2.pt());
          } else {
            ptDiffPM_before_->Fill(mu2.pt() - mu.pt());
          }
        }
      }
    }
    if (cosmic)
      return;
  
    LogTrace("") << "> Z rejection: muons above " << ptThrForZ1_ << " [GeV]: " << nmuonsForZ1;
    LogTrace("") << "> Z rejection: muons above " << ptThrForZ2_ << " [GeV]: " << nmuonsForZ2;
  
    // MET
    Handle<View<MET> > metCollection;
    if (!ev.getByToken(metToken_, metCollection)) {
      // LogWarning("") << ">>> MET collection does not exist !!!";
      return;
    }
    const MET& met = metCollection->at(0);
    double met_et = met.pt();
    LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met.px() << ", " << met.py() << " [GeV]";
    met_before_->Fill(met_et);
  
    // Vertices in the event
    Handle<View<reco::Vertex> > vertexCollection;
    if (!ev.getByToken(vertexTag_, vertexCollection)) {
      LogError("") << ">>> Vertex collection does not exist !!!";
      return;
    }
    unsigned int vertexCollectionSize = vertexCollection->size();
  
    int nvvertex = 0;
    for (unsigned int i = 0; i < vertexCollectionSize; i++) {
      const Vertex& vertex = vertexCollection->at(i);
      if (vertex.isValid())
        nvvertex++;
    }
  
    npvs_before_->Fill(nvvertex);
  
    // bool trigger_fired = false;
    //Handle<TriggerResults> triggerResults;
    // if (!ev.getByToken(trigTag_, triggerResults)) {
    // LogWarning("") << ">>> TRIGGER collection does not exist !!!";
    // return;
    // }
    // const edm::TriggerNames& trigNames = ev.triggerNames(*triggerResults);
    //  LogWarning("")<<"Loop over triggers";
  
    //HLTConfigProvider const&  hltConfigProvider = hltPrescaleProvider_.hltConfigProvider();
  
    /*  change faulty logic of triggering
    for (unsigned int i=0; i<triggerResults->size(); i++)
    {
            const std::string trigName = trigNames.triggerName(i);
  
            bool found=false;
            for(unsigned int index=0; index<trigPathNames_.size() && found==false;
    index++) {
                 size_t trigPath = trigName.find(trigPathNames_[index]); // 0 if
    found, pos if not
                 if (trigPath==0) found=true;
            }
            if(!found) {continue;}
  
            bool prescaled=false;
            for (unsigned int ps= 0; ps<  hltConfigProvider.prescaleSize();
    ps++){
                const unsigned int prescaleValue =
    hltConfigProvider.prescaleValue(ps, trigName) ;
                if (prescaleValue != 1) prescaled =true;
            }
  
            if( triggerResults->accept(i) && !prescaled){   trigger_fired=true;}
                      // LogWarning("")<<"TrigNo: "<<i<<"  "<<found<<"
    "<<trigName<<" ---> FIRED";}
    }
    */
  
    // get the prescale set for this event
    const int prescaleSet = hltPrescaleProvider_.prescaleSet(ev, iSet);
    if (prescaleSet == -1) {
      LogTrace("") << "Failed to determine prescaleSet\n";
      // std::cout << "Failed to determine prescaleSet. Check the GlobalTag in
      // cfg\n";
      return;
    }
  
    // for (unsigned int i = 0;
    //    (i < triggerResults->size()) && (trigger_fired == false); i++) {
    // skip trigger, if it did not fire
    //if (!triggerResults->accept(i)) continue;
  
    // skip trigger, if it is not on our list
    //bool found = false;
    //const std::string trigName = trigNames.triggerName(i);
    //for (unsigned int index = 0;
    //   index < trigPathNames_.size() && found == false; index++) {
    // if (trigName.find(trigPathNames_.at(index)) == 0) found = true;
    // }
    // if (!found) continue;
  
    // skip trigger, if it is prescaled
    /* if (prescaleSet != -1) {
        if (hltConfigProvider.prescaleValue(prescaleSet, trigName) != 1)
          continue;
      } else {
        // prescaleSet is not known.
        // This branch is not needed, if prescaleSet=-1 forces to skip event
        int prescaled = 0;
        for (unsigned int ps = 0;
             !prescaled && (ps < hltConfigProvider.prescaleSize()); ++ps) {
          if (hltConfigProvider.prescaleValue(ps, trigName) != 1) {
            prescaled = 1;
          }
        }
        if (prescaled) {
          // std::cout << "trigger prescaled\n";
          continue;
        }
      }
  */
    // std::cout << "found unprescaled trigger that fired: " << trigName <<
    // "\n";
    // trigger_fired = true;
    // }
    // if (trigger_fired) std::cout << "\n\tGot Trigger\n";
  
    // trig_before_->Fill(trigger_fired);
  
    // Jet collection
    Handle<View<Jet> > jetCollection;
    if (!ev.getByToken(jetToken_, jetCollection)) {
      // LogError("") << ">>> JET collection does not exist !!!";
      return;
    }
    unsigned int jetCollectionSize = jetCollection->size();
    int njets = 0;
    int LEADJET = -1;
    double max_pt = 0;
    for (unsigned int i = 0; i < jetCollectionSize; i++) {
      const Jet& jet = jetCollection->at(i);
      double minDistance = 99999;  // This is in order to use PFJets
      for (unsigned int j = 0; j < muonCollectionSize; j++) {
        const Muon& mu = muonCollection->at(j);
        double distance =
            sqrt((mu.eta() - jet.eta()) * (mu.eta() - jet.eta()) + (mu.phi() - jet.phi()) * (mu.phi() - jet.phi()));
        if (minDistance > distance)
          minDistance = distance;
      }
      if (minDistance < 0.3)
        continue;  // 0.3 is the isolation cone around the muon
      if (jet.et() > max_pt) {
        LEADJET = i;
        max_pt = jet.et();
      }
      if (jet.et() > eJetMin_) {
        njets++;
      }
    }
  
    LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
    LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
    njets_before_->Fill(njets);
    double lead_jet_pt = -1;
    double lead_jet_eta = -100;
    if (LEADJET != -1) {
      const Jet& leadJet = jetCollection->at(LEADJET);
      leadingjet_pt_before_->Fill(leadJet.pt());
      leadingjet_eta_before_->Fill(leadJet.eta());
      lead_jet_pt = leadJet.pt();
      lead_jet_eta = leadJet.eta();
    }
    // Photon Collection
    Handle<View<Photon> > photonCollection;
    if (!ev.getByToken(phoTag_, photonCollection)) {
      // LogError("")
      return;
    }
    unsigned int ngam = 0;
  
    for (unsigned int i = 0; i < photonCollection->size(); i++) {
      const Photon& ph = photonCollection->at(i);
      double photonPt = ph.pt();
      if (photonPt > ptThrForPhoton_) {
        ngam++;
        phEta_->Fill(ph.eta());
      }
      phPt_->Fill(photonPt);
    }
    nph_->Fill(ngam);
    LogTrace("") << " >>> N photons " << ngam << std::endl;
  
    nmuons_->Fill(muonCollectionSize);
  
    // Start counting
    nall++;
  
    // Histograms per event should be done only once, so keep track of them
    // bool hlt_hist_done = false;
    // bool zhlt_hist_done = false;
    bool zjets_hist_done = false;
    bool zfullsel_hist_done = false;
    bool met_hist_done = false;
    bool njets_hist_done = false;
    bool wfullsel_hist_done = false;
  
    // Central W->mu nu selection criteria
    const int NFLAGS = 10;
    bool muon_sel[NFLAGS];
    const int NFLAGSZ = 12;
    bool zmuon_sel[NFLAGSZ];
    bool muon4Z = false;
  
    double number_of_muons = 0;
    double number_of_goodMuons = 0;
  
    for (unsigned int i = 0; i < muonCollectionSize; i++) {
      for (int j = 0; j < NFLAGS; ++j) {
        muon_sel[j] = false;
      }
  
      number_of_muons++;
  
      const Muon& mu = muonCollection->at(i);
      if (!mu.isGlobalMuon())
        continue;
      if (mu.globalTrack().isNull())
        continue;
      if (mu.innerTrack().isNull())
        continue;
  
      LogTrace("") << "> Wsel: processing muon number " << i << "...";
      reco::TrackRef gm = mu.globalTrack();
      reco::TrackRef tk = mu.innerTrack();
  
      // Pt,eta cuts
      double pt = mu.pt();
      double eta = mu.eta();
      LogTrace("") << "\t... pt, eta: " << pt << " [GeV], " << eta;
      ;
      if (pt > ptCut_)
        muon_sel[0] = true;
      if (fabs(eta) < etaCut_)
        muon_sel[1] = true;
  
      double charge = mu.charge();
  
      // d0, chi2, nhits quality cuts
      double dxy = gm->dxy(beamSpotHandle->position());
      double normalizedChi2 = gm->normalizedChi2();
      double trackerHits = tk->hitPattern().numberOfValidTrackerHits();
      int pixelHits = tk->hitPattern().numberOfValidPixelHits();
      int muonHits = gm->hitPattern().numberOfValidMuonHits();
      int nMatches = mu.numberOfMatches();
  
      LogTrace("") << "\t... dxy, normalizedChi2, trackerHits, isTrackerMuon?: " << dxy << " [cm], " << normalizedChi2
                   << ", " << trackerHits << ", " << mu.isTrackerMuon();
      if (fabs(dxy) < dxyCut_)
        muon_sel[2] = true;
  
      bool quality = true;
  
      if (normalizedChi2 > normalizedChi2Cut_)
        quality = false;
      if (trackerHits < trackerHitsCut_)
        quality = false;
      if (pixelHits < pixelHitsCut_)
        quality = false;
      if (muonHits < muonHitsCut_)
        quality = false;
      ;
      if (!mu.isTrackerMuon())
        quality = false;
      if (nMatches < nMatchesCut_)
        quality = false;
      muon_sel[3] = quality;
      if (quality)
        number_of_goodMuons++;
  
      pt_before_->Fill(pt);
      eta_before_->Fill(eta);
      dxy_before_->Fill(dxy);
      muoncharge_before_->Fill(charge);
      goodewkmuon_before_->Fill(quality);
  
      // Charge asymmetry
      // if (quality) {
      //  if (charge>0) ptPlus_before_->Fill(pt);
      //  if (charge<0) ptMinus_before_->Fill(pt);
      //}
  
      // Isolation cuts
      double isovar = mu.isolationR03().sumPt;
      if (isCombinedIso_) {
        isovar += mu.isolationR03().emEt;
        isovar += mu.isolationR03().hadEt;
      }
      if (isRelativeIso_)
        isovar /= pt;
      if (isovar < isoCut03_)
        muon_sel[4] = true;
  
      LogTrace("") << "\t... isolation value" << isovar << ", isolated? " << muon_sel[6];
      iso_before_->Fill(isovar);
  
      // HLT (not mtched to muon for the time being)
      // if (trigger_fired) muon_sel[5] = true;
  
      // For Z:
      if (pt > ptThrForZ1_ && fabs(eta) < etaCut_ && fabs(dxy) < dxyCut_ && quality && isovar < isoCut03_) {
        muon4Z = true;
      }
  
      // MET/MT cuts
      double w_et = met_et + mu.pt();
      double w_px = met.px() + mu.px();
      double w_py = met.py() + mu.py();
  
      double massT = w_et * w_et - w_px * w_px - w_py * w_py;
      massT = (massT > 0) ? sqrt(massT) : 0;
  
      LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << w_px << ", " << w_py
                   << " [GeV]";
      if (massT > mtMin_ && massT < mtMax_)
        muon_sel[5] = true;
      mt_before_->Fill(massT);
      if (met_et > metMin_ && met_et < metMax_)
        muon_sel[6] = true;
  
      // Acoplanarity cuts
      Geom::Phi<double> deltaphi(mu.phi() - atan2(met.py(), met.px()));
      double acop = deltaphi.value();
      if (acop < 0)
        acop = -acop;
      acop = M_PI - acop;
      LogTrace("") << "\t... acoplanarity: " << acop;
      if (acop < acopCut_)
        muon_sel[7] = true;
      acop_before_->Fill(acop);
  
      // Remaining flags (from global event information)
      if (nmuonsForZ1 < 1 || nmuonsForZ2 < 2)
        muon_sel[8] = true;
      if (njets <= nJetMax_)
        muon_sel[9] = true;
  
      // Collect necessary flags "per muon"
      int flags_passed = 0;
      for (int j = 0; j < NFLAGS; ++j) {
        if (muon_sel[j])
          flags_passed += 1;
      }
  
      // Do N-1 histograms now (and only once for global event quantities)
      if (flags_passed >= (NFLAGS - 1)) {
        if (!muon_sel[0] || flags_passed == NFLAGS)
          pt_after_->Fill(pt);
        if (!muon_sel[1] || flags_passed == NFLAGS)
          eta_after_->Fill(eta);
        if (!muon_sel[2] || flags_passed == NFLAGS)
          dxy_after_->Fill(dxy);
        if (!muon_sel[3] || flags_passed == NFLAGS)
          goodewkmuon_after_->Fill(quality);
        if (!muon_sel[4] || flags_passed == NFLAGS)
          iso_after_->Fill(isovar);
        // if (!muon_sel[5] || flags_passed == NFLAGS)
        //   if (!hlt_hist_done) trig_after_->Fill(trigger_fired);
        //hlt_hist_done = true;
        if (!muon_sel[5] || flags_passed == NFLAGS)
          mt_after_->Fill(massT);
        if (!muon_sel[6] || flags_passed == NFLAGS)
          if (!met_hist_done)
            met_after_->Fill(met_et);
        met_hist_done = true;
        if (!muon_sel[7] || flags_passed == NFLAGS)
          acop_after_->Fill(acop);
        // no action here for muon_sel[8]
        if (!muon_sel[9] || flags_passed == NFLAGS) {
          if (!njets_hist_done) {
            njets_after_->Fill(njets);
            leadingjet_pt_after_->Fill(lead_jet_pt);
            leadingjet_eta_after_->Fill(lead_jet_eta);
          }
          njets_hist_done = true;
        }
        if (flags_passed == NFLAGS) {
          if (!wfullsel_hist_done) {
            npvs_after_->Fill(nvvertex);
            muoncharge_after_->Fill(charge);
            // if (charge>0) ptPlus_afterW_->Fill(pt);
            // if (charge<0) ptMinus_afterW_->Fill(pt);
          }
          wfullsel_hist_done = true;
        }
      }
  
      // The cases in which the event is rejected as a Z are considered
      // independently:
      if (muon4Z && !muon_sel[8]) {
        // Plots for 2 muons
        for (unsigned int j = i + 1; j < muonCollectionSize; j++) {
          for (int ij = 0; ij < NFLAGSZ; ++ij) {
            zmuon_sel[ij] = false;
          }
  
          for (int ji = 0; ji < 5; ++ji) {
            zmuon_sel[ji] = muon_sel[ji];
          }
  
          const Muon& mu2 = muonCollection->at(j);
          if (!mu2.isGlobalMuon())
            continue;
          if (mu2.charge() * charge != -1)
            continue;
          reco::TrackRef gm2 = mu2.globalTrack();
          reco::TrackRef tk2 = mu2.innerTrack();
          double pt2 = mu2.pt();
          if (pt2 > ptThrForZ2_)
            zmuon_sel[5] = true;
          double eta2 = mu2.eta();
          if (fabs(eta2) < etaCut_)
            zmuon_sel[6] = true;
          double dxy2 = gm2->dxy(beamSpotHandle->position());
          if (fabs(dxy2) < dxyCut_)
            zmuon_sel[7] = true;
          double normalizedChi22 = gm2->normalizedChi2();
          double trackerHits2 = tk2->hitPattern().numberOfValidTrackerHits();
          int pixelHits2 = tk2->hitPattern().numberOfValidPixelHits();
          int muonHits2 = gm2->hitPattern().numberOfValidMuonHits();
          int nMatches2 = mu2.numberOfMatches();
          bool quality2 = true;
          if (normalizedChi22 > normalizedChi2Cut_)
            quality2 = false;
          if (trackerHits2 < trackerHitsCut_)
            quality2 = false;
          if (pixelHits2 < pixelHitsCut_)
            quality2 = false;
          if (muonHits2 < muonHitsCut_)
            quality2 = false;
          if (!mu2.isTrackerMuon())
            quality2 = false;
          if (nMatches2 < nMatchesCut_)
            quality2 = false;
          zmuon_sel[8] = quality2;
          double isovar2 = mu2.isolationR03().sumPt;
          if (isCombinedIso_) {
            isovar2 += mu2.isolationR03().emEt;
            isovar2 += mu2.isolationR03().hadEt;
          }
          if (isRelativeIso_)
            isovar2 /= pt2;
          if (isovar2 < isoCut03_)
            zmuon_sel[9] = true;
          //  if (trigger_fired) zmuon_sel[10] = true;
          const math::XYZTLorentzVector ZRecoGlb(
              mu.px() + mu2.px(), mu.py() + mu2.py(), mu.pz() + mu2.pz(), mu.p() + mu2.p());
          if (ZRecoGlb.mass() > dimuonMassMin_ && ZRecoGlb.mass() < dimuonMassMax_)
            zmuon_sel[10] = true;
  
          // jet flag
          if (njets <= nJetMax_)
            zmuon_sel[11] = true;
  
          // start filling histos: N-1 plots
          int flags_passed_z = 0;
  
          for (int jj = 0; jj < NFLAGSZ; ++jj) {
            if (zmuon_sel[jj])
              ++flags_passed_z;
          }
  
          if (flags_passed_z >= (NFLAGSZ - 1)) {
            if (!zmuon_sel[0] || flags_passed_z == NFLAGSZ) {
              pt1_afterZ_->Fill(pt);
            }
            if (!zmuon_sel[1] || flags_passed_z == NFLAGSZ) {
              eta1_afterZ_->Fill(eta);
            }
            if (!zmuon_sel[2] || flags_passed_z == NFLAGSZ) {
              dxy1_afterZ_->Fill(dxy);
            }
            if (!zmuon_sel[3] || flags_passed_z == NFLAGSZ) {
              goodewkmuon1_afterZ_->Fill(quality);
            }
            if (!zmuon_sel[4] || flags_passed_z == NFLAGSZ) {
              iso1_afterZ_->Fill(isovar);
            }
            if (!zmuon_sel[5] || flags_passed_z == NFLAGSZ) {
              pt2_afterZ_->Fill(pt2);
            }
            if (!zmuon_sel[6] || flags_passed_z == NFLAGSZ) {
              eta2_afterZ_->Fill(eta2);
            }
            if (!zmuon_sel[7] || flags_passed_z == NFLAGSZ) {
              dxy2_afterZ_->Fill(dxy2);
            }
            if (!zmuon_sel[8] || flags_passed_z == NFLAGSZ) {
              goodewkmuon2_afterZ_->Fill(quality2);
            }
            if (!zmuon_sel[9] || flags_passed_z == NFLAGSZ) {
              iso2_afterZ_->Fill(isovar2);
            }
            //     if (!zmuon_sel[10] || flags_passed_z == NFLAGSZ) {
            //     if (!zhlt_hist_done) ztrig_afterZ_->Fill(trigger_fired);
            //   zhlt_hist_done = true;
            //  }
            if (!zmuon_sel[10] || flags_passed_z == NFLAGSZ) {
              dimuonmass_afterZ_->Fill(ZRecoGlb.mass());
            }
            if (!zmuon_sel[11] || flags_passed_z == NFLAGSZ) {
              if (!zjets_hist_done) {
                njets_afterZ_->Fill(njets);
                leadingjet_pt_afterZ_->Fill(lead_jet_pt);
                leadingjet_eta_afterZ_->Fill(lead_jet_eta);
              }
              zjets_hist_done = true;
            }
            if (flags_passed_z == NFLAGSZ) {
              met_afterZ_->Fill(met_et);
              if (!zfullsel_hist_done) {
                npvs_afterZ_->Fill(nvvertex);
                muoncharge_afterZ_->Fill(charge);
                if (charge > 0) {
                  // ptPlus_afterZ_->Fill(mu.pt());
                  // ptMinus_afterZ_->Fill(mu2.pt());
                  ptDiffPM_afterZ_->Fill(mu.pt() - mu2.pt());
                } else {
                  // ptPlus_afterZ_->Fill(mu2.pt());
                  // ptMinus_afterZ_->Fill(mu.pt());
                  ptDiffPM_afterZ_->Fill(mu2.pt() - mu.pt());
                }
              }
              zfullsel_hist_done = true;
            }
          }
        }
      }
    }
  
    if (zfullsel_hist_done) {
      // here was a Z candidate
      n_zselPt1thr_->Fill(nmuonsForZ1);
      n_zselPt2thr_->Fill(nmuonsForZ2);
    }
  
    // nmuons_->Fill(number_of_muons);
    // nmuons_->Fill(muonCollectionSize);
    ngoodmuons_->Fill(number_of_goodMuons);
  
    return;
  }
  
  // Local Variables:
  // show-trailing-whitespace: t
  // truncate-lines: t
  // End:
  