1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Tools/Logging.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include "Rivet/Tools/ParticleIdUtils.hh"
  9#include "Rivet/Particle.hh"
 10
 11namespace Rivet {
 12
 13
 14  class ATLAS_2013_I1243871 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    ATLAS_2013_I1243871()
 19      : Analysis("ATLAS_2013_I1243871")
 20    {    }
 21
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25      // Set up projections
 26      const FinalState fs((Cuts::etaIn(-4.5, 4.5)));
 27      declare(fs, "ALL_FS");
 28
 29      /// Get electrons from truth record
 30      IdentifiedFinalState elec_fs(Cuts::abseta < 2.47 && Cuts::pT > 25*GeV);
 31      elec_fs.acceptIdPair(PID::ELECTRON);
 32      declare(elec_fs, "ELEC_FS");
 33
 34      /// Get muons which pass the initial kinematic cuts:
 35      IdentifiedFinalState muon_fs(Cuts::abseta < 2.5 && Cuts::pT > 20*GeV);
 36      muon_fs.acceptIdPair(PID::MUON);
 37      declare(muon_fs, "MUON_FS");
 38
 39      // Final state used as input for jet-finding.
 40      // We include everything except the muons and neutrinos
 41      VetoedFinalState jet_input(fs);
 42      jet_input.vetoNeutrinos();
 43      jet_input.addVetoPairId(PID::MUON);
 44      declare(jet_input, "JET_INPUT");
 45
 46      // Get the jets
 47      FastJets jets(jet_input, JetAlg::ANTIKT, 0.4);
 48      declare(jets, "JETS");
 49
 50      // Book histograms
 51      for (size_t d = 0; d < 5; ++d) {
 52        book(_p_b_rho[d] ,d+1, 1, 1);
 53        book(_p_l_rho[d] ,d+1, 2, 1);
 54        book(_p_b_Psi[d] ,d+1, 1, 2);
 55        book(_p_l_Psi[d] ,d+1, 2, 2);
 56      }
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      /// Get the various sets of final state particles
 63      const Particles& elecFS = apply<IdentifiedFinalState>(event, "ELEC_FS").particlesByPt();
 64      const Particles& muonFS = apply<IdentifiedFinalState>(event, "MUON_FS").particlesByPt();
 65
 66      // Get all jets with pT > 7 GeV (ATLAS standard jet collection)
 67      /// @todo Why rewrite the jets collection as a vector of pointers?
 68      const Jets& jets = apply<FastJets>(event, "JETS").jetsByPt(Cuts::pT > 7*GeV);
 69      vector<const Jet*> allJets;
 70      for (const Jet& j : jets) allJets.push_back(&j);
 71
 72      // Keep any jets that pass the pt cut
 73      vector<const Jet*> pt_jets;
 74      for (const Jet* j : allJets) {
 75        /// @todo Use direct kinematics access
 76        const double pt = j->momentum().pT();
 77        const double eta = j->momentum().eta();
 78        if (pt > 25*GeV && fabs(eta) < 2.5) pt_jets.push_back(j);
 79      }
 80
 81      // Remove jets too close to an electron
 82      vector<const Jet*> good_jets;
 83      for (const Jet* j : pt_jets) {
 84        bool isElectron = 0;
 85        for (const Particle& e : elecFS) {
 86          const double elec_jet_dR = deltaR(e.momentum(), j->momentum());
 87          if (elec_jet_dR < 0.2) { isElectron = true; break; }
 88        }
 89        if (!isElectron) good_jets.push_back(j);
 90      }
 91
 92      // Classify the event type
 93      const size_t nElec = elecFS.size();
 94      const size_t nMuon = muonFS.size();
 95      bool isSemilepton = false, isDilepton = false;
 96      if (nElec == 1 && nMuon == 0) {
 97        isSemilepton = true;
 98      } else if (nElec == 0 && nMuon == 1) {
 99        isSemilepton = true;
100      } else if (nElec == 2 && nMuon == 0) {
101        if (charge(elecFS[0]) != charge(elecFS[1])) isDilepton = true;
102      } else if (nElec == 1 && nMuon == 1) {
103        if (charge(elecFS[0]) != charge(muonFS[0])) isDilepton = true;
104      } else if (nElec == 0 && nMuon == 2) {
105        if (charge(muonFS[0]) != charge(muonFS[1])) isDilepton = true;
106      }
107      const bool isGoodEvent = (isSemilepton && good_jets.size() >= 4) || (isDilepton && good_jets.size() >= 2);
108      if (!isGoodEvent) vetoEvent;
109
110      // Select b-hadrons
111      /// @todo Use built-in identification on Particle, avoid HepMC
112      vector<ConstGenParticlePtr> b_hadrons;
113      vector<ConstGenParticlePtr> allParticles = HepMCUtils::particles(event.genEvent());
114      for (size_t i = 0; i < allParticles.size(); i++) {
115        ConstGenParticlePtr p = allParticles.at(i);
116        if ( !(PID::isHadron( p->pdg_id() ) && PID::hasBottom( p->pdg_id() )) ) continue;
117        if (p->momentum().perp() < 5*GeV) continue;
118        b_hadrons.push_back(p);
119      }
120
121      // Select b-jets as those containing a b-hadron
122      /// @todo Use built-in dR < 0.3 Jet tagging, avoid HepMC
123      vector<const Jet*> b_jets;
124      for (const Jet* j : good_jets) {
125        bool isbJet = false;
126        for (ConstGenParticlePtr b : b_hadrons) {
127          /// @todo Use direct momentum accessor / delta functions
128          const FourMomentum hadron = b->momentum();
129          const double hadron_jet_dR = deltaR(j->momentum(), hadron);
130          if (hadron_jet_dR < 0.3) { isbJet = true; break; }
131        }
132        // Check if it is overlapped to any other jet
133        bool isOverlapped = false;
134        for (const Jet* k : allJets) {
135          if (j == k) continue;
136          double dRjj = deltaR(j->momentum(), k->momentum());
137          if (dRjj < 0.8) { isOverlapped = true; break; }
138        }
139        if (isbJet && !isOverlapped) b_jets.push_back(j);
140      }
141      MSG_DEBUG(b_jets.size() << " b-jets selected");
142
143
144      // Select light-jets as the pair of non-b-jets with invariant mass closest to the W mass
145      /// @todo Use built-in b-tagging (dR < 0.3 defn), avoid HepMC
146      const double nominalW = 80.4*GeV;
147      double deltaM = 500*GeV;
148      const Jet* light1 = NULL; const Jet* light2 = NULL; // NB: const Jets, not const pointers!
149      for (const Jet* i : good_jets) {
150        bool isbJet1 = false;
151        for (ConstGenParticlePtr b : b_hadrons) {
152          /// @todo Use direct momentum accessor / delta functions
153          const FourMomentum hadron = b->momentum();
154          const double hadron_jet_dR = deltaR(i->momentum(), hadron);
155          if (hadron_jet_dR < 0.3) { isbJet1 = true; break; }
156        }
157        if (isbJet1) continue;
158        for (const Jet* j : good_jets) {
159          bool isbJet2 = false;
160          for (ConstGenParticlePtr b : b_hadrons) {
161            FourMomentum hadron = b->momentum();
162            double hadron_jet_dR = deltaR(j->momentum(), hadron);
163            if (hadron_jet_dR < 0.3) { isbJet2 = true; break; }
164          }
165          if (isbJet2) continue;
166          double invMass = (i->momentum()+j->momentum()).mass();
167          if (fabs(invMass-nominalW) < deltaM){
168            deltaM = fabs(invMass - nominalW);
169            light1 = i;
170            light2 = j;
171          }
172        }
173      }
174
175      // Check that both jets are not overlapped, and populate the light jets list
176      vector<const Jet*> light_jets;
177      const bool hasGoodLight = light1 != NULL && light2 != NULL && light1 != light2;
178      if (hasGoodLight) {
179        bool isOverlap1 = false, isOverlap2 = false;
180        for (const Jet* j : allJets) {
181          if (light1 == j) continue;
182          const double dR1j = deltaR(light1->momentum(), j->momentum());
183          if (dR1j < 0.8) { isOverlap1 = true; break; }
184        }
185        for (const Jet* j : allJets) {
186          if (light2 == j) continue;
187          const double dR2j = deltaR(light2->momentum(), j->momentum());
188          if (dR2j < 0.8) { isOverlap2 = true; break; }
189        }
190        if (!isOverlap1 && !isOverlap2) {
191          light_jets.push_back(light1);
192          light_jets.push_back(light2);
193        }
194      }
195      MSG_DEBUG(light_jets.size() << " light jets selected");
196
197
198      // Calculate the jet shapes
199      /// @todo Use C++11 vector/array initialization
200      const double binWidth = 0.04; // -> 10 bins from 0.0-0.4
201      vector<double> ptEdges; ptEdges += {{ 30, 40, 50, 70, 100, 150 }};
202
203      // b-jet shapes
204      MSG_DEBUG("Filling b-jet shapes");
205      for (const Jet* bJet : b_jets) {
206        // Work out jet pT bin and skip this jet if out of range
207        const double jetPt = bJet->momentum().pT();
208        MSG_DEBUG("Jet pT = " << jetPt/GeV << " GeV");
209        if (!inRange(jetPt/GeV, 30., 150.)) continue;
210        /// @todo Use YODA bin index lookup tools
211        size_t ipt; for (ipt = 0; ipt < 5; ++ipt) if (inRange(jetPt/GeV, ptEdges[ipt], ptEdges[ipt+1])) break;
212        MSG_DEBUG("Jet pT index = " << ipt);
213
214        // Calculate jet shape
215        vector<double> rings(10, 0);
216        for (const Particle& p : bJet->particles()) {
217          const double dR = deltaR(bJet->momentum(), p.momentum());
218          const size_t idR = (size_t) floor(dR/binWidth);
219          for (size_t i = idR; i < 10; ++i) rings[i] += p.pT();
220        }
221
222        // Fill each dR bin of the histos for this jet pT
223        for (int iBin = 0; iBin < 10; ++iBin) {
224          const double rcenter = 0.02 + iBin*binWidth;
225          const double rhoval = (iBin != 0 ? (rings[iBin]-rings[iBin-1]) : rings[iBin]) / binWidth / rings[9];
226          const double psival = rings[iBin] / rings[9];
227          MSG_DEBUG(rcenter << ", " << rhoval << ", " << psival);
228          _p_b_rho[ipt]->fill(rcenter, rhoval);
229          _p_b_Psi[ipt]->fill(rcenter, psival);
230        }
231      }
232
233      // Light jet shapes
234      MSG_DEBUG("Filling light jet shapes");
235      for (const Jet* lJet : light_jets) {
236        // Work out jet pT bin and skip this jet if out of range
237        const double jetPt = lJet->momentum().pT();
238        MSG_DEBUG("Jet pT = " << jetPt/GeV << " GeV");
239        if (!inRange(jetPt/GeV, 30., 150.)) continue;
240        /// @todo Use YODA bin index lookup tools
241        size_t ipt; for (ipt = 0; ipt < 5; ++ipt) if (inRange(jetPt/GeV, ptEdges[ipt], ptEdges[ipt+1])) break;
242        MSG_DEBUG("Jet pT index = " << ipt);
243
244        // Calculate jet shape
245        vector<double> rings(10, 0);
246        for (const Particle& p : lJet->particles()) {
247          const double dR = deltaR(lJet->momentum(), p.momentum());
248          const size_t idR = (size_t) floor(dR/binWidth);
249          for (size_t i = idR; i < 10; ++i) rings[i] += p.pT();
250        }
251
252        // Fill each dR bin of the histos for this jet pT
253        for (int iBin = 0; iBin < 10; ++iBin) {
254          const double rcenter = 0.02 + iBin*binWidth;
255          const double rhoval = (iBin != 0 ? (rings[iBin]-rings[iBin-1]) : rings[iBin]) / binWidth / rings[9];
256          const double psival = rings[iBin] / rings[9];
257          _p_l_rho[ipt]->fill(rcenter, rhoval);
258          _p_l_Psi[ipt]->fill(rcenter, psival);
259        }
260      }
261
262    }
263
264    /// @todo why does this routine not have a finalize method?
265    /// not clear how you would combine different samples slices
266    /// correctly if you don't weight by cross-section
267
268
269  private:
270
271    Profile1DPtr _p_b_rho[5];
272    Profile1DPtr _p_l_rho[5];
273    Profile1DPtr _p_b_Psi[5];
274    Profile1DPtr _p_l_Psi[5];
275
276  };
277
278
279  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1243871);
280
281}