Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-08 05:12:09 

0001 /*Class MBUEandQCDValidation
0002  *  
0003  *  Class to fill dqm monitor elements from existing EDM file
0004  *
0005  */
0006 
0007 #include "MBUEandQCDValidation.h"
0008 #include "Validation/EventGenerator/interface/HepMCValidationHelper.h"
0009 
0010 #include "CLHEP/Units/defs.h"
0011 #include "CLHEP/Units/PhysicalConstants.h"
0012 #include "CLHEP/Units/SystemOfUnits.h"
0013 
0014 #include "fastjet/JetDefinition.hh"
0015 #include "fastjet/ClusterSequence.hh"
0016 #include "Validation/EventGenerator/interface/DQMHelper.h"
0017 using namespace edm;
0018 
0019 MBUEandQCDValidation::MBUEandQCDValidation(const edm::ParameterSet& iPSet)
0020     : wmanager_(iPSet, consumesCollector()),
0021       hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")),
0022       genchjetCollection_(iPSet.getParameter<edm::InputTag>("genChjetsCollection")),
0023       genjetCollection_(iPSet.getParameter<edm::InputTag>("genjetsCollection")),
0024       verbosity_(iPSet.getUntrackedParameter<unsigned int>("verbosity", 0)) {
0025   hepmcGPCollection.reserve(initSize);
0026   hepmcCharge.reserve(initSize);
0027 
0028   theCalo = new CaloCellManager(verbosity_);
0029 
0030   eneInCell.resize(CaloCellManager::nCaloCell);
0031 
0032   hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
0033   genjetCollectionToken_ = consumes<reco::GenJetCollection>(genjetCollection_);
0034   genchjetCollectionToken_ = consumes<reco::GenJetCollection>(genchjetCollection_);
0035   fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
0036 }
0037 
0038 MBUEandQCDValidation::~MBUEandQCDValidation() { delete theCalo; }
0039 
0040 void MBUEandQCDValidation::dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) {
0041   fPDGTable = c.getHandle(fPDGTableToken);
0042 }
0043 
0044 void MBUEandQCDValidation::bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) {
0045   ///Setting the DQM top directories
0046   DQMHelper dqm(&i);
0047   i.setCurrentFolder("Generator/MBUEandQCD");
0048 
0049   ///Booking the ME's
0050 
0051   // Number of analyzed events
0052   nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., " ", "Number of events");
0053 
0054   // Number of events with no forward trigger
0055   nNoFwdTrig = dqm.book1dHisto(
0056       "nNoFwdTrig", "n Events no forward trigger", 1, 0., 1., " ", "Number of Events with no Forward Trigger");
0057 
0058   // Number of Events with a single arm forward trigger
0059   nSaFwdTrig = dqm.book1dHisto("nSaFwdTrig",
0060                                "n Events single arm forward trigger",
0061                                1,
0062                                0.,
0063                                1.,
0064                                " ",
0065                                "Number of Events with Single Arm Forward Trigger");
0066 
0067   // Number of Events with b quark
0068   nbquark = dqm.book1dHisto("nbquark", "n Events with b quark", 1, 0., 1., " ", "Number of Events with b Quarks");
0069 
0070   // Number of Events with c and b quark
0071   ncandbquark = dqm.book1dHisto(
0072       "ncandbquark", "n Events with c and b quark", 1, 0., 1., "", "Number of Events with c and b Quark");
0073 
0074   // Number of Events with c and no b quark
0075   ncnobquark = dqm.book1dHisto(
0076       "ncnobquark", "n Events with c and no b quark", 1, 0., 1., "", "Number of Events with c and no b Quark");
0077 
0078   // Number of selected events for QCD-09-010
0079   nEvt1 = dqm.book1dHisto(
0080       "nEvt1", "n Events QCD-09-010", 1, 0., 1., "", "Number of Events passing the QCD-09-010 selection");
0081   // dNchdpt QCD-09-010
0082   dNchdpt1 = dqm.book1dHisto("dNchdpt1",
0083                              "dNchdpt QCD-09-010",
0084                              30,
0085                              0.,
0086                              6.,
0087                              "P_{t}^{charged tracks QCD-09-010 selection} (GeV)",
0088                              "Number of Charged Tracks");
0089   // dNchdeta QCD-09-010
0090   dNchdeta1 = dqm.book1dHisto("dNchdeta1",
0091                               "dNchdeta QCD-09-010",
0092                               10,
0093                               -2.5,
0094                               2.5,
0095                               "#eta^{charged tracks QCD-09-010 selection}",
0096                               "Number of Charged Tracks");
0097   // Number of selected events for QCD-10-001
0098 
0099   nEvt2 = dqm.book1dHisto(
0100       "nEvt2", "n Events QCD-10-001", 1, 0., 1., "", "Number of Events passing the QCD-10-001 selection");
0101   // Leading track pt QCD-10-001
0102   leadTrackpt = dqm.book1dHisto("leadTrackpt",
0103                                 "leading track pt QCD-10-001",
0104                                 200,
0105                                 0.,
0106                                 100.,
0107                                 "P_{t}^{lead track QCD-10-001 selection} (GeV)",
0108                                 "Number of Events");
0109   // Leading track eta QCD-10-001
0110   leadTracketa = dqm.book1dHisto("leadTracketa",
0111                                  "leading track eta QCD-10-001",
0112                                  50.,
0113                                  -2.5,
0114                                  2.5,
0115                                  "#eta^{lead track QCD-10-001 selection}",
0116                                  "Number of Events");
0117   // transverse charged particle density vs leading track pt
0118   nChaDenLpt = i.bookProfile("nChaDenLpt", "charged density vs leading pt", 200, 0., 100., 0., 100., " ");
0119   // transverse charged particle density vs leading track pt
0120   sptDenLpt = i.bookProfile("sptDenLpt", "sum pt density vs leading pt", 200, 0., 100., 0., 300., " ");
0121   // dNchdpt QCD-10-001 transverse
0122   dNchdpt2 = dqm.book1dHisto("dNchdpt2",
0123                              "dNchdpt QCD-10-001",
0124                              200,
0125                              0.,
0126                              100.,
0127                              "P_{t}^{charged tracks QCD-10-001 selection} (GeV)",
0128                              "Number of Charged Tracks");
0129   // dNchdeta QCD-10-001 transverse
0130   dNchdeta2 = dqm.book1dHisto("dNchdeta2",
0131                               "dNchdeta QCD-10-001",
0132                               50,
0133                               -2.5,
0134                               2.5,
0135                               "#eta^{charged tracks QCD-10-001 selection}",
0136                               "Number of Charged Tracks");
0137   // nCha QCD-10-001 transverse
0138   nCha = dqm.book1dHisto(
0139       "nCha", "n charged QCD-10-001", 100, 0., 100., "N^{charged tracks QCD-10-001 selection}", "Number of Events");
0140   // dNchdSpt transverse
0141   dNchdSpt = dqm.book1dHisto("dNchdSpt",
0142                              "dNchdSpt QCD-10-001",
0143                              300,
0144                              0.,
0145                              300.,
0146                              "P_{t}^{charged trackes in transverse region QCD-10-001selection} (GeV)",
0147                              "Number of Charged Tracks");
0148   // dNchdphi
0149   dNchdphi = i.bookProfile("dNchdphi", "dNchdphi QCD-10-001", nphiBin, -180., 180., 0., 30., " ");
0150   // dSptdphi
0151   dSptdphi = i.bookProfile("dSptdphi", "dSptdphi QCD-10-001", nphiBin, -180., 180., 0., 30., " ");
0152 
0153   // number of charged jets QCD-10-001
0154   nChj = dqm.book1dHisto(
0155       "nChj", "n charged jets QCD-10-001", 30, 0, 30., "N^{charged jets QCD-10-001 selection}", "Number of Events");
0156   // dNchjdeta QCD-10-001
0157   dNchjdeta = dqm.book1dHisto("dNchjdeta",
0158                               "dNchjdeta QCD-10-001",
0159                               50,
0160                               -2.5,
0161                               2.5,
0162                               "#eta^{charged jets QCD-10-001 selection}",
0163                               "Number of charged Jets");
0164   // dNchjdpt QCD-10-001
0165   dNchjdpt = dqm.book1dHisto("dNchjdpt",
0166                              "dNchjdpt QCD-10-001",
0167                              100,
0168                              0.,
0169                              100.,
0170                              "P_{t}^{charged jets QCD-10-001 selection}",
0171                              "Number of charged Jets");
0172   // leading charged jet pt QCD-10-001
0173   leadChjpt = dqm.book1dHisto("leadChjpt",
0174                               "leadChjpt QCD-10-001",
0175                               100,
0176                               0.,
0177                               100.,
0178                               "P_{t}^{lead charged jet QCD-10-001 selection}",
0179                               "Number of charged Jets");
0180   // leading charged jet eta QCD-10-001
0181   leadChjeta = dqm.book1dHisto("leadChjeta",
0182                                "leadChjeta QCD-10-001",
0183                                50,
0184                                -2.5,
0185                                2.5,
0186                                "#eta^{lead charged jet QCD-10-001 selection}",
0187                                "Number of charged Jets");
0188   // (pt1+pt2)/ptot
0189   pt1pt2optotch = i.bookProfile("pt1pt2optotch", "sum 2 leading jets over ptot", 50, 0., 100., 0., 1., " ");
0190 
0191   // particle rates in tracker acceptance
0192   nPPbar = dqm.book1dHisto(
0193       "nPPbar", "nPPbar QCD-10-001", 30, 0., 30., "N_{p/#bar{p}}^{QCD-10-001 selection}", "Number of p/#bar{p}");
0194   nKpm = dqm.book1dHisto(
0195       "nKpm", "nKpm QCD-10-001", 30, 0., 30., "N_{K^{#pm}}^{QCD-10-001 selection}", "Number of K^{#pm}");
0196   nK0s = dqm.book1dHisto("nK0s", "nK0s QCD-10-001", 30, 0., 30., "N_{K^{0}}^{QCD-10-001 selection}", "Number of K^{0}");
0197   nL0 = dqm.book1dHisto(
0198       "nL0", "nL0 QCD-10-001", 30, 0., 30., "N_{#Lambda^{0}}^{QCD-10-001 selection}", "Number of #Lambda^{0}");
0199   nXim = dqm.book1dHisto("nXim", "nXim QCD-10-001", 30, 0., 30., "N_{#Xi}^{QCD-10-001 selection}", "Number of #Xi");
0200   nOmega = dqm.book1dHisto(
0201       "nOmega", "nOmega QCD-10-001", 30, 0., 30., "N_{#Omega^{#pm}}^{QCD-10-001 selection}", "Number of #Omega^{#pm}");
0202 
0203   pPPbar = dqm.book1dHisto("pPPbar",
0204                            "Log10(pt) PPbar QCD-10-001",
0205                            25,
0206                            -2.,
0207                            3.,
0208                            "log_{10}(P_{t}^{p/#bar{p} QCD-10-001 selection}) (log_{10}(GeV))",
0209                            "Number of p/#bar{p}");
0210   pKpm = dqm.book1dHisto("pKpm",
0211                          "Log10(pt) Kpm QCD-10-001",
0212                          25,
0213                          -2.,
0214                          3.,
0215                          "log_{10}(P_{t}^{K^{#pm} QCD-10-001 selection}) (log_{10}(GeV))",
0216                          "Number of K^{#pm}");
0217   pK0s = dqm.book1dHisto("pK0s",
0218                          "Log10(pt) K0s QCD-10-001",
0219                          25,
0220                          -2.,
0221                          3.,
0222                          "log_{10}(P_{t}^{K^{0} QCD-10-001 selection}) (log_{10}(GeV))",
0223                          "Number of K^{0}");
0224   pL0 = dqm.book1dHisto("pL0",
0225                         "Log10(pt) L0 QCD-10-001",
0226                         25,
0227                         -2.,
0228                         3.,
0229                         "log_{10}(P_{t}^{#Lambda^{0} QCD-10-001 selection}) (log_{10}(GeV))",
0230                         "Number of #Lambda^{0}");
0231   pXim = dqm.book1dHisto("pXim",
0232                          "Log10(pt) Xim QCD-10-001",
0233                          25,
0234                          -2.,
0235                          3.,
0236                          "log_{10}(P_{t}^{#Xi^{#pm} QCD-10-001 selection}) (log_{10}(GeV))",
0237                          "Number of #Xi");
0238   pOmega = dqm.book1dHisto("pOmega",
0239                            "Log10(pt) Omega QCD-10-001",
0240                            25,
0241                            -2.,
0242                            3.,
0243                            "log_{10}(P_{t}^{#Omega^{#pm} QCD-10-001 selection}) (log_{10}(GeV))",
0244                            "Number of #Omega^{#pm}");
0245 
0246   // neutral rate in the barrel + HF acceptance
0247   nNNbar = dqm.book1dHisto(
0248       "nNNbar", "nNNbar QCD-10-001", 30, 0., 30., "N_{n/#bar{n}}^{QCD-10-001 selection}", "Number of Events");
0249   nGamma = dqm.book1dHisto(
0250       "nGamma", "nGamma QCD-10-001", 50, 0., 200., "N_{#gamma}^{QCD-10-001 selection}", "Number of Events");
0251 
0252   pNNbar = dqm.book1dHisto("pNNbar",
0253                            "Log10(pt) NNbar QCD-10-001",
0254                            25,
0255                            -2.,
0256                            3.,
0257                            "log_{10}(P_{t}^{n/#bar{n} QCD-10-001 selection}) (log_{10}(GeV))",
0258                            "Number of n/#bar{n}");
0259   pGamma = dqm.book1dHisto("pGamma",
0260                            "Log10(pt) Gamma QCD-10-001",
0261                            25,
0262                            -2.,
0263                            3.,
0264                            "log_{10}(P_{t}^{#gamma QCD-10-001 selection}) (log_{10}(GeV))",
0265                            "Number of #gamma");
0266 
0267   // highest pt electron spectrum
0268   elePt = dqm.book1dHisto("elePt",
0269                           "highest pt electron Log10(pt)",
0270                           30,
0271                           -2.,
0272                           4.,
0273                           "log_{10}(P_{t}^{highest e} (log_{10}(GeV))",
0274                           "Number of Events");
0275 
0276   // highest pt muon spectrum
0277   muoPt = dqm.book1dHisto("muoPt",
0278                           "highest pt muon Log10(pt)",
0279                           30,
0280                           -2.,
0281                           4.,
0282                           "log_{10}(P_{t}^{highest #mu}  (log_{10}(GeV))",
0283                           "Number of Events");
0284 
0285   // number of selected di-jet events
0286   nDijet = dqm.book1dHisto(
0287       "nDijet", "n Dijet Events", 1, 0., 1., " ", "Number of Events Passing Di-jet JME-10-001 Selection");
0288   // number of jets
0289   nj = dqm.book1dHisto("nj", "n jets ", 30, 0, 30., "N_{jets}^{JME-10-001}", "Number of Events");
0290   // dNjdeta
0291   dNjdeta = dqm.book1dHisto("dNjdeta", "dNjdeta ", 50, -5., 5., "#eta_{jets}^{JME-10-001 selection}", "Number of Jets");
0292   // dNjdpt
0293   dNjdpt = dqm.book1dHisto("dNjdpt", "dNjdpt ", 60, 0., 300., "P_{t}^{jets JME-10-001 selection}", "Number of Jets");
0294   // (pt1+pt2)/ptot
0295   pt1pt2optot = i.bookProfile("pt1pt2optot", "sum 2 leading jets over Et tot ", 60, 0., 300., 0., 1., " ");
0296   // pt1-pt2
0297   pt1pt2balance =
0298       dqm.book1dHisto("pt1pt2balance",
0299                       "2 leading jets pt difference ",
0300                       10,
0301                       0.,
0302                       1.,
0303                       "#frac{P_{t}^{1st jet}-P_{t}^{2nd jet}}{P_{t}^{1st jet}+P_{t}^{2nd jet}}^{JME-10-001 selection}",
0304                       "Number of Di-jet Events");
0305   // pt1 pt2 Delta phi
0306   pt1pt2Dphi = dqm.book1dHisto("pt1pt2Dphi",
0307                                "pt1 pt2 delta phi ",
0308                                nphiBin,
0309                                0.,
0310                                180.,
0311                                "#Delta#phi(jet^{1st},jet^{2nd})^{JME-10-001 selection} (rad)",
0312                                "Number of Di-jet Events");
0313   // pt1 pt2 invariant mass
0314   pt1pt2InvM = dqm.book1dHisto("pt1pt2InvM",
0315                                "pt1 pt2 invariant mass ",
0316                                60,
0317                                0.,
0318                                600.,
0319                                "M_{di-jet}^{JME-10-001 selection}",
0320                                "Number of di-jet events");
0321   // pt3 fraction
0322   pt3Frac = dqm.book1dHisto("pt3Frac",
0323                             "2 pt3 over pt1+pt2 ",
0324                             30,
0325                             0.,
0326                             1.,
0327                             "#frac{P_{t}^{3rd jet}}{P_{t}^{1st jet}+P_{t}^{2nd jet}}^{JME-10-001 selection}",
0328                             "Number of 3rd Jets");
0329   // sum of jets Et
0330   sumJEt = dqm.book1dHisto(
0331       "sumJEt", "sum Jet Et ", 60, 0., 300., "#Sigma E_{t}^{jets JME-10-001 selection}", "Number of di-jet events");
0332   // fraction of missing Et over sum of jets Et
0333   missEtosumJEt = dqm.book1dHisto("missEtosumJEt",
0334                                   "missing Et over sumJet Et ",
0335                                   30,
0336                                   0.,
0337                                   1.,
0338                                   "E_{t}^{miss}/#Sigma E_{t}^{jets JME-10-001 selection}",
0339                                   "Number of Di-jet Events");
0340   // sum of final state particle Pt
0341   sumPt = dqm.book1dHisto("sumPt",
0342                           "sum particle Pt ",
0343                           60,
0344                           0.,
0345                           600.,
0346                           "#Sigma P_{t}^{particles passing JME-10-001 selection}",
0347                           "Number of jets");
0348   // sum of final state charged particle Pt
0349   sumChPt = dqm.book1dHisto("sumChPt",
0350                             "sum charged particle Pt ",
0351                             60,
0352                             0.,
0353                             300.,
0354                             "#Sigma P_{t}^{charged particles passing JME-10-001 selection}",
0355                             "Number of Jets");
0356 
0357   //Number of selected events for the HF energy flux analysis
0358   nHFflow = dqm.book1dHisto("nHFflow",
0359                             "n HF flow events",
0360                             1,
0361                             0.,
0362                             1.,
0363                             " ",
0364                             "Number of Events passing JME-10-001, FWD-10-002 and Jet-Multiplicity selection");
0365   //Forward energy flow for MinBias BSC selection
0366   dEdetaHFmb = i.bookProfile("dEdetaHFmb",
0367                              "dEdeta HF MinBias",
0368                              (int)CaloCellManager::nForwardEta,
0369                              0,
0370                              (double)CaloCellManager::nForwardEta,
0371                              0.,
0372                              300.,
0373                              " ");
0374   //Forward energy flow for QCD dijet selection
0375   dEdetaHFdj = i.bookProfile("dEdetaHFdj",
0376                              "dEdeta HF QCD dijet",
0377                              (int)CaloCellManager::nForwardEta,
0378                              0,
0379                              (double)CaloCellManager::nForwardEta,
0380                              0.,
0381                              300.,
0382                              " ");
0383 
0384   // FWD-10-001 like diffraction analysis
0385   nHFSD = dqm.book1dHisto(
0386       "nHFSD", "n single diffraction in HF", 1, 0., 1., " ", "Number of single diffraction in HF FWD-10-001 selection");
0387   // E-pz HF-
0388   EmpzHFm = dqm.book1dHisto("EmpzHFm",
0389                             "E-pz HF- SD",
0390                             40,
0391                             0.,
0392                             200.,
0393                             "#Sigma E_{cal. cells/corrected for the longitudinal mometum}^{FWD-10-001 selection}",
0394                             "Number of Events");
0395   // Number of cells above threshold
0396   ntHFm = dqm.book1dHisto("ntHFm",
0397                           "number of HF- tower SD",
0398                           20,
0399                           0.,
0400                           20.,
0401                           " N_{cells over threshold for single diffraction in HF Towers}^{FWD-10-001 selection}",
0402                           "Number of Events");
0403   // Energy in HF-
0404   eneHFmSel = dqm.book1dHisto(
0405       "eneHFmSel", "energy in HF-", 40, 0., 200., "#Sigma E_{cal. cells}^{FWD-10-001 selection}", "Number of Events");
0406 
0407   // number of jets accepted in the 'Jet-Multiplicity' analysis
0408   _JM25njets = dqm.book1dHisto("JM25njets", "n jets", 15, 0, 15., "Number of JM25 Jets", "Number of Events");
0409   _JM25ht = dqm.book1dHisto("JM25ht", "HT", 80, 0, 800., "H_{t}^{JM25} (GeV)", "Number of Events");
0410   _JM25pt1 = dqm.book1dHisto("JM25pt1", "pt", 40, 0, 200., "P_{t}^{JM25,1st Jet} (GeV)", "Number of JM25 Jets");
0411   _JM25pt2 = dqm.book1dHisto("JM25pt2", "pt", 40, 0, 200., "P_{t}^{JM25,2nd Jet} (GeV)", "Number of JM25 Jets");
0412   _JM25pt3 = dqm.book1dHisto("JM25pt3", "pt", 40, 0, 200., "P_{t}^{JM25,3rd Jet} (GeV)", "Number of JM25 Jets");
0413   _JM25pt4 = dqm.book1dHisto("JM25pt4", "pt", 40, 0, 200., "P_{t}^{JM25,4th Jet} (GeV)", "Number of JM25 Jets");
0414 
0415   _JM80njets = dqm.book1dHisto("JM80njets", "n jets", 15, 0, 15., "Number of JM80 Jets", "Number of Events");
0416   _JM80ht = dqm.book1dHisto("JM80ht", "HT", 80, 300, 1100., "H_{t}^{JM80} (GeV", "Number of Events");
0417   _JM80pt1 = dqm.book1dHisto("JM80pt1", "pt", 40, 60, 260., "P_{t}^{JM80,1st Jet} (GeV)", "Number of JM80 Jets");
0418   _JM80pt2 = dqm.book1dHisto("JM80pt2", "pt", 40, 60, 260., "P_{t}^{JM80,2nd Jet} (GeV)", "Number of JM80 Jets");
0419   _JM80pt3 = dqm.book1dHisto("JM80pt3", "pt", 40, 60, 260., "P_{t}^{JM80,3rd Jet} (GeV)", "Number of JM80 Jets");
0420   _JM80pt4 = dqm.book1dHisto("JM80pt4", "pt", 40, 60, 260., "P_{t}^{JM80,4th Jet} (GeV)", "Number of JM80 Jets");
0421 
0422   // differential jet rates
0423   djr10 = dqm.book1dHisto("djr10",
0424                           "Differential Jet Rate 1#rightarrow0",
0425                           60,
0426                           -1.,
0427                           5.,
0428                           "log_{10}(d_{min}(n,n+1)) for n=0",
0429                           "Number of Events");
0430   djr21 = dqm.book1dHisto("djr21",
0431                           "Differential Jet Rate 2#rightarrow1",
0432                           60,
0433                           -1.,
0434                           5.,
0435                           "log_{10}(d_{min}(n,n+1)) for n=1",
0436                           "Number of Events");
0437   djr32 = dqm.book1dHisto("djr32",
0438                           "Differential Jet Rate 3#rightarrow2",
0439                           60,
0440                           -1.,
0441                           5.,
0442                           "log_{10}(d_{min}(n,n+1)) for n=2",
0443                           "Number of Events");
0444   djr43 = dqm.book1dHisto("djr43",
0445                           "Differential Jet Rate 4#rightarrow3",
0446                           60,
0447                           -1.,
0448                           5.,
0449                           "log_{10}(d_{min}(n,n+1)) for n=3",
0450                           "Number of Events");
0451 
0452   // sumET analysis
0453   _sumEt = dqm.book1dHisto(
0454       "sumET", "Sum of stable particles Et", 150, 0, 600., "#Sigma E_{t}^{stable particles}", "Number of Events");
0455   _sumEt1 = dqm.book1dHisto("sumET1",
0456                             "Sum of stable particles Et (eta<0.5)",
0457                             150,
0458                             0,
0459                             200.,
0460                             "#Sigma E_{t}^{stable particles (#eta<0.5)}",
0461                             "Number of Events");
0462   _sumEt2 = dqm.book1dHisto("sumET2",
0463                             "Sum of stable particles Et (0.5<eta<1.0)",
0464                             150,
0465                             0,
0466                             200.,
0467                             "#Sigma E_{t}^{stable particles (0.5<#eta<1.0)}",
0468                             "Number of Events");
0469   _sumEt3 = dqm.book1dHisto("sumET3",
0470                             "Sum of stable particles Et (1.0<eta<1.5)",
0471                             150,
0472                             0,
0473                             200.,
0474                             "#Sigma E_{t}^{stable particles (1.0<#eta<1.5)}",
0475                             "Number of Events");
0476   _sumEt4 = dqm.book1dHisto("sumET4",
0477                             "Sum of stable particles Et (1.5<eta<2.0)",
0478                             150,
0479                             0,
0480                             200.,
0481                             "#Sigma E_{t}^{stable particles (1.5<#eta<2.0)}",
0482                             "Number of Events");
0483   _sumEt5 = dqm.book1dHisto("sumET5",
0484                             "Sum of stable particles Et (2.0<eta<5.0)",
0485                             150,
0486                             0,
0487                             200.,
0488                             "#Sigma E_{t}^{stable particles (2.0<#eta<5.0)}",
0489                             "Number of Events");
0490 
0491   return;
0492 }
0493 
0494 void MBUEandQCDValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0495   ///Gathering the HepMCProduct information
0496   edm::Handle<HepMCProduct> evt;
0497   iEvent.getByToken(hepmcCollectionToken_, evt);
0498 
0499   //Get HepMC EVENT
0500   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0501 
0502   double weight = wmanager_.weight(iEvent);
0503 
0504   if (verbosity_ > 0) {
0505     myGenEvent->print();
0506   }
0507 
0508   double binW = 1.;
0509 
0510   hepmcGPCollection.clear();
0511   hepmcCharge.clear();
0512   for (unsigned int i = 0; i < eneInCell.size(); i++) {
0513     eneInCell[i] = 0.;
0514   }
0515 
0516   nEvt->Fill(0.5, weight);
0517 
0518   //Looping through HepMC::GenParticle collection to search for status 1 particles
0519   double charge = 0.;
0520   unsigned int nb = 0;
0521   unsigned int nc = 0;
0522   for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin();
0523        iter != myGenEvent->particles_end();
0524        ++iter) {
0525     if (std::fabs((*iter)->pdg_id()) == 4) {
0526       nc++;
0527     }
0528     if (std::fabs((*iter)->pdg_id()) == 5) {
0529       nb++;
0530     }
0531     if ((*iter)->status() == 1) {
0532       hepmcGPCollection.push_back(*iter);
0533       const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID((*iter)->pdg_id()));
0534       if (PData == nullptr) {
0535         charge = -999.;
0536       } else
0537         charge = PData->charge();
0538 
0539       hepmcCharge.push_back(charge);
0540 
0541       if (verbosity_ > 0) {
0542         std::cout << "HepMC " << std::setw(14) << std::fixed << (*iter)->barcode() << std::setw(14) << std::fixed
0543                   << (*iter)->pdg_id() << std::setw(14) << std::fixed << (*iter)->momentum().perp() << std::setw(14)
0544                   << std::fixed << (*iter)->momentum().eta() << std::setw(14) << std::fixed << (*iter)->momentum().phi()
0545                   << std::endl;
0546       }
0547     }
0548   }
0549 
0550   int nBSCp = 0;
0551   int nBSCm = 0;
0552   double eneHFp = 0.;
0553   double eneHFm = 0.;
0554   int nChapt05 = 0;
0555   int nChaVtx = 0;
0556   for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
0557     if (!isNeutrino(i)) {
0558       // BSC trigger
0559 
0560       if (hepmcGPCollection[i]->momentum().eta() > 3.23 && hepmcGPCollection[i]->momentum().eta() < 4.65) {
0561         nBSCp++;
0562       }
0563       if (hepmcGPCollection[i]->momentum().eta() < -3.23 && hepmcGPCollection[i]->momentum().eta() > -4.65) {
0564         nBSCm++;
0565       }
0566 
0567       // number of charged particles in different selections
0568 
0569       if (std::fabs(hepmcGPCollection[i]->momentum().eta()) < 2.5 && hepmcGPCollection[i]->momentum().perp() > 0.5 &&
0570           isCharged(i)) {
0571         nChapt05++;
0572       }
0573       if (std::fabs(hepmcGPCollection[i]->momentum().eta()) < 2.5 && hepmcGPCollection[i]->momentum().perp() > 0.1 &&
0574           isCharged(i)) {
0575         nChaVtx++;
0576       }
0577       unsigned int theIndex = theCalo->getCellIndexFromAngle(hepmcGPCollection[i]->momentum().eta(),
0578                                                              hepmcGPCollection[i]->momentum().phi());
0579       if (theIndex < CaloCellManager::nCaloCell)
0580         eneInCell[theIndex] += hepmcGPCollection[i]->momentum().rho();
0581     }
0582   }
0583 
0584   // Forward calorimeters energy
0585 
0586   for (unsigned int icell = CaloCellManager::nBarrelCell + CaloCellManager::nEndcapCell;
0587        icell < CaloCellManager::nCaloCell;
0588        icell++) {
0589     if (theCalo->getCellFromIndex(icell)->getEtaMin() < 0.) {
0590       eneHFm += eneInCell[icell];
0591     } else {
0592       eneHFp += eneInCell[icell];
0593     }
0594   }
0595 
0596   // QCD-09-010 selection
0597   bool sel1 = false;
0598   if ((nBSCp > 0 || nBSCm > 0) && eneHFp >= 3. && eneHFm >= 3.) {
0599     sel1 = true;
0600   }
0601 
0602   // QCD-10-001 selection
0603   bool sel2 = false;
0604   if ((nBSCp > 0 || nBSCm > 0) && nChaVtx >= 3 && nChapt05 > 1) {
0605     sel2 = true;
0606   }
0607 
0608   // no forward trigger selection
0609   bool sel3 = false;
0610   if (nBSCp == 0 && nBSCm == 0) {
0611     sel3 = true;
0612   }
0613 
0614   // single arm forward trigger selection
0615   bool sel4 = false;
0616   if ((nBSCp > 0 && nBSCm == 0) || (nBSCm > 0 && nBSCp == 0)) {
0617     sel4 = true;
0618   }
0619 
0620   // BSC selection
0621   bool sel5 = false;
0622   if (nBSCp > 0 && nBSCm > 0) {
0623     sel5 = true;
0624   }
0625 
0626   // basic JME-10-001, FWD-10-002 and Jet-Multiplicity selection
0627   bool sel6 = false;
0628   if (sel5 && nChaVtx > 3) {
0629     sel6 = true;
0630   }
0631 
0632   // FWD-10-001 selection
0633   bool sel7 = false;
0634   if (nChaVtx >= 3 && nBSCm > 0 && eneHFp < 8.) {
0635     sel7 = true;
0636   }
0637 
0638   // Fill selection histograms
0639   if (sel1)
0640     nEvt1->Fill(0.5, weight);
0641   if (sel2)
0642     nEvt2->Fill(0.5, weight);
0643   if (sel3)
0644     nNoFwdTrig->Fill(0.5, weight);
0645   if (sel4)
0646     nSaFwdTrig->Fill(0.5, weight);
0647   if (sel6)
0648     nHFflow->Fill(0.5, weight);
0649   if (sel7)
0650     nHFSD->Fill(0.5, weight);
0651 
0652   if (nb > 0)
0653     nbquark->Fill(0.5, weight);
0654   if (nb > 0 && nc > 0)
0655     ncandbquark->Fill(0.5, weight);
0656   if (nb == 0 && nc > 0)
0657     ncnobquark->Fill(0.5, weight);
0658 
0659   // track analyses
0660   double ptMax = 0.;
0661   unsigned int iMax = 0;
0662   double ptot = 0.;
0663   unsigned int ppbar = 0;
0664   unsigned int nnbar = 0;
0665   unsigned int kpm = 0;
0666   unsigned int k0s = 0;
0667   unsigned int l0 = 0;
0668   unsigned int gamma = 0;
0669   unsigned int xim = 0;
0670   unsigned int omega = 0;
0671   unsigned int ele = 0;
0672   unsigned int muo = 0;
0673   unsigned int eleMax = 0;
0674   unsigned int muoMax = 0;
0675 
0676   std::vector<double> hfMB(CaloCellManager::nForwardEta, 0);
0677   std::vector<double> hfDJ(CaloCellManager::nForwardEta, 0);
0678 
0679   for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
0680     double eta = hepmcGPCollection[i]->momentum().eta();
0681     double pt = hepmcGPCollection[i]->momentum().perp();
0682     int pdgId = hepmcGPCollection[i]->pdg_id();
0683     if (isCharged(i) && std::fabs(eta) < 2.5) {
0684       if (sel1) {
0685         // QCD-09-010
0686         binW = dNchdpt1->getTH1()->GetBinWidth(1);
0687         dNchdpt1->Fill(pt, 1. / binW);  // weight to account for the pt bin width
0688         binW = dNchdeta1->getTH1()->GetBinWidth(1);
0689         dNchdeta1->Fill(eta, 1. / binW);  // weight to account for the eta bin width
0690       }
0691       // search for the leading track QCD-10-001
0692       if (sel2) {
0693         if (pt > ptMax) {
0694           ptMax = pt;
0695           iMax = i;
0696         }
0697         ptot += pt;
0698 
0699         // identified charged particle
0700         if (std::abs(pdgId) == 2212) {
0701           ppbar++;
0702           pPPbar->Fill(std::log10(pt), weight);
0703         } else if (std::abs(pdgId) == 321) {
0704           kpm++;
0705           pKpm->Fill(std::log10(pt), weight);
0706         } else if (std::abs(pdgId) == 3312) {
0707           xim++;
0708           pXim->Fill(std::log10(pt), weight);
0709         } else if (std::abs(pdgId) == 3334) {
0710           omega++;
0711           pOmega->Fill(std::log10(pt), weight);
0712         } else if (std::abs(pdgId) == 11) {
0713           ele++;
0714           eleMax = i;
0715         } else if (std::abs(pdgId) == 13) {
0716           muo++;
0717           muoMax = i;
0718         }
0719       }
0720     } else if (sel2 && isNeutral(i) && std::fabs(eta) < 2.5) {
0721       if (std::abs(pdgId) == 310) {
0722         k0s++;
0723         pK0s->Fill(std::log10(pt), weight);
0724       } else if (std::abs(pdgId) == 3122) {
0725         l0++;
0726         pL0->Fill(std::log10(pt), weight);
0727       }
0728     } else if (sel2 && isNeutral(i) && std::fabs(eta) < 5.19) {
0729       if (std::abs(pdgId) == 2112) {
0730         nnbar++;
0731         pNNbar->Fill(std::log10(pt), weight);
0732       } else if (std::abs(pdgId) == 22) {
0733         gamma++;
0734         pGamma->Fill(std::log10(pt), weight);
0735       }
0736     }
0737     unsigned int iBin = getHFbin(eta);
0738     if (sel6 && !isNeutrino(i) && iBin < CaloCellManager::nForwardEta) {
0739       hfMB[iBin] += hepmcGPCollection[i]->momentum().rho();
0740     }
0741   }
0742   nPPbar->Fill(ppbar, weight);
0743   nNNbar->Fill(nnbar, weight);
0744   nKpm->Fill(kpm, weight);
0745   nK0s->Fill(k0s, weight);
0746   nL0->Fill(l0, weight);
0747   nXim->Fill(xim, weight);
0748   nOmega->Fill(omega, weight);
0749   nGamma->Fill(gamma, weight);
0750 
0751   if (ele > 0)
0752     elePt->Fill(std::log10(hepmcGPCollection[eleMax]->momentum().perp()), weight);
0753   if (muo > 0)
0754     muoPt->Fill(std::log10(hepmcGPCollection[muoMax]->momentum().perp()), weight);
0755 
0756   leadTrackpt->Fill(hepmcGPCollection[iMax]->momentum().perp(), weight);
0757   leadTracketa->Fill(hepmcGPCollection[iMax]->momentum().eta(), weight);
0758 
0759   std::vector<double> theEtaRanges(theCalo->getEtaRanges());
0760 
0761   for (unsigned int i = 0; i < CaloCellManager::nForwardEta; i++) {
0762     binW = theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i + 1] -
0763            theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i];
0764     dEdetaHFmb->Fill(i + 0.5, hfMB[i] / binW);
0765   }
0766 
0767   // FWD-10-001
0768 
0769   if (sel7) {
0770     double empz = 0.;
0771     unsigned int nCellOvTh = 0;
0772     double threshold = 0.;
0773 
0774     for (unsigned int icell = 0; icell < eneInCell.size(); icell++) {
0775       if (theCalo->getCellFromIndex(icell)->getSubSys() != CaloCellId::Forward) {
0776         threshold = 3.;
0777       } else {
0778         threshold = 4.;
0779       }
0780 
0781       if (eneInCell[icell] > threshold) {
0782         if (theCalo->getCellFromIndex(icell)->getSubSys() == CaloCellId::Forward) {
0783           nCellOvTh++;
0784         }
0785         empz += eneInCell[icell] * (1. - std::cos(theCalo->getCellFromIndex(icell)->getThetaCell()));
0786       }
0787     }
0788 
0789     EmpzHFm->Fill(empz, weight);
0790     ntHFm->Fill(nCellOvTh, weight);
0791     eneHFmSel->Fill(eneHFm, weight);
0792   }
0793 
0794   // QCD-10-001
0795   double phiMax = hepmcGPCollection[iMax]->momentum().phi();
0796   std::vector<unsigned int> nchvsphi(nphiBin, 0);
0797   std::vector<double> sptvsphi(nphiBin, 0.);
0798   unsigned int nChaTra = 0;
0799   double sptTra = 0.;
0800 
0801   double binPhiW = 360. / nphiBin;
0802   if (sel2) {
0803     for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
0804       if (isCharged(i) && std::fabs(hepmcGPCollection[i]->momentum().eta()) < 2.) {
0805         double thePhi = (hepmcGPCollection[i]->momentum().phi() - phiMax) / CLHEP::degree;
0806         if (thePhi < -180.) {
0807           thePhi += 360.;
0808         } else if (thePhi > 180.) {
0809           thePhi -= 360.;
0810         }
0811         unsigned int thePhiBin = (int)((thePhi + 180.) / binPhiW);
0812         if (thePhiBin == nphiBin) {
0813           thePhiBin -= 1;
0814         }
0815         nchvsphi[thePhiBin]++;
0816         sptvsphi[thePhiBin] += hepmcGPCollection[i]->momentum().perp();
0817         // analysis in the transverse region
0818         if (std::fabs(thePhi) > 60. && std::fabs(thePhi) < 120.) {
0819           nChaTra++;
0820           sptTra += hepmcGPCollection[i]->momentum().perp();
0821           binW = dNchdpt2->getTH1()->GetBinWidth(1);
0822           dNchdpt2->Fill(hepmcGPCollection[i]->momentum().perp(), 1. / binW);  // weight to account for the pt bin width
0823           binW = dNchdeta2->getTH1()->GetBinWidth(1);
0824           // weight to account for the eta bin width
0825           dNchdeta2->Fill(hepmcGPCollection[i]->momentum().eta(), 1. / binW);
0826         }
0827       }
0828     }
0829     nCha->Fill(nChaTra, weight);
0830     dNchdSpt->Fill(sptTra, 1.);
0831     //how do one apply weights to a profile? MonitorElement doesn't allow to
0832     nChaDenLpt->Fill(hepmcGPCollection[iMax]->momentum().perp(), nChaTra / 4. / CLHEP::twopi);
0833     sptDenLpt->Fill(hepmcGPCollection[iMax]->momentum().perp(), sptTra / 4. / CLHEP::twopi);
0834     for (unsigned int i = 0; i < nphiBin; i++) {
0835       double thisPhi = -180. + (i + 0.5) * binPhiW;
0836       dNchdphi->Fill(thisPhi, nchvsphi[i] / binPhiW / 4.);  // density in phi and eta
0837       dSptdphi->Fill(thisPhi, sptvsphi[i] / binPhiW / 4.);  // density in phi and eta
0838     }
0839   }
0840 
0841   // Gather information in the charged GenJet collection
0842   edm::Handle<reco::GenJetCollection> genChJets;
0843   iEvent.getByToken(genchjetCollectionToken_, genChJets);
0844 
0845   unsigned int nJets = 0;
0846   double pt1 = 0.;
0847   double pt2 = 0.;
0848   reco::GenJetCollection::const_iterator ij1 = genChJets->begin();
0849   reco::GenJetCollection::const_iterator ij2 = genChJets->begin();
0850   if (sel2) {
0851     for (reco::GenJetCollection::const_iterator iter = genChJets->begin(); iter != genChJets->end(); ++iter) {
0852       double eta = (*iter).eta();
0853       double pt = (*iter).pt();
0854       if (verbosity_ > 0) {
0855         std::cout << "GenJet " << std::setw(14) << std::fixed << (*iter).pt() << std::setw(14) << std::fixed
0856                   << (*iter).eta() << std::setw(14) << std::fixed << (*iter).phi() << std::endl;
0857       }
0858       if (std::fabs(eta) < 2.) {
0859         nJets++;
0860         binW = dNchjdeta->getTH1()->GetBinWidth(1);
0861         dNchjdeta->Fill(eta, 1. / binW);
0862         binW = dNchjdpt->getTH1()->GetBinWidth(1);
0863         dNchjdpt->Fill(pt, 1. / binW);
0864         if (pt >= pt1) {
0865           pt1 = pt;
0866           ij1 = iter;
0867         }
0868         if (pt < pt1 && pt >= pt2) {
0869           pt2 = pt;
0870           ij2 = iter;
0871         }
0872       }
0873     }
0874 
0875     nChj->Fill(nJets, weight);
0876     if (nJets > 0 && ij1 != genChJets->end()) {
0877       leadChjpt->Fill(pt1, weight);
0878       leadChjeta->Fill((*ij1).eta(), weight);
0879       if (nJets > 1 && ij2 != genChJets->end()) {
0880         pt1pt2optotch->Fill(pt1 + pt2, (pt1 + pt2) / ptot);
0881       }
0882     }
0883   }
0884 
0885   // Gather information in the GenJet collection
0886   edm::Handle<reco::GenJetCollection> genJets;
0887   iEvent.getByToken(genjetCollectionToken_, genJets);
0888 
0889   nJets = 0;
0890   pt1 = 0.;
0891   pt2 = 0.;
0892   double pt3 = 0.;
0893 
0894   // needed for Jet-Multiplicity Analysis
0895   int jm25njets = 0;
0896   double jm25HT = 0.;
0897   double jm25pt1 = 0.;
0898   double jm25pt2 = 0.;
0899   double jm25pt3 = 0.;
0900   double jm25pt4 = 0.;
0901 
0902   int jm80njets = 0;
0903   double jm80HT = 0.;
0904   double jm80pt1 = 0.;
0905   double jm80pt2 = 0.;
0906   double jm80pt3 = 0.;
0907   double jm80pt4 = 0.;
0908 
0909   reco::GenJetCollection::const_iterator ij3 = genJets->begin();
0910   if (sel6) {
0911     for (reco::GenJetCollection::const_iterator iter = genJets->begin(); iter != genJets->end(); ++iter) {
0912       double eta = (*iter).eta();
0913       double pt = (*iter).pt();
0914       if (verbosity_ > 0) {
0915         std::cout << "GenJet " << std::setw(14) << std::fixed << (*iter).pt() << std::setw(14) << std::fixed
0916                   << (*iter).eta() << std::setw(14) << std::fixed << (*iter).phi() << std::endl;
0917       }
0918       if (std::fabs(eta) < 5.) {
0919         nJets++;
0920         if (pt >= pt1) {
0921           pt1 = pt;
0922           ij1 = iter;
0923         }
0924         if (pt < pt1 && pt >= pt2) {
0925           pt2 = pt;
0926           ij2 = iter;
0927         }
0928         if (pt < pt2 && pt >= pt3) {
0929           pt3 = pt;
0930           ij3 = iter;
0931         }
0932       }
0933 
0934       // find variables for Jet-Multiplicity Analysis
0935       if (fabs(iter->eta()) < 3. && iter->pt() > 25.) {
0936         jm25njets++;
0937         jm25HT += iter->pt();
0938         if (iter->pt() > jm25pt1) {
0939           jm25pt4 = jm25pt3;
0940           jm25pt3 = jm25pt2;
0941           jm25pt2 = jm25pt1;
0942           jm25pt1 = iter->pt();
0943         } else if (iter->pt() > jm25pt2) {
0944           jm25pt4 = jm25pt3;
0945           jm25pt3 = jm25pt2;
0946           jm25pt2 = iter->pt();
0947         } else if (iter->pt() > jm25pt3) {
0948           jm25pt4 = jm25pt3;
0949           jm25pt3 = iter->pt();
0950         } else if (iter->pt() > jm25pt4) {
0951           jm25pt4 = iter->pt();
0952         }
0953         // even harder jets...
0954         if (iter->pt() > 80.) {
0955           jm80njets++;
0956           jm80HT += iter->pt();
0957           if (iter->pt() > jm80pt1) {
0958             jm80pt4 = jm80pt3;
0959             jm80pt3 = jm80pt2;
0960             jm80pt2 = jm80pt1;
0961             jm80pt1 = iter->pt();
0962           } else if (iter->pt() > jm80pt2) {
0963             jm80pt4 = jm80pt3;
0964             jm80pt3 = jm80pt2;
0965             jm80pt2 = iter->pt();
0966           } else if (iter->pt() > jm80pt3) {
0967             jm80pt4 = jm80pt3;
0968             jm80pt3 = iter->pt();
0969           } else if (iter->pt() > jm80pt4) {
0970             jm80pt4 = iter->pt();
0971           }
0972         }
0973       }
0974     }
0975     if (jm25njets > 3) {
0976       _JM25njets->Fill(jm25njets, weight);
0977       _JM25ht->Fill(jm25HT, weight);
0978       _JM25pt1->Fill(jm25pt1, weight);
0979       _JM25pt2->Fill(jm25pt2, weight);
0980       _JM25pt3->Fill(jm25pt3, weight);
0981       _JM25pt4->Fill(jm25pt4, weight);
0982     }
0983     if (jm80njets > 3) {
0984       _JM80njets->Fill(jm80njets, weight);
0985       _JM80ht->Fill(jm80HT, weight);
0986       _JM80pt1->Fill(jm80pt1, weight);
0987       _JM80pt2->Fill(jm80pt2, weight);
0988       _JM80pt3->Fill(jm80pt3, weight);
0989       _JM80pt4->Fill(jm80pt4, weight);
0990     }
0991 
0992     // select a di-jet event JME-10-001 variant
0993     double sumJetEt = 0;
0994     double sumPartPt = 0.;
0995     double sumChPartPt = 0.;
0996     double jpx = 0;
0997     double jpy = 0;
0998     if (nJets >= 2 && ij1 != genJets->end() && ij2 != genJets->end()) {
0999       if ((*ij1).pt() > 25. && (*ij1).pt() > 25.) {
1000         double deltaPhi = std::fabs((*ij1).phi() - (*ij2).phi()) / CLHEP::degree;
1001         if (deltaPhi > 180.)
1002           deltaPhi = 360. - deltaPhi;
1003         pt1pt2Dphi->Fill(deltaPhi, weight);
1004         if (std::fabs(deltaPhi) > 2.5 * CLHEP::degree) {
1005           nDijet->Fill(0.5, weight);
1006 
1007           for (unsigned int i = 0; i < hepmcGPCollection.size(); i++) {
1008             double eta = hepmcGPCollection[i]->momentum().eta();
1009             unsigned int iBin = getHFbin(eta);
1010             if (!isNeutrino(i) && iBin < CaloCellManager::nForwardEta) {
1011               hfDJ[iBin] += hepmcGPCollection[i]->momentum().rho();
1012             }
1013             if (!isNeutrino(i) && std::fabs(eta) < 5.) {
1014               sumPartPt += hepmcGPCollection[i]->momentum().perp();
1015               if (isCharged(i)) {
1016                 sumChPartPt += hepmcGPCollection[i]->momentum().perp();
1017               }
1018             }
1019           }
1020           for (unsigned int i = 0; i < CaloCellManager::nForwardEta; i++) {
1021             binW = theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i + 1] -
1022                    theEtaRanges[CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + i];
1023             dEdetaHFdj->Fill(i + 0.5, hfDJ[i] / binW);
1024           }
1025 
1026           double invMass = (*ij1).energy() * (*ij2).energy() - (*ij1).px() * (*ij2).px() - (*ij1).py() * (*ij2).py() -
1027                            (*ij1).pz() * (*ij2).pz();
1028           invMass = std::sqrt(invMass);
1029           pt1pt2InvM->Fill(invMass, weight);
1030 
1031           sumPt->Fill(sumPartPt, weight);
1032           sumChPt->Fill(sumChPartPt, weight);
1033 
1034           unsigned int nSelJets = 0;
1035           for (reco::GenJetCollection::const_iterator iter = genJets->begin(); iter != genJets->end(); ++iter) {
1036             double pt = (*iter).pt();
1037             double eta = (*iter).eta();
1038             if (std::fabs(eta) < 5.) {
1039               nSelJets++;
1040               binW = dNjdeta->getTH1()->GetBinWidth(1);
1041               dNjdeta->Fill(eta, 1. / binW * weight);
1042               binW = dNjdpt->getTH1()->GetBinWidth(1);
1043               dNjdpt->Fill(pt, 1. / binW * weight);
1044               sumJetEt += (*iter).pt();
1045               jpx += (*iter).px();
1046               jpy += (*iter).py();
1047             }
1048           }
1049 
1050           nj->Fill(nSelJets, weight);
1051           double mEt = std::sqrt(jpx * jpx + jpy * jpy);
1052           sumJEt->Fill(sumJetEt, weight);
1053           missEtosumJEt->Fill(mEt / sumJetEt, weight);
1054 
1055           if (nSelJets >= 3) {
1056             pt3Frac->Fill((*ij3).pt() / (pt1 + pt2), weight);
1057           }
1058 
1059           pt1pt2optot->Fill(pt1 + pt2, (pt1 + pt2) / sumJetEt);
1060           pt1pt2balance->Fill((pt1 - pt2) / (pt1 + pt2), weight);
1061         }
1062       }
1063     }
1064   }
1065 
1066   //compute differential jet rates
1067   std::vector<const HepMC::GenParticle*> qcdActivity;
1068   HepMCValidationHelper::removeIsolatedLeptons(myGenEvent, 0.2, 3., qcdActivity);
1069   //HepMCValidationHelper::allStatus1(myGenEvent, qcdActivity);
1070   //fill PseudoJets to use fastjet
1071   std::vector<fastjet::PseudoJet> vecs;
1072   int counterUser = 1;
1073   std::vector<const HepMC::GenParticle*>::const_iterator iqcdact;
1074   for (iqcdact = qcdActivity.begin(); iqcdact != qcdActivity.end(); ++iqcdact) {
1075     const HepMC::FourVector& fmom = (*iqcdact)->momentum();
1076     fastjet::PseudoJet pseudoJet(fmom.px(), fmom.py(), fmom.pz(), fmom.e());
1077     pseudoJet.set_user_index(counterUser);
1078     vecs.push_back(pseudoJet);
1079     ++counterUser;
1080   }
1081   //compute jets
1082   fastjet::ClusterSequence cseq(vecs, fastjet::JetDefinition(fastjet::kt_algorithm, 1., fastjet::E_scheme));
1083   //access the cluster sequence and get the relevant info
1084   djr10->Fill(std::log10(sqrt(cseq.exclusive_dmerge(0))), weight);
1085   djr21->Fill(std::log10(sqrt(cseq.exclusive_dmerge(1))), weight);
1086   djr32->Fill(std::log10(sqrt(cseq.exclusive_dmerge(2))), weight);
1087   djr43->Fill(std::log10(sqrt(cseq.exclusive_dmerge(3))), weight);
1088 
1089   // compute sumEt for all stable particles
1090   std::vector<const HepMC::GenParticle*> allStable;
1091   HepMCValidationHelper::allStatus1(myGenEvent, allStable);
1092 
1093   double sumEt = 0.;
1094   double sumEt1 = 0.;
1095   double sumEt2 = 0.;
1096   double sumEt3 = 0.;
1097   double sumEt4 = 0.;
1098   double sumEt5 = 0.;
1099 
1100   for (std::vector<const HepMC::GenParticle*>::const_iterator iter = allStable.begin(); iter != allStable.end();
1101        ++iter) {
1102     double thisEta = fabs((*iter)->momentum().eta());
1103 
1104     if (thisEta < 5.) {
1105       const HepMC::FourVector mom = (*iter)->momentum();
1106       double px = mom.px();
1107       double py = mom.py();
1108       double pz = mom.pz();
1109       double E = mom.e();
1110       double thisSumEt = (sqrt(px * px + py * py) * E / sqrt(px * px + py * py + pz * pz));
1111       sumEt += thisSumEt;
1112       if (thisEta < 1.0)
1113         sumEt1 += thisSumEt;
1114       else if (thisEta < 2.0)
1115         sumEt2 += thisSumEt;
1116       else if (thisEta < 3.0)
1117         sumEt3 += thisSumEt;
1118       else if (thisEta < 4.0)
1119         sumEt4 += thisSumEt;
1120       else
1121         sumEt5 += thisSumEt;
1122     }
1123   }
1124 
1125   if (sumEt > 0.)
1126     _sumEt->Fill(sumEt, weight);
1127   if (sumEt1 > 0.)
1128     _sumEt1->Fill(sumEt1, weight);
1129   if (sumEt2 > 0.)
1130     _sumEt2->Fill(sumEt2, weight);
1131   if (sumEt3 > 0.)
1132     _sumEt3->Fill(sumEt3, weight);
1133   if (sumEt4 > 0.)
1134     _sumEt4->Fill(sumEt4, weight);
1135   if (sumEt5 > 0.)
1136     _sumEt5->Fill(sumEt5, weight);
1137 
1138   delete myGenEvent;
1139 }  //analyze
1140 
1141 bool MBUEandQCDValidation::isCharged(unsigned int i) {
1142   bool status = false;
1143   if (hepmcGPCollection.size() < i + 1) {
1144     return status;
1145   } else {
1146     status = (hepmcCharge[i] != 0. && hepmcCharge[i] != -999.);
1147   }
1148   return status;
1149 }
1150 
1151 bool MBUEandQCDValidation::isNeutral(unsigned int i) {
1152   bool status = false;
1153   int pdgId = std::abs(hepmcGPCollection[i]->pdg_id());
1154   if (hepmcGPCollection.size() < i + 1) {
1155     return status;
1156   } else {
1157     status = (hepmcCharge[i] == 0. && pdgId != 12 && pdgId != 14 && pdgId != 16);
1158   }
1159   return status;
1160 }
1161 
1162 bool MBUEandQCDValidation::isNeutrino(unsigned int i) {
1163   bool status = false;
1164   int pdgId = std::abs(hepmcGPCollection[i]->pdg_id());
1165   if (hepmcGPCollection.size() < i + 1) {
1166     return status;
1167   } else {
1168     status = (pdgId == 12 || pdgId == 14 || pdgId == 16);
1169   }
1170   return status;
1171 }
1172 
1173 unsigned int MBUEandQCDValidation::getHFbin(double eta) {
1174   unsigned int iBin = 999;
1175 
1176   std::vector<double> theEtaRanges(theCalo->getEtaRanges());
1177 
1178   for (unsigned int i = CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta;
1179        i < CaloCellManager::nBarrelEta + CaloCellManager::nEndcapEta + CaloCellManager::nForwardEta;
1180        i++) {
1181     if (std::fabs(eta) >= theEtaRanges[i] && std::fabs(eta) < theEtaRanges[i + 1]) {
1182       iBin = i - CaloCellManager::nBarrelEta - CaloCellManager::nEndcapEta;
1183     }
1184   }
1185 
1186   return iBin;
1187 }
1188 
1189 const unsigned int MBUEandQCDValidation::nphiBin = 36;
1190 const unsigned int MBUEandQCDValidation::initSize = 1000;
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!