text: /*
   *  See header file for a description of this class.
   *
   *  \author Valentina Gori, University of Firenze
   */
  
  #include "DQM/Physics/src/EwkDQM.h"
  
  #include <vector>
  #include <string>
  #include <cmath>
  
  #include "DataFormats/Candidate/interface/Candidate.h"
  #include "FWCore/Framework/interface/Event.h"
  #include "FWCore/Framework/interface/EventSetup.h"
  #include "FWCore/ParameterSet/interface/ParameterSet.h"
  
  #include "DQMServices/Core/interface/DQMStore.h"
  #include "FWCore/ServiceRegistry/interface/Service.h"
  
  #include "DataFormats/Common/interface/Handle.h"
  
  #include "FWCore/MessageLogger/interface/MessageLogger.h"
  #include "FWCore/Common/interface/TriggerNames.h"
  
  // Physics Objects
  #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
  #include "DataFormats/MuonReco/interface/Muon.h"
  #include "DataFormats/JetReco/interface/Jet.h"
  #include "DataFormats/METReco/interface/MET.h"
  #include "DataFormats/TrackReco/interface/Track.h"
  #include "DataFormats/TrackReco/interface/TrackFwd.h"
  #include "DataFormats/VertexReco/interface/Vertex.h"
  #include "DataFormats/Math/interface/LorentzVector.h"
  
  #include "TLorentzVector.h"
  
  using namespace std;
  using namespace edm;
  using namespace reco;
  
  // EwkDQM::EwkDQM(const ParameterSet& parameters) {
  EwkDQM::EwkDQM(const ParameterSet& parameters) {
    eJetMin_ = parameters.getUntrackedParameter<double>("EJetMin", 999999.);
  
    // riguardare questa sintassi
    // Get parameters from configuration file
    thePFJetCollectionLabel_ = parameters.getParameter<InputTag>("PFJetCollection");
    theCaloMETCollectionLabel_ = parameters.getParameter<InputTag>("caloMETCollection");
    theTriggerResultsCollection_ = parameters.getParameter<InputTag>("triggerResultsCollection");
  
    theElecTriggerPathToPass_ = parameters.getParameter<std::vector<string> >("elecTriggerPathToPass");
    theMuonTriggerPathToPass_ = parameters.getParameter<std::vector<string> >("muonTriggerPathToPass");
    theTriggerResultsToken_ =
        consumes<edm::TriggerResults>(parameters.getParameter<InputTag>("triggerResultsCollection"));
    theMuonCollectionLabel_ = consumes<reco::MuonCollection>(parameters.getParameter<InputTag>("muonCollection"));
    theElectronCollectionLabel_ =
        consumes<reco::GsfElectronCollection>(parameters.getParameter<InputTag>("electronCollection"));
    thePFJetCollectionToken_ = consumes<edm::View<reco::Jet> >(parameters.getParameter<InputTag>("PFJetCollection"));
    theCaloMETCollectionToken_ = consumes<edm::View<reco::MET> >(parameters.getParameter<InputTag>("caloMETCollection"));
    theVertexToken_ = consumes<reco::VertexCollection>(parameters.getParameter<InputTag>("vertexCollection"));
  
    // just to initialize
    isValidHltConfig_ = false;
  
    h_vertex_number = nullptr;
    h_vertex_chi2 = nullptr;
    h_vertex_numTrks = nullptr;
    h_vertex_sumTrks = nullptr;
    h_vertex_d0 = nullptr;
  
    h_jet_count = nullptr;
    h_jet_et = nullptr;
    h_jet_pt = nullptr;
    h_jet_eta = nullptr;
    h_jet_phi = nullptr;
    h_jet2_et = nullptr;
    // h_jet2_pt = 0;
    h_jet2_eta = nullptr;
    h_jet2_phi = nullptr;
  
    h_e1_et = nullptr;
    h_e2_et = nullptr;
    h_e1_eta = nullptr;
    h_e2_eta = nullptr;
    h_e1_phi = nullptr;
    h_e2_phi = nullptr;
  
    h_m1_pt = nullptr;
    h_m2_pt = nullptr;
    h_m1_eta = nullptr;
    h_m2_eta = nullptr;
    h_m1_phi = nullptr;
    h_m2_phi = nullptr;
  
    // h_t1_et = 0;
    // h_t1_eta = 0;
    // h_t1_phi = 0;
  
    h_met = nullptr;
    h_met_phi = nullptr;
  
    h_e_invWMass = nullptr;
    h_m_invWMass = nullptr;
    h_mumu_invMass = nullptr;
    h_ee_invMass = nullptr;
  }
  
  EwkDQM::~EwkDQM() {}
  
  void EwkDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
    ibooker.setCurrentFolder("Physics/EwkDQM");
  
    char chtitle[256] = "";
    const size_t title_s = sizeof(chtitle);
  
    logTraceName = "EwkAnalyzer";
  
    LogTrace(logTraceName) << "Parameters initialization";
  
    const float pi = 4 * atan(1);
  
    // Keep the number of plots and number of bins to a minimum!
  
    h_vertex_number = ibooker.book1D("vertex_number", "Number of event vertices in collection", 10, -0.5, 9.5);
    h_vertex_chi2 = ibooker.book1D("vertex_chi2", "Event Vertex #chi^{2}/n.d.o.f.", 20, 0.0, 2.0);
    h_vertex_numTrks = ibooker.book1D("vertex_numTrks", "Event Vertex, number of tracks", 20, -0.5, 59.5);
    h_vertex_sumTrks = ibooker.book1D("vertex_sumTrks", "Event Vertex, sum of track pt", 20, 0.0, 100.0);
    h_vertex_d0 = ibooker.book1D("vertex_d0", "Event Vertex d0", 20, 0.0, 0.05);
  
    h_jet_count = ibooker.book1D("jet_count", chtitle, 8, -0.5, 7.5);
  
    snprintf(
        chtitle, title_s, "Leading jet E_{T} (from %s);E_{T}(1^{st} jet) (GeV)", thePFJetCollectionLabel_.label().data());
    h_jet_et = ibooker.book1D("jet_et", chtitle, 20, 0., 200.0);
  
    snprintf(chtitle,
             title_s,
             "Leading jet p_{T} (from %s);p_{T}(1^{st} jet) (GeV/c)",
             thePFJetCollectionLabel_.label().data());
    h_jet_pt = ibooker.book1D("jet_pt", chtitle, 20, 0., 200.0);
  
    snprintf(chtitle, title_s, "Leading jet #eta (from %s); #eta (1^{st} jet)", thePFJetCollectionLabel_.label().data());
    h_jet_eta = ibooker.book1D("jet_eta", chtitle, 20, -10., 10.0);
  
    snprintf(chtitle, title_s, "Leading jet #phi (from %s); #phi(1^{st} jet)", thePFJetCollectionLabel_.label().data());
    h_jet_phi = ibooker.book1D("jet_phi", chtitle, 22, -1.1 * pi, 1.1 * pi);
  
    snprintf(chtitle,
             title_s,
             "2^{nd} leading jet E_{T} (from %s);E_{T}(2^{nd} jet) (GeV)",
             thePFJetCollectionLabel_.label().data());
    h_jet2_et = ibooker.book1D("jet2_et", chtitle, 20, 0., 200.0);
  
    snprintf(chtitle,
             title_s,
             "2^{nd} leading jet #eta (from %s); #eta (2^{nd} jet)",
             thePFJetCollectionLabel_.label().data());
    h_jet2_eta = ibooker.book1D("jet2_eta", chtitle, 20, -10., 10.0);
  
    snprintf(
        chtitle, title_s, "2^{nd} leading jet #phi (from %s); #phi(2^{nd} jet)", thePFJetCollectionLabel_.label().data());
    h_jet2_phi = ibooker.book1D("jet2_phi", chtitle, 22, -1.1 * pi, 1.1 * pi);
  
    h_e1_et = ibooker.book1D("e1_et", "E_{T} of Leading Electron;E_{T} (GeV)", 20, 0.0, 100.0);
    h_e2_et = ibooker.book1D("e2_et", "E_{T} of Second Electron;E_{T} (GeV)", 20, 0.0, 100.0);
    h_e1_eta = ibooker.book1D("e1_eta", "#eta of Leading Electron;#eta", 20, -4.0, 4.0);
  
    h_e2_eta = ibooker.book1D("e2_eta", "#eta of Second Electron;#eta", 20, -4.0, 4.0);
  
    h_e1_phi = ibooker.book1D("e1_phi", "#phi of Leading Electron;#phi", 22, -1.1 * pi, 1.1 * pi);
    h_e2_phi = ibooker.book1D("e2_phi", "#phi of Second Electron;#phi", 22, -1.1 * pi, 1.1 * pi);
    h_m1_pt = ibooker.book1D("m1_pt", "p_{T} of Leading Muon;p_{T}(1^{st} #mu) (GeV)", 20, 0.0, 100.0);
    h_m2_pt = ibooker.book1D("m2_pt", "p_{T} of Second Muon;p_{T}(2^{nd} #mu) (GeV)", 20, 0.0, 100.0);
    h_m1_eta = ibooker.book1D("m1_eta", "#eta of Leading Muon;#eta(1^{st} #mu)", 20, -4.0, 4.0);
    h_m2_eta = ibooker.book1D("m2_eta", "#eta of Second Muon;#eta(2^{nd} #mu)", 20, -4.0, 4.0);
    h_m1_phi = ibooker.book1D(
        "m1_phi", "#phi of Leading Muon;#phi(1^{st} #mu)", 20, (-1. - 1. / 10.) * pi, (1. + 1. / 10.) * pi);
    h_m2_phi =
        ibooker.book1D("m2_phi", "#phi of Second Muon;#phi(2^{nd} #mu)", 20, (-1. - 1. / 10.) * pi, (1. + 1. / 10.) * pi);
  
    snprintf(chtitle, title_s, "Missing E_{T} (%s); GeV", theCaloMETCollectionLabel_.label().data());
    h_met = ibooker.book1D("met", chtitle, 20, 0.0, 100);
    h_met_phi =
        ibooker.book1D("met_phi", "Missing E_{T} #phi;#phi(MET)", 22, (-1. - 1. / 10.) * pi, (1. + 1. / 10.) * pi);
  
    h_e_invWMass = ibooker.book1D("we_invWMass", "W-> e #nu Transverse Mass;M_{T} (GeV)", 20, 0.0, 140.0);
    h_m_invWMass = ibooker.book1D("wm_invWMass", "W-> #mu #nu Transverse Mass;M_{T} (GeV)", 20, 0.0, 140.0);
    h_mumu_invMass = ibooker.book1D("z_mm_invMass", "#mu#mu Invariant Mass;InvMass (GeV)", 20, 40.0, 140.0);
    h_ee_invMass = ibooker.book1D("z_ee_invMass", "ee Invariant Mass;InvMass (Gev)", 20, 40.0, 140.0);
  }
  
  ///
  void EwkDQM::dqmBeginRun(const edm::Run& theRun, const edm::EventSetup& theSetup) {
    // passed as parameter to HLTConfigProvider::init(), not yet used
    bool isConfigChanged = false;
  
    // isValidHltConfig_ used to short-circuit analyze() in case of problems
    const std::string hltProcessName(theTriggerResultsCollection_.process());
    isValidHltConfig_ = hltConfigProvider_.init(theRun, theSetup, hltProcessName, isConfigChanged);
  }
  
  void EwkDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
    // short-circuit if hlt problems
    if (!isValidHltConfig_)
      return;
  
    LogTrace(logTraceName) << "Analysis of event # ";
    // Did it pass certain HLT path?
    Handle<TriggerResults> HLTresults;
    iEvent.getByToken(theTriggerResultsToken_, HLTresults);
    if (!HLTresults.isValid())
      return;
  
    const edm::TriggerNames& trigNames = iEvent.triggerNames(*HLTresults);
  
    // a temporary, until we have a list of triggers of interest
    std::vector<std::string> eleTrigPathNames;
    std::vector<std::string> muTrigPathNames;
  
    // eleTrigPathNames.push_back(theElecTriggerPathToPass_);
    // muTrigPathNames.push_back(theMuonTriggerPathToPass_);
    // end of temporary
  
    bool passed_electron_HLT = false;
    bool passed_muon_HLT = false;
    for (unsigned int i = 0; i < HLTresults->size(); i++) {
      const std::string& trigName = trigNames.triggerName(i);
      // check if triggerName matches electronPath
      for (unsigned int index = 0; index < theElecTriggerPathToPass_.size() && !passed_electron_HLT; index++) {
        // 0 if found, pos if not
        size_t trigPath = trigName.find(theElecTriggerPathToPass_[index]);
        if (trigPath == 0) {
          //      cout << "MuonTrigger passed (=trigName): " << trigName <<endl;
          passed_electron_HLT = HLTresults->accept(i);
        }
      }
      // check if triggerName matches muonPath
      for (unsigned int index = 0; index < theMuonTriggerPathToPass_.size() && !passed_muon_HLT; index++) {
        // 0 if found, pos if not
        size_t trigPath = trigName.find(theMuonTriggerPathToPass_[index]);
        if (trigPath == 0) {
          //      cout << "MuonTrigger passed (=trigName): " << trigName <<endl;
          passed_muon_HLT = HLTresults->accept(i);
        }
      }
    }
  
    // we are interested in events with a valid electron or muon
    if (!(passed_electron_HLT || passed_muon_HLT))
      return;
  
    ////////////////////////////////////////////////////////////////////////////////
    // Vertex information
    Handle<VertexCollection> vertexHandle;
    iEvent.getByToken(theVertexToken_, vertexHandle);
    if (!vertexHandle.isValid())
      return;
    VertexCollection vertexCollection = *(vertexHandle.product());
    VertexCollection::const_iterator v = vertexCollection.begin();
    int vertex_number = vertexCollection.size();
    double vertex_chi2 = v->normalizedChi2();  // v->chi2();
    double vertex_d0 = sqrt(v->x() * v->x() + v->y() * v->y());
    double vertex_numTrks = v->tracksSize();
    double vertex_sumTrks = 0.0;
    // std::cout << "vertex_d0=" << vertex_d0 << "\n";
    // double vertex_ndof    = v->ndof();cout << "ndof="<<vertex_ndof<<endl;
    for (Vertex::trackRef_iterator vertex_curTrack = v->tracks_begin(); vertex_curTrack != v->tracks_end();
         vertex_curTrack++)
      vertex_sumTrks += (*vertex_curTrack)->pt();
  
    ////////////////////////////////////////////////////////////////////////////////
    // Missing ET
    Handle<View<MET> > caloMETCollection;
    iEvent.getByToken(theCaloMETCollectionToken_, caloMETCollection);
    if (!caloMETCollection.isValid())
      return;
    float missing_et = caloMETCollection->begin()->et();
    float met_phi = caloMETCollection->begin()->phi();
  
    ////////////////////////////////////////////////////////////////////////////////
    // grab "gaussian sum fitting" electrons
    Handle<GsfElectronCollection> electronCollection;
    iEvent.getByToken(theElectronCollectionLabel_, electronCollection);
    if (!electronCollection.isValid())
      return;
  
    // Find the highest and 2nd highest electron
    float electron_et = -8.0;
    float electron_eta = -8.0;
    float electron_phi = -8.0;
    float electron2_et = -9.0;
    float electron2_eta = -9.0;
    float electron2_phi = -9.0;
    float ee_invMass = -9.0;
    TLorentzVector e1, e2;
  
    // If it passed electron HLT and the collection was found, find electrons near
    // Z mass
    if (passed_electron_HLT) {
      for (reco::GsfElectronCollection::const_iterator recoElectron = electronCollection->begin();
           recoElectron != electronCollection->end();
           recoElectron++) {
        // Require electron to pass some basic cuts
        if (recoElectron->et() < 20 || fabs(recoElectron->eta()) > 2.5)
          continue;
  
        // Tighter electron cuts
        if (recoElectron->deltaPhiSuperClusterTrackAtVtx() > 0.58 ||
            recoElectron->deltaEtaSuperClusterTrackAtVtx() > 0.01 || recoElectron->sigmaIetaIeta() > 0.027)
          continue;
  
        if (recoElectron->et() > electron_et) {
          electron2_et = electron_et;  // 2nd highest gets values from current highest
          electron2_eta = electron_eta;
          electron2_phi = electron_phi;
          electron_et = recoElectron->et();  // 1st highest gets values from new highest
          electron_eta = recoElectron->eta();
          electron_phi = recoElectron->phi();
          e1 = TLorentzVector(recoElectron->momentum().x(),
                              recoElectron->momentum().y(),
                              recoElectron->momentum().z(),
                              recoElectron->p());
        } else if (recoElectron->et() > electron2_et) {
          electron2_et = recoElectron->et();
          electron2_eta = recoElectron->eta();
          electron2_phi = recoElectron->phi();
          e2 = TLorentzVector(recoElectron->momentum().x(),
                              recoElectron->momentum().y(),
                              recoElectron->momentum().z(),
                              recoElectron->p());
        }
      }  // end of loop over electrons
      if (electron2_et > 0.0) {
        TLorentzVector pair = e1 + e2;
        ee_invMass = pair.M();
      }
    }  // end of "are electrons valid"
    ////////////////////////////////////////////////////////////////////////////////
  
    ////////////////////////////////////////////////////////////////////////////////
    // Take the STA muon container
    Handle<MuonCollection> muonCollection;
    iEvent.getByToken(theMuonCollectionLabel_, muonCollection);
    if (!muonCollection.isValid())
      return;
  
    // Find the highest pt muons
    float mm_invMass = -9.0;
    float muon_pt = -9.0;
    float muon_eta = -9.0;
    float muon_phi = -9.0;
    float muon2_pt = -9.0;
    float muon2_eta = -9.0;
    float muon2_phi = -9.0;
    TLorentzVector m1, m2;
  
    if (passed_muon_HLT) {
      for (reco::MuonCollection::const_iterator recoMuon = muonCollection->begin(); recoMuon != muonCollection->end();
           recoMuon++) {
        // Require muon to pass some basic cuts
        if (recoMuon->pt() < 20 || !recoMuon->isGlobalMuon())
          continue;
        // Some tighter muon cuts
        if (recoMuon->globalTrack()->normalizedChi2() > 10)
          continue;
  
        if (recoMuon->pt() > muon_pt) {
          muon2_pt = muon_pt;  // 2nd highest gets values from current highest
          muon2_eta = muon_eta;
          muon2_phi = muon_phi;
          muon_pt = recoMuon->pt();  // 1st highest gets values from new highest
          muon_eta = recoMuon->eta();
          muon_phi = recoMuon->phi();
          m1 =
              TLorentzVector(recoMuon->momentum().x(), recoMuon->momentum().y(), recoMuon->momentum().z(), recoMuon->p());
        } else if (recoMuon->pt() > muon2_pt) {
          muon2_pt = recoMuon->pt();
          muon2_eta = recoMuon->eta();
          muon2_phi = recoMuon->phi();
          m2 =
              TLorentzVector(recoMuon->momentum().x(), recoMuon->momentum().y(), recoMuon->momentum().z(), recoMuon->p());
        }
      }
    }
    if (muon2_pt > 0.0) {
      TLorentzVector pair = m1 + m2;
      mm_invMass = pair.M();
    }
    ////////////////////////////////////////////////////////////////////////////////
  
    ////////////////////////////////////////////////////////////////////////////////
    // Find the highest et jet
  
    //  Handle<CaloJetCollection> caloJetCollection;
    Handle<View<Jet> > PFJetCollection;
    iEvent.getByToken(thePFJetCollectionToken_, PFJetCollection);
    if (!PFJetCollection.isValid())
      return;
  
    unsigned int muonCollectionSize = muonCollection->size();
    // unsigned int jetCollectionSize = jetCollection->size();
    unsigned int PFJetCollectionSize = PFJetCollection->size();
    int jet_count = 0;
    // int LEADJET=-1;  double max_pt=0;
  
    float jet_et = -80.0;
    float jet_pt = -80.0;   // prova
    float jet_eta = -80.0;  // now USED
    float jet_phi = -80.0;  // now USED
    float jet2_et = -90.0;
    float jet2_eta = -90.0;  // now USED
    float jet2_phi = -90.0;  // now USED
    //  for (CaloJetCollection::const_iterator i_calojet =
    // caloJetCollection->begin();
    //       i_calojet != caloJetCollection->end(); i_calojet++) {
    //  for (PFJetCollection::const_iterator i_pfjet = PFJetCollection->begin();
    //       i_pfjet != PFJetCollection->end(); i_pfjet++) {
    //  float jet_current_et = i_calojet->et();
    //  float jet_current_et = i_pfjet->et();            // e` identico a jet.et()
    //    jet_count++;
  
    // cleaning: va messo prima del riempimento dell'istogramma // This is in
    // order to use PFJets
    for (unsigned int i = 0; i < PFJetCollectionSize; i++) {
      const Jet& jet = PFJetCollection->at(i);
      // la classe "jet" viene definita qui!!!
      double minDistance = 99999;
      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
      // se la distanza muone-cono del jet e` minore di 0.3, passo avanti e non
      // conteggio il mio jet
  
      // If it overlaps with ELECTRON, it is not a jet
      if (electron_et > 0.0 && fabs(jet.eta() - electron_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron_phi) < 0.2)
        continue;
      if (electron2_et > 0.0 && fabs(jet.eta() - electron2_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron2_phi) < 0.2)
        continue;
  
      // provo a cambiare la parte degli elettroni in modo simmetrico alla parte
      // per i muoni
  
      // ...
      // ...
  
      // if it has too low Et, throw away
      if (jet.et() < eJetMin_)
        continue;
      jet_count++;
  
      // ovvero: incrementa jet_count se:
      //   - non c'e un muone entro 0.3 di distanza dal cono del jet;
      //   - se il jet non si sovrappone ad un elettrone;
      //   - se l'energia trasversa e` maggiore della soglia impostata (15?)
  
      // if(jet.et()>max_pt) { LEADJET=i; max_pt=jet.et();}
      // se l'energia del jet e` maggiore di max_pt, diventa "i"
      // l'indice del jet piu` energetico e max_pt la sua energia
  
      // riguardare questo!!!
      // fino ad ora, jet_et era inizializzato a -8.0
      if (jet.et() > jet_et) {
        jet2_et = jet_et;  // 2nd highest jet gets et from current highest
        // perche` prende l'energia del primo jet??
        jet2_eta = jet_eta;  // now USED
        jet2_phi = jet_phi;  // now USED
        // jet_et = i_calojet->et(); // current highest jet gets
        // et from the new highest
        jet_et = jet.et();  // current highest jet gets et from the new highest
        // ah, ok! lo riaggiorna solo dopo!
        jet_pt = jet.pt();                          // e` il pT del leading jet
        jet_eta = jet.eta();                        // now USED
        jet_phi = jet.phi() * (Geom::pi() / 180.);  // now USED
      } else if (jet.et() > jet2_et) {
        //      jet2_et  = i_calojet->et();
        jet2_et = jet.et();
        //      jet2_eta = i_calojet->eta();  // UNUSED
        //      jet2_phi = i_calojet->phi();  // UNUSED
        jet2_eta = jet.eta();  // now USED
        jet2_phi = jet.phi();  // now USED
      }
      // questo elseif funziona
    }
    ////////////////////////////////////////////////////////////////////////////////
  
    ////////////////////////////////////////////////////////////////////////////////
    //                 Fill Histograms //
    ////////////////////////////////////////////////////////////////////////////////
  
    bool fill_e1 = false;
    bool fill_e2 = false;
    bool fill_m1 = false;
    bool fill_m2 = false;
    bool fill_met = false;
  
    // Was Z->ee found?
    if (ee_invMass > 0.0) {
      h_ee_invMass->Fill(ee_invMass);
      fill_e1 = true;
      fill_e2 = true;
    }
  
    // Was Z->mu mu found?
    if (mm_invMass > 0.0) {
      h_mumu_invMass->Fill(mm_invMass);
      fill_m1 = true;
      fill_m2 = true;
      h_jet2_et->Fill(jet2_et);
    }
  
    // Was W->e nu found?
    if (electron_et > 0.0 && missing_et > 20.0) {
      float dphiW = fabs(met_phi - electron_phi);
      float W_mt_e = sqrt(2 * missing_et * electron_et * (1 - cos(dphiW)));
      h_e_invWMass->Fill(W_mt_e);
      fill_e1 = true;
      fill_met = true;
    }
  
    // Was W->mu nu found?
    if (muon_pt > 0.0 && missing_et > 20.0) {
      float dphiW = fabs(met_phi - muon_phi);
      float W_mt_m = sqrt(2 * missing_et * muon_pt * (1 - cos(dphiW)));
      h_m_invWMass->Fill(W_mt_m);
      fill_m1 = true;
      fill_met = true;
    }
  
    if (jet_et > -10.0) {
      h_jet_et->Fill(jet_et);
      h_jet_count->Fill(jet_count);
    }
  
    if (jet_pt > 0.) {
      h_jet_pt->Fill(jet_pt);
    }
  
    if (jet_eta > -50.) {
      h_jet_eta->Fill(jet_eta);
    }
  
    if (jet_phi > -10.) {
      h_jet_phi->Fill(jet_phi);
    }
  
    if (jet2_et > -10.0) {
      h_jet2_et->Fill(jet2_et);
    }
  
    // if (jet2_pt>0.) {
    //  h_jet2_pt   ->Fill(jet2_pt);
    // }
  
    if (jet2_eta > -50.) {
      h_jet2_eta->Fill(jet2_eta);
    }
  
    if (jet2_phi > -10.) {
      h_jet2_phi->Fill(jet2_phi);
    }
  
    if (fill_e1 || fill_m1) {
      h_vertex_number->Fill(vertex_number);
      h_vertex_chi2->Fill(vertex_chi2);
      h_vertex_d0->Fill(vertex_d0);
      h_vertex_numTrks->Fill(vertex_numTrks);
      h_vertex_sumTrks->Fill(vertex_sumTrks);
    }
  
    if (fill_e1) {
      h_e1_et->Fill(electron_et);
      h_e1_eta->Fill(electron_eta);
      h_e1_phi->Fill(electron_phi);
    }
    if (fill_e2) {
      h_e2_et->Fill(electron2_et);
      h_e2_eta->Fill(electron2_eta);
      h_e2_phi->Fill(electron2_phi);
    }
    if (fill_m1) {
      h_m1_pt->Fill(muon_pt);
      h_m1_eta->Fill(muon_eta);
      h_m1_phi->Fill(muon_phi);
    }
    if (fill_m2) {
      h_m2_pt->Fill(muon2_pt);
      h_m2_eta->Fill(muon2_eta);
      h_m2_phi->Fill(muon2_phi);
    }
    if (fill_met) {
      h_met->Fill(missing_et);
      h_met_phi->Fill(met_phi);
    }
    ////////////////////////////////////////////////////////////////////////////////
  }
  
  // This always returns only a positive deltaPhi
  double EwkDQM::calcDeltaPhi(double phi1, double phi2) {
    double deltaPhi = phi1 - phi2;
  
    if (deltaPhi < 0)
      deltaPhi = -deltaPhi;
  
    if (deltaPhi > 3.1415926)
      deltaPhi = 2 * 3.1415926 - deltaPhi;
  
    return deltaPhi;
  }
  
  // Local Variables:
  // show-trailing-whitespace: t
  // truncate-lines: t
  // End:
  