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