EwkMuDQM.cc

Go to the documentation of this file.

#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))
{
  isValidHltConfig_ = false;

  theDbe = Service<DQMStore>().operator->();
  theDbe->setCurrentFolder("Physics/EwkMuDQM");
  init_histograms();

}

void EwkMuDQM::beginRun(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_ = hltConfigProvider_.init( iRun, iSet, "HLT", isConfigChanged );

}

void EwkMuDQM::beginJob() {

}

void EwkMuDQM::init_histograms() {

  char chtitle[256] = "";

  pt_before_ = theDbe->book1D("PT_BEFORECUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);
  pt_after_ = theDbe->book1D("PT_AFTERWCUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);

  eta_before_ = theDbe->book1D("ETA_BEFORECUTS","Muon pseudo-rapidity",50,-2.5,2.5);
  eta_after_ = theDbe->book1D("ETA_AFTERWCUTS","Muon pseudo-rapidity",50,-2.5,2.5);

  dxy_before_ = theDbe->book1D("DXY_BEFORECUTS","Muon transverse distance to beam spot [cm]",100,-0.5,0.5);
  dxy_after_ = theDbe->book1D("DXY_AFTERWCUTS","Muon transverse distance to beam spot [cm]",100,-0.5,0.5);

  goodewkmuon_before_ = theDbe->book1D("GOODEWKMUON_BEFORECUTS","Quality-muon flag",2,-0.5,1.5);
  goodewkmuon_after_ = theDbe->book1D("GOODEWKMUON_AFTERWCUTS","Quality-muon flag",2,-0.5,1.5);

  if (isRelativeIso_) {
    if (isCombinedIso_) {
      iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Relative (combined) isolation variable",100, 0., 1.);
      iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Relative (combined) isolation variable",100, 0., 1.);
    } else {
      iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Relative (tracker) isolation variable",100, 0., 1.);
      iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Relative (tracker) isolation variable",100, 0., 1.);
    }
  } else {
    if (isCombinedIso_) {
      iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
      iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
    } else {
      iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
      iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
    }
  }

  trig_before_ = theDbe->book1D("TRIG_BEFORECUTS","Trigger response (boolean of muon triggers)",2,-0.5,1.5);
  trig_after_ = theDbe->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_ = theDbe->book1D("MT_BEFORECUTS",chtitle,150,0.,300.);
  mt_after_ = theDbe->book1D("MT_AFTERWCUTS",chtitle,150,0.,300.);

  snprintf(chtitle, 255, "Missing transverse energy (%s) [GeV]", metTag_.label().data());
  met_before_ = theDbe->book1D("MET_BEFORECUTS",chtitle,100,0.,200.);
  met_after_ = theDbe->book1D("MET_AFTERWCUTS",chtitle,100,0.,200.);
  met_afterZ_ = theDbe->book1D("MET_AFTERZCUTS",chtitle,100,0.,200.);

  snprintf(chtitle, 255, "MU-MET (%s) acoplanarity", metTag_.label().data());
  acop_before_ = theDbe->book1D("ACOP_BEFORECUTS",chtitle,50,0.,M_PI);
  acop_after_ = theDbe->book1D("ACOP_AFTERWCUTS",chtitle,50,0.,M_PI);

  snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ1_);
  n_zselPt1thr_ = theDbe->book1D("NZSELPT1THR",chtitle,10,-0.5,9.5);
  snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ2_);
  n_zselPt2thr_ = theDbe->book1D("NZSELPT2THR",chtitle,10,-0.5,9.5);

  snprintf(chtitle, 255, "Number of jets (%s) above %.2f GeV", jetTag_.label().data(), eJetMin_);
  njets_before_ = theDbe->book1D("NJETS_BEFORECUTS",chtitle,16,-0.5,15.5);
  njets_after_ = theDbe->book1D("NJETS_AFTERWCUTS",chtitle,16,-0.5,15.5);
  njets_afterZ_ = theDbe->book1D("NJETS_AFTERZCUTS",chtitle,16,-0.5,15.5);

  leadingjet_pt_before_ = theDbe->book1D("LEADINGJET_PT_BEFORECUTS","Leading Jet transverse momentum",300,0.,300.);
  leadingjet_pt_after_ = theDbe->book1D("LEADINGJET_PT_AFTERWCUTS","Leading Jet transverse momentum",300,0.,300.);
  leadingjet_pt_afterZ_ = theDbe->book1D("LEADINGJET_PT_AFTERZCUTS","Leading Jet transverse momentum",300,0.,300.);

  leadingjet_eta_before_ = theDbe->book1D("LEADINGJET_ETA_BEFORECUTS","Leading Jet pseudo-rapidity",50,-2.5,2.5);
  leadingjet_eta_after_ = theDbe->book1D("LEADINGJET_ETA_AFTERWCUTS","Leading Jet pseudo-rapidity",50,-2.5,2.5);
  leadingjet_eta_afterZ_ = theDbe->book1D("LEADINGJET_ETA_AFTERZCUTS","Leading Jet pseudo-rapidity",50,-2.5,2.5);

  //ptPlus_before_ = theDbe->book1D("PTPLUS_BEFORE_CUTS","Muon+ transverse momentum before cuts [GeV]",100,0.,100.);
  //ptMinus_before_ = theDbe->book1D("PTMINUS_BEFORE_CUTS","Muon- transverse momentum before cuts [GeV]",100,0.,100.);
  //ptPlus_afterW_ = theDbe->book1D("PTPLUS_AFTERW_CUTS","Muon+ transverse momentum after W cuts [GeV]",100,0.,100.);
  //ptMinus_afterW_ = theDbe->book1D("PTMINUS_AFTERW_CUTS","Muon- transverse momentum after W cuts [GeV]",100,0.,100.);
  //ptPlus_afterZ_ = theDbe->book1D("PTPLUS_AFTERZ_CUTS","Muon+ transverse momentum after Z cuts [GeV]",100,0.,100.);
  //ptMinus_afterZ_ = theDbe->book1D("PTMINUS_AFTERZ_CUTS","Muon- transverse momentum after Z cuts [GeV]",100,0.,100.);
  ptDiffPM_before_ = theDbe->book1D("PTDIFFPM_BEFORE_CUTS","pt(Muon+)-pt(Muon-) after Z cuts [GeV]",200,-100.,100.);
  ptDiffPM_afterZ_ = theDbe->book1D("PTDIFFPM_AFTERZ_CUTS","pt(Muon+)-pt(Muon-) after Z cuts [GeV]",200,-100.,100.);

  pt1_afterZ_ = theDbe->book1D("PT1_AFTERZCUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);
  eta1_afterZ_ = theDbe->book1D("ETA1_AFTERZCUTS","Muon pseudo-rapidity",50,-2.5,2.5);
  dxy1_afterZ_ = theDbe->book1D("DXY1_AFTERZCUTS","Muon transverse distance to beam spot [cm]",100,-0.5,0.5);
  goodewkmuon1_afterZ_ = theDbe->book1D("GOODEWKMUON1_AFTERZCUTS","Quality-muon flag",2,-0.5,1.5);

  if (isRelativeIso_) {
    if (isCombinedIso_) {
      iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Relative (combined) isolation variable",100, 0., 1.);
      iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Relative (combined) isolation variable",100, 0., 1.);
    } else {
      iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Relative (tracker) isolation variable",100, 0., 1.);
      iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Relative (tracker) isolation variable",100, 0., 1.);
    }
  } else {
    if (isCombinedIso_) {
      iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
      iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
    } else {
      iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
      iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
    }
  }

  pt2_afterZ_ = theDbe->book1D("PT2_AFTERZCUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);
  eta2_afterZ_ = theDbe->book1D("ETA2_AFTERZCUTS","Muon pseudo-rapidity",50,-2.5,2.5);
  dxy2_afterZ_ = theDbe->book1D("DXY2_AFTERZCUTS","Muon transverse distance to beam spot [cm]",100,-0.5,0.5);
  goodewkmuon2_afterZ_ = theDbe->book1D("GOODEWKMUON2_AFTERZCUTS","Quality-muon flag",2,-0.5,1.5);
  ztrig_afterZ_ = theDbe->book1D("ZTRIG_AFTERZCUTS","Trigger response (boolean of muon triggers)",2,-0.5,1.5);
  dimuonmass_before_= theDbe->book1D("DIMUONMASS_BEFORECUTS","DiMuonMass (2 globals)",100,0,200);
  dimuonmass_afterZ_= theDbe->book1D("DIMUONMASS_AFTERZCUTS","DiMuonMass (2 globals)",100,0,200);
  npvs_before_ = theDbe->book1D("NPVs_BEFORECUTS","Number of Valid Primary Vertices",51,-0.5,50.5);
  npvs_after_ = theDbe->book1D("NPVs_AFTERWCUTS","Number of Valid Primary Vertices",51,-0.5,50.5);
  npvs_afterZ_ = theDbe->book1D("NPVs_AFTERZCUTS","Number of Valid Primary Vertices",51,-0.5,50.5);
  muoncharge_before_ = theDbe->book1D("MUONCHARGE_BEFORECUTS","Muon Charge",3,-1.5,1.5);
  muoncharge_after_ = theDbe->book1D("MUONCHARGE_AFTERWCUTS","Muon Charge",3,-1.5,1.5);
  muoncharge_afterZ_ = theDbe->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_ = theDbe->book1D("NMuons","Number of muons in the event",10,-0.5,9.5);
  ngoodmuons_ = theDbe->book1D("NGoodMuons","Number of muons passing the quality criteria",10,-0.5,9.5);

  nph_ = theDbe->book1D("nph","Number of photons in the event",20,0.,20.);
  //npfph_ = theDbe->book1D("npfph","Number of PF photons in the event",20,0.,20.);
  phPt_ = theDbe->book1D("phPt","Photon transverse momentum [GeV]",100,0.,1000.);
  //pfphPt_ = theDbe->book1D("pfphPt","PF Photon transverse momentum [GeV]",1000,0.,1000.);
  snprintf(chtitle, 255, "Photon pseudorapidity (pT>%4.1f)",ptThrForPhoton_);
  phEta_ = theDbe->book1D("phEta",chtitle,100,-2.5,2.5);
  //pfphEta_ = theDbe->book1D("pfphEta","PF Photon pseudorapidity",100,-2.5,2.5);

}


void EwkMuDQM::endJob() {
}

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";

      /*  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=hltConfigProvider_.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 = 11;
      bool muon_sel[NFLAGS];
      const int NFLAGSZ = 13;
      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 && trigger_fired && 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[6] = true;
            mt_before_->Fill(massT);
            if (met_et>metMin_ && met_et<metMax_) muon_sel[7] = 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[8] = true;
            acop_before_->Fill(acop);

            // Remaining flags (from global event information)
            if (nmuonsForZ1<1 || nmuonsForZ2<2) muon_sel[9] = true;
            if (njets<=nJetMax_) muon_sel[10] = 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[6] || flags_passed==NFLAGS)
                        mt_after_->Fill(massT);
                  if (!muon_sel[7] || flags_passed==NFLAGS)
                        if (!met_hist_done) met_after_->Fill(met_et);
                        met_hist_done = true;
                  if (!muon_sel[8] || flags_passed==NFLAGS)
                        acop_after_->Fill(acop);
          // no action here for muon_sel[9]
                  if (!muon_sel[10] || 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[9]){

                   // 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[11] = true;

                    //jet flag
                    if (njets <=nJetMax_) zmuon_sel[12] = 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[11] || flags_passed_z==NFLAGSZ) {dimuonmass_afterZ_->Fill(ZRecoGlb.mass()); }
                       if (!zmuon_sel[12] || 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: