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