1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// Drell-Yan differential cross section measurement @ 13 TeV
 11  class CMS_2018_I1711625 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2018_I1711625);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // prompt final state electrons
 25      const PromptFinalState el_pfs = PromptFinalState(Cuts::abspid == PID::ELECTRON);
 26      declare(el_pfs, "PromptFinalStateElectrons");
 27
 28      // prompt final state muons
 29      const PromptFinalState mu_pfs = PromptFinalState(Cuts::abspid == PID::MUON);
 30      declare(mu_pfs, "PromptFinalStateMuons");
 31
 32      // dressed leptons
 33      const FinalState photon_fs = FinalState(Cuts::abspid == PID::PHOTON);
 34
 35      const LeptonFinder mu_dressed(mu_pfs, photon_fs, 0.1, Cuts::open());
 36      declare(mu_dressed, "DressedMuons");
 37
 38      book(_h_massMuMu, 3, 1, 1); /// muon channel result in full-phase space @ dressed level
 39      book(_h_massMuMuFiducial, 5, 1, 1); /// muon channel result in fiducial region @ post-FSR level
 40      book(_h_massEEFiducial, 6, 1, 1); /// electron channel result in fiducial region @ post-FSR level
 41    }
 42
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46
 47      const LeptonFinder muons_dressed = apply<LeptonFinder>(event, "DressedMuons");
 48      bool filled_mu = FillHistogram_DressedLepton(muons_dressed, 13);
 49      if ( filled_mu ) {
 50
 51        const PromptFinalState muons_PFS = apply<PromptFinalState>(event, "PromptFinalStateMuons");
 52        FillHistogram_PFSLepton(muons_PFS, 13);
 53      }
 54      else { // electron channel
 55
 56        const PromptFinalState electrons_PFS = apply<PromptFinalState>(event, "PromptFinalStateElectrons");
 57        FillHistogram_PFSLepton(electrons_PFS, 11);
 58      }
 59    }
 60
 61
 62    /// Normalise histograms etc., after the run
 63    void finalize() {
 64
 65      scale(_h_massMuMu, crossSection()/picobarn/sumOfWeights()); /// norm to cross section
 66      scale(_h_massMuMuFiducial, crossSection()/picobarn/sumOfWeights()); /// norm to cross section
 67      scale(_h_massEEFiducial, crossSection()/picobarn/sumOfWeights()); /// norm to cross section
 68
 69    }
 70
 71    /// @}
 72
 73
 74  private:
 75
 76
 77    /// @name Histograms
 78    /// @{
 79    Histo1DPtr _h_massMuMu;
 80    Histo1DPtr _h_massMuMuFiducial;
 81    Histo1DPtr _h_massEEFiducial;
 82    /// @}
 83
 84
 85    // select two opposite sign leptons with highest pT & fill the histogram for full-phase space diff. x-section
 86    bool FillHistogram_DressedLepton(LeptonFinder leptons_dressed, int leptonID) {
 87      bool filled = false;
 88
 89      vector< DressedLepton > vec_dressedLepByPt = leptons_dressed.dressedLeptons();
 90
 91      int nLepton_dressed = (int)vec_dressedLepByPt.size();
 92      if ( nLepton_dressed >= 2 ) {
 93        int index_lepton1 = -1;
 94        int index_lepton2 = -1;
 95        FindDressedLeptonPair_HighestPt(vec_dressedLepByPt, index_lepton1, index_lepton2);
 96
 97        if ( index_lepton1 != -1 && index_lepton2 != -1 ) {
 98          DressedLepton lepton1_dressed = vec_dressedLepByPt[index_lepton1];
 99          DressedLepton lepton2_dressed = vec_dressedLepByPt[index_lepton2];
100
101          const FourMomentum pVec_diLep = lepton1_dressed.mom() + lepton2_dressed.mom();
102          double mass = pVec_diLep.mass();
103
104          // // fill histograms
105          if ( leptonID == 13 ) _h_massMuMu->fill(mass/GeV);
106
107          filled = true;
108        }
109      } // end of if( nLepton_dressed >= 2 )
110
111      return filled;
112    }
113
114
115    void FindDressedLeptonPair_HighestPt(DressedLeptons& vec_dressedLepByPt, int& index_lepton1, int& index_lepton2) {
116      // 1st lepton: lepton with highest pT
117      int nLepton_dressed = int(vec_dressedLepByPt.size());
118      for (int i=0; i<nLepton_dressed; ++i) {
119        auto& lepton = vec_dressedLepByPt[i]; // decreasing order of pT
120        if ( lepton.isLepton() ) {
121          index_lepton1 = i;
122          break;
123        }
124      }
125
126      // if no lepton is found in the leptons_dressed
127      if ( index_lepton1 < 0 ) {
128        index_lepton1 = -1;
129        index_lepton2 = -1;
130        return;
131      }
132
133      // 2nd lepton: lepton with highest-pT among the leptons with the opposite sign with 1st lepton
134      int pdgID_lepton1 = vec_dressedLepByPt[index_lepton1].pid();
135      for (int i=index_lepton1+1; i<nLepton_dressed; ++i) { // starting after lepton1
136        auto& lepton = vec_dressedLepByPt[i];
137        if ( lepton.isLepton() && lepton.pid() == (-1)*pdgID_lepton1 ) {
138          index_lepton2 = i;
139          break;
140        }
141      }
142    }
143
144
145    // select two opposite sign leptons with highest pT & fill the histogram for the fiducial diff. x-section
146    void FillHistogram_PFSLepton(PromptFinalState leptons_PFS, int leptonID) {
147      vector< Particle > vec_PFSLepByPt = leptons_PFS.particlesByPt();
148
149      int nLepton_PFS = int(vec_PFSLepByPt.size());
150      if ( nLepton_PFS >= 2 ) {
151        int index_lepton1 = -1;
152        int index_lepton2 = -1;
153        FindPFSLeptonPair_HighestPtWithinAcc(vec_PFSLepByPt, leptonID, index_lepton1, index_lepton2);
154        if ( index_lepton1 != -1 && index_lepton2 != -1 ) {
155          Particle lepton1_PFS = leptons_PFS.particlesByPt()[index_lepton1];
156          Particle lepton2_PFS = leptons_PFS.particlesByPt()[index_lepton2];
157
158          const FourMomentum pVec_diLep = lepton1_PFS.mom() + lepton2_PFS.mom();
159          double mass = pVec_diLep.mass();
160
161          if      ( leptonID == 11 ) _h_massEEFiducial->fill(mass/GeV);
162          else if ( leptonID == 13 ) _h_massMuMuFiducial->fill(mass/GeV);
163        }
164      } // end of if ( nLepton_PFS >= 2 )
165    }
166
167
168    void FindPFSLeptonPair_HighestPtWithinAcc(vector<Particle>& vec_PFSLepByPt, int pdgID, int& index_lepton1, int& index_lepton2) {
169      double pTCut_lead = 0;
170      if ( pdgID == 11 ) pTCut_lead = 30.0;
171      if ( pdgID == 13 ) pTCut_lead = 22.0;
172      double pTCut_sub = 10.0; // same for both channel
173
174      double etaCut_lead = 0;
175      if ( pdgID == 11 ) etaCut_lead = 2.5;
176      if ( pdgID == 13 ) etaCut_lead = 2.4;
177      double etaCut_sub = etaCut_lead;
178
179      int nLepton = int(vec_PFSLepByPt.size());
180      for (int i=0; i<nLepton; ++i) {
181        auto& lepton = vec_PFSLepByPt[i];
182        if ( lepton.isLepton() && lepton.pT() > pTCut_lead && lepton.abseta() < etaCut_lead ) {
183          if ( pdgID == 11 ) { // electron channel: check ECAL gap
184            if ( !(lepton.abseta() > 1.4442 && lepton.abseta() < 1.566) ) {
185              index_lepton1 = i;
186              break;
187            }
188          }
189          else { // muon channel
190            index_lepton1 = i;
191            break;
192          }
193        }
194      } // end of lepton iteration
195
196      // if no lepton is found in the leptons_PFS
197      if ( index_lepton1 < 0 ) {
198        index_lepton1 = -1;
199        index_lepton2 = -1;
200        return;
201      }
202
203      // 2nd lepton: lepton with highest-pT among the leptons with the opposite sign with 1st lepton
204      int pdgID_lepton1 = vec_PFSLepByPt[index_lepton1].pid();
205      for (int i=index_lepton1+1; i<nLepton; ++i) { // starting after lepton1
206        auto& lepton = vec_PFSLepByPt[i];
207        if ( lepton.isLepton() && lepton.pid() == (-1)*pdgID_lepton1 &&
208             lepton.pT() > pTCut_sub && lepton.abseta() < etaCut_sub ) {
209          if ( pdgID == 11 ) { // electron channel: check ECAL gap
210            if ( !(lepton.abseta() > 1.4442 && lepton.abseta() < 1.566) ) {
211              index_lepton2 = i;
212              break;
213            }
214          }
215          else { // muon channel
216            index_lepton2 = i;
217            break;
218          }
219        }
220      } // end of lepton iteration
221    }
222  };
223
224
225  RIVET_DECLARE_PLUGIN(CMS_2018_I1711625);
226
227
228}