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
| // -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
namespace Rivet {
class CMS_2015_I1356998 : public Analysis {
public:
CMS_2015_I1356998()
: Analysis("CMS_2015_I1356998"), edge(4.7)
{ }
void init() {
declare(FinalState(),"FS");
book(_h_noCASTORtag ,1, 1, 1);
book(_h_CASTORtag ,2, 1, 1);
book(_h_centralGap ,3, 1, 1);
book(_h_sigmaVis ,4, 1, 1);
book(_h_maxFwdGap ,5, 1, 1);
}
void analyze(const Event& event) {
const double weight = 1.0;
const FinalState& fs = apply<FinalState>(event, "FS");
// A vector containing a lot of eta values
vector<double> detparticles;
detparticles.push_back(-edge);
for (const Particle& p : fs.particles(Cuts::pT > 0.2*GeV && Cuts::abseta<edge, cmpMomByEta) ) {
detparticles.push_back(p.momentum().eta());
}
detparticles.push_back(edge);
// Find maximum gap size
vector <double>::iterator iter;
vector<double> detgaps;
for (iter = detparticles.begin()+1; iter != detparticles.end(); ++iter) {
const double detgap = *iter - *(iter-1);
detgaps.push_back(detgap);
}
double detgapbwd = detgaps.front();
double detgapfwd = detgaps.back();
double detfmax = max(detgapbwd, detgapfwd);
// Fill rapidity gap histo
if (detfmax != 2*edge ) {
_h_maxFwdGap->fill(detfmax, weight);
}
// Everything that follows has to do with the cross-section measurements
if (fs.size() < 2) vetoEvent;
// Gap center calculations
const Particles particlesByRapidity = fs.particles(cmpMomByRap); //ByRapidity();
vector<double> gaps;
vector<double> midpoints;
for (size_t ip = 1; ip < particlesByRapidity.size(); ++ip) {
const Particle& p1 = particlesByRapidity[ip-1];
const Particle& p2 = particlesByRapidity[ip];
const double gap = p2.momentum().rapidity() - p1.momentum().rapidity();
const double mid = (p2.momentum().rapidity() + p1.momentum().rapidity()) / 2.;
gaps.push_back(gap);
midpoints.push_back(mid);
}
int imid = std::distance(gaps.begin(), max_element(gaps.begin(), gaps.end()));
double gapcenter = midpoints[imid];
// Calculations for cross-sections
FourMomentum MxFourVector(0.,0.,0.,0.);
FourMomentum MyFourVector(0.,0.,0.,0.);
for(const Particle& p : fs.particles(cmpMomByEta)) {
if (p.momentum().rapidity() > gapcenter) {
MxFourVector += p.momentum();
}
else {
MyFourVector += p.momentum();
}
}
double Mx = MxFourVector.mass();
double My = MyFourVector.mass();
const double xix = (Mx*Mx)/(sqrtS()/GeV * sqrtS()/GeV);
if (log10(My) < 0.5) {
_h_noCASTORtag->fill(log10(xix), weight);
if (log10(xix) > -5.5 && log10(xix) < -2.5) _h_sigmaVis->fill(0.5, weight);
}
else if (log10(My) < 1.1) {
_h_CASTORtag->fill(log10(xix), weight);
if (log10(xix) > -5.5 && log10(xix) < -2.5) _h_sigmaVis->fill(1.5, weight);
}
// Central gap x-section
double xigen = (Mx*Mx) * (My*My) / (sqrtS()/GeV * sqrtS()/GeV * 0.93827 * 0.93827); // Proton masses...
double dy0 = -log(xigen);
if (dy0 > 3.) {
if (log10(My) > 1.1 && log10(Mx) > 1.1) {
_h_centralGap->fill(dy0, weight);
_h_sigmaVis->fill(2.5, weight);
}
}
}
void finalize() {
double xs = crossSection()/millibarn/sumOfWeights();
scale(_h_noCASTORtag, xs);
scale(_h_CASTORtag , xs);
scale(_h_centralGap , xs);
scale(_h_sigmaVis , xs);
scale(_h_maxFwdGap , xs);
}
private:
Histo1DPtr _h_noCASTORtag;
Histo1DPtr _h_CASTORtag;
Histo1DPtr _h_centralGap;
Histo1DPtr _h_sigmaVis;
Histo1DPtr _h_maxFwdGap;
double edge;
};
RIVET_DECLARE_PLUGIN(CMS_2015_I1356998);
}
|