1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// Double-parton scattering in inclusive production of four jets with low pT in 13 TeV pp
 11  class CMS_2021_I1932460 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2021_I1932460);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25
 26      // the basic final-state projection:
 27      // all final-state particles within
 28      // the given eta acceptance: 4.7 (jet range) + 0.4 (cone size)
 29      const FinalState fs(Cuts::abseta < 5.1);
 30
 31      // the final-state particles declared above are clustered using FastJet with
 32      // the anti-kT algorithm and a jet-radius parameter 0.4
 33      // muons and neutrinos are excluded from the clustering
 34      FastJets jetfs(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 35      declare(jetfs, "jets");
 36
 37      // Book histograms
 38      book(_h["JetPt1"], 1, 1, 1);
 39      book(_h["JetPt2"], 2, 1, 1);
 40      book(_h["JetPt3"], 3, 1, 1);
 41      book(_h["JetPt4"], 4, 1, 1);
 42      book(_h["JetEta1"], 5, 1, 1);
 43      book(_h["JetEta2"], 6, 1, 1);
 44      book(_h["JetEta3"], 7, 1, 1);
 45      book(_h["JetEta4"], 8, 1, 1);
 46
 47      book(_h["DeltaPhiSoft_binNorm"], 9, 1, 1);
 48      book(_h["DeltaPhi3_binNorm"], 10, 1, 1);
 49      book(_h["DeltaY_binNorm"], 11, 1, 1);
 50      book(_h["DeltaPhiY_binNorm"], 12, 1, 1);
 51      book(_h["DeltaPtSoft_binNorm"], 13, 1, 1);
 52      book(_h["DeltaS_binNorm"], 14, 1, 1);
 53
 54      book(_h["DeltaPhiSoft"], 45, 1, 1);
 55      book(_h["DeltaPhi3"], 46, 1, 1);
 56      book(_h["DeltaY"], 47, 1, 1);
 57      book(_h["DeltaPhiY"], 48, 1, 1);
 58      book(_h["DeltaPtSoft"], 49, 1, 1);
 59      book(_h["DeltaS"], 50, 1, 1);
 60
 61    }
 62
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      // retrieve clustered jets, sorted by pT, with a minimum pT cut 10 GeV and eta range 4.7 (similar to PFJet collection)
 68      Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::abseta < 4.7 && Cuts::pT > 10*GeV);
 69
 70      // fill only if there are at least 4 jets
 71      if (jets.size() < 4) vetoEvent;
 72
 73      double pt0 = jets[0].pt();
 74      double pt1 = jets[1].pt();
 75      double pt2 = jets[2].pt();
 76      double pt3 = jets[3].pt();
 77
 78      // pt selection
 79      if (pt0 < 35.0 || pt1 < 30.0 || pt2 < 25.0 || pt3 < 20.0) vetoEvent;
 80
 81      double phi0 = jets[0].phi();
 82      double phi1 = jets[1].phi();
 83      double phi2 = jets[2].phi();
 84      double phi3 = jets[3].phi();
 85
 86      _h["JetPt1"]->fill(pt0);
 87      _h["JetPt2"]->fill(pt1);
 88      _h["JetPt3"]->fill(pt2);
 89      _h["JetPt4"]->fill(pt3);
 90
 91      _h["JetEta1"]->fill(jets[0].eta());
 92      _h["JetEta2"]->fill(jets[1].eta());
 93      _h["JetEta3"]->fill(jets[2].eta());
 94      _h["JetEta4"]->fill(jets[3].eta());
 95
 96      // delta phi and eta of the 2 soft jets
 97      _h["DeltaPhiSoft"]->fill(abs(deltaPhi(phi2, phi3)));
 98      _h["DeltaPhiSoft_binNorm"]->fill(abs(deltaPhi(phi2, phi3)));
 99
100      // delta pt between soft jets
101      double DptSoft = sqrt(pow(pt2*cos(phi2) + pt3*cos(phi3), 2) + pow(pt2*sin(phi2) + pt3*sin(phi3), 2))/(pt2 + pt3);
102      _h["DeltaPtSoft"]->fill(DptSoft);
103      _h["DeltaPtSoft_binNorm"]->fill(DptSoft);
104
105      // delta S
106      if (pt0 > 50.0 && pt1 > 30.0 && pt2 > 30.0 && pt3 > 30.0) {
107        double phiH = atan2(pt0*sin(phi0) + pt1*sin(phi1) , pt0*cos(phi0) + pt1*cos(phi1));
108        double phiS = atan2(pt2*sin(phi2) + pt3*sin(phi3) , pt2*cos(phi2) + pt3*cos(phi3));
109        double DS = abs(deltaPhi(phiH, phiS));
110        _h["DeltaS"]->fill(DS);
111        _h["DeltaS_binNorm"]->fill(DS);
112      }
113
114      // delta Y: most remote jets in rapidity, find min & max eta
115      double mineta = 99999;
116      double maxeta = -99999;
117      int minetapos = -1;
118      int maxetapos = -1;
119
120      for (int i = 0; i < 4; ++i) {
121        if (jets[i].eta() < mineta) {
122          mineta = jets[i].eta();
123          minetapos = i;
124        }
125        if (jets[i].eta() > maxeta) {
126          maxeta = jets[i].eta();
127          maxetapos = i;
128        }
129      }
130
131      _h["DeltaY"]->fill(abs(jets[minetapos].eta() - jets[maxetapos].eta()));
132      _h["DeltaY_binNorm"]->fill(abs(jets[minetapos].eta() - jets[maxetapos].eta()));
133
134      // Delta phi Y: azimuthal angle between most remote jets in eta
135      _h["DeltaPhiY"]->fill(abs(deltaPhi(jets[minetapos].phi(), jets[maxetapos].phi())));
136      _h["DeltaPhiY_binNorm"]->fill(abs(deltaPhi(jets[minetapos].phi(), jets[maxetapos].phi())));
137
138      // delta phi3
139      double minphi3 = 999;
140      for (int iphi1 = 0; iphi1 < 4; ++iphi1) {
141        for (int iphi2 = 0; iphi2 < 4; ++iphi2) {
142          for (int iphi3 = 0; iphi3 < 4; ++iphi3) {
143            if ( !(iphi1 == iphi2 || iphi2 == iphi3 || iphi1 == iphi3) ) {
144              double temp_phi1= jets[iphi1].phi();
145              double temp_phi2= jets[iphi2].phi();
146              double temp_phi3= jets[iphi3].phi();
147              double temp_minphi3 = abs(deltaPhi(temp_phi1, temp_phi2)) + abs(deltaPhi(temp_phi2, temp_phi3));
148              if (temp_minphi3 < minphi3) minphi3 = temp_minphi3;
149            }
150          }
151        }
152      }
153
154      _h["DeltaPhi3"]->fill(minphi3);
155      _h["DeltaPhi3_binNorm"]->fill(minphi3);
156
157    }
158
159
160    /// Normalise histograms etc., after the run
161    void finalize() {
162
163      scale(_h["JetPt1"], crossSection()/picobarn/sumOfWeights()); // norm to cross section
164      scale(_h["JetPt2"], crossSection()/picobarn/sumOfWeights());
165      scale(_h["JetPt3"], crossSection()/picobarn/sumOfWeights());
166      scale(_h["JetPt4"], crossSection()/picobarn/sumOfWeights());
167      scale(_h["JetEta1"], crossSection()/picobarn/sumOfWeights());
168      scale(_h["JetEta2"], crossSection()/picobarn/sumOfWeights());
169      scale(_h["JetEta3"], crossSection()/picobarn/sumOfWeights());
170      scale(_h["JetEta4"], crossSection()/picobarn/sumOfWeights());
171      scale(_h["DeltaPhiSoft"], crossSection()/picobarn/sumOfWeights());
172      scale(_h["DeltaPhi3"], crossSection()/picobarn/sumOfWeights());
173      scale(_h["DeltaY"], crossSection()/picobarn/sumOfWeights());
174      scale(_h["DeltaPhiY"], crossSection()/picobarn/sumOfWeights());
175      scale(_h["DeltaPtSoft"], crossSection()/picobarn/sumOfWeights());
176      scale(_h["DeltaS"], crossSection()/picobarn/sumOfWeights());
177
178      // create bin normalised histograms
179
180      // Correct for binwidths: Rivet automatically normalises histograms to binwidth when plotting AFTER the normalisation executed here.
181      // So we must calculate an extra correction here so that finally our bin-normalised histograms end up around 1 as in the paper.
182      // in YODA bin index starts from 0
183
184      // For DeltaY, which has variable binwidths we need to do following steps
185      // divide histograms by binwidth
186      for (unsigned int i = 1; i < _h["DeltaY_binNorm"]->numBins()+1; ++i) {
187        _h["DeltaY_binNorm"]->bin(i).scaleW(1.0/_h["DeltaY_binNorm"]->bin(i).xWidth());
188      }
189
190      // normalise to average of first 4 bins
191      scale(_h["DeltaY_binNorm"], 1.0/(_h["DeltaY_binNorm"]->integralRange(1,4)/4.0));
192
193      // multiply again with binwidth
194      for (unsigned int i = 1; i < _h["DeltaY_binNorm"]->numBins()+1; ++i) {
195        _h["DeltaY_binNorm"]->bin(i).scaleW(_h["DeltaY_binNorm"]->bin(i).xWidth());
196      }
197
198      // DeltaPhiSoft and DeltaPhi3 histograms have uniform binwidths, so multiply with first binwidth is sufficient
199      scale(_h["DeltaPhiSoft_binNorm"], _h["DeltaPhiSoft_binNorm"]->bin(1).xWidth()/(_h["DeltaPhiSoft_binNorm"]->integralRange(1,5)/5.0));
200      scale(_h["DeltaPhi3_binNorm"], _h["DeltaPhi3_binNorm"]->bin(1).xWidth()/(_h["DeltaPhi3_binNorm"]->integralRange(1,4)/4.0));
201
202      // DeltaPhiY, DeltaPtSoft and DeltaS are normalised to last bin
203      scale(_h["DeltaPhiY_binNorm"], _h["DeltaPhiY_binNorm"]->bin(12).xWidth()/_h["DeltaPhiY_binNorm"]->bin(12).sumW() );
204      scale(_h["DeltaPtSoft_binNorm"], _h["DeltaPtSoft_binNorm"]->bin(8).xWidth()/_h["DeltaPtSoft_binNorm"]->bin(8).sumW() );
205      scale(_h["DeltaS_binNorm"], _h["DeltaS_binNorm"]->bin(7).xWidth()/_h["DeltaS_binNorm"]->bin(7).sumW() );
206    }
207
208    /// @}
209
210    map<string, Histo1DPtr> _h;
211
212  };
213
214
215  RIVET_DECLARE_PLUGIN(CMS_2021_I1932460);
216
217}