diff --git a/analyses/pluginCMS/CMS_2017_I1499471.info b/analyses/pluginCMS/CMS_2017_I1499471.info --- a/analyses/pluginCMS/CMS_2017_I1499471.info +++ b/analyses/pluginCMS/CMS_2017_I1499471.info @@ -1,55 +1,55 @@ Name: CMS_2017_I1499471 Year: 2017 Summary: Measurements of the associated production of a Z boson and b jets in pp collisions at $\sqrt{s} = 8$ TeV Experiment: CMS Collider: LHC -InspireID: I1499471 +InspireID: 1499471 Status: VALIDATED Authors: - Fabio Cossutti (fabio.cossutti@ts.infn.it) References: - DOI:10.1140/epjc/s10052-017-5140-y - arXiv:1611.06507 - CMS-SMP-14-010 RunInfo: 'Run MC generators with Z decaying leptonically into both electrons and muons at 8 TeV CoM energy. For jet-flavour inclusive samples, order of 20 million unweighted events can give a reasonable global comparison, but precision in the high jet multiplicity region/high jet pt may require substantially larger samples or statistical enhancement of high jet multiplicities.' NumEvents: 20000000 NeedCrossSection: yes Beams: [p+, p+] Energies: [8000] Description: 'Measurements of the associated production of a Z boson with at least one jet originating from a b quark in proton-proton collisions at $\sqrt{s} = 8$ TeV are presented. Differential cross sections are measured with data collected by the CMS experiment corresponding to an integrated luminosity of 19.8 $\text{fb}^{-1}$. Z bosons are reconstructed through their decays to electrons and muons. Cross sections are measured as a function of observables characterizing the kinematics of the b jet and the Z boson. Ratios of differential cross sections for the associated production with at least one b jet to the associated production with any jet are also presented. The production of a Z boson with at least two b jets is investigated, and differential cross sections are measured for the dijet system. Results are compared to theoretical predictions, testing two different flavour schemes for the choice of initial-state partons. Note of the authors: the Z boson selection and b-tagging algorithms used in the code are based on standard Rivet code and do not exactly coincide with those used in the paper. They anyway reproduce the paper results at the percent level or better for most of the bins.' BibKey: Khachatryan:2016iob BibTeX: '@article{Khachatryan:2016iob, author = "Khachatryan, Vardan and others", title = "{Measurements of the associated production of a Z boson and b jets in pp collisions at $\sqrt{s}$ = 8 TeV}", collaboration = "CMS", journal = "Submitted to: Eur. Phys. J. C", year = "2016", eprint = "1611.06507", archivePrefix = "arXiv", primaryClass = "hep-ex", reportNumber = "CMS-SMP-14-010, CERN-EP-2016-254", SLACcitation = "%%CITATION = ARXIV:1611.06507;%%" }' diff --git a/analyses/pluginCMS/CMS_2018_I1690148.cc b/analyses/pluginCMS/CMS_2018_I1690148.cc --- a/analyses/pluginCMS/CMS_2018_I1690148.cc +++ b/analyses/pluginCMS/CMS_2018_I1690148.cc @@ -1,439 +1,453 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/ChargedLeptons.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/InvMassFinalState.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Math/MatrixN.hh" -#include "Rivet/Math/MatrixDiag.hh" #include "Rivet/Tools/fjcontrib/Nsubjettiness.hh" #include "Rivet/Tools/fjcontrib/EnergyCorrelator.hh" namespace Rivet { /// Measurement of jet substructure observables in $t\bar{t}$ events from $pp$ collisions at 13~TeV class CMS_2018_I1690148 : public Analysis { public: enum Reconstruction { CHARGED=0, ALL=1 }; enum Observable { MULT=0, PTDS=1, GA_LHA=2, GA_WIDTH=3, GA_THRUST=4, ECC=5, ZG=6, ZGDR=7, NSD=8, TAU21=9, TAU32=10, TAU43=11, C1_00=12, C1_02=13, C1_05=14, C1_10=15, C1_20=16, C2_00=17, C2_02=18, C2_05=19, C2_10=20, C2_20=21, C3_00=22, C3_02=23, C3_05=24, C3_10=25, C3_20=26, M2_B1=27, N2_B1=28, N3_B1=29, M2_B2=30, N2_B2=31, N3_B2=32 }; enum Flavor { INCL=0, BOTTOM=1, QUARK=2, GLUON=3 }; /// Minimal constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2018_I1690148); /// @name Analysis methods //@{ /// Set up projections and book histograms void init() { // Cuts particle_cut = (Cuts::abseta < 5.0 && Cuts::pT > 0.*GeV); lepton_cut = (Cuts::abseta < 2.4 && Cuts::pT > 15.*GeV); jet_cut = (Cuts::abseta < 2.4 && Cuts::pT > 30.*GeV); // Generic final state FinalState fs(particle_cut); // Dressed leptons ChargedLeptons charged_leptons(fs); IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); PromptFinalState prompt_leptons(charged_leptons); prompt_leptons.acceptMuonDecays(true); prompt_leptons.acceptTauDecays(true); PromptFinalState prompt_photons(photons); prompt_photons.acceptMuonDecays(true); prompt_photons.acceptTauDecays(true); // NB. useDecayPhotons=true allows for photons with tau ancestor; photons from hadrons are vetoed by the PromptFinalState; DressedLeptons dressed_leptons(prompt_photons, prompt_leptons, 0.1, lepton_cut, true); addProjection(dressed_leptons, "DressedLeptons"); // Projection for jets VetoedFinalState fsForJets(fs); fsForJets.addVetoOnThisFinalState(dressed_leptons); addProjection(FastJets(fsForJets, FastJets::ANTIKT, 0.4, JetAlg::ALL_MUONS, JetAlg::NO_INVISIBLES), "Jets"); // Booking of histograms int d = 0; for (int r = 0; r < 2; ++r) { // reconstruction (charged, all) for (int o = 0; o < 33; ++o) { // observable d += 1; for (int f = 0; f < 4; ++f) { // flavor - char buffer [11]; + char buffer [20]; sprintf(buffer, "d%02d-x01-y%02d", d, f+1); _h[r][o][f] = bookHisto1D(buffer); } } } } void analyze(const Event& event) { const double weight = event.weight(); // select ttbar -> lepton+jets const vector& leptons = applyProjection(event, "DressedLeptons").dressedLeptons(); int nsel_leptons = 0; for (const DressedLepton& lepton : leptons) { if (lepton.pt() > 26.) nsel_leptons += 1; else vetoEvent; // found veto lepton } if (nsel_leptons != 1) vetoEvent; const Jets all_jets = applyProjection(event, "Jets").jetsByPt(jet_cut); if (all_jets.size() < 4) vetoEvent; // categorize jets int nsel_bjets = 0; int nsel_wjets = 0; Jets jets[4]; for (const Jet& jet : all_jets) { // check for jet-lepton overlap -> do not consider for selection if (deltaR(jet, leptons[0]) < 0.4) continue; bool overlap = false; bool w_jet = false; for (const Jet& jet2 : all_jets) { if (jet.momentum() == jet2.momentum()) continue; // check for jet-jet overlap -> do not consider for analysis if (deltaR(jet, jet2) < 0.8) overlap = true; // check for W candidate if (jet.bTagged() or jet2.bTagged()) continue; FourMomentum w_cand = jet.momentum() + jet2.momentum(); if (abs(w_cand.mass() - 80.4) < 15.) w_jet = true; } // count jets for event selection if (jet.bTagged()) nsel_bjets += 1; if (w_jet) nsel_wjets += 1; // jets for analysis if (jet.abseta() > 2. or overlap) continue; if (jet.bTagged()) { jets[BOTTOM].push_back(jet); } else if (w_jet) { jets[QUARK].push_back(jet); } else { jets[GLUON].push_back(jet); } } if (nsel_bjets != 2) vetoEvent; if (nsel_wjets < 2) vetoEvent; // substructure analysis // no loop over incl jets -> more loc but faster for (int f = 1; f < 4; ++f) { for (const Jet& jet : jets[f]) { // apply cuts on constituents vector particles[2]; for (const Particle& p : jet.particles(Cuts::pT > 1.*GeV)) { particles[ALL].push_back( PseudoJet(p.px(), p.py(), p.pz(), p.energy()) ); if (p.charge3() != 0) particles[CHARGED].push_back( PseudoJet(p.px(), p.py(), p.pz(), p.energy()) ); } if (particles[CHARGED].size() == 0) continue; // recluster with C/A and anti-kt+WTA PseudoJet ca_jet[2]; JetDefinition ca_def(fastjet::cambridge_algorithm, fastjet::JetDefinition::max_allowable_R); ClusterSequence ca_charged(particles[CHARGED], ca_def); ClusterSequence ca_all(particles[ALL], ca_def); ca_jet[CHARGED] = ca_charged.exclusive_jets(1)[0]; ca_jet[ALL] = ca_all.exclusive_jets(1)[0]; PseudoJet akwta_jet[2]; JetDefinition akwta_def(fastjet::antikt_algorithm, fastjet::JetDefinition::max_allowable_R, fastjet::RecombinationScheme::WTA_pt_scheme); ClusterSequence akwta_charged(particles[CHARGED], akwta_def); ClusterSequence akwta_all(particles[ALL], akwta_def); akwta_jet[CHARGED] = akwta_charged.exclusive_jets(1)[0]; akwta_jet[ALL] = akwta_all.exclusive_jets(1)[0]; // calculate observables for (int r = 0; r < 2; ++r) { int mult = akwta_jet[r].constituents().size(); // generalized angularities _h[r][MULT][INCL]->fill(mult, weight); _h[r][MULT][f]->fill(mult, weight); if (mult > 1) { double ptds = getPtDs(akwta_jet[r]); double ga_lha = calcGA(0.5, 1., akwta_jet[r]); double ga_width = calcGA(1., 1., akwta_jet[r]); double ga_thrust = calcGA(2., 1., akwta_jet[r]); _h[r][PTDS][INCL]->fill(ptds, weight); _h[r][PTDS][f]->fill(ptds, weight); _h[r][GA_LHA][INCL]->fill(ga_lha, weight); _h[r][GA_LHA][f]->fill(ga_lha, weight); _h[r][GA_WIDTH][INCL]->fill(ga_width, weight); _h[r][GA_WIDTH][f]->fill(ga_width, weight); _h[r][GA_THRUST][INCL]->fill(ga_thrust, weight); _h[r][GA_THRUST][f]->fill(ga_thrust, weight); } // eccentricity if (mult > 3) { double ecc = getEcc(akwta_jet[r]); _h[r][ECC][INCL]->fill(ecc, weight); _h[r][ECC][f]->fill(ecc, weight); } // N-subjettiness if (mult > 2) { double tau21 = getTau(2, 1, ca_jet[r]); _h[r][TAU21][INCL]->fill(tau21, weight); _h[r][TAU21][f]->fill(tau21, weight); } if (mult > 3) { double tau32 = getTau(3, 2, ca_jet[r]); _h[r][TAU32][INCL]->fill(tau32, weight); _h[r][TAU32][f]->fill(tau32, weight); } if (mult > 4) { double tau43 = getTau(4, 3, ca_jet[r]); _h[r][TAU43][INCL]->fill(tau43, weight); _h[r][TAU43][f]->fill(tau43, weight); } // soft drop if (mult > 1) { vector sd_results = getZg(ca_jet[r]); if (sd_results[0] > 0.) { _h[r][ZG][INCL]->fill(sd_results[0], weight); _h[r][ZG][f]->fill(sd_results[0], weight); _h[r][ZGDR][INCL]->fill(sd_results[1], weight); _h[r][ZGDR][f]->fill(sd_results[1], weight); } } int nsd = getNSD(0.007, -1., ca_jet[r]); _h[r][NSD][INCL]->fill(nsd, weight); _h[r][NSD][f]->fill(nsd, weight); // C-series energy correlation ratios if (mult > 1) { double cn_00 = getC(1, 0.0, ca_jet[r]); double cn_02 = getC(1, 0.2, ca_jet[r]); double cn_05 = getC(1, 0.5, ca_jet[r]); double cn_10 = getC(1, 1.0, ca_jet[r]); double cn_20 = getC(1, 2.0, ca_jet[r]); _h[r][C1_00][INCL]->fill(cn_00, weight); _h[r][C1_02][INCL]->fill(cn_02, weight); _h[r][C1_05][INCL]->fill(cn_05, weight); _h[r][C1_10][INCL]->fill(cn_10, weight); _h[r][C1_20][INCL]->fill(cn_20, weight); _h[r][C1_00][f]->fill(cn_00, weight); _h[r][C1_02][f]->fill(cn_02, weight); _h[r][C1_05][f]->fill(cn_05, weight); _h[r][C1_10][f]->fill(cn_10, weight); _h[r][C1_20][f]->fill(cn_20, weight); } if (mult > 2) { double cn_00 = getC(2, 0.0, ca_jet[r]); double cn_02 = getC(2, 0.2, ca_jet[r]); double cn_05 = getC(2, 0.5, ca_jet[r]); double cn_10 = getC(2, 1.0, ca_jet[r]); double cn_20 = getC(2, 2.0, ca_jet[r]); _h[r][C2_00][INCL]->fill(cn_00, weight); _h[r][C2_02][INCL]->fill(cn_02, weight); _h[r][C2_05][INCL]->fill(cn_05, weight); _h[r][C2_10][INCL]->fill(cn_10, weight); _h[r][C2_20][INCL]->fill(cn_20, weight); _h[r][C2_00][f]->fill(cn_00, weight); _h[r][C2_02][f]->fill(cn_02, weight); _h[r][C2_05][f]->fill(cn_05, weight); _h[r][C2_10][f]->fill(cn_10, weight); _h[r][C2_20][f]->fill(cn_20, weight); } if (mult > 3) { double cn_00 = getC(3, 0.0, ca_jet[r]); double cn_02 = getC(3, 0.2, ca_jet[r]); double cn_05 = getC(3, 0.5, ca_jet[r]); double cn_10 = getC(3, 1.0, ca_jet[r]); double cn_20 = getC(3, 2.0, ca_jet[r]); _h[r][C3_00][INCL]->fill(cn_00, weight); _h[r][C3_02][INCL]->fill(cn_02, weight); _h[r][C3_05][INCL]->fill(cn_05, weight); _h[r][C3_10][INCL]->fill(cn_10, weight); _h[r][C3_20][INCL]->fill(cn_20, weight); _h[r][C3_00][f]->fill(cn_00, weight); _h[r][C3_02][f]->fill(cn_02, weight); _h[r][C3_05][f]->fill(cn_05, weight); _h[r][C3_10][f]->fill(cn_10, weight); _h[r][C3_20][f]->fill(cn_20, weight); } // M/N-series energy correlation ratios if (mult > 2) { double m2_b1 = getM(2, 1., ca_jet[r]); double m2_b2 = getM(2, 2., ca_jet[r]); double n2_b1 = getN(2, 1., ca_jet[r]); double n2_b2 = getN(2, 2., ca_jet[r]); _h[r][M2_B1][INCL]->fill(m2_b1, weight); _h[r][M2_B2][INCL]->fill(m2_b2, weight); _h[r][N2_B1][INCL]->fill(n2_b1, weight); _h[r][N2_B2][INCL]->fill(n2_b2, weight); _h[r][M2_B1][f]->fill(m2_b1, weight); _h[r][M2_B2][f]->fill(m2_b2, weight); _h[r][N2_B1][f]->fill(n2_b1, weight); _h[r][N2_B2][f]->fill(n2_b2, weight); } if (mult > 3) { double n3_b1 = getN(3, 1., ca_jet[r]); double n3_b2 = getN(3, 2., ca_jet[r]); _h[r][N3_B1][INCL]->fill(n3_b1, weight); _h[r][N3_B2][INCL]->fill(n3_b2, weight); _h[r][N3_B1][f]->fill(n3_b1, weight); _h[r][N3_B2][f]->fill(n3_b2, weight); } } } } } void finalize() { for (int r = 0; r < 2; ++r) { // reconstruction (charged, all) for (int o = 0; o < 33; ++o) { // observable for (int f = 0; f < 4; ++f) { // flavor normalize(_h[r][o][f], 1.0, false); } } } } //@} private: double deltaR(PseudoJet j1, PseudoJet j2) { double deta = j1.eta() - j2.eta(); double dphi = j1.delta_phi_to(j2); return sqrt(deta*deta + dphi*dphi); } double getPtDs(PseudoJet jet) { double mult = jet.constituents().size(); double sumpt = 0.; // would be jet.pt() in WTA scheme but better keep it generic double sumpt2 = 0.; for (auto p : jet.constituents()) { sumpt += p.pt(); sumpt2 += pow(p.pt(), 2); } double ptd = sumpt2/pow(sumpt,2); return max(0., sqrt((ptd-1./mult) * mult/(mult-1.))); } double calcGA(double beta, double kappa, PseudoJet jet) { double sumpt = 0.; for (const auto& p : jet.constituents()) { sumpt += p.pt(); } double ga = 0.; for (auto p : jet.constituents()) { ga += pow(p.pt()/sumpt, kappa) * pow(deltaR(jet, p)/0.4, beta); } return ga; } double getEcc(PseudoJet jet) { // Covariance matrix Matrix<2> M; for (const auto& p : jet.constituents()) { Matrix<2> MPart; MPart.set(0, 0, (p.eta() - jet.eta()) * (p.eta() - jet.eta())); MPart.set(0, 1, (p.eta() - jet.eta()) * mapAngleMPiToPi(p.phi() - jet.phi())); MPart.set(1, 0, mapAngleMPiToPi(p.phi() - jet.phi()) * (p.eta() - jet.eta())); MPart.set(1, 1, mapAngleMPiToPi(p.phi() - jet.phi()) * mapAngleMPiToPi(p.phi() - jet.phi())); M += MPart * p.e(); } // Calculate eccentricity from eigenvalues - const EigenSystem<2> eigen = diagonalize(M); - return 1. - eigen.getEigenValues()[1]/eigen.getEigenValues()[0]; + // Check that the matrix is symmetric. + const bool isSymm = M.isSymm(); + if (!isSymm) { + MSG_ERROR("Error: energy tensor not symmetric:"); + MSG_ERROR("[0,1] vs. [1,0]: " << M.get(0,1) << ", " << M.get(1,0)); + } + // If not symmetric, something's wrong (we made sure the error msg appeared first). + + assert(isSymm); + const double a = M.get(0,0); + const double b = M.get(1,1); + const double c = M.get(1,0); + + const double l1 = 0.5*(a+b+sqrt( (a-b)*(a-b) + 4 *c*c)); + const double l2 = 0.5*(a+b-sqrt( (a-b)*(a-b) + 4 *c*c)); + + return 1. - l2/l1; } double getTau(int N, int M, PseudoJet jet) { fjcontrib::Nsubjettiness::NsubjettinessRatio tau_ratio(N, M, fjcontrib::Nsubjettiness::OnePass_WTA_KT_Axes(), fjcontrib::Nsubjettiness::NormalizedMeasure(1.0, 0.4)); return tau_ratio(jet); } vector getZg(PseudoJet jet) { PseudoJet jet0 = jet; PseudoJet jet1, jet2; double zg = 0.; while (zg < 0.1 && jet0.has_parents(jet1, jet2)) { zg = jet2.pt()/jet0.pt(); jet0 = jet1; } if (zg < 0.1) return {-1., -1.}; return {zg, jet1.delta_R(jet2)}; } int getNSD(double zcut, double beta, PseudoJet jet) { PseudoJet jet0 = jet; PseudoJet jet1, jet2; int nsd = 0.; double zg = 0.; while (jet0.has_parents(jet1, jet2)) { zg = jet2.pt()/jet0.pt(); if (zg > zcut * pow(jet1.delta_R(jet2)/0.4, beta)) nsd += 1; jet0 = jet1; } return nsd; } double getC(int N, double beta, PseudoJet jet) { fjcontrib::EnergyCorrelatorDoubleRatio C(N, beta); return C(jet); } double getM(int N, double beta, PseudoJet jet) { fjcontrib::EnergyCorrelatorMseries CM(N, beta); return CM(jet); } double getN(int N, double beta, PseudoJet jet) { fjcontrib::EnergyCorrelatorNseries CN(N, beta); return CN(jet); } private: // @name Histogram data members //@{ Cut particle_cut, lepton_cut, jet_cut; Histo1DPtr _h[2][33][4]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2018_I1690148); }