1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/IdentifiedFinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  8#include "Rivet/Projections/UnstableParticles.hh"
  9
 10namespace Rivet {
 11
 12
 13  class ATLAS_2012_I1082009 : public Analysis {
 14  public:
 15
 16    /// @name Constructors etc.
 17    //@{
 18
 19    /// Constructor
 20    ATLAS_2012_I1082009()
 21      : Analysis("ATLAS_2012_I1082009")
 22    {    }
 23
 24    //@}
 25
 26
 27  public:
 28
 29    /// @name Analysis methods
 30    //@{
 31
 32    /// Book histograms and initialise projections before the run
 33    void init() {
 34
 35      // Input for the jets: No neutrinos, no muons
 36      VetoedFinalState veto;
 37      veto.addVetoPairId(PID::MUON);
 38      veto.vetoNeutrinos();
 39      FastJets jets(veto, FastJets::ANTIKT, 0.6);
 40      declare(jets, "jets");
 41      // unstable final-state for D*
 42      declare(UnstableParticles(), "UFS");
 43
 44      book(_weight25_30, "_weight_25_30");
 45      book(_weight30_40, "_weight_30_40");
 46      book(_weight40_50, "_weight_40_50");
 47      book(_weight50_60, "_weight_50_60");
 48      book(_weight60_70, "_weight_60_70");
 49      book(_weight25_70, "_weight_25_70");
 50
 51      book(_h_pt25_30 , 8,1,1);
 52      book(_h_pt30_40 , 9,1,1);
 53      book(_h_pt40_50 ,10,1,1);
 54      book(_h_pt50_60 ,11,1,1);
 55      book(_h_pt60_70 ,12,1,1);
 56      book(_h_pt25_70 ,13,1,1);
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      // get the jets
 63      Jets jets;
 64      for (const Jet& jet : apply<FastJets>(event, "jets").jetsByPt(25.0*GeV)) {
 65        if ( jet.abseta() < 2.5 ) jets.push_back(jet);
 66      }
 67      // get the D* mesons
 68      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 69      Particles Dstar;
 70      for (const Particle& p : ufs.particles()) {
 71        const int id = p.abspid();
 72        if(id==413) Dstar.push_back(p);
 73      }
 74
 75      // loop over the jobs
 76      for (const Jet& jet : jets ) {
 77        double perp = jet.perp();
 78        bool found = false;
 79        double z(0.);
 80        if(perp<25.||perp>70.) continue;
 81        for(const Particle & p : Dstar) {
 82          if(p.perp()<7.5) continue;
 83          if(deltaR(p, jet.momentum())<0.6) {
 84            Vector3 axis = jet.p3().unit();
 85            z = axis.dot(p.p3())/jet.E();
 86            if(z<0.3) continue;
 87            found = true;
 88            break;
 89          }
 90        }
 91        _weight25_70->fill();
 92        if(found) _h_pt25_70->fill(z);
 93        if(perp>=25.&&perp<30.) {
 94          _weight25_30->fill();
 95          if(found) _h_pt25_30->fill(z);
 96        }
 97        else if(perp>=30.&&perp<40.) {
 98          _weight30_40->fill();
 99          if(found) _h_pt30_40->fill(z);
100        }
101        else if(perp>=40.&&perp<50.) {
102          _weight40_50->fill();
103          if(found) _h_pt40_50->fill(z);
104        }
105        else if(perp>=50.&&perp<60.) {
106          _weight50_60->fill();
107          if(found) _h_pt50_60->fill(z);
108        }
109        else if(perp>=60.&&perp<70.) {
110          _weight60_70->fill();
111          if(found) _h_pt60_70->fill(z);
112        }
113      }
114    }
115
116
117    /// Normalise histograms etc., after the run
118    void finalize() {
119      scale(_h_pt25_30,1./ *_weight25_30);
120      scale(_h_pt30_40,1./ *_weight30_40);
121      scale(_h_pt40_50,1./ *_weight40_50);
122      scale(_h_pt50_60,1./ *_weight50_60);
123      scale(_h_pt60_70,1./ *_weight60_70);
124      scale(_h_pt25_70,1./ *_weight25_70);
125    }
126
127    //@}
128
129
130  private:
131
132    /// @name Histograms
133    //@{
134    CounterPtr _weight25_30,_weight30_40,_weight40_50;
135    CounterPtr _weight50_60,_weight60_70,_weight25_70;
136
137    Histo1DPtr _h_pt25_30;
138    Histo1DPtr _h_pt30_40;
139    Histo1DPtr _h_pt40_50;
140    Histo1DPtr _h_pt50_60;
141    Histo1DPtr _h_pt60_70;
142    Histo1DPtr _h_pt25_70;
143    //@}
144
145  };
146
147  // The hook for the plugin system
148  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1082009);
149
150}