1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/DirectFinalState.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9#include "Rivet/Projections/InvisibleFinalState.hh"
 10
 11namespace Rivet {
 12
 13  /// @brief ATLAS ZZ --> 4l and 2l2v
 14  class ATLAS_2016_I1494075 : public Analysis {
 15  public:
 16    /// Constructor
 17    //
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1494075);
 19    /// @name Analysis methods
 20    /// @{
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23        _mode = 0;
 24        if (getOption("LMODE") == "2L2NU") _mode = 2;
 25        if (getOption("LMODE") == "4L") _mode = 1;
 26
 27        PromptFinalState prompt_photons(Cuts::abspid == PID::PHOTON);
 28        PromptFinalState prompt_ele(Cuts::abspid == PID::ELECTRON);
 29        PromptFinalState prompt_mu(Cuts::abspid == PID::MUON);
 30
 31        // Wide lepton cuts which cover both channels and are used for the jet veto.
 32        Cut dressedele_cuts = (Cuts::abseta < 4.9) && (Cuts::pT > 7*GeV);
 33        Cut dressedmu_cuts = (Cuts::abseta < 2.7) && (Cuts::pT > 7*GeV);
 34        const LeptonFinder dressedelectrons(prompt_ele, prompt_photons, 0.1, dressedele_cuts);
 35        const LeptonFinder dressedmuons(prompt_mu, prompt_photons, 0.1, dressedmu_cuts);
 36
 37        declare(dressedelectrons, "electrons");
 38        declare(dressedmuons, "muons");
 39
 40        VisibleFinalState vfs;
 41        VetoedFinalState jetinput(vfs);
 42        jetinput.addVetoOnThisFinalState(dressedmuons);
 43
 44        if (_mode != 1)  declare(InvisibleFinalState(OnlyPrompt::YES), "MET");
 45
 46        FastJets fastjets(jetinput, JetAlg::ANTIKT, 0.4);
 47        declare (fastjets, "Jets");
 48
 49        // ZZ to four leptons channel
 50        book(_h["leading_ll_pt"], 2, 1, 1);
 51        book(_h["Njets"], 3, 1, 1);
 52        book(_h["leading_ll_phi"], 4, 1, 1);
 53        book(_h["ZZ_rapidity"], 5, 1, 1);
 54        // ZZ to lvlv channel
 55        book(_h["dilepton_pt"], 6, 1, 1);
 56        book(_h["llphi_lvchannel"], 7, 1, 1);
 57        book(_h["mzz_lvchannel"], 8, 1, 1);
 58
 59    }
 60
 61    struct Zstate : public ParticlePair {
 62      Zstate() { }
 63      Zstate(ParticlePair _particlepair) : ParticlePair(_particlepair) { }
 64      FourMomentum mom() const { return first.momentum() + second.momentum();}
 65      double rapid() const { return ((first.momentum() + second.momentum()).rapidity()); }
 66      double dphi() const { return deltaPhi(first.ptvec(), second.ptvec()); }
 67    };
 68
 69
 70    /// Perform the per-event analysis
 71    void analyze(const Event& event) {
 72
 73      // find out how many good jets
 74      //Jets jets = apply<FastJets>(event, "Jets").jetsByPt();
 75      Particles cand_e = apply<LeptonFinder>(event, "electrons").particlesByPt();
 76      Particles cand_mu = apply<LeptonFinder>(event, "muons").particlesByPt();
 77
 78      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::abseta < 4.5 && Cuts::pT>25*GeV);
 79      idiscardIfAnyDeltaRLess(jets, cand_e, 0.3);
 80
 81      //llvv channel
 82      if (_mode != 1) {
 83
 84        // jet veto
 85        if (jets.empty()){
 86
 87          Particles selected_pair;
 88          Vector3 met_vec;
 89
 90          const FinalState& metfs = apply<InvisibleFinalState>(event, "MET");
 91          for (const Particle& p : metfs.particles())  met_vec += p.mom().perpVec();
 92
 93          for ( const Particle& mu : cand_mu ) {
 94            if (mu.pT() > 25*GeV && mu.abseta() < 2.5 ) {
 95              selected_pair.push_back(mu);
 96            }
 97          }
 98          for ( const Particle& e : cand_e ) {
 99            if (e.pT() > 25*GeV && e.abseta() < 2.5 ) {
100              selected_pair.push_back(e);
101            }
102          }
103
104          // selections on pair.
105          if ((selected_pair.size() == 2) // exactly two leptons
106              && (selected_pair[0].abspid() == selected_pair[1].abspid()) //same flavour
107          && (deltaR(selected_pair[0],selected_pair[1]) > 0.3) // deltaR > 0.3
108          && (selected_pair[0].pid() * selected_pair[1].pid() < 0)) {  // opposite sign
109
110            //ensure the first one is leading lepton
111            if (selected_pair[0].momentum().pT() < selected_pair[1].momentum().pT()) {
112              std::swap(selected_pair[0], selected_pair[1]);
113            }
114
115            //ZZ four momentum, three momentum if from invisble final states
116            const FourMomentum Z_1_mom = selected_pair[0].momentum() + selected_pair[1].momentum();
117
118            const double axial_Etmiss = -1.0*met_vec.mod()*cos(deltaPhi(met_vec, Z_1_mom.ptvec()));
119
120            double pT_balance = fabs( (met_vec.mod() - Z_1_mom.pT()) /Z_1_mom.pT() );
121            if (axial_Etmiss > 90*GeV && pT_balance < 0.4 && inRange(Z_1_mom.mass(), 76*GeV, 106*GeV) ) {
122
123              double mz_pdg2 = 91.1876*91.1876*GeV*GeV;
124              // transverse mass
125              double mTrans = sqrt(
126                 sqr((sqrt(Z_1_mom.pT()*Z_1_mom.pT() + mz_pdg2) + sqrt(met_vec.mod2() + mz_pdg2)))
127                 - (Z_1_mom.ptvec() + met_vec.perpVec()).mod2()
128              );
129              _h["mzz_lvchannel"]->fill(mTrans/GeV);
130              _h["llphi_lvchannel"]->fill(deltaPhi(selected_pair[0].momentum().ptvec(), selected_pair[1].momentum().ptvec()));
131              _h["dilepton_pt"]->fill(Z_1_mom.pT()/GeV);
132
133            }
134          }
135        }
136      }
137
138      //for llll
139      if (_mode != 2) {
140
141        ///////////
142        // Insert selected muons then electrons into the lepton 4l final state
143        ///////////
144        DressedLeptons leptonsFS_sel4l;
145        leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), cand_mu.begin(), cand_mu.end() );
146        leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), cand_e.begin(), cand_e.end() );
147
148        ////////////
149        // Cut dR>0.2 between all leptons
150        Particles n_parts;
151        for (const DressedLepton& l1 : leptonsFS_sel4l) {
152          bool isolated = true;
153          for (DressedLepton& l2 : leptonsFS_sel4l){
154                  const double fourL_dR = deltaR(l1, l2);
155                  if (fourL_dR < 0.2 && !isSame(l1, l2)) {
156              isolated = false;
157              break;
158            }
159          }
160          if (isolated) n_parts.push_back(l1);
161        }
162
163        double totalCharge = 0;
164        for (const Particle& p : n_parts) totalCharge += p.pid();
165
166        if (n_parts.size() == 4 && totalCharge == 0 ) {
167
168          Zstate lead_Z, sub_Z;
169          identifyZstates(lead_Z, sub_Z, n_parts);
170          if (lead_Z.mom().pT() < sub_Z.mom().pT()) {
171            std::swap(lead_Z, sub_Z);
172          }
173
174          DressedLeptons lepton4l;
175          lepton4l.insert( lepton4l.end(), leptonsFS_sel4l.begin(), leptonsFS_sel4l.end() );
176          std::sort(lepton4l.begin(), lepton4l.end(), [](const DressedLepton& l1, const DressedLepton& l2) {
177            return (l1.abseta() > l2.abseta());
178          });
179          if (lead_Z.first.abspid() == 11 && sub_Z.first.abspid() == 11) {
180            if (lepton4l[1].abseta() > 2.5 || lepton4l[2].abseta() > 2.5 || lepton4l[3].abseta() > 2.5) vetoEvent;
181          }
182          else if (lead_Z.first.abspid() != sub_Z.first.abspid()) {
183            if (lead_Z.first.abspid() == 11) {
184               if (std::min(lead_Z.first.abseta(), lead_Z.second.abseta()) > 2.5)  vetoEvent;
185            }
186            if (lead_Z.first.abspid() == 13) {
187               if (std::min(sub_Z.first.abseta(), sub_Z.second.abseta()) > 2.5)  vetoEvent;
188            }
189          }
190
191          double m_Z1      = lead_Z.mom().mass();
192          double m_Z2      = sub_Z.mom().mass();
193          double lead_Z_Pt = lead_Z.mom().pT();
194          double lead_dPhi  = lead_Z.dphi();
195
196          double ZZ_rap    = fabs(lead_Z.rapid() - sub_Z.rapid());
197
198          //Z mass selections
199          if ( inRange(m_Z1, 66*GeV, 116*GeV) && inRange(m_Z2, 66*GeV, 116*GeV) ) {
200            _h["leading_ll_pt"]->fill(lead_Z_Pt/GeV);
201            _h["leading_ll_phi"]->fill(lead_dPhi);
202            _h["ZZ_rapidity"]->fill(ZZ_rap);
203            _h["Njets"]->fill(jets.size());
204          }
205        }
206      }
207    };
208
209
210
211    /// Normalise histograms etc., after the run
212    void finalize() {
213      // histo1D is divided by bin width when converting to scatter2D so no need of further normalisation for this one.
214      const double sf  = crossSectionPerEvent()/femtobarn;
215      const double sf2 = crossSectionPerEvent()/picobarn;
216      // 4l is divided by branching ratio to make a ZZ cross section
217      const double br = (3.3632 + 3.3662)/100.;
218      scale(_h["leading_ll_pt"],  sf/sqr(br));
219      scale(_h["Njets"],          sf/sqr(br));
220      scale(_h["leading_ll_phi"], sf2/sqr(br));
221      scale(_h["ZZ_rapidity"],    sf2/sqr(br));
222      // llvv is cross section for a single flavour. The analysis assumes we have run on e and mu,
223      // so divides by 2.
224      scale(_h["dilepton_pt"], sf/2.0);
225      scale(_h["mzz_lvchannel"], sf/2.0);
226      scale(_h["llphi_lvchannel"],sf/2.0);
227
228    }
229
230  private:
231
232    void identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& n_parts);
233    const double Zmass = 91.1876*GeV; // GeV
234    map<string, Histo1DPtr> _h;
235    size_t _mode;
236
237  };
238
239
240  void ATLAS_2016_I1494075::identifyZstates(Zstate& Z1, Zstate& Z2, const Particles& n_parts){
241
242
243    // first find the lepton types
244    Particles part_pos_el, part_neg_el, part_pos_mu, part_neg_mu;
245    for (const Particle& l : n_parts) {
246      if (l.abspid() == PID::ELECTRON) {
247        if (l.pid() < 0) part_neg_el.push_back(l);
248        if (l.pid() > 0) part_pos_el.push_back(l);
249      }
250      else if (l.abspid() == PID::MUON) {
251        if (l.pid() < 0) part_neg_mu.push_back(l);
252        if (l.pid() > 0) part_pos_mu.push_back(l);
253      }
254    }
255
256    //4e/4mu channel, pairing ambiguity
257    if (part_neg_el.size() == 2 || part_neg_mu.size() == 2) {
258      Zstate Zcand1, Zcand2, Zcand3, Zcand4;
259      if (part_neg_el.size() == 2) {
260        Zcand1 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[0] ) );
261        Zcand2 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[1] ) );
262        Zcand3 = Zstate( ParticlePair( part_neg_el[1], part_pos_el[0] ) );
263        Zcand4 = Zstate( ParticlePair( part_neg_el[1], part_pos_el[1] ) );
264      } else {
265        Zcand1 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[0] ) );
266        Zcand2 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[1] ) );
267        Zcand3 = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[0] ) );
268        Zcand4 = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[1] ) );
269      }
270
271      // pairing should be |1 + 4| and |2 + 3| in mass order
272      double V1, V2;
273      V1 = fabs( Zcand1.mom().mass() - Zmass ) + fabs( Zcand4.mom().mass() - Zmass);
274      V2 = fabs( Zcand2.mom().mass() - Zmass ) + fabs( Zcand3.mom().mass() - Zmass);
275
276      if (V1 > V2) {
277        Z1 = Zcand2;
278        Z2 = Zcand3;
279      }
280      else {
281        Z1 = Zcand1;
282        Z2 = Zcand4;
283      }
284      //2e2mu
285    }
286    else if (part_neg_el.size() == 1 && part_neg_mu.size() == 1) {
287      Z1 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[0] ) );
288      Z2 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[0] ) );
289    }
290
291  }
292
293  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1494075);
294}