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 high-mass multi-jet events
 11  class CDF_1996_I418504 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_1996_I418504);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      /// Initialise and register projections here
 24      const FinalState fs(Cuts::abseta < 4.2);
 25      FastJets fj(fs, JetAlg::CDFJETCLU, 0.7);
 26      declare(fj, "Jets");
 27
 28      // Smear energy and mass with the 10% uncertainty quoted in the paper
 29      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()); });
 30      declare(sj_E, "SmearedJets");
 31
 32      /// Book histograms here, e.g.:
 33      book(_h_3_mNJ ,1, 1, 1);
 34      book(_h_3_X3 ,2, 1, 1);
 35      book(_h_3_X4 ,3, 1, 1);
 36      book(_h_3_costheta3 ,8, 1, 1);
 37      book(_h_3_psi3 ,9, 1, 1);
 38      book(_h_3_f3 ,14, 1, 1);
 39      book(_h_3_f4 ,14, 1, 2);
 40      book(_h_3_f5 ,14, 1, 3);
 41
 42      book(_h_4_mNJ ,1, 1, 2);
 43      book(_h_4_X3 ,4, 1, 1);
 44      book(_h_4_X4 ,5, 1, 1);
 45      book(_h_4_costheta3 ,10, 1, 1);
 46      book(_h_4_psi3 ,11, 1, 1);
 47      book(_h_4_f3 ,15, 1, 1);
 48      book(_h_4_f4 ,15, 1, 2);
 49      book(_h_4_f5 ,15, 1, 3);
 50      book(_h_4_XA ,17, 1, 1);
 51      book(_h_4_psiAB ,19, 1, 1);
 52      book(_h_4_fA ,21, 1, 1);
 53      book(_h_4_fB ,21, 1, 2);
 54
 55      book(_h_5_mNJ ,1, 1, 3);
 56      book(_h_5_X3 ,6, 1, 1);
 57      book(_h_5_X4 ,7, 1, 1);
 58      book(_h_5_costheta3 ,12, 1, 1);
 59      book(_h_5_psi3 ,13, 1, 1);
 60      book(_h_5_f3 ,16, 1, 1);
 61      book(_h_5_f4 ,16, 1, 2);
 62      book(_h_5_f5 ,16, 1, 3);
 63      book(_h_5_XA ,18, 1, 1);
 64      book(_h_5_XC ,18, 1, 2);
 65      book(_h_5_psiAB ,20, 1, 1);
 66      book(_h_5_psiCD ,20, 1, 2);
 67      book(_h_5_fA ,22, 1, 1);
 68      book(_h_5_fB ,23, 1, 1);
 69      book(_h_5_fC ,24, 1, 1);
 70      book(_h_5_fD ,25, 1, 1);
 71    }
 72
 73
 74    void analyze(const Event& event) {
 75      Jets jets;
 76      FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0);
 77      for (const Jet& jet : apply<JetFinder>(event, "SmearedJets").jets(Cuts::Et > 20.0*GeV, cmpMomByEt)) {
 78        bool separated = true;
 79        for (const Jet& ref : jets) {
 80          if (deltaR(jet, ref) < 0.9) {
 81            separated = false;
 82            break;
 83          }
 84        }
 85        if (!separated) continue;
 86        jets.push_back(jet);
 87        jetsystem += jet.momentum();
 88        if (jets.size() >= 5) break;
 89      }
 90
 91      if (jets.size() > 4) {
 92        _fiveJetAnalysis(jets);
 93        jets.resize(4);
 94      }
 95      if (jets.size() > 3) {
 96        _fourJetAnalysis(jets);
 97        jets.resize(3);
 98      }
 99      if (jets.size() > 2) {
100        _threeJetAnalysis(jets);
101      }
102    }
103
104
105    void _threeJetAnalysis(const Jets& jets) {
106      MSG_DEBUG("3 jet analysis");
107
108      double sumEt = 0.0;
109      FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0);
110      for (const Jet& jet : jets) {
111        sumEt += jet.Et();
112        jetsystem += jet.momentum();
113      }
114      if (sumEt < 420.0*GeV) return;
115
116      const double m3J = _safeMass(jetsystem);
117      if (m3J < 600*GeV) return;
118
119      const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec());
120      vector<FourMomentum> jets3;
121      for (Jet jet : jets) {
122        jets3.push_back(cms_boost.transform(jet.momentum()));
123      }
124      isortBy(jets3, cmpMomByE);
125      FourMomentum p3(jets3[0]), p4(jets3[1]), p5(jets3[2]);
126
127      FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m3J, jetsystem.rapidity()));
128      double costheta3 = pAV.p3().unit().dot(p3.p3().unit());
129      if (fabs(costheta3) > 0.6) return;
130
131      double X3 = 2.0*p3.E()/m3J;
132      if (X3 > 0.9) return;
133
134      const double X4 = 2.0*p4.E()/m3J;
135      const double psi3 = _psi(p3, pAV, p4, p5);
136      const double f3 = _safeMass(p3)/m3J;
137      const double f4 = _safeMass(p4)/m3J;
138      const double f5 = _safeMass(p5)/m3J;
139
140      _h_3_mNJ->fill(m3J);
141      _h_3_X3->fill(X3);
142      _h_3_X4->fill(X4);
143      _h_3_costheta3->fill(costheta3);
144      _h_3_psi3->fill(psi3);
145      _h_3_f3->fill(f3);
146      _h_3_f4->fill(f4);
147      _h_3_f5->fill(f5);
148
149    }
150
151
152    void _fourJetAnalysis(const Jets& jets) {
153      MSG_DEBUG("4 jet analysis");
154
155      double sumEt=0.0;
156      FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0);
157      for (const Jet& jet : jets) {
158        sumEt+=jet.Et();
159        jetsystem+=jet.momentum();
160      }
161      if (sumEt < 420.0*GeV) return;
162
163      const double m4J = _safeMass(jetsystem);
164      if (m4J < 650*GeV) return;
165
166      const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec());
167      vector<FourMomentum> jets4;
168      for (Jet jet : jets) {
169        jets4.push_back(cms_boost.transform(jet.momentum()));
170      }
171      isortBy(jets4, cmpMomByE);
172
173      FourMomentum pA, pB;
174      vector<FourMomentum> jets3(_reduce(jets4, pA, pB));
175      isortBy(jets3, cmpMomByE);
176      FourMomentum p3(jets3[0]);
177      FourMomentum p4(jets3[1]);
178      FourMomentum p5(jets3[2]);
179
180      FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m4J, jetsystem.rapidity()));
181      double costheta3=pAV.p3().unit().dot(p3.p3().unit());
182      if (fabs(costheta3)>0.8) {
183        return;
184      }
185
186      const double X3 = 2.0*p3.E()/m4J;
187      if (X3>0.9) {
188        return;
189      }
190
191      // fill histograms
192      const double X4 = 2.0*p4.E()/m4J;
193      const double psi3 = _psi(p3, pAV, p4, p5);
194      const double f3 = _safeMass(p3)/m4J;
195      const double f4 = _safeMass(p4)/m4J;
196      const double f5 = _safeMass(p5)/m4J;
197      const double fA = _safeMass(pA)/m4J;
198      const double fB = _safeMass(pB)/m4J;
199      const double XA = pA.E()/(pA.E()+pB.E());
200      const double psiAB = _psi(pA, pB, pA+pB, pAV);
201
202      _h_4_mNJ->fill(m4J);
203      _h_4_X3->fill(X3);
204      _h_4_X4->fill(X4);
205      _h_4_costheta3->fill(costheta3);
206      _h_4_psi3->fill(psi3);
207      _h_4_f3->fill(f3);
208      _h_4_f4->fill(f4);
209      _h_4_f5->fill(f5);
210      _h_4_XA->fill(XA);
211      _h_4_psiAB->fill(psiAB);
212      _h_4_fA->fill(fA);
213      _h_4_fB->fill(fB);
214    }
215
216
217    void _fiveJetAnalysis(const Jets& jets) {
218      MSG_DEBUG("5 jet analysis");
219
220      double sumEt=0.0;
221      FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0);
222      for (const Jet& jet : jets) {
223        sumEt+=jet.Et();
224        jetsystem+=jet.momentum();
225      }
226      if (sumEt < 420.0*GeV) return;
227
228      const double m5J = _safeMass(jetsystem);
229      if (m5J < 750*GeV) return;
230
231      const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec());
232      vector<FourMomentum> jets5;
233      for (Jet jet : jets) {
234        jets5.push_back(cms_boost.transform(jet.momentum()));
235      }
236      isortBy(jets5, cmpMomByE);
237
238      FourMomentum pC, pD;
239      vector<FourMomentum> jets4(_reduce(jets5, pC, pD));
240      isortBy(jets4, cmpMomByE);
241
242      FourMomentum pA, pB;
243      vector<FourMomentum> jets3(_reduce(jets4, pA, pB));
244      isortBy(jets3, cmpMomByE);
245      FourMomentum p3(jets3[0]);
246      FourMomentum p4(jets3[1]);
247      FourMomentum p5(jets3[2]);
248
249      // fill histograms
250      FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m5J, jetsystem.rapidity()));
251      const double costheta3 = pAV.p3().unit().dot(p3.p3().unit());
252      const double X3 = 2.0*p3.E()/m5J;
253      const double X4 = 2.0*p4.E()/m5J;
254      const double psi3 = _psi(p3, pAV, p4, p5);
255      const double f3 = _safeMass(p3)/m5J;
256      const double f4 = _safeMass(p4)/m5J;
257      const double f5 = _safeMass(p5)/m5J;
258      const double fA = _safeMass(pA)/m5J;
259      const double fB = _safeMass(pB)/m5J;
260      const double XA = pA.E()/(pA.E()+pB.E());
261      const double psiAB = _psi(pA, pB, pA+pB, pAV);
262      const double fC = _safeMass(pC)/m5J;
263      const double fD = _safeMass(pD)/m5J;
264      const double XC = pC.E()/(pC.E()+pD.E());
265      const double psiCD = _psi(pC, pD, pC+pD, pAV);
266
267      _h_5_mNJ->fill(m5J);
268      _h_5_X3->fill(X3);
269      _h_5_X4->fill(X4);
270      _h_5_costheta3->fill(costheta3);
271      _h_5_psi3->fill(psi3);
272      _h_5_f3->fill(f3);
273      _h_5_f4->fill(f4);
274      _h_5_f5->fill(f5);
275      _h_5_XA->fill(XA);
276      _h_5_psiAB->fill(psiAB);
277      _h_5_fA->fill(fA);
278      _h_5_fB->fill(fB);
279      _h_5_XC->fill(XC);
280      _h_5_psiCD->fill(psiCD);
281      _h_5_fC->fill(fC);
282      _h_5_fD->fill(fD);
283    }
284
285
286    /// Normalise histograms etc., after the run
287    void finalize() {
288
289      /// Normalise, scale and otherwise manipulate histograms here
290      normalize(_h_3_mNJ);
291      normalize(_h_3_X3);
292      normalize(_h_3_X4);
293      normalize(_h_3_costheta3);
294      normalize(_h_3_psi3);
295      normalize(_h_3_f3);
296      normalize(_h_3_f4);
297      normalize(_h_3_f5);
298
299      normalize(_h_4_mNJ);
300      normalize(_h_4_X3);
301      normalize(_h_4_X4);
302      normalize(_h_4_costheta3);
303      normalize(_h_4_psi3);
304      normalize(_h_4_f3);
305      normalize(_h_4_f4);
306      normalize(_h_4_f5);
307      normalize(_h_4_XA);
308      normalize(_h_4_psiAB);
309      normalize(_h_4_fA);
310      normalize(_h_4_fB);
311
312      normalize(_h_5_mNJ);
313      normalize(_h_5_X3);
314      normalize(_h_5_X4);
315      normalize(_h_5_costheta3);
316      normalize(_h_5_psi3);
317      normalize(_h_5_f3);
318      normalize(_h_5_f4);
319      normalize(_h_5_f5);
320      normalize(_h_5_XA);
321      normalize(_h_5_XC);
322      normalize(_h_5_psiAB);
323      normalize(_h_5_psiCD);
324      normalize(_h_5_fA);
325      normalize(_h_5_fB);
326      normalize(_h_5_fC);
327      normalize(_h_5_fD);
328
329    }
330
331    /// @}
332
333
334  private:
335
336    vector<FourMomentum> _reduce(const vector<FourMomentum>& jets, FourMomentum& combined1, FourMomentum& combined2) {
337      double minMass2 = 1e9;
338      size_t idx1(jets.size()), idx2(jets.size());
339      for (size_t i=0; i<jets.size(); ++i) {
340        for (size_t j=i+1; j<jets.size(); ++j) {
341          double mass2 = FourMomentum(jets[i]+jets[j]).mass2();
342          if (mass2<minMass2) {
343            idx1=i;
344            idx2=j;
345          }
346        }
347      }
348      vector<FourMomentum> newjets;
349      for (size_t i=0; i<jets.size(); ++i) {
350        if (i!=idx1 && i!=idx2) newjets.push_back(jets[i]);
351      }
352      newjets.push_back(jets[idx1]+jets[idx2]);
353      combined1 = jets[idx1];
354      combined2 = jets[idx2];
355      return newjets;
356    }
357
358
359    FourMomentum _avg_beam_in_lab(const double& m, const double& y) {
360      const double mt = m/2.0;
361      FourMomentum beam1(mt, 0, 0, mt);
362      FourMomentum beam2(mt, 0, 0, -mt);
363      if (fabs(y)>1e-3) {
364        FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y));
365        const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(boostvec.betaVec()).inverse();
366        beam1 = cms_boost.transform(beam1);
367        beam2 = cms_boost.transform(beam2);
368      }
369      return (beam1.E() > beam2.E()) ? beam1-beam2 : beam2-beam1;
370    }
371
372
373    double _psi(const FourMomentum& p1, const FourMomentum& p2,
374                const FourMomentum& p3, const FourMomentum& p4) {
375      Vector3 p1xp2 = p1.p3().cross(p2.p3());
376      Vector3 p3xp4 = p3.p3().cross(p4.p3());
377      return mapAngle0ToPi(acos(p1xp2.unit().dot(p3xp4.unit())));
378    }
379
380
381    double _safeMass(const FourMomentum& p) {
382      double mass2=p.mass2();
383      if (mass2>0.0) return sqrt(mass2);
384      else if (mass2<-1.0e-5) {
385        MSG_WARNING("m2 = " << m2 << ". Assuming m2=0.");
386        return 0.0;
387      }
388      else return 0.0;
389    }
390
391
392  private:
393
394    /// @name Histograms
395    /// @{
396    Histo1DPtr _h_3_mNJ;
397    Histo1DPtr _h_3_X3;
398    Histo1DPtr _h_3_X4;
399    Histo1DPtr _h_3_costheta3;
400    Histo1DPtr _h_3_psi3;
401    Histo1DPtr _h_3_f3;
402    Histo1DPtr _h_3_f4;
403    Histo1DPtr _h_3_f5;
404
405    Histo1DPtr _h_4_mNJ;
406    Histo1DPtr _h_4_X3;
407    Histo1DPtr _h_4_X4;
408    Histo1DPtr _h_4_costheta3;
409    Histo1DPtr _h_4_psi3;
410    Histo1DPtr _h_4_f3;
411    Histo1DPtr _h_4_f4;
412    Histo1DPtr _h_4_f5;
413    Histo1DPtr _h_4_XA;
414    Histo1DPtr _h_4_psiAB;
415    Histo1DPtr _h_4_fA;
416    Histo1DPtr _h_4_fB;
417
418    Histo1DPtr _h_5_mNJ;
419    Histo1DPtr _h_5_X3;
420    Histo1DPtr _h_5_X4;
421    Histo1DPtr _h_5_costheta3;
422    Histo1DPtr _h_5_psi3;
423    Histo1DPtr _h_5_f3;
424    Histo1DPtr _h_5_f4;
425    Histo1DPtr _h_5_f5;
426    Histo1DPtr _h_5_XA;
427    Histo1DPtr _h_5_XC;
428    Histo1DPtr _h_5_psiAB;
429    Histo1DPtr _h_5_psiCD;
430    Histo1DPtr _h_5_fA;
431    Histo1DPtr _h_5_fB;
432    Histo1DPtr _h_5_fC;
433    Histo1DPtr _h_5_fD;
434    /// @}
435
436  };
437
438
439
440  RIVET_DECLARE_ALIASED_PLUGIN(CDF_1996_I418504, CDF_1996_S3349578);
441
442}