Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:31 

0001 //
0002 // Jet Analyzer class for heavy ion jets. for DQM jet analysis monitoring
0003 // For CMSSW_7_4_X, especially reading background subtracted jets
0004 // author: Raghav Kunnawalkam Elayavalli, Mohammed Zakaria (co Author)
0005 //         Jan 12th 2015
0006 //         Rutgers University, email: raghav.k.e at CERN dot CH
0007 //         UIC, email: mzakaria @ CERN dot CH
0008 
0009 #include "DQMOffline/JetMET/interface/JetAnalyzer_HeavyIons.h"
0010 
0011 using namespace edm;
0012 using namespace reco;
0013 using namespace std;
0014 
0015 // declare the constructors:
0016 
0017 JetAnalyzer_HeavyIons::JetAnalyzer_HeavyIons(const edm::ParameterSet &iConfig)
0018     : mInputCollection(iConfig.getParameter<edm::InputTag>("src")),
0019       mInputVtxCollection(iConfig.getUntrackedParameter<edm::InputTag>("srcVtx", edm::InputTag("hiSelectedVertex"))),
0020       mInputPFCandCollection(iConfig.getParameter<edm::InputTag>("PFcands")),
0021       mInputCsCandCollection(iConfig.exists("CScands") ? iConfig.getParameter<edm::InputTag>("CScands")
0022                                                        : edm::InputTag()),
0023       mOutputFile(iConfig.getUntrackedParameter<std::string>("OutputFile", "")),
0024       JetType(iConfig.getUntrackedParameter<std::string>("JetType")),
0025       UEAlgo(iConfig.getUntrackedParameter<std::string>("UEAlgo")),
0026       mRecoJetPtThreshold(iConfig.getParameter<double>("recoJetPtThreshold")),
0027       mReverseEnergyFractionThreshold(iConfig.getParameter<double>("reverseEnergyFractionThreshold")),
0028       mRThreshold(iConfig.getParameter<double>("RThreshold")),
0029       JetCorrectionService(iConfig.getParameter<std::string>("JetCorrections")) {
0030   std::string inputCollectionLabel(mInputCollection.label());
0031 
0032   isCaloJet = (std::string("calo") == JetType);
0033   isJPTJet = (std::string("jpt") == JetType);
0034   isPFJet = (std::string("pf") == JetType);
0035 
0036   //consumes
0037   pvToken_ = consumes<std::vector<reco::Vertex>>(edm::InputTag("offlinePrimaryVertices"));
0038   caloTowersToken_ = consumes<CaloTowerCollection>(edm::InputTag("towerMaker"));
0039   if (isCaloJet)
0040     caloJetsToken_ = consumes<reco::CaloJetCollection>(mInputCollection);
0041   if (isJPTJet)
0042     jptJetsToken_ = consumes<reco::JPTJetCollection>(mInputCollection);
0043   if (isPFJet) {
0044     if (std::string("Pu") == UEAlgo)
0045       basicJetsToken_ = consumes<reco::BasicJetCollection>(mInputCollection);
0046     if (std::string("Cs") == UEAlgo)
0047       pfJetsToken_ = consumes<reco::PFJetCollection>(mInputCollection);
0048   }
0049 
0050   pfCandToken_ = consumes<reco::PFCandidateCollection>(mInputPFCandCollection);
0051   csCandToken_ = mayConsume<reco::PFCandidateCollection>(mInputCsCandCollection);
0052   pfCandViewToken_ = consumes<reco::CandidateView>(mInputPFCandCollection);
0053   caloCandViewToken_ = consumes<reco::CandidateView>(edm::InputTag("towerMaker"));
0054 
0055   centralityTag_ = iConfig.getParameter<InputTag>("centralitycollection");
0056   centralityToken = consumes<reco::Centrality>(centralityTag_);
0057 
0058   centralityBinToken = mayConsume<int>(iConfig.exists("centralitybincollection")
0059                                            ? iConfig.getParameter<edm::InputTag>("centralitybincollection")
0060                                            : edm::InputTag());
0061 
0062   hiVertexToken_ = consumes<std::vector<reco::Vertex>>(mInputVtxCollection);
0063 
0064   etaToken_ = mayConsume<std::vector<double>>(iConfig.exists("etaMap") ? iConfig.getParameter<edm::InputTag>("etaMap")
0065                                                                        : edm::InputTag());
0066   rhoToken_ = mayConsume<std::vector<double>>(iConfig.exists("rho") ? iConfig.getParameter<edm::InputTag>("rho")
0067                                                                     : edm::InputTag());
0068   rhomToken_ = mayConsume<std::vector<double>>(iConfig.exists("rhom") ? iConfig.getParameter<edm::InputTag>("rhom")
0069                                                                       : edm::InputTag());
0070 
0071   // need to initialize the PF cand histograms : which are also event variables
0072   if (isPFJet) {
0073     mNPFpart = nullptr;
0074     mPFPt = nullptr;
0075     mPFEta = nullptr;
0076     mPFPhi = nullptr;
0077 
0078     mSumPFPt = nullptr;
0079     mSumPFPt_eta = nullptr;
0080     mSumSquaredPFPt = nullptr;
0081     mSumSquaredPFPt_eta = nullptr;
0082     mSumPFPt_HF = nullptr;
0083 
0084     mPFDeltaR = nullptr;
0085     mPFDeltaR_Scaled_R = nullptr;
0086 
0087     for (int ieta = 0; ieta < etaBins_; ieta++) {
0088       mSumPFPtEtaDep[ieta] = nullptr;
0089     }
0090 
0091     //cs-specific histograms
0092     mRhoDist_vsEta = nullptr;
0093     mRhoMDist_vsEta = nullptr;
0094     mRhoDist_vsPt = nullptr;
0095     mRhoMDist_vsPt = nullptr;
0096 
0097     rhoEtaRange = nullptr;
0098     for (int ieta = 0; ieta < etaBins_; ieta++) {
0099       mCSCandpT_vsPt[ieta] = nullptr;
0100       mRhoDist_vsCent[ieta] = nullptr;
0101       mRhoMDist_vsCent[ieta] = nullptr;
0102       for (int ipt = 0; ipt < ptBins_; ipt++) {
0103         mSubtractedEFrac[ipt][ieta] = nullptr;
0104         mSubtractedE[ipt][ieta] = nullptr;
0105       }
0106     }
0107 
0108     mPFCandpT_vs_eta_Unknown = nullptr;        // pf id 0
0109     mPFCandpT_vs_eta_ChargedHadron = nullptr;  // pf id - 1
0110     mPFCandpT_vs_eta_electron = nullptr;       // pf id - 2
0111     mPFCandpT_vs_eta_muon = nullptr;           // pf id - 3
0112     mPFCandpT_vs_eta_photon = nullptr;         // pf id - 4
0113     mPFCandpT_vs_eta_NeutralHadron = nullptr;  // pf id - 5
0114     mPFCandpT_vs_eta_HadE_inHF = nullptr;      // pf id - 6
0115     mPFCandpT_vs_eta_EME_inHF = nullptr;       // pf id - 7
0116 
0117     mPFCandpT_Barrel_Unknown = nullptr;        // pf id 0
0118     mPFCandpT_Barrel_ChargedHadron = nullptr;  // pf id - 1
0119     mPFCandpT_Barrel_electron = nullptr;       // pf id - 2
0120     mPFCandpT_Barrel_muon = nullptr;           // pf id - 3
0121     mPFCandpT_Barrel_photon = nullptr;         // pf id - 4
0122     mPFCandpT_Barrel_NeutralHadron = nullptr;  // pf id - 5
0123     mPFCandpT_Barrel_HadE_inHF = nullptr;      // pf id - 6
0124     mPFCandpT_Barrel_EME_inHF = nullptr;       // pf id - 7
0125 
0126     mPFCandpT_Endcap_Unknown = nullptr;        // pf id 0
0127     mPFCandpT_Endcap_ChargedHadron = nullptr;  // pf id - 1
0128     mPFCandpT_Endcap_electron = nullptr;       // pf id - 2
0129     mPFCandpT_Endcap_muon = nullptr;           // pf id - 3
0130     mPFCandpT_Endcap_photon = nullptr;         // pf id - 4
0131     mPFCandpT_Endcap_NeutralHadron = nullptr;  // pf id - 5
0132     mPFCandpT_Endcap_HadE_inHF = nullptr;      // pf id - 6
0133     mPFCandpT_Endcap_EME_inHF = nullptr;       // pf id - 7
0134 
0135     mPFCandpT_Forward_Unknown = nullptr;        // pf id 0
0136     mPFCandpT_Forward_ChargedHadron = nullptr;  // pf id - 1
0137     mPFCandpT_Forward_electron = nullptr;       // pf id - 2
0138     mPFCandpT_Forward_muon = nullptr;           // pf id - 3
0139     mPFCandpT_Forward_photon = nullptr;         // pf id - 4
0140     mPFCandpT_Forward_NeutralHadron = nullptr;  // pf id - 5
0141     mPFCandpT_Forward_HadE_inHF = nullptr;      // pf id - 6
0142     mPFCandpT_Forward_EME_inHF = nullptr;       // pf id - 7
0143   }
0144   if (isCaloJet) {
0145     mNCalopart = nullptr;
0146     mCaloPt = nullptr;
0147     mCaloEta = nullptr;
0148     mCaloPhi = nullptr;
0149 
0150     mSumCaloPt = nullptr;
0151     mSumCaloPt_eta = nullptr;
0152     mSumSquaredCaloPt = nullptr;
0153     mSumSquaredCaloPt_eta = nullptr;
0154     mSumCaloPt_HF = nullptr;
0155 
0156     for (int ieta = 0; ieta < etaBins_; ieta++) {
0157       mSumCaloPtEtaDep[ieta] = nullptr;
0158     }
0159   }
0160 
0161   mSumpt = nullptr;
0162 
0163   // Events variables
0164   mNvtx = nullptr;
0165   mHF = nullptr;
0166 
0167   // added Jan 12th 2015
0168 
0169   // Jet parameters
0170   mEta = nullptr;
0171   mPhi = nullptr;
0172   mEnergy = nullptr;
0173   mP = nullptr;
0174   mPt = nullptr;
0175   mMass = nullptr;
0176   mConstituents = nullptr;
0177   mJetArea = nullptr;
0178   mjetpileup = nullptr;
0179   mNJets_40 = nullptr;
0180   mNJets = nullptr;
0181 }
0182 
0183 void JetAnalyzer_HeavyIons::bookHistograms(DQMStore::IBooker &ibooker,
0184                                            edm::Run const &iRun,
0185                                            edm::EventSetup const &iSetup) {
0186   ibooker.setCurrentFolder("JetMET/HIJetValidation/" + mInputCollection.label());
0187 
0188   TH2F *h2D_etabins_vs_pt2 =
0189       new TH2F("h2D_etabins_vs_pt2", ";#eta;sum p_{T}^{2}", etaBins_, edge_pseudorapidity, 10000, 0, 10000);
0190   TH2F *h2D_etabins_vs_pt =
0191       new TH2F("h2D_etabins_vs_pt", ";#eta;sum p_{T}", etaBins_, edge_pseudorapidity, 500, 0, 500);
0192   TH2F *h2D_pfcand_etabins_vs_pt =
0193       new TH2F("h2D_etabins_vs_pt", ";#eta;sum p_{T}", etaBins_, edge_pseudorapidity, 300, 0, 300);
0194 
0195   const int nHihfBins = 100;
0196   const double hihfBins[nHihfBins + 1] = {
0197       0,           11.282305,   11.82962,    12.344717,   13.029054,   13.698554,   14.36821,    15.140326,
0198       15.845786,   16.684441,   17.449186,   18.364939,   19.247023,   20.448898,   21.776642,   22.870239,
0199       24.405788,   26.366919,   28.340206,   30.661842,   33.657627,   36.656773,   40.028049,   44.274784,
0200       48.583706,   52.981358,   56.860199,   61.559853,   66.663689,   72.768196,   78.265915,   84.744431,
0201       92.483459,   100.281021,  108.646576,  117.023911,  125.901093,  135.224899,  147.046875,  159.864258,
0202       171.06015,   184.76535,   197.687103,  212.873535,  229.276413,  245.175369,  262.498322,  280.54599,
0203       299.570801,  317.188446,  336.99881,   357.960144,  374.725922,  400.638367,  426.062103,  453.07251,
0204       483.99704,   517.556396,  549.421143,  578.050781,  608.358643,  640.940979,  680.361755,  719.215027,
0205       757.798645,  793.882385,  839.83728,   887.268127,  931.233276,  980.856689,  1023.191833, 1080.281494,
0206       1138.363892, 1191.303345, 1251.439453, 1305.288818, 1368.290894, 1433.700684, 1501.597412, 1557.918335,
0207       1625.636475, 1695.08374,  1761.771484, 1848.941162, 1938.178345, 2027.55603,  2127.364014, 2226.186523,
0208       2315.188965, 2399.225342, 2501.608643, 2611.077881, 2726.316162, 2848.74707,  2972.975342, 3096.565674,
0209       3219.530762, 3361.178223, 3568.028564, 3765.690186, 50000};
0210 
0211   TH2F *h2D_etabins_forRho =
0212       new TH2F("etabinsForRho", "#rho vs. #eta;#eta;#rho", etaBins_, edge_pseudorapidity, 500, 0, 300);
0213   TH2F *h2D_ptBins_forRho = new TH2F("ptBinsForRho", "#rho vs. p_{T};p_{T};#rho", 300, 0, 300, 500, 0, 300);
0214   TH2F *h2D_centBins_forRho = new TH2F("centBinsForRho", "dummy;HIHF;#rho", nHihfBins, hihfBins, 500, 0, 300);
0215 
0216   TH2F *h2D_etabins_forRhoM =
0217       new TH2F("etabinsForRho", "#rho_{M} vs. #eta;#eta;#rho_{M}", etaBins_, edge_pseudorapidity, 100, 0, 1.5);
0218   TH2F *h2D_ptBins_forRhoM = new TH2F("ptBinsForRho", "#rho_{M} vs. p_{T};p_{T};#rho_{M}", 300, 0, 300, 100, 0, 1.5);
0219   TH2F *h2D_centBins_forRhoM = new TH2F("centBinsForRho", "dummy;HIHF;#rho_{M}", nHihfBins, hihfBins, 100, 0, 1.5);
0220 
0221   if (isPFJet) {
0222     mNPFpart = ibooker.book1D("NPFpart", "No of particle flow candidates", 1000, 0, 1000);
0223     mPFPt = ibooker.book1D("PFPt", "PF candidate p_{T}", 10000, -500, 500);
0224     mPFEta = ibooker.book1D("PFEta", "PF candidate #eta", 120, -6, 6);
0225     mPFPhi = ibooker.book1D("PFPhi", "PF candidate #phi", 70, -3.5, 3.5);
0226 
0227     mPFDeltaR = ibooker.book1D("PFDeltaR", "PF candidate DeltaR", 100, 0, 4);  //MZ
0228     mPFDeltaR_Scaled_R =
0229         ibooker.book1D("PFDeltaR_Scaled_R", "PF candidate DeltaR Divided by DeltaR square", 100, 0, 4);  //MZ
0230 
0231     mSumPFPt = ibooker.book1D("SumPFPt", "Sum of initial PF p_{T}", 1000, -10000, 10000);
0232     mSumPFPt_eta = ibooker.book2D("SumPFPt_etaBins", h2D_etabins_vs_pt);
0233 
0234     mSumSquaredPFPt = ibooker.book1D("SumSquaredPFPt", "Sum of initial PF p_{T} squared", 10000, 0, 10000);
0235     mSumSquaredPFPt_eta = ibooker.book2D("SumSquaredPFPt_etaBins", h2D_etabins_vs_pt2);
0236 
0237     mSumPFPt_HF = ibooker.book2D(
0238         "SumPFPt_HF", "HF energy (y axis) vs Sum initial PF p_{T} (x axis)", 1000, -1000, 1000, 1000, 0, 10000);
0239 
0240     for (int ieta = 0; ieta < etaBins_; ieta++) {
0241       int range = 1000;
0242       if (ieta < 2 || etaBins_ - ieta <= 2)
0243         range = 500;
0244       const char *lc = edge_pseudorapidity[ieta] < 0 ? "n" : "p";
0245       const char *rc = edge_pseudorapidity[ieta + 1] < 0 ? "n" : "p";
0246       std::string histoName =
0247           Form("mSumCaloPt_%s%.3g_%s%.3g", lc, abs(edge_pseudorapidity[ieta]), rc, abs(edge_pseudorapidity[ieta + 1]));
0248       for (int id = 0; id < 2; id++) {
0249         if (histoName.find('.') != std::string::npos) {
0250           histoName.replace(histoName.find('.'), 1, "p");
0251         }
0252       }
0253       mSumPFPtEtaDep[ieta] = ibooker.book1D(
0254           histoName.c_str(),
0255           Form("Sum PFPt in the eta range %.3g to %.3g", edge_pseudorapidity[ieta], edge_pseudorapidity[ieta + 1]),
0256           500,
0257           0,
0258           range);
0259     }
0260 
0261     if (std::string("Cs") == UEAlgo) {
0262       mRhoDist_vsEta = ibooker.book2D("rhoDist_vsEta", h2D_etabins_forRho);
0263       mRhoMDist_vsEta = ibooker.book2D("rhoMDist_vsEta", h2D_etabins_forRhoM);
0264       mRhoDist_vsPt = ibooker.book2D("rhoDist_vsPt", h2D_ptBins_forRho);
0265       mRhoMDist_vsPt = ibooker.book2D("rhoMDist_vsPt", h2D_ptBins_forRhoM);
0266 
0267       //this is kind of a janky way to fill the eta since i can't get it from the edm::Event here... - kjung
0268       rhoEtaRange = ibooker.book1D("rhoEtaRange", "", 500, -5.5, 5.5);
0269       for (int ieta = 0; ieta < etaBins_; ieta++) {
0270         mCSCandpT_vsPt[ieta] =
0271             ibooker.book1D(Form("csCandPt_etaBin%d", ieta), "CS candidate pt, eta-by-eta", 150, 0, 300);
0272 
0273         const char *lc = edge_pseudorapidity[ieta] < 0 ? "n" : "p";
0274         const char *rc = edge_pseudorapidity[ieta + 1] < 0 ? "n" : "p";
0275         std::string histoName = Form(
0276             "Dist_vsCent_%s%.3g_%s%.3g", lc, abs(edge_pseudorapidity[ieta]), rc, abs(edge_pseudorapidity[ieta + 1]));
0277         for (int id = 0; id < 2; id++) {
0278           if (histoName.find('.') != std::string::npos) {
0279             histoName.replace(histoName.find('.'), 1, "p");
0280           }
0281         }
0282         std::string rhoName = "rho";
0283         rhoName.append(histoName);
0284         h2D_centBins_forRho->SetTitle(Form(
0285             "#rho vs. HIHF in the range %.3g < #eta < %.3g", edge_pseudorapidity[ieta], edge_pseudorapidity[ieta + 1]));
0286         mRhoDist_vsCent[ieta] = ibooker.book2D(rhoName.c_str(), h2D_centBins_forRho);
0287         std::string rhoMName = "rhoM";
0288         rhoMName.append(histoName);
0289         h2D_centBins_forRhoM->SetTitle(Form("#rho_{M} vs. HIHF in the range %.3g < #eta < %.3g",
0290                                             edge_pseudorapidity[ieta],
0291                                             edge_pseudorapidity[ieta + 1]));
0292         mRhoMDist_vsCent[ieta] = ibooker.book2D(rhoMName.c_str(), h2D_centBins_forRhoM);
0293         for (int ipt = 0; ipt < ptBins_; ipt++) {
0294           mSubtractedEFrac[ipt][ieta] =
0295               ibooker.book1D(Form("subtractedEFrac_JetPt%d_to_%d_etaBin%d", ptBin[ipt], ptBin[ipt + 1], ieta),
0296                              "subtracted fraction of CS jet",
0297                              50,
0298                              0,
0299                              1);
0300           mSubtractedE[ipt][ieta] =
0301               ibooker.book1D(Form("subtractedE_JetPt%d_to_%d_etaBin%d", ptBin[ipt], ptBin[ipt + 1], ieta),
0302                              "subtracted total of CS jet",
0303                              300,
0304                              0,
0305                              300);
0306         }
0307         mCSCand_corrPFcand[ieta] = ibooker.book2D(
0308             Form("csCandCorrPF%d", ieta), "CS to PF candidate correlation, eta-by-eta", 300, 0, 300, 300, 0, 300);
0309       }
0310     }
0311 
0312     mPFCandpT_vs_eta_Unknown = ibooker.book2D("PF_cand_X_unknown", h2D_pfcand_etabins_vs_pt);         // pf id 0
0313     mPFCandpT_vs_eta_ChargedHadron = ibooker.book2D("PF_cand_chargedHad", h2D_pfcand_etabins_vs_pt);  // pf id - 1
0314     mPFCandpT_vs_eta_electron = ibooker.book2D("PF_cand_electron", h2D_pfcand_etabins_vs_pt);         // pf id - 2
0315     mPFCandpT_vs_eta_muon = ibooker.book2D("PF_cand_muon", h2D_pfcand_etabins_vs_pt);                 // pf id - 3
0316     mPFCandpT_vs_eta_photon = ibooker.book2D("PF_cand_photon", h2D_pfcand_etabins_vs_pt);             // pf id - 4
0317     mPFCandpT_vs_eta_NeutralHadron = ibooker.book2D("PF_cand_neutralHad", h2D_pfcand_etabins_vs_pt);  // pf id - 5
0318     mPFCandpT_vs_eta_HadE_inHF = ibooker.book2D("PF_cand_HadEner_inHF", h2D_pfcand_etabins_vs_pt);    // pf id - 6
0319     mPFCandpT_vs_eta_EME_inHF = ibooker.book2D("PF_cand_EMEner_inHF", h2D_pfcand_etabins_vs_pt);      // pf id - 7
0320 
0321     mPFCandpT_Barrel_Unknown = ibooker.book1D("mPFCandpT_Barrel_Unknown",
0322                                               Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0323                                               300,
0324                                               0,
0325                                               300);  // pf id  - 0
0326     mPFCandpT_Barrel_ChargedHadron = ibooker.book1D("mPFCandpT_Barrel_ChargedHadron",
0327                                                     Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0328                                                     300,
0329                                                     0,
0330                                                     300);  // pf id - 1
0331     mPFCandpT_Barrel_electron = ibooker.book1D("mPFCandpT_Barrel_electron",
0332                                                Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0333                                                300,
0334                                                0,
0335                                                300);  // pf id - 2
0336     mPFCandpT_Barrel_muon = ibooker.book1D("mPFCandpT_Barrel_muon",
0337                                            Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0338                                            300,
0339                                            0,
0340                                            300);  // pf id - 3
0341     mPFCandpT_Barrel_photon = ibooker.book1D("mPFCandpT_Barrel_photon",
0342                                              Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0343                                              300,
0344                                              0,
0345                                              300);  // pf id - 4
0346     mPFCandpT_Barrel_NeutralHadron = ibooker.book1D("mPFCandpT_Barrel_NeutralHadron",
0347                                                     Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0348                                                     300,
0349                                                     0,
0350                                                     300);  // pf id - 5
0351     mPFCandpT_Barrel_HadE_inHF = ibooker.book1D("mPFCandpT_Barrel_HadE_inHF",
0352                                                 Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0353                                                 300,
0354                                                 0,
0355                                                 300);  // pf id - 6
0356     mPFCandpT_Barrel_EME_inHF = ibooker.book1D("mPFCandpT_Barrel_EME_inHF",
0357                                                Form(";PF candidate p_{T}, |#eta|<%2.2f; counts", BarrelEta),
0358                                                300,
0359                                                0,
0360                                                300);  // pf id - 7
0361 
0362     mPFCandpT_Endcap_Unknown =
0363         ibooker.book1D("mPFCandpT_Endcap_Unknown",
0364                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0365                        300,
0366                        0,
0367                        300);  // pf id - 0
0368     mPFCandpT_Endcap_ChargedHadron =
0369         ibooker.book1D("mPFCandpT_Endcap_ChargedHadron",
0370                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0371                        300,
0372                        0,
0373                        300);  // pf id - 1
0374     mPFCandpT_Endcap_electron =
0375         ibooker.book1D("mPFCandpT_Endcap_electron",
0376                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0377                        300,
0378                        0,
0379                        300);  // pf id - 2
0380     mPFCandpT_Endcap_muon =
0381         ibooker.book1D("mPFCandpT_Endcap_muon",
0382                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0383                        300,
0384                        0,
0385                        300);  // pf id - 3
0386     mPFCandpT_Endcap_photon =
0387         ibooker.book1D("mPFCandpT_Endcap_photon",
0388                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0389                        300,
0390                        0,
0391                        300);  // pf id - 4
0392     mPFCandpT_Endcap_NeutralHadron =
0393         ibooker.book1D("mPFCandpT_Endcap_NeutralHadron",
0394                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0395                        300,
0396                        0,
0397                        300);  // pf id - 5
0398     mPFCandpT_Endcap_HadE_inHF =
0399         ibooker.book1D("mPFCandpT_Endcap_HadE_inHF",
0400                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0401                        300,
0402                        0,
0403                        300);  // pf id - 6
0404     mPFCandpT_Endcap_EME_inHF =
0405         ibooker.book1D("mPFCandpT_Endcap_EME_inHF",
0406                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", BarrelEta, EndcapEta),
0407                        300,
0408                        0,
0409                        300);  // pf id - 7
0410 
0411     mPFCandpT_Forward_Unknown =
0412         ibooker.book1D("mPFCandpT_Forward_Unknown",
0413                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0414                        300,
0415                        0,
0416                        300);  // pf id - 0
0417     mPFCandpT_Forward_ChargedHadron =
0418         ibooker.book1D("mPFCandpT_Forward_ChargedHadron",
0419                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0420                        300,
0421                        0,
0422                        300);  // pf id - 1
0423     mPFCandpT_Forward_electron =
0424         ibooker.book1D("mPFCandpT_Forward_electron",
0425                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0426                        300,
0427                        0,
0428                        300);  // pf id - 2
0429     mPFCandpT_Forward_muon =
0430         ibooker.book1D("mPFCandpT_Forward_muon",
0431                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0432                        300,
0433                        0,
0434                        300);  // pf id - 3
0435     mPFCandpT_Forward_photon =
0436         ibooker.book1D("mPFCandpT_Forward_photon",
0437                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0438                        300,
0439                        0,
0440                        300);  // pf id - 4
0441     mPFCandpT_Forward_NeutralHadron =
0442         ibooker.book1D("mPFCandpT_Forward_NeutralHadron",
0443                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0444                        300,
0445                        0,
0446                        300);  // pf id - 5
0447     mPFCandpT_Forward_HadE_inHF =
0448         ibooker.book1D("mPFCandpT_Forward_HadE_inHF",
0449                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0450                        300,
0451                        0,
0452                        300);  // pf id - 6
0453     mPFCandpT_Forward_EME_inHF =
0454         ibooker.book1D("mPFCandpT_Forward_EME_inHF",
0455                        Form(";PF candidate p_{T}, %2.2f<|#eta|<%2.2f; counts", EndcapEta, ForwardEta),
0456                        300,
0457                        0,
0458                        300);  // pf id - 7
0459   }
0460 
0461   if (isCaloJet) {
0462     mNCalopart = ibooker.book1D("NCalopart", "No of particle flow candidates", 1000, 0, 10000);
0463     mCaloPt = ibooker.book1D("CaloPt", "Calo candidate p_{T}", 1000, -5000, 5000);
0464     mCaloEta = ibooker.book1D("CaloEta", "Calo candidate #eta", 120, -6, 6);
0465     mCaloPhi = ibooker.book1D("CaloPhi", "Calo candidate #phi", 70, -3.5, 3.5);
0466 
0467     mSumCaloPt = ibooker.book1D("SumCaloPt", "Sum Calo p_{T}", 1000, -10000, 10000);
0468     mSumCaloPt_eta = ibooker.book2D("SumCaloPt_etaBins", h2D_etabins_vs_pt);
0469 
0470     mSumSquaredCaloPt = ibooker.book1D("SumSquaredCaloPt", "Sum of initial Calo tower p_{T} squared", 10000, 0, 10000);
0471     mSumSquaredCaloPt_eta = ibooker.book2D("SumSquaredCaloPt_etaBins", h2D_etabins_vs_pt2);
0472 
0473     mSumCaloPt_HF =
0474         ibooker.book2D("SumCaloPt_HF", "HF Energy (y axis) vs Sum Calo tower p_{T}", 1000, -1000, 1000, 1000, 0, 10000);
0475 
0476     for (int ieta = 0; ieta < etaBins_; ieta++) {
0477       int range = 1000;
0478       if (ieta < 2 || etaBins_ - ieta <= 2)
0479         range = 5000;
0480       const char *lc = edge_pseudorapidity[ieta] < 0 ? "n" : "p";
0481       const char *rc = edge_pseudorapidity[ieta + 1] < 0 ? "n" : "p";
0482       std::string histoName =
0483           Form("mSumCaloPt_%s%.3g_%s%.3g", lc, abs(edge_pseudorapidity[ieta]), rc, abs(edge_pseudorapidity[ieta + 1]));
0484       for (int id = 0; id < 2; id++) {
0485         if (histoName.find('.') != std::string::npos) {
0486           histoName.replace(histoName.find('.'), 1, "p");
0487         }
0488       }
0489       mSumCaloPtEtaDep[ieta] = ibooker.book1D(histoName.c_str(),
0490                                               Form("Sum Calo tower Pt in the eta range %.3g to %.3g",
0491                                                    edge_pseudorapidity[ieta],
0492                                                    edge_pseudorapidity[ieta + 1]),
0493                                               1000,
0494                                               -1 * range,
0495                                               range);
0496     }
0497   }
0498 
0499   // particle flow variables histograms
0500   mSumpt = ibooker.book1D("SumpT", "Sum p_{T} of all the PF candidates per event", 1000, 0, 10000);
0501 
0502   // Event variables
0503   mNvtx = ibooker.book1D("Nvtx", "number of vertices", 60, 0, 60);
0504   mHF = ibooker.book1D("HF", "HF energy distribution", 1000, 0, 10000);
0505 
0506   // Jet parameters
0507   mEta = ibooker.book1D("Eta", "Eta", 120, -6, 6);
0508   mPhi = ibooker.book1D("Phi", "Phi", 70, -3.5, 3.5);
0509   mPt = ibooker.book1D("Pt", "Pt", 1000, 0, 500);
0510   mP = ibooker.book1D("P", "P", 100, 0, 1000);
0511   mEnergy = ibooker.book1D("Energy", "Energy", 100, 0, 1000);
0512   mMass = ibooker.book1D("Mass", "Mass", 100, 0, 200);
0513   mConstituents = ibooker.book1D("Constituents", "Constituents", 100, 0, 100);
0514   mJetArea = ibooker.book1D("JetArea", "JetArea", 100, 0, 4);
0515   mjetpileup = ibooker.book1D("jetPileUp", "jetPileUp", 100, 0, 150);
0516   mNJets_40 = ibooker.book1D("NJets_pt_greater_40", "NJets pT > 40 GeV", 50, 0, 50);
0517   mNJets = ibooker.book1D("NJets", "NJets", 100, 0, 100);
0518 
0519   if (mOutputFile.empty())
0520     LogInfo("OutputInfo") << " Histograms will NOT be saved";
0521   else
0522     LogInfo("OutputInfo") << " Histograms will be saved to file:" << mOutputFile;
0523 
0524   delete h2D_etabins_vs_pt2;
0525   delete h2D_etabins_vs_pt;
0526   delete h2D_pfcand_etabins_vs_pt;
0527   delete h2D_etabins_forRho;
0528   delete h2D_ptBins_forRho;
0529   delete h2D_centBins_forRho;
0530   delete h2D_centBins_forRhoM;
0531 }
0532 
0533 //------------------------------------------------------------------------------
0534 // ~JetAnalyzer_HeavyIons
0535 //------------------------------------------------------------------------------
0536 JetAnalyzer_HeavyIons::~JetAnalyzer_HeavyIons() {}
0537 
0538 //------------------------------------------------------------------------------
0539 // analyze
0540 //------------------------------------------------------------------------------
0541 void JetAnalyzer_HeavyIons::analyze(const edm::Event &mEvent, const edm::EventSetup &mSetup) {
0542   // switch(mEvent.id().event() == 15296770)
0543   // case 1:
0544   //   break;
0545 
0546   // Get the primary vertices
0547   edm::Handle<vector<reco::Vertex>> pvHandle;
0548   mEvent.getByToken(pvToken_, pvHandle);
0549   reco::Vertex::Point vtx(0, 0, 0);
0550   edm::Handle<reco::VertexCollection> vtxs;
0551 
0552   mEvent.getByToken(hiVertexToken_, vtxs);
0553   int greatestvtx = 0;
0554   int nVertex = vtxs->size();
0555 
0556   for (unsigned int i = 0; i < vtxs->size(); ++i) {
0557     unsigned int daughter = (*vtxs)[i].tracksSize();
0558     if (daughter > (*vtxs)[greatestvtx].tracksSize())
0559       greatestvtx = i;
0560   }
0561 
0562   if (nVertex <= 0) {
0563     vtx = reco::Vertex::Point(0, 0, 0);
0564   }
0565   vtx = (*vtxs)[greatestvtx].position();
0566 
0567   int nGoodVertices = 0;
0568 
0569   if (pvHandle.isValid()) {
0570     for (unsigned i = 0; i < pvHandle->size(); i++) {
0571       if ((*pvHandle)[i].ndof() > 4 && (fabs((*pvHandle)[i].z()) <= 24) && (fabs((*pvHandle)[i].position().rho()) <= 2))
0572         nGoodVertices++;
0573     }
0574   }
0575 
0576   mNvtx->Fill(nGoodVertices);
0577 
0578   // Get the Jet collection
0579   //----------------------------------------------------------------------------
0580 
0581   std::vector<Jet> recoJets;
0582   recoJets.clear();
0583 
0584   edm::Handle<CaloJetCollection> caloJets;
0585   edm::Handle<JPTJetCollection> jptJets;
0586   edm::Handle<PFJetCollection> pfJets;
0587   edm::Handle<BasicJetCollection> basicJets;
0588 
0589   // Get the Particle flow candidates and the Voronoi variables
0590   edm::Handle<reco::PFCandidateCollection> pfCandidates;
0591   edm::Handle<reco::PFCandidateCollection> csCandidates;
0592   edm::Handle<CaloTowerCollection> caloCandidates;
0593   edm::Handle<reco::CandidateView> pfcandidates_;
0594   edm::Handle<reco::CandidateView> calocandidates_;
0595 
0596   //Get the new CS stuff
0597   edm::Handle<std::vector<double>> etaRanges;
0598   edm::Handle<std::vector<double>> rho;
0599   edm::Handle<std::vector<double>> rhom;
0600   if (std::string("Cs") == UEAlgo) {
0601     mEvent.getByToken(etaToken_, etaRanges);
0602     mEvent.getByToken(rhoToken_, rho);
0603     mEvent.getByToken(rhomToken_, rhom);
0604     const int rhoSize = (int)etaRanges->size();
0605     double rhoRange[rhoSize];
0606     for (int irho = 0; irho < rhoSize; irho++) {
0607       rhoRange[irho] = etaRanges->at(irho);
0608     }
0609     double yaxisLimits[501];
0610     for (int ibin = 0; ibin < 501; ibin++)
0611       yaxisLimits[ibin] = ibin * 2;
0612     if (mRhoDist_vsEta->getNbinsX() != rhoSize - 1) {
0613       mRhoDist_vsEta->getTH2F()->SetBins(
0614           rhoSize - 1, const_cast<double *>(rhoRange), 500, const_cast<double *>(yaxisLimits));
0615       mRhoMDist_vsEta->getTH2F()->SetBins(
0616           rhoSize - 1, const_cast<double *>(rhoRange), 500, const_cast<double *>(yaxisLimits));
0617     }
0618   }
0619 
0620   if (isCaloJet)
0621     mEvent.getByToken(caloJetsToken_, caloJets);
0622   if (isJPTJet)
0623     mEvent.getByToken(jptJetsToken_, jptJets);
0624   if (isPFJet) {
0625     if (std::string("Pu") == UEAlgo)
0626       mEvent.getByToken(basicJetsToken_, basicJets);
0627     if (std::string("Cs") == UEAlgo)
0628       mEvent.getByToken(pfJetsToken_, pfJets);
0629     if (std::string("Vs") == UEAlgo)
0630       return;  //avoid running Vs jets
0631   }
0632 
0633   mEvent.getByToken(pfCandToken_, pfCandidates);
0634   mEvent.getByToken(pfCandViewToken_, pfcandidates_);
0635   if (std::string("Cs") == UEAlgo)
0636     mEvent.getByToken(csCandToken_, csCandidates);
0637 
0638   mEvent.getByToken(caloTowersToken_, caloCandidates);
0639   mEvent.getByToken(caloCandViewToken_, calocandidates_);
0640 
0641   // get the centrality
0642   edm::Handle<reco::Centrality> cent;
0643   mEvent.getByToken(centralityToken, cent);  //_centralitytag comes from the cfg
0644 
0645   if (!cent.isValid())
0646     return;
0647 
0648   mHF->Fill(cent->EtHFtowerSum());
0649   Float_t HF_energy = cent->EtHFtowerSum();
0650 
0651   //for later when centrality gets added to RelVal
0652   //edm::Handle<int> cbin;
0653   //mEvent.getByToken(centralityBinToken, cbin);
0654 
0655   /*int hibin = -999;
0656   if(cbin.isValid()){
0657     hibin = *cbin;
0658   }*/
0659 
0660   const reco::PFCandidateCollection *pfCandidateColl = pfCandidates.product();
0661 
0662   Int_t NPFpart = 0;
0663   Int_t NCaloTower = 0;
0664   Float_t pfPt = 0;
0665   Float_t pfEta = 0;
0666   Float_t pfPhi = 0;
0667   Int_t pfID = 0;
0668   Float_t pfDeltaR = 0;
0669   Float_t caloPt = 0;
0670   Float_t caloEta = 0;
0671   Float_t caloPhi = 0;
0672   Float_t SumPt_value = 0;
0673 
0674   vector<vector<float>> numbers;
0675   vector<float> tempVector;
0676   numbers.clear();
0677   tempVector.clear();
0678 
0679   if (isCaloJet) {
0680     Float_t SumCaloPt[etaBins_];
0681     Float_t SumSquaredCaloPt[etaBins_];
0682 
0683     // Need to set up histograms to get the RMS values for each pT bin
0684     TH1F *hSumCaloPt[nedge_pseudorapidity - 1];
0685 
0686     for (int i = 0; i < etaBins_; ++i) {
0687       SumCaloPt[i] = 0;
0688       SumSquaredCaloPt[i] = 0;
0689       hSumCaloPt[i] = new TH1F(Form("hSumCaloPt_%d", i), "", 10000, -10000, 10000);
0690     }
0691 
0692     for (unsigned icand = 0; icand < caloCandidates->size(); icand++) {
0693       const CaloTower &tower = (*caloCandidates)[icand];
0694       reco::CandidateViewRef ref(calocandidates_, icand);
0695       //10 is tower pT min
0696       if (tower.p4(vtx).Et() < 0.1)
0697         continue;
0698 
0699       NCaloTower++;
0700 
0701       caloPt = tower.p4(vtx).Et();
0702       caloEta = tower.p4(vtx).Eta();
0703       caloPhi = tower.p4(vtx).Phi();
0704 
0705       for (size_t k = 0; k < nedge_pseudorapidity - 1; k++) {
0706         if (caloEta >= edge_pseudorapidity[k] && caloEta < edge_pseudorapidity[k + 1]) {
0707           SumCaloPt[k] = SumCaloPt[k] + caloPt;
0708           SumSquaredCaloPt[k] = SumSquaredCaloPt[k] + caloPt * caloPt;
0709           break;
0710         }  // eta selection statement
0711 
0712       }  // eta bin loop
0713 
0714       SumPt_value = SumPt_value + caloPt;
0715 
0716       mCaloPt->Fill(caloPt);
0717       mCaloEta->Fill(caloEta);
0718       mCaloPhi->Fill(caloPhi);
0719 
0720     }  // calo tower candidate  loop
0721 
0722     for (int k = 0; k < nedge_pseudorapidity - 1; k++) {
0723       hSumCaloPt[k]->Fill(SumCaloPt[k]);
0724 
0725     }  // eta bin loop
0726 
0727     Float_t Evt_SumCaloPt = 0;
0728     Float_t Evt_SumSquaredCaloPt = 0;
0729 
0730     for (int ieta = 0; ieta < etaBins_; ieta++) {
0731       mSumCaloPtEtaDep[ieta]->Fill(SumCaloPt[ieta]);
0732     }
0733 
0734     for (size_t k = 0; k < nedge_pseudorapidity - 1; k++) {
0735       Evt_SumCaloPt = Evt_SumCaloPt + SumCaloPt[k];
0736       mSumCaloPt_eta->Fill(edge_pseudorapidity[k], SumCaloPt[k]);
0737 
0738       Evt_SumSquaredCaloPt = Evt_SumSquaredCaloPt + SumSquaredCaloPt[k];
0739       mSumSquaredCaloPt_eta->Fill(edge_pseudorapidity[k], hSumCaloPt[k]->GetRMS(1));
0740 
0741       delete hSumCaloPt[k];
0742 
0743     }  // eta bin loop
0744 
0745     mSumCaloPt->Fill(Evt_SumCaloPt);
0746     mSumCaloPt_HF->Fill(Evt_SumCaloPt, HF_energy);
0747 
0748     mSumSquaredCaloPt->Fill(Evt_SumSquaredCaloPt);
0749 
0750     mNCalopart->Fill(NCaloTower);
0751     mSumpt->Fill(SumPt_value);
0752 
0753   }  // is calo jet
0754 
0755   if (isPFJet) {
0756     Float_t SumPFPt[etaBins_];
0757 
0758     Float_t SumSquaredPFPt[etaBins_];
0759 
0760     // Need to set up histograms to get the RMS values for each pT bin
0761     TH1F *hSumPFPt[nedge_pseudorapidity - 1];
0762 
0763     for (int i = 0; i < etaBins_; i++) {
0764       SumPFPt[i] = 0;
0765       SumSquaredPFPt[i] = 0;
0766 
0767       hSumPFPt[i] = new TH1F(Form("hSumPFPt_%d", i), "", 10000, -10000, 10000);
0768     }
0769 
0770     vector<vector<float>> PF_Space(1, vector<float>(3));
0771 
0772     if (std::string("Cs") == UEAlgo) {
0773       const reco::PFCandidateCollection *csCandidateColl = csCandidates.product();
0774 
0775       for (unsigned iCScand = 0; iCScand < csCandidateColl->size(); iCScand++) {
0776         assert(csCandidateColl->size() <= pfCandidateColl->size());
0777         const reco::PFCandidate csCandidate = csCandidateColl->at(iCScand);
0778         const reco::PFCandidate pfCandidate = pfCandidateColl->at(iCScand);
0779         int ieta = 0;
0780         while (csCandidate.eta() > edge_pseudorapidity[ieta] && ieta < etaBins_ - 1)
0781           ieta++;
0782         mCSCandpT_vsPt[ieta]->Fill(csCandidate.pt());
0783         mCSCand_corrPFcand[ieta]->Fill(csCandidate.pt(), pfCandidate.pt());
0784       }
0785     }
0786 
0787     for (unsigned icand = 0; icand < pfCandidateColl->size(); icand++) {
0788       const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
0789       reco::CandidateViewRef ref(pfcandidates_, icand);
0790       if (pfCandidate.pt() < 5)
0791         continue;
0792 
0793       NPFpart++;
0794       pfPt = pfCandidate.pt();
0795       pfEta = pfCandidate.eta();
0796       pfPhi = pfCandidate.phi();
0797       pfID = pfCandidate.particleId();
0798 
0799       bool isBarrel = false;
0800       bool isEndcap = false;
0801       bool isForward = false;
0802 
0803       if (fabs(pfEta) < BarrelEta)
0804         isBarrel = true;
0805       if (fabs(pfEta) >= BarrelEta && fabs(pfEta) < EndcapEta)
0806         isEndcap = true;
0807       if (fabs(pfEta) >= EndcapEta && fabs(pfEta) < ForwardEta)
0808         isForward = true;
0809 
0810       switch (pfID) {
0811         case 0:
0812           mPFCandpT_vs_eta_Unknown->Fill(pfPt, pfEta);
0813           if (isBarrel)
0814             mPFCandpT_Barrel_Unknown->Fill(pfPt);
0815           if (isEndcap)
0816             mPFCandpT_Endcap_Unknown->Fill(pfPt);
0817           if (isForward)
0818             mPFCandpT_Forward_Unknown->Fill(pfPt);
0819           break;
0820         case 1:
0821           mPFCandpT_vs_eta_ChargedHadron->Fill(pfPt, pfEta);
0822           if (isBarrel)
0823             mPFCandpT_Barrel_ChargedHadron->Fill(pfPt);
0824           if (isEndcap)
0825             mPFCandpT_Endcap_ChargedHadron->Fill(pfPt);
0826           if (isForward)
0827             mPFCandpT_Forward_ChargedHadron->Fill(pfPt);
0828           break;
0829         case 2:
0830           mPFCandpT_vs_eta_electron->Fill(pfPt, pfEta);
0831           if (isBarrel)
0832             mPFCandpT_Barrel_electron->Fill(pfPt);
0833           if (isEndcap)
0834             mPFCandpT_Endcap_electron->Fill(pfPt);
0835           if (isForward)
0836             mPFCandpT_Forward_electron->Fill(pfPt);
0837           break;
0838         case 3:
0839           mPFCandpT_vs_eta_muon->Fill(pfPt, pfEta);
0840           if (isBarrel)
0841             mPFCandpT_Barrel_muon->Fill(pfPt);
0842           if (isEndcap)
0843             mPFCandpT_Endcap_muon->Fill(pfPt);
0844           if (isForward)
0845             mPFCandpT_Forward_muon->Fill(pfPt);
0846           break;
0847         case 4:
0848           mPFCandpT_vs_eta_photon->Fill(pfPt, pfEta);
0849           if (isBarrel)
0850             mPFCandpT_Barrel_photon->Fill(pfPt);
0851           if (isEndcap)
0852             mPFCandpT_Endcap_photon->Fill(pfPt);
0853           if (isForward)
0854             mPFCandpT_Forward_photon->Fill(pfPt);
0855           break;
0856         case 5:
0857           mPFCandpT_vs_eta_NeutralHadron->Fill(pfPt, pfEta);
0858           if (isBarrel)
0859             mPFCandpT_Barrel_NeutralHadron->Fill(pfPt);
0860           if (isEndcap)
0861             mPFCandpT_Endcap_NeutralHadron->Fill(pfPt);
0862           if (isForward)
0863             mPFCandpT_Forward_NeutralHadron->Fill(pfPt);
0864           break;
0865         case 6:
0866           mPFCandpT_vs_eta_HadE_inHF->Fill(pfPt, pfEta);
0867           if (isBarrel)
0868             mPFCandpT_Barrel_HadE_inHF->Fill(pfPt);
0869           if (isEndcap)
0870             mPFCandpT_Endcap_HadE_inHF->Fill(pfPt);
0871           if (isForward)
0872             mPFCandpT_Forward_HadE_inHF->Fill(pfPt);
0873           break;
0874         case 7:
0875           mPFCandpT_vs_eta_EME_inHF->Fill(pfPt, pfEta);
0876           if (isBarrel)
0877             mPFCandpT_Barrel_EME_inHF->Fill(pfPt);
0878           if (isEndcap)
0879             mPFCandpT_Endcap_EME_inHF->Fill(pfPt);
0880           if (isForward)
0881             mPFCandpT_Forward_EME_inHF->Fill(pfPt);
0882           break;
0883       }
0884 
0885       //Fill 2d vector matrix
0886       tempVector.push_back(pfPt);
0887       tempVector.push_back(pfEta);
0888       tempVector.push_back(pfPhi);
0889 
0890       numbers.push_back(tempVector);
0891       tempVector.clear();
0892 
0893       for (size_t k = 0; k < nedge_pseudorapidity - 1; k++) {
0894         if (pfEta >= edge_pseudorapidity[k] && pfEta < edge_pseudorapidity[k + 1]) {
0895           SumPFPt[k] = SumPFPt[k] + pfPt;
0896           SumSquaredPFPt[k] = SumSquaredPFPt[k] + pfPt * pfPt;
0897           break;
0898         }  // eta selection statement
0899 
0900       }  // eta bin loop
0901 
0902       SumPt_value = SumPt_value + pfPt;
0903 
0904       mPFPt->Fill(pfPt);
0905       mPFEta->Fill(pfEta);
0906       mPFPhi->Fill(pfPhi);
0907 
0908     }  // pf candidate loop
0909 
0910     for (int k = 0; k < nedge_pseudorapidity - 1; k++) {
0911       hSumPFPt[k]->Fill(SumPFPt[k]);
0912 
0913     }  // eta bin loop
0914 
0915     Float_t Evt_SumPFPt = 0;
0916     Float_t Evt_SumSquaredPFPt = 0;
0917 
0918     for (int ieta = 0; ieta < etaBins_; ieta++) {
0919       mSumPFPtEtaDep[ieta]->Fill(SumPFPt[ieta]);
0920     }
0921 
0922     for (size_t k = 0; k < nedge_pseudorapidity - 1; k++) {
0923       Evt_SumPFPt = Evt_SumPFPt + SumPFPt[k];
0924       mSumPFPt_eta->Fill(edge_pseudorapidity[k], SumPFPt[k]);
0925 
0926       Evt_SumSquaredPFPt = Evt_SumSquaredPFPt + SumSquaredPFPt[k];
0927       mSumSquaredPFPt_eta->Fill(edge_pseudorapidity[k], hSumPFPt[k]->GetRMS(1));
0928 
0929       delete hSumPFPt[k];
0930 
0931     }  // eta bin loop
0932 
0933     mSumPFPt->Fill(Evt_SumPFPt);
0934     mSumPFPt_HF->Fill(Evt_SumPFPt, HF_energy);
0935 
0936     mSumSquaredPFPt->Fill(Evt_SumSquaredPFPt);
0937 
0938     mNPFpart->Fill(NPFpart);
0939     mSumpt->Fill(SumPt_value);
0940   }
0941 
0942   if (isCaloJet) {
0943     for (unsigned ijet = 0; ijet < caloJets->size(); ijet++) {
0944       recoJets.push_back((*caloJets)[ijet]);
0945     }
0946   }
0947 
0948   if (isJPTJet) {
0949     for (unsigned ijet = 0; ijet < jptJets->size(); ijet++)
0950       recoJets.push_back((*jptJets)[ijet]);
0951   }
0952 
0953   if (isPFJet) {
0954     if (std::string("Pu") == UEAlgo) {
0955       for (unsigned ijet = 0; ijet < basicJets->size(); ijet++) {
0956         recoJets.push_back((*basicJets)[ijet]);
0957       }
0958     }
0959     if (std::string("Cs") == UEAlgo) {
0960       for (unsigned ijet = 0; ijet < pfJets->size(); ijet++) {
0961         recoJets.push_back((*pfJets)[ijet]);
0962       }
0963     }
0964   }
0965 
0966   if (isCaloJet && !caloJets.isValid()) {
0967     return;
0968   }
0969   if (isJPTJet && !jptJets.isValid()) {
0970     return;
0971   }
0972   if (isPFJet) {
0973     if (std::string("Pu") == UEAlgo) {
0974       if (!basicJets.isValid())
0975         return;
0976     }
0977     if (std::string("Cs") == UEAlgo) {
0978       if (!pfJets.isValid())
0979         return;
0980     }
0981     if (std::string("Vs") == UEAlgo) {
0982       return;
0983     }
0984   }
0985 
0986   int nJet_40 = 0;
0987 
0988   mNJets->Fill(recoJets.size());
0989 
0990   for (unsigned ijet = 0; ijet < recoJets.size(); ijet++) {
0991     if (recoJets[ijet].pt() > mRecoJetPtThreshold) {
0992       //counting forward and barrel jets
0993       // get an idea of no of jets with pT>40 GeV
0994       if (recoJets[ijet].pt() > 40)
0995         nJet_40++;
0996       if (mEta)
0997         mEta->Fill(recoJets[ijet].eta());
0998       if (mjetpileup)
0999         mjetpileup->Fill(recoJets[ijet].pileup());
1000       if (mJetArea)
1001         mJetArea->Fill(recoJets[ijet].jetArea());
1002       if (mPhi)
1003         mPhi->Fill(recoJets[ijet].phi());
1004       if (mEnergy)
1005         mEnergy->Fill(recoJets[ijet].energy());
1006       if (mP)
1007         mP->Fill(recoJets[ijet].p());
1008       if (mPt)
1009         mPt->Fill(recoJets[ijet].pt());
1010       if (mMass)
1011         mMass->Fill(recoJets[ijet].mass());
1012       if (mConstituents)
1013         mConstituents->Fill(recoJets[ijet].nConstituents());
1014 
1015       if (std::string("Cs") == UEAlgo) {
1016         int ipt = 0, ieta = 0;
1017         while (recoJets[ijet].pt() > ptBin[ipt + 1] && ipt < ptBins_ - 1)
1018           ipt++;
1019         while (recoJets[ijet].eta() > etaRanges->at(ieta + 1) && ieta < (int)(rho->size() - 1))
1020           ieta++;
1021         mSubtractedEFrac[ipt][ieta]->Fill((double)recoJets[ijet].pileup() / (double)recoJets[ijet].energy());
1022         mSubtractedE[ipt][ieta]->Fill(recoJets[ijet].pileup());
1023 
1024         for (unsigned irho = 0; irho < rho->size(); irho++) {
1025           mRhoDist_vsEta->Fill(recoJets[ijet].eta(), rho->at(irho));
1026           mRhoMDist_vsEta->Fill(recoJets[ijet].eta(), rhom->at(irho));
1027           mRhoDist_vsPt->Fill(recoJets[ijet].pt(), rho->at(irho));
1028           mRhoMDist_vsPt->Fill(recoJets[ijet].pt(), rhom->at(irho));
1029           mRhoDist_vsCent[ieta]->Fill(HF_energy, rho->at(irho));
1030           mRhoMDist_vsCent[ieta]->Fill(HF_energy, rhom->at(irho));
1031         }
1032       }
1033 
1034       for (size_t iii = 0; iii < numbers.size(); iii++) {
1035         pfDeltaR = sqrt((numbers[iii][2] - recoJets[ijet].phi()) * (numbers[iii][2] - recoJets[ijet].phi()) +
1036                         (numbers[iii][1] - recoJets[ijet].eta()) * (numbers[iii][1] - recoJets[ijet].eta()));  //MZ
1037 
1038         mPFDeltaR->Fill(pfDeltaR);                                  //MZ
1039         mPFDeltaR_Scaled_R->Fill(pfDeltaR, 1. / pow(pfDeltaR, 2));  //MZ
1040       }
1041     }
1042   }
1043   if (mNJets_40)
1044     mNJets_40->Fill(nJet_40);
1045 
1046   numbers.clear();
1047 }
[ 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!