1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/DileptonFinder.hh"
  8#include "fastjet/JetDefinition.hh"
  9#include "fastjet/ClusterSequence.hh"
 10#include "fastjet/tools/Recluster.hh"
 11#include "fastjet/contrib/SoftDrop.hh"
 12#include <algorithm>
 13
 14namespace Rivet {
 15
 16  // This analysis uses all the FJ stuff...
 17  using namespace fastjet;
 18
 19
 20  /// Routine for QG substructure analysis
 21  class CMS_2021_I1920187 : public Analysis {
 22  public:
 23
 24
 25    /// Constructor
 26    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2021_I1920187);
 27
 28
 29    /// Book histograms and initialise projections before the run
 30    void init() {
 31
 32      _mode = 0;
 33      if ( getOption("MODE") == "DIJET" ) _mode = 1;
 34      else if ( getOption("MODE") == "ZJET" ) _mode = 2;
 35      else {
 36        MSG_WARNING("Mode not specified in CMS_2021_I1920187, using DIJET");
 37        _mode = 1;
 38      }
 39
 40      // Initialise and register projections
 41      FinalState fs(Cuts::abseta < 5);
 42      // Z-jet
 43      if (_mode == 2) {
 44        // for the muons
 45        double mu_pt = 26.;
 46        double mz_min = (90-20);
 47        double mz_max = (90+20);
 48        double eta_max = 2.4;
 49        DileptonFinder zfinder(fs, 91.2*GeV, 0.1, Cuts::pT > mu_pt*GeV  && Cuts::abseta < eta_max &&
 50                                   Cuts::abspid == PID::MUON, Cuts::massIn(mz_min*GeV, mz_max*GeV));
 51        declare(zfinder, "DileptonFinder");
 52
 53        eta_max = 2.4;
 54        FinalState fs_muons(Cuts::abseta < eta_max);
 55        IdentifiedFinalState muons_noCut(fs_muons, {PID::MUON, PID::ANTIMUON});
 56        declare(muons_noCut, "MUONS_NOCUT");
 57        // Particles for the jets
 58        VetoedFinalState jet_input(fs);
 59        jet_input.vetoNeutrinos();
 60        jet_input.addVetoOnThisFinalState(getProjection<DileptonFinder>("DileptonFinder"));
 61        declare(jet_input, "JET_INPUT");
 62        _ptBinsGen = { 50, 65, 88, 120, 150, 186, 254, 326, 408, 1500};
 63      }
 64      // dijet
 65      else {
 66        // Particles for the jets
 67        VetoedFinalState jet_input(fs);
 68        jet_input.vetoNeutrinos();
 69        declare(jet_input, "JET_INPUT");
 70        _ptBinsGen = {50, 65, 88, 120, 150, 186, 254, 326, 408, 481, 614, 800, 1000, 4000};
 71      }
 72      // Book histograms
 73      // resize vectors appropriately
 74      size_t nHistsRadii = _jetRadii.size();
 75      size_t nHistsLambda = _lambdaVars.size();
 76      size_t nHistsPt = _ptBinsGen.size()-1;
 77
 78      // Z jet
 79      if (_mode == 2) {
 80        _h_zpj.resize(nHistsRadii, vector<vector<Histo1DPtr> >(nHistsLambda, vector<Histo1DPtr>(nHistsPt)));
 81        _h_zpj_groomed.resize(nHistsRadii, vector<vector<Histo1DPtr> >(nHistsLambda, vector<Histo1DPtr>(nHistsPt)));
 82
 83        // Now book histos
 84        // remember 1-indexed
 85
 86        // yoda plot naming scheme
 87        // --------------------------------------------------------------------------
 88        // channel (ak4/8 radii [000, 100] + zjet / dijet cen / fwd [00, 10, 20] + groomed versions [0, 1])
 89        // lambda variable; neutral+charged & charged-only are treated separately [0000,..,4000]
 90        // pT bin [00000,..,120000]
 91        for (size_t radiusInd=0; radiusInd < _jetRadii.size(); radiusInd++) {
 92          for (size_t lambdaInd=0; lambdaInd < _lambdaVars.size(); lambdaInd++) {
 93            for (size_t ptInd=0; ptInd < _ptBinsGen.size()-1; ptInd++) {
 94              book(_h_zpj[radiusInd][lambdaInd][ptInd], _hepdata_index.at(0+100*radiusInd+1000*lambdaInd+10000*ptInd), 1, 1);
 95              book(_h_zpj_groomed[radiusInd][lambdaInd][ptInd], _hepdata_index.at(10+100*radiusInd+1000*lambdaInd+10000*ptInd), 1, 1);
 96            }
 97          }
 98        }
 99      }
100      // di jet
101      else {
102        _h_dijet_cen.resize(nHistsRadii, vector<vector<Histo1DPtr> >(nHistsLambda, vector<Histo1DPtr>(nHistsPt)));
103        _h_dijet_cen_groomed.resize(nHistsRadii, vector<vector<Histo1DPtr> >(nHistsLambda, vector<Histo1DPtr>(nHistsPt)));
104        _h_dijet_fwd.resize(nHistsRadii, vector<vector<Histo1DPtr> >(nHistsLambda, vector<Histo1DPtr>(nHistsPt)));
105        _h_dijet_fwd_groomed.resize(nHistsRadii, vector<vector<Histo1DPtr> >(nHistsLambda, vector<Histo1DPtr>(nHistsPt)));
106
107        // Now book histos
108        // remember 1-indexed
109
110        // yoda plot naming scheme
111        // --------------------------------------------------------------------------
112        // channel (ak4/8 radii [000, 100] + zjet / dijet cen / fwd [00, 10, 20] + groomed versions [0, 1])
113        // lambda variable; neutral+charged & charged-only are treated separately [0000,..,4000]
114        // pT bin [00000,..,120000]
115        for (size_t radiusInd=0; radiusInd < _jetRadii.size(); radiusInd++) {
116          for (size_t lambdaInd=0; lambdaInd < _lambdaVars.size(); lambdaInd++) {
117            for (size_t ptInd=0; ptInd < _ptBinsGen.size()-1; ptInd++) {
118              book(_h_dijet_cen[radiusInd][lambdaInd][ptInd], _hepdata_index.at(1+100*radiusInd+1000*lambdaInd+10000*ptInd), 1, 1);
119              book(_h_dijet_cen_groomed[radiusInd][lambdaInd][ptInd], _hepdata_index.at(11+100*radiusInd+1000*lambdaInd+10000*ptInd), 1, 1);
120              book(_h_dijet_fwd[radiusInd][lambdaInd][ptInd], _hepdata_index.at(2+100*radiusInd+1000*lambdaInd+10000*ptInd), 1, 1);
121              book(_h_dijet_fwd_groomed[radiusInd][lambdaInd][ptInd], _hepdata_index.at(12+100*radiusInd+1000*lambdaInd+10000*ptInd), 1, 1);
122            }
123          }
124        }
125      }
126    }
127
128    // Get index of largest bin smaller than value in vector
129    // e.g. what you'd need when binning a continuous variable
130    size_t getBinIndex(float value, const vector<float>& bins) {
131      auto itr = std::lower_bound(bins.begin(), bins.end(), value);
132      return itr - bins.begin() - 1;
133    }
134
135    /// Perform the per-event analysis
136    void analyze(const Event& event) {
137      // Convert Particles into PseudoJets for clustering
138      const VetoedFinalState& fs = apply<VetoedFinalState>(event, "JET_INPUT");
139      const Particles& fsParticles = fs.particles();
140      vector<PseudoJet> particles;
141      particles.reserve(fsParticles.size());
142      for (size_t iFS=0; iFS<fsParticles.size(); iFS++){
143        PseudoJet p = fsParticles[iFS].pseudojet();
144        p.set_user_index(fsParticles[iFS].isCharged()); // for later reference to charge
145        particles.push_back(p);
146      }
147
148      // z jet
149      if (_mode == 2) {
150        for (size_t radiusInd=0; radiusInd < _jetRadii.size(); radiusInd++) {
151          float jetRadius = _jetRadii.at(radiusInd);
152
153          JetDefinition jet_def(antikt_algorithm, jetRadius);
154          vector<PseudoJet> jets = (SelectorPtMin(15))(jet_def(particles));
155
156          const FinalState& muons = apply<IdentifiedFinalState>(event, "MUONS_NOCUT");
157          if (muons.size() >= 2) {
158            Particle muon1 = muons.particlesByPt()[0];
159            Particle muon2 = muons.particlesByPt()[1];
160            FourMomentum z = muon1.momentum() + muon2.momentum();
161            if (jets.size() > 0) {
162              PseudoJet jet1 = jets[0];
163            }
164          }
165
166          // Reconstruct Z
167          const DileptonFinder& zfinder = apply<DileptonFinder>(event, "DileptonFinder");
168          if (zfinder.bosons().size() < 1) continue;
169
170          const Particle & z = zfinder.bosons()[0];
171          double zpt = z.pt();
172
173          // Now do selection criteria
174          bool passZpJ = false;
175          if (jets.size() < 1) continue;
176          PseudoJet jet1 = jets[0];
177          float jet1pt = jet1.pt();
178          float asym = fabs((jet1pt - zpt) / (jet1pt+zpt));
179          float dphi = Rivet::deltaPhi(jet1.phi(), z.phi());
180          passZpJ = ((fabs(jet1.rapidity()) < 1.7) && (zpt > 30) && (asym < 0.3) && (dphi > 2.0));
181
182          if (!passZpJ) continue;
183
184          // Now calculate lambda variables and fill hists
185
186          // Simplify life - ignore this jet if it is below 1st hist pt range
187          // Note that we don't apply it to the original jet pt cut - since
188          // we have phase space where one jet is > 50, and one < 50
189          if (jet1pt < _ptBinsGen[0]) continue;
190          // ignore jet if beyond the last bin
191          if (jet1pt > _ptBinsGen.back()) continue;
192
193          // Need to use original, ungroomed jet pT to bin
194          size_t ptBinInd = getBinIndex(jet1pt, _ptBinsGen);
195
196          // UNGROOMED VERSION
197          // -------------------------------------------------------------------
198          vector<PseudoJet> chargedParticles;
199          for (size_t iC=0; iC<jet1.constituents().size(); iC++)
200            if (jet1.constituents()[iC].user_index())
201              chargedParticles.push_back(jet1.constituents()[iC]);
202          vector<PseudoJet> chargedJets = jet_def(chargedParticles);
203
204          // Fill hists for each lambda variable
205          for (size_t lambdaInd=0; lambdaInd < _lambdaVars.size(); lambdaInd++) {
206            const LambdaVar & thisLambdaVar = _lambdaVars[lambdaInd];
207            Angularity angularity(thisLambdaVar.beta, jetRadius, thisLambdaVar.kappa, thisLambdaVar.constitCut);
208            float val = -1;
209            if (thisLambdaVar.isCharged)
210              val = (chargedJets.size()>0) ? angularity(chargedJets[0]) : -1;
211            else
212              val = angularity(jet1);
213            if (val<0) continue;
214            _h_zpj[radiusInd][lambdaInd][ptBinInd]->fill(val);
215          }
216
217          // GROOMED VERSION
218          // -------------------------------------------------------------------
219          // Get groomed jet
220          fastjet::contrib::SoftDrop sd(0, 0.1, jetRadius);
221          PseudoJet groomedJet = sd(jet1);
222          PseudoJet groomedJetCharged;
223          if (chargedJets.size()>0)
224            groomedJetCharged= sd(chargedJets[0]);
225
226          // Fill hists for each lambda variable
227          for (size_t lambdaInd=0; lambdaInd < _lambdaVars.size(); lambdaInd++) {
228            const LambdaVar & thisLambdaVar = _lambdaVars[lambdaInd];
229            Angularity angularity(thisLambdaVar.beta, jetRadius, thisLambdaVar.kappa, thisLambdaVar.constitCut);
230            float val = -1;
231            if (thisLambdaVar.isCharged)
232              val = (chargedJets.size()>0) ? angularity(groomedJetCharged) : -1;
233            else
234              val = angularity(groomedJet);
235            if (val<0) continue;
236            _h_zpj_groomed[radiusInd][lambdaInd][ptBinInd]->fill(val);
237          }
238        } // end loop over jet radii
239      }
240      // di jet
241      else {
242        for (size_t radiusInd=0; radiusInd < _jetRadii.size(); radiusInd++) {
243          float jetRadius = _jetRadii.at(radiusInd);
244
245          JetDefinition jet_def(antikt_algorithm, jetRadius);
246          vector<PseudoJet> jets = (SelectorNHardest(2) * SelectorPtMin(15))(jet_def(particles));
247
248          bool passDijet = false;
249          if (jets.size() < 2) continue;
250          const auto & jet1 = jets.at(0);
251          const auto & jet2 = jets.at(1);
252          float jet1pt = jet1.pt();
253          float jet2pt = jet2.pt();
254          float asym = (jet1pt - jet2pt) / (jet1pt+jet2pt);
255          float dphi = Rivet::deltaPhi(jet1.phi(), jet2.phi());
256          passDijet = ((fabs(jet1.rapidity()) < 1.7) && (fabs(jet2.rapidity()) < 1.7) && (asym < 0.3) && (dphi > 2.0));
257
258          if (!passDijet) continue;
259
260          // Sort by increasing absolute rapidity
261          vector<PseudoJet> dijets = {jet1, jet2};
262          std::sort(dijets.begin(), dijets.end(),
263              [] (const PseudoJet & A, const PseudoJet & B)
264              { return fabs(A.rapidity()) < fabs(B.rapidity()); }
265              );
266
267          for (size_t iJ=0; iJ<dijets.size(); iJ++) {
268            bool isCentral = (iJ == 0);
269            PseudoJet & jetItr = dijets[iJ];
270
271            // Simplify life - ignore this jet if it is below 1st hist pt range
272            // Note that we don't apply it to the original jet pt cut - since
273            // we have phase space where one jet is > 50, and one < 50
274            if (jetItr.pt() < _ptBinsGen[0]) continue;
275            // ignore jet if beyond the last bin
276            if (jetItr.pt() > _ptBinsGen.back()) continue;
277
278            // Need to use original, ungroomed jet pT to bin
279            size_t ptBinInd = getBinIndex(jetItr.pt(), _ptBinsGen);
280
281            // UNGROOMED VERSION
282            // -------------------------------------------------------------------
283            vector<PseudoJet> chargedParticles;
284            for (size_t iC=0; iC<jetItr.constituents().size(); iC++)
285              if (jetItr.constituents()[iC].user_index())
286          chargedParticles.push_back(jetItr.constituents()[iC]);
287            vector<PseudoJet> chargedJets = jet_def(chargedParticles);
288
289            // Fill hists for each lambda variable
290            for (size_t lambdaInd=0; lambdaInd < _lambdaVars.size(); lambdaInd++) {
291              const LambdaVar & thisLambdaVar = _lambdaVars[lambdaInd];
292              Angularity angularity(thisLambdaVar.beta, jetRadius, thisLambdaVar.kappa, thisLambdaVar.constitCut);
293              float val = -1;
294              if (thisLambdaVar.isCharged)
295                val = (chargedJets.size()>0) ? angularity(chargedJets[0]) : -1;
296              else
297                val = angularity(jetItr);
298              if (val<0) continue;
299
300              if (isCentral) {
301                _h_dijet_cen[radiusInd][lambdaInd][ptBinInd]->fill(val);
302              } else {
303                _h_dijet_fwd[radiusInd][lambdaInd][ptBinInd]->fill(val);
304              }
305            }
306
307            // GROOMED VERSION
308            // -------------------------------------------------------------------
309            // Get groomed jet
310            fastjet::contrib::SoftDrop sd(0, 0.1, jetRadius);
311            PseudoJet groomedJet = sd(jetItr);
312            PseudoJet groomedJetCharged;
313            if (chargedJets.size()>0)
314              groomedJetCharged= sd(chargedJets[0]);
315
316            // Fill hists for each lambda variable
317            for (size_t lambdaInd=0; lambdaInd < _lambdaVars.size(); lambdaInd++) {
318              const LambdaVar & thisLambdaVar = _lambdaVars[lambdaInd];
319              Angularity angularity(thisLambdaVar.beta, jetRadius, thisLambdaVar.kappa, thisLambdaVar.constitCut);
320              float val = -1;
321              if (thisLambdaVar.isCharged)
322                val = (chargedJets.size()>0) ? angularity(groomedJetCharged) : -1;
323              else
324                val = angularity(groomedJet);
325              if (val<0) continue;
326
327              if (isCentral) {
328                _h_dijet_cen_groomed[radiusInd][lambdaInd][ptBinInd]->fill(val);
329              } else {
330                _h_dijet_fwd_groomed[radiusInd][lambdaInd][ptBinInd]->fill(val);
331              }
332            }
333
334          } // end loop over dijets
335        } // end loop over jet radii
336      }
337    } // end analyze() function
338
339    void finalize() {
340      for (auto& hvec3 : { _h_dijet_cen, _h_dijet_cen_groomed,
341                          _h_dijet_fwd, _h_dijet_fwd_groomed,
342                          _h_zpj, _h_zpj_groomed }) {
343        for (auto& hvec2 : hvec3) {
344          for (auto& hvec : hvec2) {
345            normalize(hvec, 1.0, false);
346          }
347        }
348      }
349    }
350
351    /// @class Angularity
352    /// Definition of angularity
353    ///
354    class Angularity : public fastjet::FunctionOfPseudoJet<double> {
355    public:
356
357      /// ctor
358      Angularity(double alpha, double jet_radius,
359                 double kappa=1.0, Selector constitCut=SelectorPtMin(0.))
360        : _alpha(alpha), _radius(jet_radius), _kappa(kappa), _constitCut(constitCut) {}
361
362      /// computation of the angularity itself
363      double result(const PseudoJet& jet) const{
364        // get the jet constituents
365        vector<PseudoJet> constits = jet.constituents();
366
367        // get the reference axis
368        PseudoJet reference_axis = _get_reference_axis(jet);
369
370        // do the actual coputation
371        double numerator = 0.0, denominator = 0.0;
372        for (const auto &c : constits){
373          if (!_constitCut.pass(c)) continue;
374          double pt = c.pt();
375          // Note: better compute (dist^2)^(alpha/2) to avoid an extra square root
376          numerator   += pow(pt, _kappa) * pow(c.squared_distance(reference_axis), 0.5*_alpha);
377          denominator += pt;
378        }
379        if (denominator == 0) return -1;
380        // the formula is only correct for the the typical angularities which satisfy either kappa==1 or alpha==0.
381        else return numerator/(pow(denominator, _kappa)*pow(_radius, _alpha));
382      }
383
384    protected:
385
386      PseudoJet _get_reference_axis(const PseudoJet& jet) const{
387        if (_alpha>1) return jet;
388
389        Recluster recluster(JetDefinition(antikt_algorithm, JetDefinition::max_allowable_R, WTA_pt_scheme));
390        return recluster(jet);
391      }
392
393      double _alpha, _radius, _kappa;
394      Selector _constitCut;
395    };
396
397
398    /// Lightweight class to hold info about Lambda variable
399    class LambdaVar {
400    public:
401      LambdaVar(const std::string & name_, float kappa_, float beta_, bool isCharged_, Selector constitCut_):
402        name(name_),
403        kappa(kappa_),
404        beta(beta_),
405        isCharged(isCharged_),
406        constitCut(constitCut_)
407      {}
408
409      string name;
410      float kappa;
411      float beta;
412      bool isCharged;
413      Selector constitCut;
414    };
415
416
417    // Order matters here
418    const vector<float> _jetRadii = {0.4, 0.8};
419
420    // This order is important! index in vector used to create YODA plot name
421    // Must match that in extracRivetPlotsDijet.py
422    const vector<LambdaVar> _lambdaVars = {
423      LambdaVar("jet_multiplicity", 0, 0, false, SelectorPtMin(1.)),
424      LambdaVar("jet_pTD", 2, 0, false, SelectorPtMin(0.)),
425      LambdaVar("jet_LHA", 1, 0.5, false, SelectorPtMin(0.)),
426      LambdaVar("jet_width", 1, 1, false, SelectorPtMin(0.)),
427      LambdaVar("jet_thrust", 1, 2, false, SelectorPtMin(0.)),
428      LambdaVar("jet_multiplicity_charged", 0, 0, true, SelectorPtMin(1.)),
429      LambdaVar("jet_pTD_charged", 2, 0, true, SelectorPtMin(0.)),
430      LambdaVar("jet_LHA_charged", 1, 0.5, true, SelectorPtMin(0.)),
431      LambdaVar("jet_width_charged", 1, 1, true, SelectorPtMin(0.)),
432      LambdaVar("jet_thrust_charged", 1, 2, true, SelectorPtMin(0.)),
433    };
434
435    vector<float> _ptBinsGen;
436
437    const map<size_t, size_t> _hepdata_index {
438      {30000,1},{30001,2},{31000,3},{31001,4},{34000,5},{34001,6},
439      {33000,7},{33001,8},{32000,9},{32001,10},{82000,11},
440      {122001,12},{32100,13},{32101,14},{37000,15},{37001,16},
441      {32010,17},{32011,18},{2000,139},{12000,140},{22000,141},
442      {42000,142},{52000,143},{62000,144},{72000,145},{2001,155},
443      {12001,156},{22001,157},{42001,158},{52001,159},{62001,160},
444      {72001,161},{82001,162},{92001,163},{102001,164},{112001,165},
445      {2002,179},{12002,180},{22002,181},{32002,182},{42002,183},
446      {52002,184},{62002,185},{72002,186},{82002,187},{92002,188},
447      {102002,189},{112002,190},{122002,191},{7000,192},{17000,193},
448      {27000,194},{47000,195},{57000,196},{67000,197},{77000,198},
449      {87000,199},{7001,209},{17001,210},{27001,211},{47001,212},
450      {57001,213},{67001,214},{77001,215},{87001,216},{97001,217},
451      {107001,218},{117001,219},{127001,220},{7002,234},{17002,235},
452      {27002,236},{37002,237},{47002,238},{57002,239},{67002,240},
453      {77002,241},{87002,242},{97002,243},{107002,244},{117002,245},
454      {127002,246},{5000,247},{15000,248},{25000,249},{35000,250},
455      {45000,251},{55000,252},{65000,253},{75000,254},{85000,255},
456      {5001,265},{15001,266},{25001,267},{35001,268},{45001,269},
457      {55001,270},{65001,271},{75001,272},{85001,273},{95001,274},
458      {105001,275},{115001,276},{125001,277},{5002,291},{15002,292},
459      {25002,293},{35002,294},{45002,295},{55002,296},{65002,297},
460      {75002,298},{85002,299},{95002,300},{105002,301},{115002,302},
461      {125002,303},{6000,304},{16000,305},{26000,306},{36000,307},
462      {46000,308},{56000,309},{66000,310},{76000,311},{86000,312},
463      {6001,322},{16001,323},{26001,324},{36001,325},{46001,326},
464      {56001,327},{66001,328},{76001,329},{86001,330},{96001,331},
465      {106001,332},{116001,333},{126001,334},{6002,348},{16002,349},
466      {26002,350},{36002,351},{46002,352},{56002,353},{66002,354},
467      {76002,355},{86002,356},{96002,357},{106002,358},{116002,359},
468      {126002,360},{9000,361},{19000,362},{29000,363},{39000,364},
469      {49000,365},{59000,366},{69000,367},{79000,368},{89000,369},
470      {9001,379},{19001,380},{29001,381},{39001,382},{49001,383},
471      {59001,384},{69001,385},{79001,386},{89001,387},{99001,388},
472      {109001,389},{119001,390},{129001,391},{9002,405},{19002,406},
473      {29002,407},{39002,408},{49002,409},{59002,410},{69002,411},
474      {79002,412},{89002,413},{99002,414},{109002,415},{119002,416},
475      {129002,417},{8000,418},{18000,419},{28000,420},{38000,421},
476      {48000,422},{58000,423},{68000,424},{78000,425},{88000,426},
477      {8001,436},{18001,437},{28001,438},{38001,439},{48001,440},
478      {58001,441},{68001,442},{78001,443},{88001,444},{98001,445},
479      {108001,446},{118001,447},{128001,448},{8002,462},{18002,463},
480      {28002,464},{38002,465},{48002,466},{58002,467},{68002,468},
481      {78002,469},{88002,470},{98002,471},{108002,472},{118002,473},
482      {128002,474},{0,475},{10000,476},{20000,477},{40000,478},
483      {50000,479},{60000,480},{70000,481},{80000,482},{1,492},
484      {10001,493},{20001,494},{40001,495},{50001,496},{60001,497},
485      {70001,498},{80001,499},{90001,500},{100001,501},{110001,502},
486      {120001,503},{2,517},{10002,518},{20002,519},{30002,520},
487      {40002,521},{50002,522},{60002,523},{70002,524},{80002,525},
488      {90002,526},{100002,527},{110002,528},{120002,529},{1000,530},
489      {11000,531},{21000,532},{41000,533},{51000,534},{61000,535},
490      {71000,536},{81000,537},{1001,547},{11001,548},{21001,549},
491      {41001,550},{51001,551},{61001,552},{71001,553},{81001,554},
492      {91001,555},{101001,556},{111001,557},{121001,558},{1002,572},
493      {11002,573},{21002,574},{31002,575},{41002,576},{51002,577},
494      {61002,578},{71002,579},{81002,580},{91002,581},{101002,582},
495      {111002,583},{121002,584},{4000,585},{14000,586},{24000,587},
496      {44000,588},{54000,589},{64000,590},{74000,591},{84000,592},
497      {4001,602},{14001,603},{24001,604},{44001,605},{54001,606},
498      {64001,607},{74001,608},{84001,609},{94001,610},{104001,611},
499      {114001,612},{124001,613},{4002,627},{14002,628},{24002,629},
500      {34002,630},{44002,631},{54002,632},{64002,633},{74002,634},
501      {84002,635},{94002,636},{104002,637},{114002,638},{124002,639},
502      {3000,640},{13000,641},{23000,642},{43000,643},{53000,644},
503      {63000,645},{73000,646},{83000,647},{3001,657},{13001,658},
504      {23001,659},{43001,660},{53001,661},{63001,662},{73001,663},
505      {83001,664},{93001,665},{103001,666},{113001,667},{123001,668},
506      {3002,682},{13002,683},{23002,684},{33002,685},{43002,686},
507      {53002,687},{63002,688},{73002,689},{83002,690},{93002,691},
508      {103002,692},{113002,693},{123002,694},{2010,695},{12010,696},
509      {22010,697},{42010,698},{52010,699},{62010,700},{72010,701},
510      {82010,702},{2011,712},{12011,713},{22011,714},{42011,715},
511      {52011,716},{62011,717},{72011,718},{82011,719},{92011,720},
512      {102011,721},{112011,722},{122011,723},{2012,737},{12012,738},
513      {22012,739},{32012,740},{42012,741},{52012,742},{62012,743},
514      {72012,744},{82012,745},{92012,746},{102012,747},{112012,748},
515      {122012,749},{7010,750},{17010,751},{27010,752},{37010,753},
516      {47010,754},{57010,755},{67010,756},{77010,757},{87010,758},
517      {7011,768},{17011,769},{27011,770},{37011,771},{47011,772},
518      {57011,773},{67011,774},{77011,775},{87011,776},{97011,777},
519      {107011,778},{117011,779},{127011,780},{7012,794},{17012,795},
520      {27012,796},{37012,797},{47012,798},{57012,799},{67012,800},
521      {77012,801},{87012,802},{97012,803},{107012,804},{117012,805},
522      {127012,806},{5010,807},{15010,808},{25010,809},{35010,810},
523      {45010,811},{55010,812},{65010,813},{75010,814},{85010,815},
524      {5011,825},{15011,826},{25011,827},{35011,828},{45011,829},
525      {55011,830},{65011,831},{75011,832},{85011,833},{95011,834},
526      {105011,835},{115011,836},{125011,837},{5012,851},{15012,852},
527      {25012,853},{35012,854},{45012,855},{55012,856},{65012,857},
528      {75012,858},{85012,859},{95012,860},{105012,861},{115012,862},
529      {125012,863},{6010,864},{16010,865},{26010,866},{36010,867},
530      {46010,868},{56010,869},{66010,870},{76010,871},{86010,872},
531      {6011,882},{16011,883},{26011,884},{36011,885},{46011,886},
532      {56011,887},{66011,888},{76011,889},{86011,890},{96011,891},
533      {106011,892},{116011,893},{126011,894},{6012,908},{16012,909},
534      {26012,910},{36012,911},{46012,912},{56012,913},{66012,914},
535      {76012,915},{86012,916},{96012,917},{106012,918},{116012,919},
536      {126012,920},{9010,921},{19010,922},{29010,923},{39010,924},
537      {49010,925},{59010,926},{69010,927},{79010,928},{89010,929},
538      {9011,939},{19011,940},{29011,941},{39011,942},{49011,943},
539      {59011,944},{69011,945},{79011,946},{89011,947},{99011,948},
540      {109011,949},{119011,950},{129011,951},{9012,965},{19012,966},
541      {29012,967},{39012,968},{49012,969},{59012,970},{69012,971},
542      {79012,972},{89012,973},{99012,974},{109012,975},{119012,976},
543      {129012,977},{8010,978},{18010,979},{28010,980},{38010,981},
544      {48010,982},{58010,983},{68010,984},{78010,985},{88010,986},
545      {8011,996},{18011,997},{28011,998},{38011,999},{48011,1000},
546      {58011,1001},{68011,1002},{78011,1003},{88011,1004},{98011,1005},
547      {108011,1006},{118011,1007},{128011,1008},{8012,1022},
548      {18012,1023},{28012,1024},{38012,1025},{48012,1026},
549      {58012,1027},{68012,1028},{78012,1029},{88012,1030},
550      {98012,1031},{108012,1032},{118012,1033},{128012,1034},
551      {10,1035},{10010,1036},{20010,1037},{30010,1038},{40010,1039},
552      {50010,1040},{60010,1041},{70010,1042},{80010,1043},{11,1053},
553      {10011,1054},{20011,1055},{30011,1056},{40011,1057},{50011,1058},
554      {60011,1059},{70011,1060},{80011,1061},{90011,1062},{100011,1063},
555      {110011,1064},{120011,1065},{12,1079},{10012,1080},{20012,1081},
556      {30012,1082},{40012,1083},{50012,1084},{60012,1085},{70012,1086},
557      {80012,1087},{90012,1088},{100012,1089},{110012,1090},
558      {120012,1091},{1010,1092},{11010,1093},{21010,1094},{31010,1095},
559      {41010,1096},{51010,1097},{61010,1098},{71010,1099},{81010,1100},
560      {1011,1110},{11011,1111},{21011,1112},{31011,1113},{41011,1114},
561      {51011,1115},{61011,1116},{71011,1117},{81011,1118},{91011,1119},
562      {101011,1120},{111011,1121},{121011,1122},{1012,1136},{11012,1137},
563      {21012,1138},{31012,1139},{41012,1140},{51012,1141},{61012,1142},
564      {71012,1143},{81012,1144},{91012,1145},{101012,1146},{111012,1147},
565      {121012,1148},{4010,1149},{14010,1150},{24010,1151},{34010,1152},
566      {44010,1153},{54010,1154},{64010,1155},{74010,1156},{84010,1157},
567      {4011,1167},{14011,1168},{24011,1169},{34011,1170},{44011,1171},
568      {54011,1172},{64011,1173},{74011,1174},{84011,1175},{94011,1176},
569      {104011,1177},{114011,1178},{124011,1179},{4012,1193},{14012,1194},
570      {24012,1195},{34012,1196},{44012,1197},{54012,1198},{64012,1199},
571      {74012,1200},{84012,1201},{94012,1202},{104012,1203},{114012,1204},
572      {124012,1205},{3010,1206},{13010,1207},{23010,1208},{33010,1209},
573      {43010,1210},{53010,1211},{63010,1212},{73010,1213},{83010,1214},
574      {3011,1224},{13011,1225},{23011,1226},{33011,1227},{43011,1228},
575      {53011,1229},{63011,1230},{73011,1231},{83011,1232},{93011,1233},
576      {103011,1234},{113011,1235},{123011,1236},{3012,1250},{13012,1251},
577      {23012,1252},{33012,1253},{43012,1254},{53012,1255},{63012,1256},
578      {73012,1257},{83012,1258},{93012,1259},{103012,1260},{113012,1261},
579      {123012,1262},{2100,1263},{12100,1264},{22100,1265},{42100,1266},
580      {52100,1267},{62100,1268},{72100,1269},{82100,1270},{2101,1271},
581      {12101,1272},{22101,1273},{42101,1274},{52101,1275},{62101,1276},
582      {72101,1277},{82101,1278},{92101,1279},{102101,1280},{112101,1281},
583      {122101,1282},{2102,1283},{12102,1284},{22102,1285},{32102,1286},
584      {42102,1287},{52102,1288},{62102,1289},{72102,1290},{82102,1291},
585      {92102,1292},{102102,1293},{112102,1294},{122102,1295},{7100,1296},
586      {17100,1297},{27100,1298},{37100,1299},{47100,1300},{57100,1301},
587      {67100,1302},{77100,1303},{87100,1304},{7101,1314},{17101,1315},
588      {27101,1316},{37101,1317},{47101,1318},{57101,1319},{67101,1320},
589      {77101,1321},{87101,1322},{97101,1323},{107101,1324},{117101,1325},
590      {127101,1326},{7102,1327},{17102,1328},{27102,1329},{37102,1330},
591      {47102,1331},{57102,1332},{67102,1333},{77102,1334},{87102,1335},
592      {97102,1336},{107102,1337},{117102,1338},{127102,1339},{5100,1340},
593      {15100,1341},{25100,1342},{35100,1343},{45100,1344},{55100,1345},
594      {65100,1346},{75100,1347},{85100,1348},{5101,1358},{15101,1359},
595      {25101,1360},{35101,1361},{45101,1362},{55101,1363},{65101,1364},
596      {75101,1365},{85101,1366},{95101,1367},{105101,1368},{115101,1369},
597      {125101,1370},{5102,1371},{15102,1372},{25102,1373},{35102,1374},
598      {45102,1375},{55102,1376},{65102,1377},{75102,1378},{85102,1379},
599      {95102,1380},{105102,1381},{115102,1382},{125102,1383},{6100,1384},
600      {16100,1385},{26100,1386},{36100,1387},{46100,1388},{56100,1389},
601      {66100,1390},{76100,1391},{86100,1392},{6101,1402},{16101,1403},
602      {26101,1404},{36101,1405},{46101,1406},{56101,1407},{66101,1408},
603      {76101,1409},{86101,1410},{96101,1411},{106101,1412},{116101,1413},
604      {126101,1414},{6102,1415},{16102,1416},{26102,1417},{36102,1418},
605      {46102,1419},{56102,1420},{66102,1421},{76102,1422},{86102,1423},
606      {96102,1424},{106102,1425},{116102,1426},{126102,1427},{9100,1428},
607      {19100,1429},{29100,1430},{39100,1431},{49100,1432},{59100,1433},
608      {69100,1434},{79100,1435},{89100,1436},{9101,1446},{19101,1447},
609      {29101,1448},{39101,1449},{49101,1450},{59101,1451},{69101,1452},
610      {79101,1453},{89101,1454},{99101,1455},{109101,1456},{119101,1457},
611      {129101,1458},{9102,1459},{19102,1460},{29102,1461},{39102,1462},
612      {49102,1463},{59102,1464},{69102,1465},{79102,1466},{89102,1467},
613      {99102,1468},{109102,1469},{119102,1470},{129102,1471},{8100,1472},
614      {18100,1473},{28100,1474},{38100,1475},{48100,1476},{58100,1477},
615      {68100,1478},{78100,1479},{88100,1480},{8101,1490},{18101,1491},
616      {28101,1492},{38101,1493},{48101,1494},{58101,1495},{68101,1496},
617      {78101,1497},{88101,1498},{98101,1499},{108101,1500},{118101,1501},
618      {128101,1502},{8102,1503},{18102,1504},{28102,1505},{38102,1506},
619      {48102,1507},{58102,1508},{68102,1509},{78102,1510},{88102,1511},
620      {98102,1512},{108102,1513},{118102,1514},{128102,1515},{100,1516},
621      {10100,1517},{20100,1518},{30100,1519},{40100,1520},{50100,1521},
622      {60100,1522},{70100,1523},{80100,1524},{101,1525},{10101,1526},
623      {20101,1527},{30101,1528},{40101,1529},{50101,1530},{60101,1531},
624      {70101,1532},{80101,1533},{90101,1534},{100101,1535},{110101,1536},
625      {120101,1537},{102,1538},{10102,1539},{20102,1540},{30102,1541},
626      {40102,1542},{50102,1543},{60102,1544},{70102,1545},{80102,1546},
627      {90102,1547},{100102,1548},{110102,1549},{120102,1550},{1100,1551},
628      {11100,1552},{21100,1553},{31100,1554},{41100,1555},{51100,1556},
629      {61100,1557},{71100,1558},{81100,1559},{1101,1560},{11101,1561},
630      {21101,1562},{31101,1563},{41101,1564},{51101,1565},{61101,1566},
631      {71101,1567},{81101,1568},{91101,1569},{101101,1570},{111101,1571},
632      {121101,1572},{1102,1573},{11102,1574},{21102,1575},{31102,1576},
633      {41102,1577},{51102,1578},{61102,1579},{71102,1580},{81102,1581},
634      {91102,1582},{101102,1583},{111102,1584},{121102,1585},{4100,1586},
635      {14100,1587},{24100,1588},{34100,1589},{44100,1590},{54100,1591},
636      {64100,1592},{74100,1593},{84100,1594},{4101,1595},{14101,1596},
637      {24101,1597},{34101,1598},{44101,1599},{54101,1600},{64101,1601},
638      {74101,1602},{84101,1603},{94101,1604},{104101,1605},{114101,1606},
639      {124101,1607},{4102,1608},{14102,1609},{24102,1610},{34102,1611},
640      {44102,1612},{54102,1613},{64102,1614},{74102,1615},{84102,1616},
641      {94102,1617},{104102,1618},{114102,1619},{124102,1620},{3100,1621},
642      {13100,1622},{23100,1623},{33100,1624},{43100,1625},{53100,1626},
643      {63100,1627},{73100,1628},{83100,1629},{3101,1630},{13101,1631},
644      {23101,1632},{33101,1633},{43101,1634},{53101,1635},{63101,1636},
645      {73101,1637},{83101,1638},{93101,1639},{103101,1640},{113101,1641},
646      {123101,1642},{3102,1643},{13102,1644},{23102,1645},{33102,1646},
647      {43102,1647},{53102,1648},{63102,1649},{73102,1650},{83102,1651},
648      {93102,1652},{103102,1653},{113102,1654},{123102,1655},{2110,1656},
649      {12110,1657},{22110,1658},{32110,1659},{42110,1660},{52110,1661},
650      {62110,1662},{72110,1663},{82110,1664},{2111,1674},{12111,1675},
651      {22111,1676},{32111,1677},{42111,1678},{52111,1679},{62111,1680},
652      {72111,1681},{82111,1682},{92111,1683},{102111,1684},{112111,1685},
653      {122111,1686},{2112,1687},{12112,1688},{22112,1689},{32112,1690},
654      {42112,1691},{52112,1692},{62112,1693},{72112,1694},{82112,1695},
655      {92112,1696},{102112,1697},{112112,1698},{122112,1699},{7110,1700},
656      {17110,1701},{27110,1702},{37110,1703},{47110,1704},{57110,1705},
657      {67110,1706},{77110,1707},{87110,1708},{7111,1718},{17111,1719},
658      {27111,1720},{37111,1721},{47111,1722},{57111,1723},{67111,1724},
659      {77111,1725},{87111,1726},{97111,1727},{107111,1728},{117111,1729},
660      {127111,1730},{7112,1731},{17112,1732},{27112,1733},{37112,1734},
661      {47112,1735},{57112,1736},{67112,1737},{77112,1738},{87112,1739},
662      {97112,1740},{107112,1741},{117112,1742},{127112,1743},{5110,1744},
663      {15110,1745},{25110,1746},{35110,1747},{45110,1748},{55110,1749},
664      {65110,1750},{75110,1751},{85110,1752},{5111,1762},{15111,1763},
665      {25111,1764},{35111,1765},{45111,1766},{55111,1767},{65111,1768},
666      {75111,1769},{85111,1770},{95111,1771},{105111,1772},{115111,1773},
667      {125111,1774},{5112,1775},{15112,1776},{25112,1777},{35112,1778},
668      {45112,1779},{55112,1780},{65112,1781},{75112,1782},{85112,1783},
669      {95112,1784},{105112,1785},{115112,1786},{125112,1787},{6110,1788},
670      {16110,1789},{26110,1790},{36110,1791},{46110,1792},{56110,1793},
671      {66110,1794},{76110,1795},{86110,1796},{6111,1806},{16111,1807},
672      {26111,1808},{36111,1809},{46111,1810},{56111,1811},{66111,1812},
673      {76111,1813},{86111,1814},{96111,1815},{106111,1816},{116111,1817},
674      {126111,1818},{6112,1819},{16112,1820},{26112,1821},{36112,1822},
675      {46112,1823},{56112,1824},{66112,1825},{76112,1826},{86112,1827},
676      {96112,1828},{106112,1829},{116112,1830},{126112,1831},{9110,1832},
677      {19110,1833},{29110,1834},{39110,1835},{49110,1836},{59110,1837},
678      {69110,1838},{79110,1839},{89110,1840},{9111,1850},{19111,1851},
679      {29111,1852},{39111,1853},{49111,1854},{59111,1855},{69111,1856},
680      {79111,1857},{89111,1858},{99111,1859},{109111,1860},{119111,1861},
681      {129111,1862},{9112,1863},{19112,1864},{29112,1865},{39112,1866},
682      {49112,1867},{59112,1868},{69112,1869},{79112,1870},{89112,1871},
683      {99112,1872},{109112,1873},{119112,1874},{129112,1875},{8110,1876},
684      {18110,1877},{28110,1878},{38110,1879},{48110,1880},{58110,1881},
685      {68110,1882},{78110,1883},{88110,1884},{8111,1894},{18111,1895},
686      {28111,1896},{38111,1897},{48111,1898},{58111,1899},{68111,1900},
687      {78111,1901},{88111,1902},{98111,1903},{108111,1904},{118111,1905},
688      {128111,1906},{8112,1907},{18112,1908},{28112,1909},{38112,1910},
689      {48112,1911},{58112,1912},{68112,1913},{78112,1914},{88112,1915},
690      {98112,1916},{108112,1917},{118112,1918},{128112,1919},{110,1920},
691      {10110,1921},{20110,1922},{30110,1923},{40110,1924},{50110,1925},
692      {60110,1926},{70110,1927},{80110,1928},{111,1938},{10111,1939},
693      {20111,1940},{30111,1941},{40111,1942},{50111,1943},{60111,1944},
694      {70111,1945},{80111,1946},{90111,1947},{100111,1948},{110111,1949},
695      {120111,1950},{112,1951},{10112,1952},{20112,1953},{30112,1954},
696      {40112,1955},{50112,1956},{60112,1957},{70112,1958},{80112,1959},
697      {90112,1960},{100112,1961},{110112,1962},{120112,1963},{1110,1964},
698      {11110,1965},{21110,1966},{31110,1967},{41110,1968},{51110,1969},
699      {61110,1970},{71110,1971},{81110,1972},{1111,1982},{11111,1983},
700      {21111,1984},{31111,1985},{41111,1986},{51111,1987},{61111,1988},
701      {71111,1989},{81111,1990},{91111,1991},{101111,1992},{111111,1993},
702      {121111,1994},{1112,1995},{11112,1996},{21112,1997},{31112,1998},
703      {41112,1999},{51112,2000},{61112,2001},{71112,2002},{81112,2003},
704      {91112,2004},{101112,2005},{111112,2006},{121112,2007},{4110,2008},
705      {14110,2009},{24110,2010},{34110,2011},{44110,2012},{54110,2013},
706      {64110,2014},{74110,2015},{84110,2016},{4111,2026},{14111,2027},
707      {24111,2028},{34111,2029},{44111,2030},{54111,2031},{64111,2032},
708      {74111,2033},{84111,2034},{94111,2035},{104111,2036},{114111,2037},
709      {124111,2038},{4112,2039},{14112,2040},{24112,2041},{34112,2042},
710      {44112,2043},{54112,2044},{64112,2045},{74112,2046},{84112,2047},
711      {94112,2048},{104112,2049},{114112,2050},{124112,2051},{3110,2052},
712      {13110,2053},{23110,2054},{33110,2055},{43110,2056},{53110,2057},
713      {63110,2058},{73110,2059},{83110,2060},{3111,2070},{13111,2071},
714      {23111,2072},{33111,2073},{43111,2074},{53111,2075},{63111,2076},
715      {73111,2077},{83111,2078},{93111,2079},{103111,2080},{113111,2081},
716      {123111,2082},{3112,2083},{13112,2084},{23112,2085},{33112,2086},
717      {43112,2087},{53112,2088},{63112,2089},{73112,2090},{83112,2091},
718      {93112,2092},{103112,2093},{113112,2094},{123112,2095} };
719
720    // mode for the analysis
721    unsigned int _mode;
722
723    // 3D vector: [jet radius][lambda variable][pt bin]
724    // since each pt bin has its own normalised distribution
725    vector<vector<vector<Histo1DPtr> > > _h_dijet_cen,
726                                         _h_dijet_cen_groomed,
727                                         _h_dijet_fwd,
728                                         _h_dijet_fwd_groomed;
729
730    // 3D vector: [jet radius][lambda variable][pt bin]
731    // since each pt bin has its own normalised distribution
732    vector<vector<vector<Histo1DPtr> > > _h_zpj, _h_zpj_groomed;
733
734  };
735
736
737
738  RIVET_DECLARE_PLUGIN(CMS_2021_I1920187);
739
740}