MBUEandQCDValidation.cc

Go to the documentation of this file.

/*Class MBUEandQCDValidation
 *  
 *  Class to fill dqm monitor elements from existing EDM file
 *
 */

#include "MBUEandQCDValidation.h"
#include "Validation/EventGenerator/interface/HepMCValidationHelper.h"

#include "CLHEP/Units/defs.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include "CLHEP/Units/SystemOfUnits.h"

#include "fastjet/JetDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "Validation/EventGenerator/interface/DQMHelper.h"
using namespace edm;

MBUEandQCDValidation::MBUEandQCDValidation(const edm::ParameterSet& iPSet)
    : wmanager_(iPSet, consumesCollector()),
      hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")),
      genchjetCollection_(iPSet.getParameter<edm::InputTag>("genChjetsCollection")),
      genjetCollection_(iPSet.getParameter<edm::InputTag>("genjetsCollection")),
      verbosity_(iPSet.getUntrackedParameter<unsigned int>("verbosity", 0)) {
  hepmcGPCollection.reserve(initSize);
  hepmcCharge.reserve(initSize);

  theCalo = new CaloCellManager(verbosity_);

  eneInCell.resize(CaloCellManager::nCaloCell);

  hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
  genjetCollectionToken_ = consumes<reco::GenJetCollection>(genjetCollection_);
  genchjetCollectionToken_ = consumes<reco::GenJetCollection>(genchjetCollection_);
  fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
}

MBUEandQCDValidation::~MBUEandQCDValidation() { delete theCalo; }

void MBUEandQCDValidation::dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) {
  fPDGTable = c.getHandle(fPDGTableToken);
}

void MBUEandQCDValidation::bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) {
  DQMHelper dqm(&i);
  i.setCurrentFolder("Generator/MBUEandQCD");


  // Number of analyzed events
  nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., " ", "Number of events");

  // Number of events with no forward trigger
  nNoFwdTrig = dqm.book1dHisto(
      "nNoFwdTrig", "n Events no forward trigger", 1, 0., 1., " ", "Number of Events with no Forward Trigger");

  // Number of Events with a single arm forward trigger
  nSaFwdTrig = dqm.book1dHisto("nSaFwdTrig",
                               "n Events single arm forward trigger",
                               1,
                               0.,
                               1.,
                               " ",
                               "Number of Events with Single Arm Forward Trigger");

  // Number of Events with b quark
  nbquark = dqm.book1dHisto("nbquark", "n Events with b quark", 1, 0., 1., " ", "Number of Events with b Quarks");

  // Number of Events with c and b quark
  ncandbquark = dqm.book1dHisto(
      "ncandbquark", "n Events with c and b quark", 1, 0., 1., "", "Number of Events with c and b Quark");

  // Number of Events with c and no b quark
  ncnobquark = dqm.book1dHisto(
      "ncnobquark", "n Events with c and no b quark", 1, 0., 1., "", "Number of Events with c and no b Quark");

  // Number of selected events for QCD-09-010
  nEvt1 = dqm.book1dHisto(
      "nEvt1", "n Events QCD-09-010", 1, 0., 1., "", "Number of Events passing the QCD-09-010 selection");
  // dNchdpt QCD-09-010
  dNchdpt1 = dqm.book1dHisto("dNchdpt1",
                             "dNchdpt QCD-09-010",
                             30,
                             0.,
                             6.,
                             "P_{t}^{charged tracks QCD-09-010 selection} (GeV)",
                             "Number of Charged Tracks");
  // dNchdeta QCD-09-010
  dNchdeta1 = dqm.book1dHisto("dNchdeta1",
                              "dNchdeta QCD-09-010",
                              10,
                              -2.5,
                              2.5,
                              "#eta^{charged tracks QCD-09-010 selection}",
                              "Number of Charged Tracks");
  // Number of selected events for QCD-10-001

  nEvt2 = dqm.book1dHisto(
      "nEvt2", "n Events QCD-10-001", 1, 0., 1., "", "Number of Events passing the QCD-10-001 selection");
  // Leading track pt QCD-10-001
  leadTrackpt = dqm.book1dHisto("leadTrackpt",
                                "leading track pt QCD-10-001",
                                200,
                                0.,
                                100.,
                                "P_{t}^{lead track QCD-10-001 selection} (GeV)",
                                "Number of Events");
  // Leading track eta QCD-10-001
  leadTracketa = dqm.book1dHisto("leadTracketa",
                                 "leading track eta QCD-10-001",
                                 50.,
                                 -2.5,
                                 2.5,
                                 "#eta^{lead track QCD-10-001 selection}",
                                 "Number of Events");
  // transverse charged particle density vs leading track pt
  nChaDenLpt = i.bookProfile("nChaDenLpt", "charged density vs leading pt", 200, 0., 100., 0., 100., " ");
  // transverse charged particle density vs leading track pt
  sptDenLpt = i.bookProfile("sptDenLpt", "sum pt density vs leading pt", 200, 0., 100., 0., 300., " ");
  // dNchdpt QCD-10-001 transverse
  dNchdpt2 = dqm.book1dHisto("dNchdpt2",
                             "dNchdpt QCD-10-001",
                             200,
                             0.,
                             100.,
                             "P_{t}^{charged tracks QCD-10-001 selection} (GeV)",
                             "Number of Charged Tracks");
  // dNchdeta QCD-10-001 transverse
  dNchdeta2 = dqm.book1dHisto("dNchdeta2",
                              "dNchdeta QCD-10-001",
                              50,
                              -2.5,
                              2.5,
                              "#eta^{charged tracks QCD-10-001 selection}",
                              "Number of Charged Tracks");
  // nCha QCD-10-001 transverse
  nCha = dqm.book1dHisto(
      "nCha", "n charged QCD-10-001", 100, 0., 100., "N^{charged tracks QCD-10-001 selection}", "Number of Events");
  // dNchdSpt transverse
  dNchdSpt = dqm.book1dHisto("dNchdSpt",
                             "dNchdSpt QCD-10-001",
                             300,
                             0.,
                             300.,
                             "P_{t}^{charged trackes in transverse region QCD-10-001selection} (GeV)",
                             "Number of Charged Tracks");
  // dNchdphi
  dNchdphi = i.bookProfile("dNchdphi", "dNchdphi QCD-10-001", nphiBin, -180., 180., 0., 30., " ");
  // dSptdphi
  dSptdphi = i.bookProfile("dSptdphi", "dSptdphi QCD-10-001", nphiBin, -180., 180., 0., 30., " ");

  // number of charged jets QCD-10-001
  nChj = dqm.book1dHisto(
      "nChj", "n charged jets QCD-10-001", 30, 0, 30., "N^{charged jets QCD-10-001 selection}", "Number of Events");
  // dNchjdeta QCD-10-001
  dNchjdeta = dqm.book1dHisto("dNchjdeta",
                              "dNchjdeta QCD-10-001",
                              50,
                              -2.5,
                              2.5,
                              "#eta^{charged jets QCD-10-001 selection}",
                              "Number of charged Jets");
  // dNchjdpt QCD-10-001
  dNchjdpt = dqm.book1dHisto("dNchjdpt",
                             "dNchjdpt QCD-10-001",
                             100,
                             0.,
                             100.,
                             "P_{t}^{charged jets QCD-10-001 selection}",
                             "Number of charged Jets");
  // leading charged jet pt QCD-10-001
  leadChjpt = dqm.book1dHisto("leadChjpt",
                              "leadChjpt QCD-10-001",
                              100,
                              0.,
                              100.,
                              "P_{t}^{lead charged jet QCD-10-001 selection}",
                              "Number of charged Jets");
  // leading charged jet eta QCD-10-001
  leadChjeta = dqm.book1dHisto("leadChjeta",
                               "leadChjeta QCD-10-001",
                               50,
                               -2.5,
                               2.5,
                               "#eta^{lead charged jet QCD-10-001 selection}",
                               "Number of charged Jets");
  // (pt1+pt2)/ptot
  pt1pt2optotch = i.bookProfile("pt1pt2optotch", "sum 2 leading jets over ptot", 50, 0., 100., 0., 1., " ");

  // particle rates in tracker acceptance
  nPPbar = dqm.book1dHisto(
      "nPPbar", "nPPbar QCD-10-001", 30, 0., 30., "N_{p/#bar{p}}^{QCD-10-001 selection}", "Number of p/#bar{p}");
  nKpm = dqm.book1dHisto(
      "nKpm", "nKpm QCD-10-001", 30, 0., 30., "N_{K^{#pm}}^{QCD-10-001 selection}", "Number of K^{#pm}");
  nK0s = dqm.book1dHisto("nK0s", "nK0s QCD-10-001", 30, 0., 30., "N_{K^{0}}^{QCD-10-001 selection}", "Number of K^{0}");
  nL0 = dqm.book1dHisto(
      "nL0", "nL0 QCD-10-001", 30, 0., 30., "N_{#Lambda^{0}}^{QCD-10-001 selection}", "Number of #Lambda^{0}");
  nXim = dqm.book1dHisto("nXim", "nXim QCD-10-001", 30, 0., 30., "N_{#Xi}^{QCD-10-001 selection}", "Number of #Xi");
  nOmega = dqm.book1dHisto(
      "nOmega", "nOmega QCD-10-001", 30, 0., 30., "N_{#Omega^{#pm}}^{QCD-10-001 selection}", "Number of #Omega^{#pm}");

  pPPbar = dqm.book1dHisto("pPPbar",
                           "Log10(pt) PPbar QCD-10-001",
                           25,
                           -2.,
                           3.,
                           "log_{10}(P_{t}^{p/#bar{p} QCD-10-001 selection}) (log_{10}(GeV))",
                           "Number of p/#bar{p}");
  pKpm = dqm.book1dHisto("pKpm",
                         "Log10(pt) Kpm QCD-10-001",
                         25,
                         -2.,
                         3.,
                         "log_{10}(P_{t}^{K^{#pm} QCD-10-001 selection}) (log_{10}(GeV))",
                         "Number of K^{#pm}");
  pK0s = dqm.book1dHisto("pK0s",
                         "Log10(pt) K0s QCD-10-001",
                         25,
                         -2.,
                         3.,
                         "log_{10}(P_{t}^{K^{0} QCD-10-001 selection}) (log_{10}(GeV))",
                         "Number of K^{0}");
  pL0 = dqm.book1dHisto("pL0",
                        "Log10(pt) L0 QCD-10-001",
                        25,
                        -2.,
                        3.,
                        "log_{10}(P_{t}^{#Lambda^{0} QCD-10-001 selection}) (log_{10}(GeV))",
                        "Number of #Lambda^{0}");
  pXim = dqm.book1dHisto("pXim",
                         "Log10(pt) Xim QCD-10-001",
                         25,
                         -2.,
                         3.,
                         "log_{10}(P_{t}^{#Xi^{#pm} QCD-10-001 selection}) (log_{10}(GeV))",
                         "Number of #Xi");
  pOmega = dqm.book1dHisto("pOmega",
                           "Log10(pt) Omega QCD-10-001",
                           25,
                           -2.,
                           3.,
                           "log_{10}(P_{t}^{#Omega^{#pm} QCD-10-001 selection}) (log_{10}(GeV))",
                           "Number of #Omega^{#pm}");

  // neutral rate in the barrel + HF acceptance
  nNNbar = dqm.book1dHisto(
      "nNNbar", "nNNbar QCD-10-001", 30, 0., 30., "N_{n/#bar{n}}^{QCD-10-001 selection}", "Number of Events");
  nGamma = dqm.book1dHisto(
      "nGamma", "nGamma QCD-10-001", 50, 0., 200., "N_{#gamma}^{QCD-10-001 selection}", "Number of Events");

  pNNbar = dqm.book1dHisto("pNNbar",
                           "Log10(pt) NNbar QCD-10-001",
                           25,
                           -2.,
                           3.,
                           "log_{10}(P_{t}^{n/#bar{n} QCD-10-001 selection}) (log_{10}(GeV))",
                           "Number of n/#bar{n}");
  pGamma = dqm.book1dHisto("pGamma",
                           "Log10(pt) Gamma QCD-10-001",
                           25,
                           -2.,
                           3.,
                           "log_{10}(P_{t}^{#gamma QCD-10-001 selection}) (log_{10}(GeV))",
                           "Number of #gamma");

  // highest pt electron spectrum
  elePt = dqm.book1dHisto("elePt",
                          "highest pt electron Log10(pt)",
                          30,
                          -2.,
                          4.,
                          "log_{10}(P_{t}^{highest e} (log_{10}(GeV))",
                          "Number of Events");

  // highest pt muon spectrum
  muoPt = dqm.book1dHisto("muoPt",
                          "highest pt muon Log10(pt)",
                          30,
                          -2.,
                          4.,
                          "log_{10}(P_{t}^{highest #mu}  (log_{10}(GeV))",
                          "Number of Events");

  // number of selected di-jet events
  nDijet = dqm.book1dHisto(
      "nDijet", "n Dijet Events", 1, 0., 1., " ", "Number of Events Passing Di-jet JME-10-001 Selection");
  // number of jets
  nj = dqm.book1dHisto("nj", "n jets ", 30, 0, 30., "N_{jets}^{JME-10-001}", "Number of Events");
  // dNjdeta
  dNjdeta = dqm.book1dHisto("dNjdeta", "dNjdeta ", 50, -5., 5., "#eta_{jets}^{JME-10-001 selection}", "Number of Jets");
  // dNjdpt
  dNjdpt = dqm.book1dHisto("dNjdpt", "dNjdpt ", 60, 0., 300., "P_{t}^{jets JME-10-001 selection}", "Number of Jets");
  // (pt1+pt2)/ptot
  pt1pt2optot = i.bookProfile("pt1pt2optot", "sum 2 leading jets over Et tot ", 60, 0., 300., 0., 1., " ");
  // pt1-pt2
  pt1pt2balance =
      dqm.book1dHisto("pt1pt2balance",
                      "2 leading jets pt difference ",
                      10,
                      0.,
                      1.,
                      "#frac{P_{t}^{1st jet}-P_{t}^{2nd jet}}{P_{t}^{1st jet}+P_{t}^{2nd jet}}^{JME-10-001 selection}",
                      "Number of Di-jet Events");
  // pt1 pt2 Delta phi
  pt1pt2Dphi = dqm.book1dHisto("pt1pt2Dphi",
                               "pt1 pt2 delta phi ",
                               nphiBin,
                               0.,
                               180.,
                               "#Delta#phi(jet^{1st},jet^{2nd})^{JME-10-001 selection} (rad)",
                               "Number of Di-jet Events");
  // pt1 pt2 invariant mass
  pt1pt2InvM = dqm.book1dHisto("pt1pt2InvM",
                               "pt1 pt2 invariant mass ",
                               60,
                               0.,
                               600.,
                               "M_{di-jet}^{JME-10-001 selection}",
                               "Number of di-jet events");
  // pt3 fraction
  pt3Frac = dqm.book1dHisto("pt3Frac",
                            "2 pt3 over pt1+pt2 ",
                            30,
                            0.,
                            1.,
                            "#frac{P_{t}^{3rd jet}}{P_{t}^{1st jet}+P_{t}^{2nd jet}}^{JME-10-001 selection}",
                            "Number of 3rd Jets");
  // sum of jets Et
  sumJEt = dqm.book1dHisto(
      "sumJEt", "sum Jet Et ", 60, 0., 300., "#Sigma E_{t}^{jets JME-10-001 selection}", "Number of di-jet events");
  // fraction of missing Et over sum of jets Et
  missEtosumJEt = dqm.book1dHisto("missEtosumJEt",
                                  "missing Et over sumJet Et ",
                                  30,
                                  0.,
                                  1.,
                                  "E_{t}^{miss}/#Sigma E_{t}^{jets JME-10-001 selection}",
                                  "Number of Di-jet Events");
  // sum of final state particle Pt
  sumPt = dqm.book1dHisto("sumPt",
                          "sum particle Pt ",
                          60,
                          0.,
                          600.,
                          "#Sigma P_{t}^{particles passing JME-10-001 selection}",
                          "Number of jets");
  // sum of final state charged particle Pt
  sumChPt = dqm.book1dHisto("sumChPt",
                            "sum charged particle Pt ",
                            60,
                            0.,
                            300.,
                            "#Sigma P_{t}^{charged particles passing JME-10-001 selection}",
                            "Number of Jets");

  //Number of selected events for the HF energy flux analysis
  nHFflow = dqm.book1dHisto("nHFflow",
                            "n HF flow events",
                            1,
                            0.,
                            1.,
                            " ",
                            "Number of Events passing JME-10-001, FWD-10-002 and Jet-Multiplicity selection");
  //Forward energy flow for MinBias BSC selection
  dEdetaHFmb = i.bookProfile("dEdetaHFmb",
                             "dEdeta HF MinBias",
                             (int)CaloCellManager::nForwardEta,
                             0,
                             (double)CaloCellManager::nForwardEta,
                             0.,
                             300.,
                             " ");
  //Forward energy flow for QCD dijet selection
  dEdetaHFdj = i.bookProfile("dEdetaHFdj",
                             "dEdeta HF QCD dijet",
                             (int)CaloCellManager::nForwardEta,
                             0,
                             (double)CaloCellManager::nForwardEta,
                             0.,
                             300.,
                             " ");

  // FWD-10-001 like diffraction analysis
  nHFSD = dqm.book1dHisto(
      "nHFSD", "n single diffraction in HF", 1, 0., 1., " ", "Number of single diffraction in HF FWD-10-001 selection");
  // E-pz HF-
  EmpzHFm = dqm.book1dHisto("EmpzHFm",
                            "E-pz HF- SD",
                            40,
                            0.,
                            200.,
                            "#Sigma E_{cal. cells/corrected for the longitudinal mometum}^{FWD-10-001 selection}",
                            "Number of Events");
  // Number of cells above threshold
  ntHFm = dqm.book1dHisto("ntHFm",
                          "number of HF- tower SD",
                          20,
                          0.,
                          20.,
                          " N_{cells over threshold for single diffraction in HF Towers}^{FWD-10-001 selection}",
                          "Number of Events");
  // Energy in HF-
  eneHFmSel = dqm.book1dHisto(
      "eneHFmSel", "energy in HF-", 40, 0., 200., "#Sigma E_{cal. cells}^{FWD-10-001 selection}", "Number of Events");

  // number of jets accepted in the 'Jet-Multiplicity' analysis
  _JM25njets = dqm.book1dHisto("JM25njets", "n jets", 15, 0, 15., "Number of JM25 Jets", "Number of Events");
  _JM25ht = dqm.book1dHisto("JM25ht", "HT", 80, 0, 800., "H_{t}^{JM25} (GeV)", "Number of Events");
  _JM25pt1 = dqm.book1dHisto("JM25pt1", "pt", 40, 0, 200., "P_{t}^{JM25,1st Jet} (GeV)", "Number of JM25 Jets");
  _JM25pt2 = dqm.book1dHisto("JM25pt2", "pt", 40, 0, 200., "P_{t}^{JM25,2nd Jet} (GeV)", "Number of JM25 Jets");
  _JM25pt3 = dqm.book1dHisto("JM25pt3", "pt", 40, 0, 200., "P_{t}^{JM25,3rd Jet} (GeV)", "Number of JM25 Jets");
  _JM25pt4 = dqm.book1dHisto("JM25pt4", "pt", 40, 0, 200., "P_{t}^{JM25,4th Jet} (GeV)", "Number of JM25 Jets");

  _JM80njets = dqm.book1dHisto("JM80njets", "n jets", 15, 0, 15., "Number of JM80 Jets", "Number of Events");
  _JM80ht = dqm.book1dHisto("JM80ht", "HT", 80, 300, 1100., "H_{t}^{JM80} (GeV", "Number of Events");
  _JM80pt1 = dqm.book1dHisto("JM80pt1", "pt", 40, 60, 260., "P_{t}^{JM80,1st Jet} (GeV)", "Number of JM80 Jets");
  _JM80pt2 = dqm.book1dHisto("JM80pt2", "pt", 40, 60, 260., "P_{t}^{JM80,2nd Jet} (GeV)", "Number of JM80 Jets");
  _JM80pt3 = dqm.book1dHisto("JM80pt3", "pt", 40, 60, 260., "P_{t}^{JM80,3rd Jet} (GeV)", "Number of JM80 Jets");
  _JM80pt4 = dqm.book1dHisto("JM80pt4", "pt", 40, 60, 260., "P_{t}^{JM80,4th Jet} (GeV)", "Number of JM80 Jets");

  // differential jet rates
  djr10 = dqm.book1dHisto("djr10",
                          "Differential Jet Rate 1#rightarrow0",
                          60,
                          -1.,
                          5.,
                          "log_{10}(d_{min}(n,n+1)) for n=0",
                          "Number of Events");
  djr21 = dqm.book1dHisto("djr21",
                          "Differential Jet Rate 2#rightarrow1",
                          60,
                          -1.,
                          5.,
                          "log_{10}(d_{min}(n,n+1)) for n=1",
                          "Number of Events");
  djr32 = dqm.book1dHisto("djr32",
                          "Differential Jet Rate 3#rightarrow2",
                          60,
                          -1.,
                          5.,
                          "log_{10}(d_{min}(n,n+1)) for n=2",
                          "Number of Events");
  djr43 = dqm.book1dHisto("djr43",
                          "Differential Jet Rate 4#rightarrow3",
                          60,
                          -1.,
                          5.,
                          "log_{10}(d_{min}(n,n+1)) for n=3",
                          "Number of Events");

  // sumET analysis
  _sumEt = dqm.book1dHisto(
      "sumET", "Sum of stable particles Et", 150, 0, 600., "#Sigma E_{t}^{stable particles}", "Number of Events");
  _sumEt1 = dqm.book1dHisto("sumET1",
                            "Sum of stable particles Et (eta<0.5)",
                            150,
                            0,
                            200.,
                            "#Sigma E_{t}^{stable particles (#eta<0.5)}",
                            "Number of Events");
  _sumEt2 = dqm.book1dHisto("sumET2",
                            "Sum of stable particles Et (0.5<eta<1.0)",
                            150,
                            0,
                            200.,
                            "#Sigma E_{t}^{stable particles (0.5<#eta<1.0)}",
                            "Number of Events");
  _sumEt3 = dqm.book1dHisto("sumET3",
                            "Sum of stable particles Et (1.0<eta<1.5)",
                            150,
                            0,
                            200.,
                            "#Sigma E_{t}^{stable particles (1.0<#eta<1.5)}",
                            "Number of Events");
  _sumEt4 = dqm.book1dHisto("sumET4",
                            "Sum of stable particles Et (1.5<eta<2.0)",
                            150,
                            0,
                            200.,
                            "#Sigma E_{t}^{stable particles (1.5<#eta<2.0)}",
                            "Number of Events");
  _sumEt5 = dqm.book1dHisto("sumET5",
                            "Sum of stable particles Et (2.0<eta<5.0)",
                            150,
                            0,
                            200.,
                            "#Sigma E_{t}^{stable particles (2.0<#eta<5.0)}",
                            "Number of Events");

  return;
}

void MBUEandQCDValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  edm::Handle<HepMCProduct> evt;
  iEvent.getByToken(hepmcCollectionToken_, evt);

  //Get HepMC EVENT
  HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));

  double weight = wmanager_.weight(iEvent);

  if (verbosity_ > 0) {
    myGenEvent->print();
  }

  double binW = 1.;

  hepmcGPCollection.clear();
  hepmcCharge.clear();
  for (unsigned int i = 0; i < eneInCell.size(); i++) {
    eneInCell[i] = 0.;
  }

  nEvt->Fill(0.5, weight);

  //Looping through HepMC::GenParticle collection to search for status 1 particles
  double charge = 0.;
  unsigned int nb = 0;
  unsigned int nc = 0;
  for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin();
       iter != myGenEvent->particles_end();
       ++iter) {
    if (std::fabs((*iter)->pdg_id()) == 4) {
      nc++;
    }
    if (std::fabs((*iter)->pdg_id()) == 5) {
      nb++;
    }
    if ((*iter)->status() == 1) {
      hepmcGPCollection.push_back(*iter);
      const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID((*iter)->pdg_id()));
      if (PData == nullptr) {
        charge = -999.;
      } else
        charge = PData->charge();

      hepmcCharge.push_back(charge);

      if (verbosity_ > 0) {
        std::cout << "HepMC " << std::setw(14) << std::fixed << (*iter)->barcode() << std::setw(14) << std::fixed
                  << (*iter)->pdg_id() << std::setw(14) << std::fixed << (*iter)->momentum().perp() << std::setw(14)
                  << std::fixed << (*iter)->momentum().eta() << std::setw(14) << std::fixed << (*iter)->momentum().phi()
                  << std::endl;
      }
    }
  }

  int nBSCp = 0;
  int nBSCm = 0;
  double eneHFp = 0.;
  double eneHFm = 0.;
  int nChapt05 = 0;
  int nChaVtx = 0;
  for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
    if (!isNeutrino(i)) {
      // BSC trigger

      if (hepmcGPCollection[i]->momentum().eta() > 3.23 && hepmcGPCollection[i]->momentum().eta() < 4.65) {
        nBSCp++;
      }
      if (hepmcGPCollection[i]->momentum().eta() < -3.23 && hepmcGPCollection[i]->momentum().eta() > -4.65) {
        nBSCm++;
      }

      // number of charged particles in different selections

      if (std::fabs(hepmcGPCollection[i]->momentum().eta()) < 2.5 && hepmcGPCollection[i]->momentum().perp() > 0.5 &&
          isCharged(i)) {
        nChapt05++;
      }
      if (std::fabs(hepmcGPCollection[i]->momentum().eta()) < 2.5 && hepmcGPCollection[i]->momentum().perp() > 0.1 &&
          isCharged(i)) {
        nChaVtx++;
      }
      unsigned int theIndex = theCalo->getCellIndexFromAngle(hepmcGPCollection[i]->momentum().eta(),
                                                             hepmcGPCollection[i]->momentum().phi());
      if (theIndex < CaloCellManager::nCaloCell)
        eneInCell[theIndex] += hepmcGPCollection[i]->momentum().rho();
    }
  }

  // Forward calorimeters energy

  for (unsigned int icell = CaloCellManager::nBarrelCell + CaloCellManager::nEndcapCell;
       icell < CaloCellManager::nCaloCell;
       icell++) {
    if (theCalo->getCellFromIndex(icell)->getEtaMin() < 0.) {
      eneHFm += eneInCell[icell];
    } else {
      eneHFp += eneInCell[icell];
    }
  }

  // QCD-09-010 selection
  bool sel1 = false;
  if ((nBSCp > 0 || nBSCm > 0) && eneHFp >= 3. && eneHFm >= 3.) {
    sel1 = true;
  }

  // QCD-10-001 selection
  bool sel2 = false;
  if ((nBSCp > 0 || nBSCm > 0) && nChaVtx >= 3 && nChapt05 > 1) {
    sel2 = true;
  }

  // no forward trigger selection
  bool sel3 = false;
  if (nBSCp == 0 && nBSCm == 0) {
    sel3 = true;
  }

  // single arm forward trigger selection
  bool sel4 = false;
  if ((nBSCp > 0 && nBSCm == 0) || (nBSCm > 0 && nBSCp == 0)) {
    sel4 = true;
  }

  // BSC selection
  bool sel5 = false;
  if (nBSCp > 0 && nBSCm > 0) {
    sel5 = true;
  }

  // basic JME-10-001, FWD-10-002 and Jet-Multiplicity selection
  bool sel6 = false;
  if (sel5 && nChaVtx > 3) {
    sel6 = true;
  }

  // FWD-10-001 selection
  bool sel7 = false;
  if (nChaVtx >= 3 && nBSCm > 0 && eneHFp < 8.) {
    sel7 = true;
  }

  // Fill selection histograms
  if (sel1)
    nEvt1->Fill(0.5, weight);
  if (sel2)
    nEvt2->Fill(0.5, weight);
  if (sel3)
    nNoFwdTrig->Fill(0.5, weight);
  if (sel4)
    nSaFwdTrig->Fill(0.5, weight);
  if (sel6)
    nHFflow->Fill(0.5, weight);
  if (sel7)
    nHFSD->Fill(0.5, weight);

  if (nb > 0)
    nbquark->Fill(0.5, weight);
  if (nb > 0 && nc > 0)
    ncandbquark->Fill(0.5, weight);
  if (nb == 0 && nc > 0)
    ncnobquark->Fill(0.5, weight);

  // track analyses
  double ptMax = 0.;
  unsigned int iMax = 0;
  double ptot = 0.;
  unsigned int ppbar = 0;
  unsigned int nnbar = 0;
  unsigned int kpm = 0;
  unsigned int k0s = 0;
  unsigned int l0 = 0;
  unsigned int gamma = 0;
  unsigned int xim = 0;
  unsigned int omega = 0;
  unsigned int ele = 0;
  unsigned int muo = 0;
  unsigned int eleMax = 0;
  unsigned int muoMax = 0;

  std::vector<double> hfMB(CaloCellManager::nForwardEta, 0);
  std::vector<double> hfDJ(CaloCellManager::nForwardEta, 0);

  for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
    double eta = hepmcGPCollection[i]->momentum().eta();
    double pt = hepmcGPCollection[i]->momentum().perp();
    int pdgId = hepmcGPCollection[i]->pdg_id();
    if (isCharged(i) && std::fabs(eta) < 2.5) {
      if (sel1) {
        // QCD-09-010
        binW = dNchdpt1->getTH1()->GetBinWidth(1);
        dNchdpt1->Fill(pt, 1. / binW);  // weight to account for the pt bin width
        binW = dNchdeta1->getTH1()->GetBinWidth(1);
        dNchdeta1->Fill(eta, 1. / binW);  // weight to account for the eta bin width
      }
      // search for the leading track QCD-10-001
      if (sel2) {
        if (pt > ptMax) {
          ptMax = pt;
          iMax = i;
        }
        ptot += pt;

        // identified charged particle
        if (std::abs(pdgId) == 2212) {
          ppbar++;
          pPPbar->Fill(std::log10(pt), weight);
        } else if (std::abs(pdgId) == 321) {
          kpm++;
          pKpm->Fill(std::log10(pt), weight);
        } else if (std::abs(pdgId) == 3312) {
          xim++;
          pXim->Fill(std::log10(pt), weight);
        } else if (std::abs(pdgId) == 3334) {
          omega++;
          pOmega->Fill(std::log10(pt), weight);
        } else if (std::abs(pdgId) == 11) {
          ele++;
          eleMax = i;
        } else if (std::abs(pdgId) == 13) {
          muo++;
          muoMax = i;
        }
      }
    } else if (sel2 && isNeutral(i) && std::fabs(eta) < 2.5) {
      if (std::abs(pdgId) == 310) {
        k0s++;
        pK0s->Fill(std::log10(pt), weight);
      } else if (std::abs(pdgId) == 3122) {
        l0++;
        pL0->Fill(std::log10(pt), weight);
      }
    } else if (sel2 && isNeutral(i) && std::fabs(eta) < 5.19) {
      if (std::abs(pdgId) == 2112) {
        nnbar++;
        pNNbar->Fill(std::log10(pt), weight);
      } else if (std::abs(pdgId) == 22) {
        gamma++;
        pGamma->Fill(std::log10(pt), weight);
      }
    }
    unsigned int iBin = getHFbin(eta);
    if (sel6 && !isNeutrino(i) && iBin < CaloCellManager::nForwardEta) {
      hfMB[iBin] += hepmcGPCollection[i]->momentum().rho();
    }
  }
  nPPbar->Fill(ppbar, weight);
  nNNbar->Fill(nnbar, weight);
  nKpm->Fill(kpm, weight);
  nK0s->Fill(k0s, weight);
  nL0->Fill(l0, weight);
  nXim->Fill(xim, weight);
  nOmega->Fill(omega, weight);
  nGamma->Fill(gamma, weight);

  if (ele > 0)
    elePt->Fill(std::log10(hepmcGPCollection[eleMax]->momentum().perp()), weight);
  if (muo > 0)
    muoPt->Fill(std::log10(hepmcGPCollection[muoMax]->momentum().perp()), weight);

  leadTrackpt->Fill(hepmcGPCollection[iMax]->momentum().perp(), weight);
  leadTracketa->Fill(hepmcGPCollection[iMax]->momentum().eta(), weight);

  std::vector<double> theEtaRanges(theCalo->getEtaRanges());

  for (unsigned int i = 0; i < CaloCellManager::nForwardEta; i++) {
    binW = theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i + 1] -
           theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i];
    dEdetaHFmb->Fill(i + 0.5, hfMB[i] / binW);
  }

  // FWD-10-001

  if (sel7) {
    double empz = 0.;
    unsigned int nCellOvTh = 0;
    double threshold = 0.;

    for (unsigned int icell = 0; icell < eneInCell.size(); icell++) {
      if (theCalo->getCellFromIndex(icell)->getSubSys() != CaloCellId::Forward) {
        threshold = 3.;
      } else {
        threshold = 4.;
      }

      if (eneInCell[icell] > threshold) {
        if (theCalo->getCellFromIndex(icell)->getSubSys() == CaloCellId::Forward) {
          nCellOvTh++;
        }
        empz += eneInCell[icell] * (1. - std::cos(theCalo->getCellFromIndex(icell)->getThetaCell()));
      }
    }

    EmpzHFm->Fill(empz, weight);
    ntHFm->Fill(nCellOvTh, weight);
    eneHFmSel->Fill(eneHFm, weight);
  }

  // QCD-10-001
  double phiMax = hepmcGPCollection[iMax]->momentum().phi();
  std::vector<unsigned int> nchvsphi(nphiBin, 0);
  std::vector<double> sptvsphi(nphiBin, 0.);
  unsigned int nChaTra = 0;
  double sptTra = 0.;

  double binPhiW = 360. / nphiBin;
  if (sel2) {
    for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
      if (isCharged(i) && std::fabs(hepmcGPCollection[i]->momentum().eta()) < 2.) {
        double thePhi = (hepmcGPCollection[i]->momentum().phi() - phiMax) / CLHEP::degree;
        if (thePhi < -180.) {
          thePhi += 360.;
        } else if (thePhi > 180.) {
          thePhi -= 360.;
        }
        unsigned int thePhiBin = (int)((thePhi + 180.) / binPhiW);
        if (thePhiBin == nphiBin) {
          thePhiBin -= 1;
        }
        nchvsphi[thePhiBin]++;
        sptvsphi[thePhiBin] += hepmcGPCollection[i]->momentum().perp();
        // analysis in the transverse region
        if (std::fabs(thePhi) > 60. && std::fabs(thePhi) < 120.) {
          nChaTra++;
          sptTra += hepmcGPCollection[i]->momentum().perp();
          binW = dNchdpt2->getTH1()->GetBinWidth(1);
          dNchdpt2->Fill(hepmcGPCollection[i]->momentum().perp(), 1. / binW);  // weight to account for the pt bin width
          binW = dNchdeta2->getTH1()->GetBinWidth(1);
          // weight to account for the eta bin width
          dNchdeta2->Fill(hepmcGPCollection[i]->momentum().eta(), 1. / binW);
        }
      }
    }
    nCha->Fill(nChaTra, weight);
    binW = dNchdSpt->getTH1()->GetBinWidth(1);
    dNchdSpt->Fill(sptTra, 1.);
    //how do one apply weights to a profile? MonitorElement doesn't allow to
    nChaDenLpt->Fill(hepmcGPCollection[iMax]->momentum().perp(), nChaTra / 4. / CLHEP::twopi);
    sptDenLpt->Fill(hepmcGPCollection[iMax]->momentum().perp(), sptTra / 4. / CLHEP::twopi);
    for (unsigned int i = 0; i < nphiBin; i++) {
      double thisPhi = -180. + (i + 0.5) * binPhiW;
      dNchdphi->Fill(thisPhi, nchvsphi[i] / binPhiW / 4.);  // density in phi and eta
      dSptdphi->Fill(thisPhi, sptvsphi[i] / binPhiW / 4.);  // density in phi and eta
    }
  }

  // Gather information in the charged GenJet collection
  edm::Handle<reco::GenJetCollection> genChJets;
  iEvent.getByToken(genchjetCollectionToken_, genChJets);

  unsigned int nJets = 0;
  double pt1 = 0.;
  double pt2 = 0.;
  reco::GenJetCollection::const_iterator ij1 = genChJets->begin();
  reco::GenJetCollection::const_iterator ij2 = genChJets->begin();
  if (sel2) {
    for (reco::GenJetCollection::const_iterator iter = genChJets->begin(); iter != genChJets->end(); ++iter) {
      double eta = (*iter).eta();
      double pt = (*iter).pt();
      if (verbosity_ > 0) {
        std::cout << "GenJet " << std::setw(14) << std::fixed << (*iter).pt() << std::setw(14) << std::fixed
                  << (*iter).eta() << std::setw(14) << std::fixed << (*iter).phi() << std::endl;
      }
      if (std::fabs(eta) < 2.) {
        nJets++;
        binW = dNchjdeta->getTH1()->GetBinWidth(1);
        dNchjdeta->Fill(eta, 1. / binW);
        binW = dNchjdpt->getTH1()->GetBinWidth(1);
        dNchjdpt->Fill(pt, 1. / binW);
        if (pt >= pt1) {
          pt1 = pt;
          ij1 = iter;
        }
        if (pt < pt1 && pt >= pt2) {
          pt2 = pt;
          ij2 = iter;
        }
      }
    }

    nChj->Fill(nJets, weight);
    if (nJets > 0 && ij1 != genChJets->end()) {
      leadChjpt->Fill(pt1, weight);
      leadChjeta->Fill((*ij1).eta(), weight);
      if (nJets > 1 && ij2 != genChJets->end()) {
        pt1pt2optotch->Fill(pt1 + pt2, (pt1 + pt2) / ptot);
      }
    }
  }

  // Gather information in the GenJet collection
  edm::Handle<reco::GenJetCollection> genJets;
  iEvent.getByToken(genjetCollectionToken_, genJets);

  nJets = 0;
  pt1 = 0.;
  pt2 = 0.;
  double pt3 = 0.;

  // needed for Jet-Multiplicity Analysis
  int jm25njets = 0;
  double jm25HT = 0.;
  double jm25pt1 = 0.;
  double jm25pt2 = 0.;
  double jm25pt3 = 0.;
  double jm25pt4 = 0.;

  int jm80njets = 0;
  double jm80HT = 0.;
  double jm80pt1 = 0.;
  double jm80pt2 = 0.;
  double jm80pt3 = 0.;
  double jm80pt4 = 0.;

  reco::GenJetCollection::const_iterator ij3 = genJets->begin();
  if (sel6) {
    for (reco::GenJetCollection::const_iterator iter = genJets->begin(); iter != genJets->end(); ++iter) {
      double eta = (*iter).eta();
      double pt = (*iter).pt();
      if (verbosity_ > 0) {
        std::cout << "GenJet " << std::setw(14) << std::fixed << (*iter).pt() << std::setw(14) << std::fixed
                  << (*iter).eta() << std::setw(14) << std::fixed << (*iter).phi() << std::endl;
      }
      if (std::fabs(eta) < 5.) {
        nJets++;
        if (pt >= pt1) {
          pt1 = pt;
          ij1 = iter;
        }
        if (pt < pt1 && pt >= pt2) {
          pt2 = pt;
          ij2 = iter;
        }
        if (pt < pt2 && pt >= pt3) {
          pt3 = pt;
          ij3 = iter;
        }
      }

      // find variables for Jet-Multiplicity Analysis
      if (fabs(iter->eta()) < 3. && iter->pt() > 25.) {
        jm25njets++;
        jm25HT += iter->pt();
        if (iter->pt() > jm25pt1) {
          jm25pt4 = jm25pt3;
          jm25pt3 = jm25pt2;
          jm25pt2 = jm25pt1;
          jm25pt1 = iter->pt();
        } else if (iter->pt() > jm25pt2) {
          jm25pt4 = jm25pt3;
          jm25pt3 = jm25pt2;
          jm25pt2 = iter->pt();
        } else if (iter->pt() > jm25pt3) {
          jm25pt4 = jm25pt3;
          jm25pt3 = iter->pt();
        } else if (iter->pt() > jm25pt4) {
          jm25pt4 = iter->pt();
        }
        // even harder jets...
        if (iter->pt() > 80.) {
          jm80njets++;
          jm80HT += iter->pt();
          if (iter->pt() > jm80pt1) {
            jm80pt4 = jm80pt3;
            jm80pt3 = jm80pt2;
            jm80pt2 = jm80pt1;
            jm80pt1 = iter->pt();
          } else if (iter->pt() > jm80pt2) {
            jm80pt4 = jm80pt3;
            jm80pt3 = jm80pt2;
            jm80pt2 = iter->pt();
          } else if (iter->pt() > jm80pt3) {
            jm80pt4 = jm80pt3;
            jm80pt3 = iter->pt();
          } else if (iter->pt() > jm80pt4) {
            jm80pt4 = iter->pt();
          }
        }
      }
    }
    if (jm25njets > 3) {
      _JM25njets->Fill(jm25njets, weight);
      _JM25ht->Fill(jm25HT, weight);
      _JM25pt1->Fill(jm25pt1, weight);
      _JM25pt2->Fill(jm25pt2, weight);
      _JM25pt3->Fill(jm25pt3, weight);
      _JM25pt4->Fill(jm25pt4, weight);
    }
    if (jm80njets > 3) {
      _JM80njets->Fill(jm80njets, weight);
      _JM80ht->Fill(jm80HT, weight);
      _JM80pt1->Fill(jm80pt1, weight);
      _JM80pt2->Fill(jm80pt2, weight);
      _JM80pt3->Fill(jm80pt3, weight);
      _JM80pt4->Fill(jm80pt4, weight);
    }

    // select a di-jet event JME-10-001 variant
    double sumJetEt = 0;
    double sumPartPt = 0.;
    double sumChPartPt = 0.;
    double jpx = 0;
    double jpy = 0;
    if (nJets >= 2 && ij1 != genJets->end() && ij2 != genJets->end()) {
      if ((*ij1).pt() > 25. && (*ij1).pt() > 25.) {
        double deltaPhi = std::fabs((*ij1).phi() - (*ij2).phi()) / CLHEP::degree;
        if (deltaPhi > 180.)
          deltaPhi = 360. - deltaPhi;
        pt1pt2Dphi->Fill(deltaPhi, weight);
        if (std::fabs(deltaPhi) > 2.5 * CLHEP::degree) {
          nDijet->Fill(0.5, weight);

          for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
            double eta = hepmcGPCollection[i]->momentum().eta();
            unsigned int iBin = getHFbin(eta);
            if (!isNeutrino(i) && iBin < CaloCellManager::nForwardEta) {
              hfDJ[iBin] += hepmcGPCollection[i]->momentum().rho();
            }
            if (!isNeutrino(i) && std::fabs(eta) < 5.) {
              sumPartPt += hepmcGPCollection[i]->momentum().perp();
              if (isCharged(i)) {
                sumChPartPt += hepmcGPCollection[i]->momentum().perp();
              }
            }
          }
          for (unsigned int i = 0; i < CaloCellManager::nForwardEta; i++) {
            binW = theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i + 1] -
                   theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i];
            dEdetaHFdj->Fill(i + 0.5, hfDJ[i] / binW);
          }

          double invMass = (*ij1).energy() * (*ij2).energy() - (*ij1).px() * (*ij2).px() - (*ij1).py() * (*ij2).py() -
                           (*ij1).pz() * (*ij2).pz();
          invMass = std::sqrt(invMass);
          pt1pt2InvM->Fill(invMass, weight);

          sumPt->Fill(sumPartPt, weight);
          sumChPt->Fill(sumChPartPt, weight);

          unsigned int nSelJets = 0;
          for (reco::GenJetCollection::const_iterator iter = genJets->begin(); iter != genJets->end(); ++iter) {
            double pt = (*iter).pt();
            double eta = (*iter).eta();
            if (std::fabs(eta) < 5.) {
              nSelJets++;
              binW = dNjdeta->getTH1()->GetBinWidth(1);
              dNjdeta->Fill(eta, 1. / binW * weight);
              binW = dNjdpt->getTH1()->GetBinWidth(1);
              dNjdpt->Fill(pt, 1. / binW * weight);
              sumJetEt += (*iter).pt();
              jpx += (*iter).px();
              jpy += (*iter).py();
            }
          }

          nj->Fill(nSelJets, weight);
          double mEt = std::sqrt(jpx * jpx + jpy * jpy);
          sumJEt->Fill(sumJetEt, weight);
          missEtosumJEt->Fill(mEt / sumJetEt, weight);

          if (nSelJets >= 3) {
            pt3Frac->Fill((*ij3).pt() / (pt1 + pt2), weight);
          }

          pt1pt2optot->Fill(pt1 + pt2, (pt1 + pt2) / sumJetEt);
          pt1pt2balance->Fill((pt1 - pt2) / (pt1 + pt2), weight);
        }
      }
    }
  }

  //compute differential jet rates
  std::vector<const HepMC::GenParticle*> qcdActivity;
  HepMCValidationHelper::removeIsolatedLeptons(myGenEvent, 0.2, 3., qcdActivity);
  //HepMCValidationHelper::allStatus1(myGenEvent, qcdActivity);
  //fill PseudoJets to use fastjet
  std::vector<fastjet::PseudoJet> vecs;
  int counterUser = 1;
  std::vector<const HepMC::GenParticle*>::const_iterator iqcdact;
  for (iqcdact = qcdActivity.begin(); iqcdact != qcdActivity.end(); ++iqcdact) {
    const HepMC::FourVector& fmom = (*iqcdact)->momentum();
    fastjet::PseudoJet pseudoJet(fmom.px(), fmom.py(), fmom.pz(), fmom.e());
    pseudoJet.set_user_index(counterUser);
    vecs.push_back(pseudoJet);
    ++counterUser;
  }
  //compute jets
  fastjet::ClusterSequence cseq(vecs, fastjet::JetDefinition(fastjet::kt_algorithm, 1., fastjet::E_scheme));
  //access the cluster sequence and get the relevant info
  djr10->Fill(std::log10(sqrt(cseq.exclusive_dmerge(0))), weight);
  djr21->Fill(std::log10(sqrt(cseq.exclusive_dmerge(1))), weight);
  djr32->Fill(std::log10(sqrt(cseq.exclusive_dmerge(2))), weight);
  djr43->Fill(std::log10(sqrt(cseq.exclusive_dmerge(3))), weight);

  // compute sumEt for all stable particles
  std::vector<const HepMC::GenParticle*> allStable;
  HepMCValidationHelper::allStatus1(myGenEvent, allStable);

  double sumEt = 0.;
  double sumEt1 = 0.;
  double sumEt2 = 0.;
  double sumEt3 = 0.;
  double sumEt4 = 0.;
  double sumEt5 = 0.;

  for (std::vector<const HepMC::GenParticle*>::const_iterator iter = allStable.begin(); iter != allStable.end();
       ++iter) {
    double thisEta = fabs((*iter)->momentum().eta());

    if (thisEta < 5.) {
      const HepMC::FourVector mom = (*iter)->momentum();
      double px = mom.px();
      double py = mom.py();
      double pz = mom.pz();
      double E = mom.e();
      double thisSumEt = (sqrt(px * px + py * py) * E / sqrt(px * px + py * py + pz * pz));
      sumEt += thisSumEt;
      if (thisEta < 1.0)
        sumEt1 += thisSumEt;
      else if (thisEta < 2.0)
        sumEt2 += thisSumEt;
      else if (thisEta < 3.0)
        sumEt3 += thisSumEt;
      else if (thisEta < 4.0)
        sumEt4 += thisSumEt;
      else
        sumEt5 += thisSumEt;
    }
  }

  if (sumEt > 0.)
    _sumEt->Fill(sumEt, weight);
  if (sumEt1 > 0.)
    _sumEt1->Fill(sumEt1, weight);
  if (sumEt2 > 0.)
    _sumEt2->Fill(sumEt2, weight);
  if (sumEt3 > 0.)
    _sumEt3->Fill(sumEt3, weight);
  if (sumEt4 > 0.)
    _sumEt4->Fill(sumEt4, weight);
  if (sumEt5 > 0.)
    _sumEt5->Fill(sumEt5, weight);

  delete myGenEvent;
}  //analyze

bool MBUEandQCDValidation::isCharged(unsigned int i) {
  bool status = false;
  if (hepmcGPCollection.size() < i + 1) {
    return status;
  } else {
    status = (hepmcCharge[i] != 0. && hepmcCharge[i] != -999.);
  }
  return status;
}

bool MBUEandQCDValidation::isNeutral(unsigned int i) {
  bool status = false;
  int pdgId = std::abs(hepmcGPCollection[i]->pdg_id());
  if (hepmcGPCollection.size() < i + 1) {
    return status;
  } else {
    status = (hepmcCharge[i] == 0. && pdgId != 12 && pdgId != 14 && pdgId != 16);
  }
  return status;
}

bool MBUEandQCDValidation::isNeutrino(unsigned int i) {
  bool status = false;
  int pdgId = std::abs(hepmcGPCollection[i]->pdg_id());
  if (hepmcGPCollection.size() < i + 1) {
    return status;
  } else {
    status = (pdgId == 12 || pdgId == 14 || pdgId == 16);
  }
  return status;
}

unsigned int MBUEandQCDValidation::getHFbin(double eta) {
  unsigned int iBin = 999;

  std::vector<double> theEtaRanges(theCalo->getEtaRanges());

  for (unsigned int i = CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta;
       i < CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + CaloCellManager::nForwardEta;
       i++) {
    if (std::fabs(eta) >= theEtaRanges[i] && std::fabs(eta) < theEtaRanges[i + 1]) {
      iBin = i - CaloCellManager::nBarrelEta - CaloCellManager::nEndcapEta;
    }
  }

  return iBin;
}

const unsigned int MBUEandQCDValidation::nphiBin = 36;
const unsigned int MBUEandQCDValidation::initSize = 1000;