1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
| // -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Projections/InvMassFinalState.hh"
namespace Rivet {
/// CMS W + 2 jet double parton scattering analysis
class CMS_2013_I1272853 : public Analysis {
public:
/// Constructor
CMS_2013_I1272853()
: Analysis("CMS_2013_I1272853") { }
/// Book histograms and initialise projections before the run
void init() {
const FinalState fs;
declare(fs, "FS");
/// @todo Use C++11 initialisation syntax
vector<PdgIdPair> vidsW;
vidsW += make_pair(PID::MUON, PID::NU_MUBAR), make_pair(PID::ANTIMUON, PID::NU_MU);
InvMassFinalState invfsW(fs, vidsW, 20*GeV, 1e6*GeV);
declare(invfsW, "INVFSW");
VetoedFinalState vfs(fs);
vfs.addVetoOnThisFinalState(invfsW);
declare(vfs, "VFS");
declare(FastJets(vfs, FastJets::ANTIKT, 0.5), "Jets");
book(_h_deltaS_eq2jet_Norm ,1,1,1);
book(_h_rel_deltaPt_eq2jet_Norm ,2,1,1);
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// Find Ws
const InvMassFinalState& invMassFinalStateW = apply<InvMassFinalState>(event, "INVFSW");
if (invMassFinalStateW.empty()) vetoEvent;
const Particles& WDecayProducts = invMassFinalStateW.particles();
if (WDecayProducts.size() < 2) vetoEvent;
// Cuts on W decay properties
const int iNU_MU = (WDecayProducts[1].abspid() == PID::NU_MU) ? 1 : 0;
const int iAN_MU = 1 - iNU_MU;
const double pt1 = WDecayProducts[iAN_MU].pT();
const double pt2 = WDecayProducts[iNU_MU].Et();
const double eta1 = WDecayProducts[iAN_MU].abseta();
const double phi1 = WDecayProducts[iAN_MU].phi();
const double phi2 = WDecayProducts[iNU_MU].phi();
const double mt = sqrt(2 * pt1 * pt2 * (1 - cos(phi1-phi2)));
if (mt < 50*GeV || pt1 < 35*GeV || eta1 > 2.1 || pt2 < 30*GeV) vetoEvent;
// Get jets and make sure there are at least two of them in |y| < 2
const FastJets& jetpro = apply<FastJets>(event, "Jets");
/// @todo Collapse this into jetpro.jetsByPt(ptGtr(20*GeV) & rapIn(2.0))
vector<FourMomentum> jets;
for (const Jet& jet : jetpro.jetsByPt(20*GeV))
if (jet.absrap() < 2.0) jets.push_back(jet.momentum());
if (jets.size() != 2) vetoEvent;
const double mupx = pt1 * cos(phi1);
const double mupy = pt1 * sin(phi1);
const double met_x = pt2 * cos(phi2);
const double met_y = pt2 * sin(phi2);
const double dpt = add_quad(jets[0].px() + jets[1].px(), jets[0].py() + jets[1].py());
const double rel_dpt = dpt / (jets[0].pT() + jets[1].pT());
const double pT2 = sqr(mupx + met_x) + sqr(mupy + met_y);
const double Px = (mupx + met_x)*(jets[0].px() + jets[1].px());
const double Py = (mupy + met_y)*(jets[0].py() + jets[1].py());
const double p1p2_mag = dpt * sqrt(pT2);
const double dS = acos((Px+Py) / p1p2_mag);
const double weight = 1.0;
_h_rel_deltaPt_eq2jet_Norm->fill(rel_dpt, weight);
_h_deltaS_eq2jet_Norm->fill(dS, weight);
}
/// Normalise histograms etc., after the run
void finalize() {
const double rel_dpt_bw = 1.0002 / 30.0;
const double dphi_bw = 3.14160 / 30.0;
normalize(_h_rel_deltaPt_eq2jet_Norm, rel_dpt_bw);
normalize(_h_deltaS_eq2jet_Norm, dphi_bw);
}
private:
Histo1DPtr _h_rel_deltaPt_eq2jet_Norm;
Histo1DPtr _h_deltaS_eq2jet_Norm;
};
RIVET_DECLARE_PLUGIN(CMS_2013_I1272853);
}
|