1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/SmearedJets.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief CDF properties of 6-jet events with large 6-jet mass
 11  class CDF_1997_I442265 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_1997_I442265);
 15
 16
 17    void init() {
 18
 19      // Find true jets
 20      const FinalState fs(Cuts::abseta < 4.2);
 21      FastJets fj(fs, JetAlg::CDFJETCLU, 0.7);
 22
 23      // Smear jet energy and mass with the 10% uncertainty quoted in the paper
 24      SmearedJets sj_E(fj, [](const Jet& jet){ return P4_SMEAR_MASS_GAUSS(P4_SMEAR_E_GAUSS(jet, 0.1*jet.E()), 0.1*jet.mass()); });
 25      declare(sj_E, "Jets");
 26
 27      book(_h_m6J ,1, 1, 1);
 28      book(_h_X3ppp ,2, 1, 1);
 29      book(_h_X4ppp ,3, 1, 1);
 30      book(_h_costheta3ppp ,4, 1, 1);
 31      book(_h_psi3ppp ,5, 1, 1);
 32      book(_h_f3ppp ,6, 1, 1);
 33      book(_h_f4ppp ,6, 1, 2);
 34      book(_h_f5ppp ,6, 1, 3);
 35      book(_h_XApp ,7, 1, 1);
 36      book(_h_XCp ,8, 1, 1);
 37      book(_h_XE ,9, 1, 1);
 38      book(_h_psiAppBpp ,10, 1, 1);
 39      book(_h_psiCpDp ,11, 1, 1);
 40      book(_h_psiEF ,12, 1, 1);
 41      book(_h_fApp ,13, 1, 1);
 42      book(_h_fBpp ,14, 1, 1);
 43      book(_h_fCp ,15, 1, 1);
 44      book(_h_fDp ,16, 1, 1);
 45      book(_h_fE ,17, 1, 1);
 46      book(_h_fF ,18, 1, 1);
 47    }
 48
 49
 50    void analyze(const Event& event) {
 51
 52      const Jets alljets = apply<JetFinder>(event, "Jets").jets(Cuts::Et > 20*GeV && Cuts::abseta < 3, cmpMomByEt);
 53
 54      Jets jets;
 55      double sumEt = 0.0;
 56      FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0);
 57      for (const Jet& jet : alljets) {
 58        double Et = jet.Et();
 59        bool separated = true;
 60        for (const Jet& ref : jets) {
 61          if (deltaR(jet, ref) < 0.9) {
 62            separated = false;
 63            break;
 64          }
 65        }
 66        if (!separated) continue;
 67        jets.push_back(jet);
 68        sumEt += Et;
 69        jetsystem += jet.momentum();
 70        if (jets.size() >= 6) break;
 71      }
 72
 73      if (jets.size() < 6) vetoEvent;
 74      if (sumEt < 320.0*GeV) vetoEvent;
 75
 76      double m6J = _safeMass(jetsystem);
 77      if (m6J < 520.0*GeV) vetoEvent;
 78
 79      if (getLog().isActive(Log::DEBUG)) {
 80        stringstream ss;
 81        ss << "Jets:\n";
 82        for (const Jet& j : jets) ss << j << "\n";
 83        MSG_DEBUG(ss.str());
 84      }
 85
 86      const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec());
 87      vector<FourMomentum> jets6;
 88      for (const Jet& jet : jets) {
 89        jets6.push_back(cms_boost.transform(jet.momentum()));
 90      }
 91      isortBy(jets6, cmpMomByE);
 92
 93      FourMomentum pE, pF;
 94      vector<FourMomentum> jets5(_reduce(jets6, pE, pF));
 95      isortBy(jets5, cmpMomByE);
 96
 97      FourMomentum pCp, pDp;
 98      vector<FourMomentum> jets4(_reduce(jets5, pCp, pDp));
 99      isortBy(jets4, cmpMomByE);
100
101      FourMomentum pApp, pBpp;
102      vector<FourMomentum> jets3(_reduce(jets4, pApp, pBpp));
103      isortBy(jets3, cmpMomByE);
104      FourMomentum p3ppp(jets3[0]);
105      FourMomentum p4ppp(jets3[1]);
106      FourMomentum p5ppp(jets3[2]);
107
108      double X3ppp = 2.0*p3ppp.E()/m6J;
109      if (X3ppp > 0.9) vetoEvent;
110
111      FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m6J, jetsystem.rapidity()));
112      double costheta3ppp = pAV.p3().unit().dot(p3ppp.p3().unit());
113      if (fabs(costheta3ppp) > 0.9) vetoEvent;
114
115      // 3-jet-system variables
116      _h_m6J->fill(m6J);
117      _h_X3ppp->fill(X3ppp);
118      _h_X4ppp->fill(2.0*p4ppp.E()/m6J);
119      _h_costheta3ppp->fill(costheta3ppp);
120      const double psi3ppp = _psi(p3ppp, pAV, p4ppp, p5ppp);
121      _h_psi3ppp->fill(psi3ppp);
122      _h_f3ppp->fill(_safeMass(p3ppp)/m6J);
123      _h_f4ppp->fill(_safeMass(p4ppp)/m6J);
124      _h_f5ppp->fill(_safeMass(p5ppp)/m6J);
125
126      // 4 -> 3 jet variables
127      _h_fApp->fill(_safeMass(pApp)/m6J);
128      _h_fBpp->fill(_safeMass(pApp)/m6J);
129      _h_XApp->fill(pApp.E()/(pApp.E()+pBpp.E()));
130      const double psiAppBpp = _psi(pApp, pBpp, pApp+pBpp, pAV);
131      _h_psiAppBpp->fill(psiAppBpp);
132
133      // 5 -> 4 jet variables
134      _h_fCp->fill(_safeMass(pCp)/m6J);
135      _h_fDp->fill(_safeMass(pDp)/m6J);
136      _h_XCp->fill(pCp.E()/(pCp.E()+pDp.E()));
137      const double psiCpDp = _psi(pCp, pDp, pCp+pDp, pAV);
138      _h_psiCpDp->fill(psiCpDp);
139
140      // 6 -> 5 jet variables
141      _h_fE->fill(_safeMass(pE)/m6J);
142      _h_fF->fill(_safeMass(pF)/m6J);
143      _h_XE->fill(pE.E()/(pE.E()+pF.E()));
144      const double psiEF = _psi(pE, pF, pE+pF, pAV);
145      _h_psiEF->fill(psiEF);
146    }
147
148
149    void finalize() {
150      normalize(_h_m6J);
151      normalize(_h_X3ppp);
152      normalize(_h_X4ppp);
153      normalize(_h_costheta3ppp);
154      normalize(_h_psi3ppp);
155      normalize(_h_f3ppp);
156      normalize(_h_f4ppp);
157      normalize(_h_f5ppp);
158      normalize(_h_XApp);
159      normalize(_h_XCp);
160      normalize(_h_XE);
161      normalize(_h_psiAppBpp);
162      normalize(_h_psiCpDp);
163      normalize(_h_psiEF);
164      normalize(_h_fApp);
165      normalize(_h_fBpp);
166      normalize(_h_fCp);
167      normalize(_h_fDp);
168      normalize(_h_fE);
169      normalize(_h_fF);
170    }
171
172
173  private:
174
175    vector<FourMomentum> _reduce(const vector<FourMomentum>& jets,
176                                 FourMomentum& combined1,
177                                 FourMomentum& combined2) {
178      double minMass2 = 1e9;
179      size_t idx1(jets.size()), idx2(jets.size());
180      for (size_t i = 0; i < jets.size(); ++i) {
181        for (size_t j = i+1; j < jets.size(); ++j) {
182          double mass2 = FourMomentum(jets[i] + jets[j]).mass2();
183          if (mass2 < minMass2) {
184            idx1 = i;
185            idx2 = j;
186          }
187        }
188      }
189      vector<FourMomentum> newjets;
190      for (size_t i = 0; i < jets.size(); ++i) {
191        if (i != idx1 && i != idx2) newjets.push_back(jets[i]);
192      }
193      newjets.push_back(jets[idx1] + jets[idx2]);
194      combined1 = jets[idx1];
195      combined2 = jets[idx2];
196      return newjets;
197    }
198
199
200    FourMomentum _avg_beam_in_lab(const double& m, const double& y) {
201      const double mt = m/2.0;
202      FourMomentum beam1(mt, 0, 0,  mt);
203      FourMomentum beam2(mt, 0, 0, -mt);
204      if (fabs(y) > 1e-3) {
205        FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y));
206        const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(boostvec.betaVec()).inverse();
207        beam1 = cms_boost.transform(beam1);
208        beam2 = cms_boost.transform(beam2);
209      }
210      return (beam1.E() > beam2.E()) ? beam1 - beam2 : beam2 - beam1;
211    }
212
213
214    double _psi(const FourMomentum& p1, const FourMomentum& p2,
215                const FourMomentum& p3, const FourMomentum& p4) {
216      Vector3 p1xp2 = p1.p3().cross(p2.p3());
217      Vector3 p3xp4 = p3.p3().cross(p4.p3());
218      return mapAngle0ToPi(acos(p1xp2.unit().dot(p3xp4.unit())));
219    }
220
221
222    double _safeMass(const FourMomentum& p) {
223      double mass2 = p.mass2();
224      if (mass2 > 0.0) return sqrt(mass2);
225      if (mass2 < -1e-5) MSG_WARNING("m2 = " << m2 << ". Assuming m2=0.");
226      return 0.0;
227    }
228
229
230  private:
231
232    Histo1DPtr _h_m6J;
233    Histo1DPtr _h_X3ppp, _h_X4ppp;
234    Histo1DPtr _h_costheta3ppp;
235    Histo1DPtr _h_psi3ppp;
236    Histo1DPtr _h_f3ppp, _h_f4ppp, _h_f5ppp;
237    Histo1DPtr _h_XApp, _h_XCp, _h_XE;
238    Histo1DPtr _h_psiAppBpp, _h_psiCpDp, _h_psiEF;
239    Histo1DPtr _h_fApp, _h_fBpp, _h_fCp, _h_fDp, _h_fE, _h_fF;
240
241  };
242
243
244
245  RIVET_DECLARE_ALIASED_PLUGIN(CDF_1997_I442265, CDF_1997_S3541940);
246
247}