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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
| // -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/WFinder.hh"
#include "Rivet/Projections/ZFinder.hh"
namespace Rivet {
/// @brief Measurement of forward W and Z boson production with jets in pp collisions at 8 TeV
class LHCB_2016_I1454404 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2016_I1454404);
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
_mode = 0;
string mode = getOption("MODE");
if (mode == "ALL" ) _mode = 0;
else if (mode == "WpJET") _mode = 1;
else if (mode == "WmJET") _mode = 2;
else if (mode == "ZJET") _mode = 3;
else if (mode == "WJET") _mode = 4;
const Cut muSel = Cuts::eta >= 2.0 && Cuts::eta <= 4.5 && Cuts::pT > 20*GeV;
// Z boson
ZFinder zfinder(FinalState(), muSel, PID::MUON, 60*GeV, 120*GeV);
declare(zfinder, "ZFinder");
// W boson
WFinder wfinder(FinalState(), muSel, PID::MUON, 0*GeV ,500*GeV ,0.);
declare(wfinder, "WFinder");
// Jet Z
FastJets jetproZ(zfinder.remainingFinalState(), FastJets::ANTIKT, 0.5);
declare(jetproZ, "JetsZ");
// Jet W
FastJets jetproW(wfinder.remainingFinalState(), FastJets::ANTIKT, 0.5);
declare(jetproW, "JetsW");
// Book histograms
/////////
if (_mode == 0 || _mode == 1 || _mode == 4) {
book(_h_wpj, 1, 1, 1);
book(_h_eta_wpj, 4, 1, 1);
book(_h_etaj_wpj, 5, 1, 1);
book(_h_ptj_wpj, 6, 1, 1);
}
/////////
if (_mode == 0 || _mode == 2 || _mode == 4) {
book(_h_wmj, 1, 1, 2);
book(_h_eta_wmj, 4, 1, 2);
book(_h_etaj_wmj, 5, 1, 2);
book(_h_ptj_wmj, 6, 1, 2);
}
/////////
if (_mode == 0 || _mode == 3) {
book(_h_zj, 1, 1, 3);
book(_h_yz_zj, 7, 1, 1);
book(_h_etaj_zj, 8, 1, 1);
book(_h_ptj_zj, 9, 1, 1);
book(_h_dphi_zj, 10, 1, 1);
}
if (_mode == 0 ){
book(_h_rwz, 2,1,1);
book(_h_rwpz, 2,1,2);
book(_h_rwmz, 2,1,3);
}
if (_mode == 0 || _mode == 4){
book(_h_rwpm, 2,1,4);
book(_h_aw, 3,1,1);
// this is a temporary histogram to construct rwz later
book(_h_wj, "_temp_wj", refData(1,1,1));
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const Cut jetSel = Cuts::eta >= 2.2 && Cuts::eta <= 4.2 && Cuts::pT > 20*GeV;
if (_mode == 0 || _mode == 3) {
//////////////////////////////////////////////////////////
///////////////ZFinder Muon //////////////////////////////
//////////////////////////////////////////////////////////
const ZFinder& zfinder = apply<ZFinder>(event, "ZFinder");
if (zfinder.bosons().size() ==1){
const Particles muon = zfinder.constituentLeptons(); //zfinder.constituents()?
const Particles Z = zfinder.bosons();
const FourMomentum Zmom = Z[0].momentum();
const Jets jetsZ = apply<FastJets>(event, "JetsZ").jetsByPt(jetSel);
const Jets cleanedJetsZ = filter_discard(jetsZ, [&](const Jet& j) {return any(muon, deltaRLess(j, 0.5)); });
if (cleanedJetsZ.size() > 0 && cleanedJetsZ.at(0).pT() > 20*GeV) {
const double yZ = Zmom.rap(); //histogram 7
const double etaj = cleanedJetsZ[0].eta(); //histogram 8
const double ptj = cleanedJetsZ[0].pT()/GeV; //histogram 9
double dphi_tmp = abs(Zmom.phi() - cleanedJetsZ[0].phi());
const double dphi = dphi_tmp < Rivet::pi ? dphi_tmp : Rivet::twopi - dphi_tmp;
_h_zj->fill(sqrtS()/GeV);
_h_dphi_zj->fill(dphi);
_h_yz_zj->fill(yZ); // boson rapidity vs diff cross section
_h_etaj_zj->fill(etaj); // jet pseudorapidity vs diff cross section
_h_ptj_zj->fill(ptj); //jet transverse momentum vs diff cross section
}
}
}
if (_mode == 0 || _mode == 1 || _mode == 2 || _mode == 4) {
//////////////////////////////////////////////////////////
///////////////WFinder Muon //////////////////////////////
//////////////////////////////////////////////////////////
const WFinder& wfinder = apply<WFinder>(event, "WFinder");
if (wfinder.bosons().size() == 1) {
const Particles Muons = wfinder.constituentLeptons();
const FourMomentum muonmom = Muons[0].momentum();
const Jets jetsW = apply<FastJets>(event, "JetsW").jetsByPt(jetSel);
const Jets cleanedJetsW = filter_discard(jetsW, [&](const Jet& j) {return any(Muons, deltaRLess(j, 0.5)); });
if (cleanedJetsW.size() > 0 && cleanedJetsW.at(0).pT() > 20*GeV) {
const double etaj = cleanedJetsW[0].eta(); //histogram 5
const double etamu = muonmom.eta(); //histogram 4
if( (_mode == 0 || _mode == 1 || _mode == 4) && Muons[0].charge() > 0) {
//fill with W related analysis
if (_mode != 1 ) _h_wj->fill(sqrtS()/GeV); // don't need this for single charge case
_h_wpj->fill(sqrtS()/GeV);
_h_eta_wpj->fill(etamu); // W+ Jet muon pseudorapidity vs diff cross section
_h_etaj_wpj->fill(etaj); // W+ Jet jet pseudorapidity vs diff cross section
_h_ptj_wpj->fill(cleanedJetsW[0].pT()/GeV); // W+ Jet jet transverse momentum vs diff cross section
}
else if( (_mode == 0 || _mode == 2 || _mode == 4 ) && Muons[0].charge() < 0) {
//fill with W related analysis
if (_mode != 2) _h_wj->fill(sqrtS()/GeV);
_h_wmj->fill(sqrtS()/GeV); // don't need this for single charge case
_h_eta_wmj->fill(etamu); // W- Jet muon pseudorapidity vs diff cross section
_h_etaj_wmj->fill(etaj); // W- Jet jet pseudorapidity vs diff cross section
_h_ptj_wmj->fill(cleanedJetsW[0].pT()/GeV); // W+ Jet jet transverse momentum vs diff cross section
}
}
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
double scalefactor = crossSection()/picobarn/sumOfWeights();
if(_mode == 0 || _mode == 1 || _mode == 4) {
scale({_h_wpj, _h_eta_wpj, _h_etaj_wpj, _h_ptj_wpj}, scalefactor);
}
if(_mode == 0 || _mode == 2 || _mode == 4) {
scale({_h_wmj, _h_eta_wmj, _h_etaj_wmj, _h_ptj_wmj}, scalefactor);
}
if(_mode == 0 || _mode == 3) {
scale({_h_zj, _h_yz_zj, _h_etaj_zj, _h_ptj_zj, _h_dphi_zj}, scalefactor);
}
if (_mode == 0 ) {
scale(_h_wj, scalefactor); // need to scale this for consistency
divide(_h_wpj, _h_zj, _h_rwpz);
divide(_h_wmj, _h_zj, _h_rwmz);
divide(_h_wj, _h_zj, _h_rwz);
}
if (_mode == 0 || _mode == 4) {
divide(_h_wpj, _h_wmj, _h_rwpm);
asymm(_h_wpj, _h_wmj, _h_aw);
}
}
protected:
Log& log = getLog(); // in Analysis or Projection
size_t _mode;
/// histograms
Histo1DPtr _h_wpj, _h_wmj, _h_wj, _h_zj;
Scatter2DPtr _h_rwz, _h_rwpz, _h_rwmz, _h_rwpm, _h_aw;
Histo1DPtr _h_eta_wpj, _h_eta_wmj, _h_etaj_wpj, _h_etaj_wmj, _h_ptj_wpj, _h_ptj_wmj;
Histo1DPtr _h_yz_zj, _h_etaj_zj, _h_ptj_zj, _h_dphi_zj;
};
RIVET_DECLARE_PLUGIN(LHCB_2016_I1454404);
}
|