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
| // -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
class ATLAS_2014_I1298811 : public Analysis {
public:
ATLAS_2014_I1298811()
: Analysis("ATLAS_2014_I1298811") { }
void init() {
// Configure projections
const FinalState fs((Cuts::etaIn(-4.8, 4.8)));
declare(fs, "FS");
const FastJets jets(fs, FastJets::ANTIKT, 0.4);
declare(jets, "Jets");
// Book histograms
for (size_t itopo = 0; itopo < 2; ++itopo) {
// Profiles
for (size_t iregion = 0; iregion < 3; ++iregion) {
book(_p_ptsumch_vs_ptlead[itopo][iregion] ,1+iregion, 1, itopo+1);
book(_p_nch_vs_ptlead[itopo][iregion] ,4+iregion, 1, itopo+1);
}
book(_p_etsum25_vs_ptlead_trans[itopo] ,7, 1, itopo+1);
book(_p_etsum48_vs_ptlead_trans[itopo] ,8, 1, itopo+1);
book(_p_chratio_vs_ptlead_trans[itopo] ,9, 1, itopo+1);
book(_p_ptmeanch_vs_ptlead_trans[itopo] ,10, 1, itopo+1);
// 1D histos
for (size_t iregion = 0; iregion < 3; ++iregion) {
for (size_t ipt = 0; ipt < 4; ++ipt) {
book(_h_ptsumch[ipt][itopo][iregion] ,13+3*ipt+iregion, 1, itopo+1);
book(_h_nch[ipt][itopo][iregion] ,25+3*ipt+iregion, 1, itopo+1);
}
}
}
book(_p_ptmeanch_vs_nch_trans[0], 11, 1, 1);
book(_p_ptmeanch_vs_nch_trans[1], 12, 1, 1);
}
void analyze(const Event& event) {
// Find the jets with pT > 20 GeV and *rapidity* within 2.8
/// @todo Use Cuts instead rather than an eta cut in the proj and a y cut after
const Jets alljets = apply<FastJets>(event, "Jets").jetsByPt(20*GeV);
Jets jets;
for (const Jet& j : alljets)
if (j.absrap() < 2.8) jets.push_back(j);
// Require at least one jet in the event
if (jets.empty()) vetoEvent;
// Identify the leading jet and its phi and pT
const FourMomentum plead = jets[0].momentum();
const double philead = plead.phi();
const double etalead = plead.eta();
const double ptlead = plead.pT();
MSG_DEBUG("Leading object: pT = " << ptlead << ", eta = " << etalead << ", phi = " << philead);
// Sum particle properties in the transverse regions
int tmpnch[2] = {0,0};
double tmpptsum[2] = {0,0};
double tmpetsum48[2] = {0,0};
double tmpetsum25[2] = {0,0};
const Particles particles = apply<FinalState>(event, "FS").particles();
for (const Particle& p : particles) {
// Only consider the transverse region(s), not toward or away
if (!inRange(deltaPhi(p.phi(), philead), PI/3.0, TWOPI/3.0)) continue;
// Work out which transverse side this particle is on
const size_t iside = (mapAngleMPiToPi(p.phi() - philead) > 0) ? 0 : 1;
MSG_TRACE(p.phi() << " vs. " << philead << ": " << iside);
// Charged or neutral particle?
const bool charged = p.charge3() != 0;
// Track observables
if (charged && fabs(p.eta()) < 2.5 && p.pT() > 500*MeV) {
tmpnch[iside] += 1;
tmpptsum[iside] += p.pT();
}
// Cluster observables
if ((charged && p.p3().mod() > 200*MeV) || (!charged && p.p3().mod() > 500*MeV)) {
tmpetsum48[iside] += p.pT();
if (fabs(p.eta()) < 2.5) tmpetsum25[iside] += p.pT();
}
}
// Construct tot/max/min counts (for trans/max/min, indexed by iregion)
const int nch[3] = { tmpnch[0] + tmpnch[1],
std::max(tmpnch[0], tmpnch[1]),
std::min(tmpnch[0], tmpnch[1]) };
const double ptsum[3] = { tmpptsum[0] + tmpptsum[1],
std::max(tmpptsum[0], tmpptsum[1]),
std::min(tmpptsum[0], tmpptsum[1]) };
const double etsum48[3] = { tmpetsum48[0] + tmpetsum48[1],
std::max(tmpetsum48[0], tmpetsum48[1]),
std::min(tmpetsum48[0], tmpetsum48[1]) };
const double etsum25[3] = { tmpetsum25[0] + tmpetsum25[1],
std::max(tmpetsum25[0], tmpetsum25[1]),
std::min(tmpetsum25[0], tmpetsum25[1]) };
//////////////////////////////////////////////////////////
// Now fill the histograms with the computed quantities
// phi sizes of each trans/max/min region (for indexing by iregion)
const double dphi[3] = { 2*PI/3.0, PI/3.0, PI/3.0 };
// Loop over inclusive jet and exclusive dijet configurations
for (size_t itopo = 0; itopo < 2; ++itopo) {
// Exit early if in the exclusive dijet iteration and the exclusive dijet cuts are not met
if (itopo == 1) {
if (jets.size() != 2) continue;
const FourMomentum psublead = jets[1].momentum();
// Delta(phi) cut
const double phisublead = psublead.phi();
if (deltaPhi(philead, phisublead) < 2.5) continue;
// pT fraction cut
const double ptsublead = psublead.pT();
if (ptsublead < 0.5*ptlead) continue;
MSG_DEBUG("Exclusive dijet event");
}
// Plot profiles and distributions which have no max/min region definition
_p_etsum25_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum25[0]/5.0/dphi[0]/GeV);
_p_etsum48_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum48[0]/9.6/dphi[0]/GeV);
if (etsum25[0] > 0) {
_p_chratio_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptsum[0]/etsum25[0]);
}
const double ptmean = safediv(ptsum[0], nch[0], -1); ///< Return -1 if div by zero
if (ptmean >= 0) {
_p_ptmeanch_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptmean/GeV);
_p_ptmeanch_vs_nch_trans[itopo]->fill(nch[0], ptmean/GeV);
}
// Plot remaining profile and 1D observables, which are defined in all 3 tot/max/min regions
for (size_t iregion = 0; iregion < 3; ++iregion) {
_p_ptsumch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, ptsum[iregion]/5.0/dphi[iregion]/GeV);
_p_nch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, nch[iregion]/5.0/dphi[iregion]);
for (size_t ipt = 0; ipt < 4; ++ipt) {
if (ipt == 1 && !inRange(ptlead/GeV, 20, 60)) continue;
if (ipt == 2 && !inRange(ptlead/GeV, 60, 210)) continue;
if (ipt == 3 && ptlead/GeV < 210) continue;
_h_ptsumch[ipt][itopo][iregion]->fill(ptsum[iregion]/5.0/dphi[iregion]/GeV);
_h_nch[ipt][itopo][iregion]->fill(nch[iregion]/5.0/dphi[iregion]);
}
}
}
}
void finalize() {
for (size_t iregion = 0; iregion < 3; ++iregion) {
for (size_t itopo = 0; itopo < 2; ++itopo) {
for (size_t ipt = 0; ipt < 4; ++ipt) {
normalize(_h_ptsumch[ipt][itopo][iregion], 1.0);
normalize(_h_nch[ipt][itopo][iregion], 1.0);
}
}
}
}
private:
/// @name Histogram arrays
//@{
Profile1DPtr _p_ptsumch_vs_ptlead[2][3];
Profile1DPtr _p_nch_vs_ptlead[2][3];
Profile1DPtr _p_ptmeanch_vs_ptlead_trans[2];
Profile1DPtr _p_etsum25_vs_ptlead_trans[2];
Profile1DPtr _p_etsum48_vs_ptlead_trans[2];
Profile1DPtr _p_chratio_vs_ptlead_trans[2];
Profile1DPtr _p_ptmeanch_vs_nch_trans[2];
Histo1DPtr _h_ptsumch[4][2][3];
Histo1DPtr _h_nch[4][2][3];
//@}
};
// The hook for the plugin system
RIVET_DECLARE_PLUGIN(ATLAS_2014_I1298811);
}
|