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
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
| // -*- C++ -*-
/**
* @name ATLAS Diffractive Gaps Rivet Analysis
* @author Tim Martin, tim.martin@cern.ch
* @version 1.0
* @date 16/01/2012
* @see http://arxiv.org/abs/1201.2808
* @note pp, sqrt(s) = 7 TeV
* @note Rapidity gap finding algorithm designed to compliment
* the ATLAS detector acceptance. Forward rapidity gap sizes
* are calculated for each event, considering all stable
* particles above pT cut values 200, 400, 600 and 800 MeV in
* turn. A forward rapidity gap is defined to be the largest
* continuous region stretching inward from either edge of the
* detector at eta = |4.9| which contains zero particles above
* pT Cut. Soft diffractive topologies are isolated at large
* gap sizes.
*
*/
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
namespace Rivet {
class ATLAS_2012_I1084540 : public Analysis {
public:
ATLAS_2012_I1084540() : Analysis("ATLAS_2012_I1084540") {}
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
//All final states. Rapidity range = ATLAS calorimetry. Lowest pT cut = 200 MeV.
const FinalState cnfs2((Cuts::etaIn(-_etaMax, _etaMax) && Cuts::pT >= 0.2 * GeV));
const FinalState cnfs4((Cuts::etaIn(-_etaMax, _etaMax) && Cuts::pT >= 0.4 * GeV));
const FinalState cnfs6((Cuts::etaIn(-_etaMax, _etaMax) && Cuts::pT >= 0.6 * GeV));
const FinalState cnfs8((Cuts::etaIn(-_etaMax, _etaMax) && Cuts::pT >= 0.8 * GeV));
declare(cnfs2, "CNFS2");
declare(cnfs4, "CNFS4");
declare(cnfs6, "CNFS6");
declare(cnfs8, "CNFS8");
_etaBinSize = (2. * _etaMax)/(double)_etaBins;
//Book histogram
book(_h_DeltaEtaF_200 ,1, 1, 1);
book(_h_DeltaEtaF_400 ,2, 1, 1);
book(_h_DeltaEtaF_600 ,3, 1, 1);
book(_h_DeltaEtaF_800 ,4, 1, 1);
}
private:
void fillMap(const FinalState& fs, bool* energyMap, double pTcut) {
// Fill true/false array by iterating over particles and compare their
// pT with pTcut
for (const Particle& p : fs.particles(cmpMomByEta)) {
int checkBin = -1;
double checkEta = -_etaMax;
while (1) {
checkEta += _etaBinSize;
++checkBin;
if (p.eta() < checkEta) {
energyMap[checkBin] = (p.pT() > pTcut * GeV);
break;
}
}
}
}
public:
/// Perform the per-event analysis
void analyze(const Event& event) {
static unsigned int event_count = 0;
++event_count;
const FinalState& fs2 = apply<FinalState>(event, "CNFS2");
const FinalState& fs4 = apply<FinalState>(event, "CNFS4");
const FinalState& fs6 = apply<FinalState>(event, "CNFS6");
const FinalState& fs8 = apply<FinalState>(event, "CNFS8");
// Set up Yes/No arrays for energy in each eta bin at each pT cut
bool energyMap_200[_etaBins];
bool energyMap_400[_etaBins];
bool energyMap_600[_etaBins];
bool energyMap_800[_etaBins];
for (int i = 0; i < _etaBins; ++i) {
energyMap_200[i] = false;
energyMap_400[i] = false;
energyMap_600[i] = false;
energyMap_800[i] = false;
}
// Veto bins based on final state particles > Cut (Where Cut = 200 - 800 MeV pT)
fillMap(fs2, energyMap_200, 0.2);
fillMap(fs4, energyMap_400, 0.4);
fillMap(fs6, energyMap_600, 0.6);
fillMap(fs8, energyMap_800, 0.8);
// Apply gap finding algorithm
// Detector layout follows...
// -Eta [Proton --- DetectorCSide --- DetectorBarrel --- DetectorASide --- Proton] +Eta
bool gapDirectionAt200 = false; //False is gap on C size, True is gap on A side.
double largestEdgeGap_200 = 0.;
double largestEdgeGap_400 = 0.;
double largestEdgeGap_600 = 0.;
double largestEdgeGap_800 = 0.;
for (int E = 200; E <= 800; E += 200) {
double EdgeGapSizeA = -1, EdgeGapSizeC = -1;
bool* energyMap = 0;
switch (E) {
case 200: energyMap = energyMap_200; break;
case 400: energyMap = energyMap_400; break;
case 600: energyMap = energyMap_600; break;
case 800: energyMap = energyMap_800; break;
}
// Look left to right
for (int a = 0; a < _etaBins; ++a) {
if (energyMap[a] == true) {
EdgeGapSizeA = (_etaBinSize * a);
break;
}
}
// And look right to left
for (int c = _etaBins-1; c >= 0; --c) {
if (energyMap[c] == true) {
EdgeGapSizeC = (2 * _etaMax) - (_etaBinSize * (c+1));
if (fuzzyEquals(EdgeGapSizeC, 4.47035e-08)) EdgeGapSizeC = 0.0;
break;
}
}
// Put your hands on your hips
// Find the largest gap
double largestEdgeGap = 0.;
if (E == 200) {
// If the 200 MeV pass, take the biggest of the two gaps. Make note of which side for higher pT cuts.
largestEdgeGap = std::max(EdgeGapSizeA,EdgeGapSizeC);
gapDirectionAt200 = (EdgeGapSizeA > EdgeGapSizeC);
} else {
// Use the direction from 200 MeV pass, most accurate measure of which side gap is on.
if (gapDirectionAt200) {
largestEdgeGap = EdgeGapSizeA;
}
else largestEdgeGap = EdgeGapSizeC;
}
// Check case of empty detector
if (largestEdgeGap < 0.0) largestEdgeGap = 2.0 * _etaMax;
// Fill bin centre
switch (E) {
case 200: _h_DeltaEtaF_200->fill(largestEdgeGap + _etaBinSize/2.); break;
case 400: _h_DeltaEtaF_400->fill(largestEdgeGap + _etaBinSize/2.); break;
case 600: _h_DeltaEtaF_600->fill(largestEdgeGap + _etaBinSize/2.); break;
case 800: _h_DeltaEtaF_800->fill(largestEdgeGap + _etaBinSize/2.); break;
}
if (E == 200) largestEdgeGap_200 = largestEdgeGap;
if (E == 400) largestEdgeGap_400 = largestEdgeGap;
if (E == 600) largestEdgeGap_600 = largestEdgeGap;
if (E == 800) largestEdgeGap_800 = largestEdgeGap;
}
// Algorithm result every 1000 events
if (event_count % 1000 == 0) {
for (int E = 200; E <= 800; E += 200) {
bool* energyMap = 0;
double largestEdgeGap = 0;
switch (E) {
case 200: energyMap = energyMap_200; largestEdgeGap = largestEdgeGap_200; break;
case 400: energyMap = energyMap_400; largestEdgeGap = largestEdgeGap_400; break;
case 600: energyMap = energyMap_600; largestEdgeGap = largestEdgeGap_600; break;
case 800: energyMap = energyMap_800; largestEdgeGap = largestEdgeGap_800; break;
}
MSG_DEBUG("Largest Forward Gap at pT Cut " << E << " MeV=" << largestEdgeGap
<< " eta, NFinalState pT > 200 in ATLAS acceptance:" << fs2.particles().size());
std::string hitPattern = "Detector HitPattern=-4.9[";
for (int a = 0; a < _etaBins; ++a) {
if (energyMap[a] == true) hitPattern += "X";
else hitPattern += "_";
}
hitPattern += "]4.9";
MSG_DEBUG(hitPattern);
std::string gapArrow = " ";
if (!gapDirectionAt200) {
int drawSpaces = (int)(_etaBins - (largestEdgeGap/_etaBinSize) + 0.5);
for (int i = 0; i < drawSpaces; ++i) gapArrow += " ";
}
int drawArrows = (int)((largestEdgeGap/_etaBinSize) + 0.5);
for (int i = 0; i < drawArrows; ++i) gapArrow += "^";
MSG_DEBUG(gapArrow);
}
}
}
/// Normalise histograms after the run, Scale to cross section
void finalize() {
MSG_DEBUG("Cross Section=" << crossSection() / millibarn << "mb, SumOfWeights=" << sumOfWeights());
scale(_h_DeltaEtaF_200, (crossSection() / millibarn)/sumOfWeights());
scale(_h_DeltaEtaF_400, (crossSection() / millibarn)/sumOfWeights());
scale(_h_DeltaEtaF_600, (crossSection() / millibarn)/sumOfWeights());
scale(_h_DeltaEtaF_800, (crossSection() / millibarn)/sumOfWeights());
}
//@}
private:
/// @name Histograms
//@{
Histo1DPtr _h_DeltaEtaF_200;
Histo1DPtr _h_DeltaEtaF_400;
Histo1DPtr _h_DeltaEtaF_600;
Histo1DPtr _h_DeltaEtaF_800;
//@}
/// @name Private variables
//@{
static constexpr int _etaBins = 49;
static constexpr double _etaMax = 4.9;
double _etaBinSize;
//@}
};
// The hook for the plugin system
RIVET_DECLARE_PLUGIN(ATLAS_2012_I1084540);
}
|