1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/HeavyHadrons.hh"
  4#include "Rivet/Projections/LeptonFinder.hh"
  5#include "Rivet/Projections/FinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/IdentifiedFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// BB to Jpsi plus mu at 8 TeV
 13  class ATLAS_2017_I1598613 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1598613);
 18
 19    struct HistoHandler {
 20      Histo1DPtr histo;
 21      Estimate1DPtr scatter;
 22      unsigned int d, x, y;
 23
 24      HistoHandler() {}
 25
 26      void fill(double value) {
 27        histo->fill(value);
 28      }
 29    };
 30
 31
 32    /// Book histograms and initialise projections before the run
 33    void init() {
 34
 35        // default to widest cut, electrons and muons.
 36        _mode = 0;
 37        if ( getOption("BMODE") == "BB" )  _mode = 1;
 38
 39      // Get the particles needed for each running mode:
 40      if (_mode == 0) {
 41        // Get photons to dress leptons
 42        FinalState photons(Cuts::abspid == PID::PHOTON);
 43        FinalState muons(Cuts::abspid == PID::MUON);
 44        Cut eta_lep = Cuts::abseta < 2.5;
 45        LeptonFinder dressedmuons(muons, photons, 0.1, eta_lep && Cuts::pT >= 6*GeV);
 46        declare(dressedmuons, "dressedmuons");
 47      } else {
 48        declare(HeavyHadrons(Cuts::absrap < 2.4 && Cuts::pT > 15.5*GeV), "BHadrons");
 49      }
 50
 51      //Book the histograms:
 52      bookHandler(_h["dR"],         1);
 53      bookHandler(_h["highpT_dR"],  4);
 54      bookHandler(_h["lowpT_dR"],   7);
 55      bookHandler(_h["dPhi"],      10);
 56      bookHandler(_h["dy"],        13);
 57      bookHandler(_h["MopT"],      16);
 58      bookHandler(_h["pToM"],      19);
 59      bookHandler(_h["pT"],        22);
 60      bookHandler(_h["M"],         25);
 61      bookHandler(_h["yboost"],    29);
 62    }
 63
 64
 65    void bookHandler(HistoHandler& handler, unsigned int id_xsec) {
 66      if (_mode) {
 67        book(handler.histo, "_aux_hist" + toString(id_xsec), refData(id_xsec, 1, 1));
 68        book(handler.scatter, id_xsec, 1, 1);
 69        handler.d = id_xsec + 1; // transfer function
 70        handler.x = 1; handler.y = 1;
 71      }
 72      else  book(handler.histo, id_xsec, 1, 1);
 73    }
 74
 75
 76    /// Perform the per-event analysis
 77    void analyze(const Event& event) {
 78
 79      if (_mode == 1) { // make the 2-b-hadron-level plots
 80        const Particles& bHadrons = apply<HeavyHadrons>(event, "BHadrons").bHadrons();
 81        if (bHadrons.size() > 1) {
 82          sortBy(bHadrons, cmpMomByPt);
 83
 84          float dphiBB = deltaPhi(bHadrons[0], bHadrons[1]);
 85          float dRBB = deltaR(bHadrons[0], bHadrons[1], RAPIDITY);
 86          float dyBB = fabs(bHadrons[0].rapidity() - bHadrons[1].rapidity());
 87          float yboostBB = 0.5*fabs(bHadrons[0].rapidity() + bHadrons[1].rapidity());
 88          FourMomentum systemBB = bHadrons[0].momentum() +  bHadrons[1].momentum();
 89          // Due to the additional particles produced in the decay,
 90          // the 3 muons carry only a fraction of the momentum of the b-hadrons,
 91          // scale down mass and pT to match 3-muon-level more closely
 92          float MBB = systemBB.mass()/1.75;
 93          float pTBB = systemBB.pT()/1.75;
 94
 95          _h["dPhi"].fill(dphiBB);
 96          _h["dy"].fill(dyBB);
 97          _h["yboost"].fill(yboostBB);
 98          _h["dR"].fill(dRBB);
 99          _h["M"].fill(MBB/GeV);
100          _h["pT"].fill(pTBB/GeV);
101          _h["MopT"].fill(MBB/pTBB);
102          _h["pToM"].fill(pTBB/MBB);
103
104          if (pTBB >= 20*GeV)  _h["highpT_dR"].fill(dRBB);
105          else                 _h["lowpT_dR"].fill(dRBB);
106        }
107      }
108
109
110      if (_mode == 0) { // the 3-muon-level analysis
111
112        // First, simply check that we have enough muons
113        const DressedLeptons muons = apply<LeptonFinder>(event, "dressedmuons").dressedLeptons();
114        if (muons.size() < 3)  vetoEvent;
115
116        // Not sure if this is going to work, but ..
117        DressedLeptons Jpsi_muons, third_muons;
118        for (const DressedLepton& mu : muons) {
119	  const Particle& baremu = mu.bareLepton();
120          if (baremu.fromBottom() && baremu.hasAncestorWith(Cuts::pid == PID::JPSI)) {
121            Jpsi_muons.push_back(mu);
122          }
123          else if (baremu.fromBottom()) {
124            third_muons.push_back(mu);
125          }
126        }
127
128        // Veto events without enough muons:
129        if (Jpsi_muons.size() < 2)  vetoEvent;
130
131        // At this point, we must always have a Jpsi. So get its 4-vector:
132        FourMomentum Jpsi = Jpsi_muons[0].momentum() + Jpsi_muons[1].momentum();
133
134        // If there is more than one J/psi, take the one closest to PDG mass,
135        // and push all the other muons back to the 3rd muon list
136        size_t mu1 = 0, mu2 = 1;
137        if (Jpsi_muons.size() > 2) {
138          for (size_t i = 0; i < Jpsi_muons.size(); ++i) {
139            for (size_t j = i; j < Jpsi_muons.size(); ++j) {
140              FourMomentum temp = Jpsi_muons[i].momentum() + Jpsi_muons[j].momentum();
141              if (fabs(temp.mass() - 3.096) < fabs(Jpsi.mass() - 3.096)) {
142                Jpsi = temp;
143                mu1 = i;
144                mu2 = j;
145              }
146            }
147          }
148
149          for (size_t i = 0; i < Jpsi_muons.size(); ++i) {
150            if (i == mu1 || i == mu2)  continue;
151            third_muons.push_back(Jpsi_muons[i]);
152          }
153        }
154
155        // There *is* more than one Jpsi:
156        if (Jpsi_muons[mu1].abseta() >= 2.3) vetoEvent;
157        if (Jpsi_muons[mu2].abseta() >= 2.3) vetoEvent;
158
159        // We should now have the full list of 3rd muons to consider. Make sure we have one:
160        if (third_muons.empty())  vetoEvent;
161
162        // Sort the third muons by pT and pick highest one
163        std::sort(third_muons.begin(), third_muons.end(), [](const DressedLepton &l1, const DressedLepton &l2) {
164          return (l1.pT() > l2.pT());
165        });
166        FourMomentum third_mu = third_muons[0].momentum();
167
168        // Finally, make the plots!
169        float dphi = deltaPhi(Jpsi, third_mu);
170        float dR = deltaR(Jpsi, third_mu, RAPIDITY);
171        float dy = fabs(Jpsi.rapidity() - third_mu.rapidity());
172        float yboost = 0.5*fabs(Jpsi.rapidity() + third_mu.rapidity());
173        FourMomentum system = Jpsi +  third_mu;
174        float M = system.mass();
175        float pT = system.pT();
176
177        _h["dPhi"].fill(dphi);
178        _h["dy"].fill(dy);
179        _h["yboost"].fill(yboost);
180        _h["dR"].fill(dR);
181        if (pT >= 20*GeV)  _h["highpT_dR"].fill(dR);
182        else  _h["lowpT_dR"].fill(dR);
183
184        _h["M"].fill(M);
185        _h["pT"].fill(pT);
186        _h["MopT"].fill(M/pT);
187        _h["pToM"].fill(pT/M);
188
189      } //< end of muon analysis.
190    }
191
192
193    /// Normalise histograms etc., after the run
194    void finalize() {
195      for (auto& hit : _h) {
196        normalize(hit.second.histo);
197        if (_mode == 1)  applyTransferFnAndNorm(hit.second);
198      }
199    }
200
201
202    void applyTransferFnAndNorm(HistoHandler &handler) { ///< @todo Pass as const reference?
203      // Load transfer function from reference data file
204      const YODA::Estimate1D& myTransferFn = refData(handler.d, handler.x, handler.y);
205      double area = 0.0;
206      for (size_t i = 1; i < handler.scatter->numBins()+1; ++i) {
207        const auto& f = myTransferFn.bin(i);
208        auto& p = handler.scatter->bin(i);
209        const auto&  b = handler.histo->bin(i);
210        double newy;
211        try {
212          newy = b.sumW();
213        } catch (const Exception&) { // LowStatsError or WeightError
214          newy = 0;
215        }
216        double newey;
217        try {
218          newey = b.errW();
219        } catch (const Exception&) { // LowStatsError or WeightError
220          newey = 0;
221        }
222        // apply transfer function here
223        newy *= f.val(); newey *= f.val();
224        double rp = safediv(newey, newy);
225        double rf = safediv(f.errAvg(), f.val());
226        newey = newy * sqrt(rp*rp + rf*rf);
227        // set new values
228        p.set(newy, newey);
229        area += newy * (p.xMax() - p.xMin());
230      }
231      if (area > 0.) { // normalise to unity
232        for (size_t i = 1; i < handler.scatter->numBins()+1; ++i) {
233          handler.scatter->bin(i).scale(1.0 / area);
234        }
235      }
236    }
237
238
239  protected:
240
241    /// Analysis-mode switch
242    size_t _mode;
243
244    /// Histograms
245    map<string, HistoHandler> _h;
246
247  };
248
249
250  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1598613);
251
252}