Rivet analyses referenceCMS_2014_I1305624Study of hadronic event-shape variables in multijet final states in $pp$ collisions at $\sqrt{s} = 7$ TeVExperiment: CMS (LHC) Inspire ID: 1305624 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
Event-shape variables, which are sensitive to perturbative and nonperturbative aspects of quantum chromodynamic (QCD) interactions, are studied in multijet events recorded in proton-proton collisions at $\sqrt{s}=7$ TeV. Events are selected with at least one jet with transverse momentum $p_T > 110$ GeV and pseudorapidity $\mid\eta\mid < 2.4$, in a data sample corresponding to integrated luminosities of up to $5 fb^{-1}$. The distributions of five event-shape variables in various leading jet $p_T$ ranges are compared to predictions from different QCD Monte Carlo event generators. Five event-shape variables are analyzed in this paper: the transverse thrust $\tau_{\perp}$, the total jet broadening $B_{tot}$, the total jet mass $\rho_{tot}$, the total transverse jet mass $\rho^T_{tot}$ and the third-jet resolution parameter $Y_{23}$. In the formulae below, $p_{T,i}$, $\eta_i$, and $\phi_i$ represent the transverse momentum, pseudorapidity, and azimuthal angle of the $i$th jet, and $\hat{n}_{T}$ is the unit vector that maximizes the sum of the projections of $\vec{p}_{T,i}$. The transverse thrust axis $\hat{n}_{T}$ and the beam form the so-called event plane. Based on the direction of $\hat{n}_{T}$, the transverse region is separated into an upper side $\cal{C}_U$, consisting of all jets with $\vec{p}_{T}\cdot\hat{n}_{T}$ $>$ 0, and a lower side $\cal{C}_L$, with $\vec{p}_{T}\cdot\hat{n}_{T}$ $<$ 0. The jet broadening and third-jet resolution variables require at least three jets, whereas the calculation of other variables requires at least two jets. The $\hat{n}_{T}$ vector is defined only up to a global sign - choosing one sign or the other has no consequence since it simply exchanges the upper and lower events regions. Transverse Thrust : The event thrust observable in the transverse plane is defined as \begin{eqnarray} \tau_{\perp} & \equiv & 1 - \max_{\hat{n}_{T}} \frac{\sum_i |\vec{p}_{T,i} \cdot \hat{n}_{T} |}{\sum_i p_{T,i}} \end{eqnarray} This variable probes the hadronization process and is sensitive to the modeling of two-jet and multijet topologies. In this paper ``multijet' refers to ``more-than-two-jet'. In the limit of a perfectly balanced two-jet event, $\tau_{\perp}$ is zero, while in isotropic multijet events it amounts to $(1-2/\pi)$. Jet Broadenings : The pseudorapidities and the azimuthal angles of the axes for the upper and lower event regions are defined by \begin{eqnarray} \eta_X & \equiv & \frac{\sum_{i\in{\cal{C}}_X} p_{T,i}\eta_i}{\sum_{i\in{\cal{C}}_X} p_{T,i}} , \\ \phi_X & \equiv & \frac{\sum_{i\in{\cal{C}}_X} p_{T,i}\phi_i}{\sum_{i\in{\cal{C}}_X} p_{T,i}} , \end{eqnarray} where $X$ refers to upper ($U$) or lower ($L$) side. From these, the jet broadening variable in each region is defined as \begin{eqnarray} B_{X} & \equiv & \frac{1}{2P_{T}} \sum_{i\in{\cal{C}}_X}p_{T,i}\sqrt{(\eta_i -\eta_X)^2 + (\phi_i - \phi_X)^2} , \end{eqnarray} where $P_{T}$ is the scalar sum of the transverse momenta of all the jets. The total jet broadening is then defined as \begin{eqnarray} B_{tot} & \equiv & B_{U} + B_{L}. \end{eqnarray} Jet Masses : The normalized squared invariant mass of the jets in the upper and lower regions of the event is defined by \begin{eqnarray} \rho_X & \equiv & \frac{M^2_X}{P^2}, \end{eqnarray} where $M_X$ is the invariant mass of the constituents of the jets in the region $X$, and $P$ is the scalar sum of the momenta of all constituents in both sides. The jet mass variable is defined as the sum of the masses in the upper and lower regions, \begin{eqnarray} \rho_{tot} & \equiv & \rho_U + \rho_L . \end{eqnarray} The corresponding jet mass in the transverse plane, $\rho^T_{tot}$, is also similarly calculated in transverse plane. Third-jet resolution parameter : The third-jet resolution parameter is defined as \begin{eqnarray} Y_{23} \equiv \frac{\mathrm{min}(p_{T,3}^2,[\mathrm{min}(p_{T,i}, p_{T,j})^2 \times (\Delta R_{ij})^2/R^2])}{P_{12}^2} , \end{eqnarray} where i, j run over all three jets, $(\Delta R_{ij})^2 = (\eta_i - \eta_j)^2 + (\phi_i - \phi_j)^2$, and $p_{T,3}$ is the transverse momentum of the third jet in the event. If there are more than three jets in the event, they are iteratively merged using the $k_T$ algorithm with a distance parameter $R = 0.6$. To compute $P_{12}$, three jets are merged into two using the procedure described above and $P_{12}$ is then defined as the scalar sum of the transverse momenta of the two remaining jets. The $Y_{23}$ variable estimates the relative strength of the $p_T$ of the third jet with respect to the other two jets. It vanishes for two-jet events, and a nonzero value of $Y_{23}$ indicates the presence of hard parton emission, which tests the parton showering model of QCD event generators. A test like this is less sensitive to the details of the underlying event (UE) and parton hadronization models than the other event-shape variables. Source code: CMS_2014_I1305624.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/FastJets.hh"
5
6namespace Rivet {
7
8
9 namespace {
10
11 /// Number of event shape variables
12 /// @todo Move into the EventShape class
13 const int NEVTVAR = 5;
14
15 /// Number of leading jet pT thresholds
16 /// @todo Move into the analysis class
17 const int NJETPTMN = 5;
18 /// Leading jet pT thresholds
19 /// @todo Move into the analysis class
20 const double LEADINGPTTHRESHOLD[NJETPTMN] = { 110.0, 170.0, 250.0, 320.0, 390.0 };
21
22
23 // Helpers for event shape calculations in hidden namespace; implementation at bottom of file
24 /// @todo Why a class? Improve/remove this junk
25 class EventShape {
26 public:
27
28 /// Constructor from vectors of four-vectors as input objects in the event to calculate the event shapes
29 EventShape(const vector<double>& px_vector, const vector<double>& py_vector, const vector<double>& pz_vector,
30 const vector<double>& e_vector, double eta_central, int irap, int nmn)
31 : _object_px(px_vector), _object_py(py_vector), _object_pz(pz_vector),
32 _object_e(e_vector), _eta_c(eta_central), _irap(irap), _nmnjet(nmn)
33 { }
34
35 /// @brief Returns the values of the five event shapes
36 ///
37 /// Event shape indices:
38 /// 0. central transverse thrust
39 /// 1. central total jet broadening
40 /// 2. central total jet mass
41 /// 3. central total transverse jet mass
42 /// 4. central three-jet resolution threshold
43 vector<double> getEventShapes() {
44 _calculate(); ///< @todo There should be some test for success/failure!!
45 return _event_shapes;
46 }
47
48 /// Returns the global thrust axis Nx, Ny, Nz=0
49 vector<double> getThrustAxis() {
50 _calculate(); ///< @todo There should be some test for success/failure!!
51 return _thrust_axis;
52 }
53
54 /// Returns the central thrust axis Nx, Ny, Nz=0
55 vector<double> getThrustAxisC() {
56 _calculate(); ///< @todo There should be some test for success/failure!!
57 return _thrust_axis_c;
58 }
59
60 // /// @brief Choice of the central region
61 // void setEtaC(double eta_central) { _eta_c = eta_central; }
62
63 // // Whether to use the rapidity y (rap==1) or the pseudorapidity eta (rap==0)
64 // void setRapType(int irap) { _irap = irap; }
65
66
67 private:
68
69 /// Calculate everything
70 int _calculate();
71
72 /// Returns the difference in phi between two vectors
73 double _delta_phi(double, double);
74
75 /// The Lorentz scalar product
76 double _lorentz_sp(const vector<double>&, const vector<double>&);
77
78 // Calculates the three-jet resolutions
79 double _three_jet_res(const vector<double>&, const vector<double>&, const vector<double>&, const vector<double>&, int);
80
81 // Calculates the thrust axis and the tau values
82 vector<double> _thrust(const vector<double>&, const vector<double>&);
83
84
85 vector<double> _object_px, _object_py, _object_pz, _object_p;
86 vector<double> _object_pt, _object_e, _object_phi, _object_eta;
87 vector<double> _event_shapes;
88 vector<double> _thrust_axis, _thrust_axis_c;
89
90 double _eta_c;
91 int _irap;
92 size_t _nmnjet;
93
94 };
95
96 }
97
98
99
100
101 class CMS_2014_I1305624 : public Analysis {
102 public:
103
104 /// Constructor
105 CMS_2014_I1305624()
106 : Analysis("CMS_2014_I1305624")
107 { }
108
109
110 /// @name Analysis methods
111
112 /// Book histograms and initialise projections before the run
113 void init() {
114 const FastJets jets(FinalState(Cuts::abseta < 2.6), JetAlg::ANTIKT, 0.5);
115 declare(jets, "Jets");
116
117 for (int ij=0; ij < NJETPTMN; ij++) {
118 book(_h_thrustc[ij] ,1, 1, ij+1);
119 book(_h_broadt[ij] ,1, 2, ij+1);
120 book(_h_tot3dmass[ij] ,1, 3, ij+1);
121 book(_h_tottrnsmass[ij] ,1, 4, ij+1);
122 book(_h_y23c[ij] ,1, 5, ij+1);
123 //
124 }
125 _needBinInit = true;
126 }
127
128
129 /// Perform the per-event analysis
130 void analyze(const Event& event) {
131
132 if (_needBinInit) {
133 for (int ij=0; ij < NJETPTMN; ij++) {
134 _alow1[ij] = _h_thrustc[ij]->xMin();
135 _alow2[ij] = _h_broadt[ij]->xMin();
136 _alow3[ij] = _h_tot3dmass[ij]->xMin();
137 _alow4[ij] = _h_tottrnsmass[ij]->xMin();
138 _alow5[ij] = _h_y23c[ij]->xMin();
139 //
140 _ahgh1[ij] = _h_thrustc[ij]->xMax();
141 _ahgh2[ij] = _h_broadt[ij]->xMax();
142 _ahgh3[ij] = _h_tot3dmass[ij]->xMax();
143 _ahgh4[ij] = _h_tottrnsmass[ij]->xMax();
144 _ahgh5[ij] = _h_y23c[ij]->xMax();
145
146 _needBinInit = false;
147 }
148 }
149
150 const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 30*GeV);
151 if (jets.size() < 2) vetoEvent;
152 if (jets[0].abseta() > 2.4 || jets[1].abseta() > 2.4) vetoEvent;
153
154 const double leadingpt = jets[0].pT();
155 if (leadingpt < 110*GeV) vetoEvent;
156
157 vector<double> jtpx, jtpy, jtpz, jten;
158 for (const Jet& j : jets) {
159 if (j.abseta() < 2.4) {
160 jtpx.push_back(j.px());
161 jtpy.push_back(j.py());
162 jtpz.push_back(j.pz());
163 jten.push_back(j.E());
164 }
165 }
166
167 EventShape eventshape(jtpx, jtpy, jtpz, jten, 2.4, 0, 2);
168 const vector<double> eventvar = eventshape.getEventShapes();
169 if (eventvar[NEVTVAR] < 0) vetoEvent; // Jets are not only one hemisphere
170
171 for (int ij = NJETPTMN-1; ij >= 0; --ij) {
172 if (leadingpt/GeV > LEADINGPTTHRESHOLD[ij]) {
173 if (inRange(eventvar[0], _alow1[ij], _ahgh1[ij])) _h_thrustc[ij]->fill(eventvar[0]);
174 if (inRange(eventvar[2], _alow3[ij], _ahgh3[ij])) _h_tot3dmass[ij]->fill(eventvar[2]);
175 if (inRange(eventvar[3], _alow4[ij], _ahgh4[ij])) _h_tottrnsmass[ij]->fill(eventvar[3]);
176 if (eventvar[NEVTVAR] >= 3) {
177 if (inRange(eventvar[1], _alow2[ij], _ahgh2[ij])) _h_broadt[ij]->fill(eventvar[1]);
178 if (inRange(eventvar[4], _alow5[ij], _ahgh5[ij])) _h_y23c[ij]->fill(eventvar[4]);
179 }
180 break;
181 }
182 }
183
184 }
185
186
187 /// Normalise histograms etc., after the run
188 void finalize() {
189 for (int ij = 0; ij < NJETPTMN; ij++) {
190 normalize(_h_thrustc[ij]);
191 normalize(_h_broadt[ij]);
192 normalize(_h_tot3dmass[ij]);
193 normalize(_h_tottrnsmass[ij]);
194 normalize(_h_y23c[ij]);
195 }
196 }
197
198
199 private:
200
201 /// @name Histograms
202 /// @{
203 Histo1DPtr _h_thrustc[NJETPTMN];
204 Histo1DPtr _h_broadt[NJETPTMN];
205 Histo1DPtr _h_tot3dmass[NJETPTMN];
206 Histo1DPtr _h_tottrnsmass[NJETPTMN];
207 Histo1DPtr _h_y23c[NJETPTMN];
208 /// @}
209
210 // Data members
211 bool _needBinInit;
212 double _alow1[NJETPTMN], _alow2[NJETPTMN], _alow3[NJETPTMN], _alow4[NJETPTMN], _alow5[NJETPTMN];
213 double _ahgh1[NJETPTMN], _ahgh2[NJETPTMN], _ahgh3[NJETPTMN], _ahgh4[NJETPTMN], _ahgh5[NJETPTMN];
214
215 };
216
217
218 RIVET_DECLARE_PLUGIN(CMS_2014_I1305624);
219
220
221
222 /////////////////////
223
224
225 namespace {
226
227 // EventShape helper class method implementations:
228
229 int EventShape::_calculate() {
230 if (!_event_shapes.empty() && !_thrust_axis.empty() && !_thrust_axis_c.empty())
231 return 1; //< return success if this appears to already have been run
232
233 const size_t length = (size_t) _object_px.size();
234
235 if (((size_t) _object_py.size() != length) ||
236 ((size_t) _object_pz.size() != length) ||
237 ((size_t) _object_e.size() != length)) {
238 /// @todo Change to exception or assert
239 // cout << "ERROR!!!! Input vectors differ in size! Change that please!" << '\n';
240 // cout<<"py_size: "<<_object_py.size()<<" ,pz_size: "<<_object_pz.size()
241 // <<" ,px_size: "<<_object_px.size()<<" ,E_size: "<<_object_e.size()<<'\n';
242 return 0;
243 }
244
245 if (!_object_p.empty()) {
246 _object_p.clear();
247 _object_pt.clear();
248 _object_eta.clear();
249 _object_phi.clear();
250 _event_shapes.clear();
251 _thrust_axis.clear();
252 _thrust_axis_c.clear();
253 }
254
255 for (size_t j = 0; j < length; j++) {
256 _object_p.push_back(0.);
257 _object_pt.push_back(0.);
258 _object_eta.push_back(0.);
259 _object_phi.push_back(0.);
260 }
261
262 for (int j = 0; j < NEVTVAR; j++) {
263 _event_shapes.push_back(-50.);
264 }
265
266 _event_shapes.push_back(double(_object_px.size())); //< WTF?
267
268 for (size_t j = 0; j < 3; j++) {
269 _thrust_axis.push_back(0.);
270 _thrust_axis_c.push_back(0.);
271 }
272
273 double theta = 0;
274
275 for (size_t k = 0; k < length; k++) {
276 _object_p[k] = sqrt(pow(_object_px[k],2) + pow(_object_py[k],2) + pow(_object_pz[k],2));
277 _object_pt[k] = sqrt(pow(_object_px[k],2) + pow(_object_py[k],2));
278 if (_object_p[k] > _object_e[k] + 1e-4) {
279 /// @todo Change to exception or assert
280 // cout << "ERROR!!! object " << k <<" has P = " << _object_p[k]
281 // << " which is bigger than E = " << _object_e[k] <<" "
282 // << _object_px[k] <<" "<< _object_py[k] <<" "
283 // << _object_pz[k] <<" of total length "<< length
284 // << '\n';
285 return 0;
286 }
287
288 //to prevent a division by zero
289 if (_irap == 0) {
290 if (fabs(_object_pz[k]) > 1e-5) {
291 theta = atan(_object_pt[k]/(_object_pz[k]));
292 } else {
293 theta = M_PI/2;
294 }
295 if (theta < 0.) theta = theta + M_PI;
296 _object_eta[k] = -log(tan(0.5*theta));
297 }
298 if (_irap == 1) {
299 if (_object_pz[k] == _object_e[k]) {
300 /// @todo Change to exception
301 // cout << "ERROR!!! object "<<k<<" has Pz "<< _object_pz[k] <<" which is equal to E = "<< _object_e[k] <<'\n';
302 return 0;
303 }
304 _object_eta[k]=0.5*log((_object_e[k]+_object_pz[k])/(_object_e[k]-_object_pz[k]));
305 }
306 if (_irap != 0 && _irap != 1) {
307 /// @todo Change to exception
308 // cout << "ERROR!!!, The choice to use the rapidity y or the pseudorapidity eta is not set correctly! Change that please!" << '\n';
309 return 0;
310 }
311 _object_phi[k] = atan2(_object_py[k], _object_px[k]);
312 }
313
314 vector<double> object_px_in, object_py_in, object_pz_in, object_pt_in, object_e_in, object_et_in, object_eta_in;
315 vector<double> object_px_out, object_py_out, object_pz_out, object_e_out, object_pt_out, object_eta_out;
316 if (!object_px_in.empty()) { //< FFS, this is impossible: it's only just been created!
317 object_px_in.clear();
318 object_py_in.clear();
319 object_pz_in.clear();
320 object_pt_in.clear();
321 object_e_in.clear();
322 object_et_in.clear();
323 object_eta_in.clear();
324 object_px_out.clear();
325 object_py_out.clear();
326 object_pz_out.clear();
327 object_pt_out.clear();
328 object_e_out.clear();
329 object_eta_out.clear();
330 }
331
332 size_t nin = 0;
333
334 for (size_t j = 0; j < length; j++) {
335 if (fabs(_object_eta[j]) < _eta_c) {
336 object_px_in.push_back(_object_px[j]);
337 object_py_in.push_back(_object_py[j]);
338 object_pz_in.push_back(_object_pz[j]);
339 object_e_in.push_back(_object_e[j]);
340 object_pt_in.push_back(sqrt(pow(_object_px[j],2)+pow(_object_py[j],2)));
341 object_et_in.push_back(sqrt((pow(_object_e[j],2)*pow(_object_pt[j],2))/(pow(_object_pt[j],2)+pow(_object_pz[j],2))));
342 object_eta_in.push_back(_object_eta[j]);
343 nin += 1;
344 } else {
345 object_px_out.push_back(_object_px[j]);
346 object_py_out.push_back(_object_py[j]);
347 object_pz_out.push_back(_object_pz[j]);
348 object_e_out.push_back(_object_e[j]);
349 object_pt_out.push_back(sqrt(pow(_object_px[j],2)+pow(_object_py[j],2)));
350 object_eta_out.push_back(_object_eta[j]);
351 }
352 }
353
354 if (object_px_in.size() != nin) {
355 /// @todo Change to exception or assert
356 cout << "ERROR!!! wrong dimension of 'in' momenta" << endl;
357 //return 0; ///< @todo Why not do this?
358 }
359 // const size_t nout = length - nin;
360
361 if (nin < _nmnjet) {
362 for (int i = 0; i < NEVTVAR; i++) {
363 _event_shapes[i] = -50.0;
364 }
365 }
366
367 _event_shapes[NEVTVAR] = nin;
368
369 if (nin >= _nmnjet) {
370 double p_sum_c = 0; //GMA
371 double pt_sum_c = 0;
372 // double eta_cw = 0;
373 // double px_sum_in = 0;
374 // double py_sum_in = 0;
375 for (size_t j = 0; j < nin; j++) {
376 pt_sum_c += object_pt_in[j];
377 p_sum_c += sqrt(pow(object_pt_in[j],2.) + pow(object_pz_in[j], 2.0)); //GMA
378 // eta_cw += object_pt_in[j]*object_eta_in[j];
379 // px_sum_in += object_px_in[j];
380 // py_sum_in += object_py_in[j];
381 }
382 // eta_cw /= pt_sum_c;
383
384 // double expTerm = 0;
385 // for (size_t j = 0; j < nout; j++) {
386 // expTerm += object_pt_out[j] * exp(-fabs(object_eta_out[j]-eta_cw));
387 // }
388 // expTerm /= pt_sum_c;
389
390 //the central global transverse thrust centrthr is calculated
391 double centrthr = 0;
392 vector<double> thrust_central = _thrust(object_px_in, object_py_in);
393
394 for (size_t l=0; l<3; l++) _thrust_axis_c[l] = thrust_central[l];
395 //the variable which gets resummed is not thrust
396 //but tau = 1 - thrust - see calculation
397 centrthr = thrust_central[3];
398 _event_shapes[0] = centrthr;
399
400 double alpha_c = atan2(_thrust_axis_c[1], _thrust_axis_c[0]);
401 //central jet masses
402 //define two jet masses in region U and D
403 double cenjm_up = 0;
404 double cenjm_down= 0;
405 double dot_product = 0;
406
407 vector<double> up_sum;
408 vector<double> down_sum;
409 for (size_t j=0; j<4;j++) {
410 up_sum.push_back(0.);
411 down_sum.push_back(0.);
412 }
413 for (size_t i=0;i<nin;i++) {
414 dot_product = object_px_in[i] * _thrust_axis_c[0] + object_py_in[i] * _thrust_axis_c[1];
415 if (dot_product >= 0) {
416 up_sum[0]+=object_px_in[i];
417 up_sum[1]+=object_py_in[i];
418 up_sum[2]+=object_pz_in[i];
419 up_sum[3]+=object_e_in[i];
420 } else {
421 down_sum[0]+=object_px_in[i];
422 down_sum[1]+=object_py_in[i];
423 down_sum[2]+=object_pz_in[i];
424 down_sum[3]+=object_e_in[i];
425 }
426 }
427 cenjm_up = _lorentz_sp(up_sum, up_sum) / pow(p_sum_c, 2.); //GMA pow(pt_sum_c,2);
428 cenjm_down = _lorentz_sp(down_sum, down_sum) / pow(p_sum_c, 2.); //GMA pow(pt_sum_c,2);
429
430 //central total jet mass centotjm
431 double centotjm=0;
432 centotjm = cenjm_up + cenjm_down;
433
434 _event_shapes[2]=centotjm;
435
436 double centrjm_up=0, centrjm_down=0;
437 vector<double> upsum;
438 vector<double> downsum;
439 for (size_t j = 0; j < 3; j++) {
440 upsum.push_back(0.);
441 downsum.push_back(0.);
442 }
443 for (size_t i = 0; i < nin; i++) {
444 dot_product = object_px_in[i]*_thrust_axis_c[0]+object_py_in[i]*_thrust_axis_c[1];
445 if (dot_product >= 0) {
446 upsum[0] += object_px_in[i];
447 upsum[1] += object_py_in[i];
448 upsum[2] += object_et_in[i];
449 } else {
450 downsum[0] += object_px_in[i];
451 downsum[1] += object_py_in[i];
452 downsum[2] += object_et_in[i];
453 }
454 }
455 centrjm_up = _lorentz_sp(upsum, upsum) / pow(pt_sum_c, 2);
456 centrjm_down = _lorentz_sp(downsum, downsum) / pow(pt_sum_c, 2);
457 double centottrjm = centrjm_up + centrjm_down;
458
459 _event_shapes[3] = centottrjm;
460
461 //central three-jet resolution threshold
462 double ceny3=0;
463 if (nin < 3) {
464 ceny3 = -1.0;
465 } else {
466 ceny3 = _three_jet_res(object_px_in, object_py_in, object_pz_in, object_e_in, _irap);
467 }
468
469 _event_shapes[4] = ceny3;
470
471 //the central jet broadenings in the up and down region
472 double cenbroad_up=0;
473 double cenbroad_down=0;
474
475 double eta_up=0;
476 size_t num_up=0;
477 double eta_down =0;
478 size_t num_down =0;
479 double phi_temp =0;
480 double phi_up_aver =0;
481 double phi_down_aver =0;
482 double pt_sum_up =0;
483 double pt_sum_down =0;
484 double dot_product_b =0;
485 vector<double> phi_up;
486 vector<double> phi_down;
487 double py_rot =0;
488 double px_rot =0;
489
490 for (size_t j = 0; j < 4; j++) {
491 up_sum.push_back(0.);
492 down_sum.push_back(0.);
493 }
494
495 for (size_t i=0;i<nin;i++) {
496 dot_product_b =sqrt(object_px_in[i]*_thrust_axis_c[0] + object_py_in[i]*_thrust_axis_c[1]);
497 if (dot_product_b>=0){
498 pt_sum_up += object_pt_in[i];
499 //rotate the coordinate system so that
500 //the central thrust axis is e_x
501 px_rot = cos(alpha_c)*object_px_in[i]+sin(alpha_c)*object_py_in[i];
502 py_rot = - sin(alpha_c)*object_px_in[i]+cos(alpha_c)*object_py_in[i];
503 //calculate the eta and phi in the rotated system
504 eta_up += object_pt_in[i]*object_eta_in[i];
505 phi_temp = atan2(py_rot,px_rot);
506
507 if(phi_temp > M_PI/2){
508 phi_temp = phi_temp - M_PI/2;
509 }
510 if (phi_temp < -M_PI/2){
511 phi_temp = phi_temp + M_PI/2;
512 }
513 phi_up.push_back(phi_temp);
514 phi_up_aver += object_pt_in[i]*phi_temp;
515 num_up += 1;
516 } else {
517 eta_down += object_pt_in[i]*object_eta_in[i];
518 pt_sum_down += object_pt_in[i];
519 px_rot = cos(alpha_c)*object_px_in[i]+sin(alpha_c)*object_py_in[i];
520 py_rot = - sin(alpha_c)*object_px_in[i]+cos(alpha_c)*object_py_in[i];
521 phi_temp = atan2(py_rot,px_rot);
522 if (phi_temp > M_PI/2) {
523 //if phi is bigger than pi/2 in the new system calculate
524 //the difference to the thrust axis
525 phi_temp = M_PI -phi_temp;
526 }
527 if (phi_temp<-M_PI/2) {
528 //if phi is smaller than
529 phi_temp = -M_PI-phi_temp;
530 }
531 phi_down.push_back(phi_temp);
532 //calculate the pt-weighted phi
533 phi_down_aver += object_pt_in[i]*phi_temp;
534 num_down += 1;
535 }
536 }
537 if (num_up!=0){
538 eta_up = eta_up/pt_sum_up;
539 phi_up_aver = phi_up_aver/pt_sum_up;
540 }
541 if (num_down!=0) {
542 eta_down = eta_down/pt_sum_down;
543 phi_down_aver = phi_down_aver/pt_sum_down;
544 }
545
546 size_t index_up=0, index_down=0;
547 for (size_t i = 0; i < nin; i++) {
548 dot_product_b = object_px_in[i]*_thrust_axis_c[0] + object_py_in[i]*_thrust_axis_c[1];
549 if (dot_product_b >= 0) {
550 //calculate the broadenings of the regions with the rotated system
551 //and the pt-weighted average of phi in the rotated system
552 cenbroad_up += object_pt_in[i]*sqrt(pow(object_eta_in[i]-eta_up, 2) +
553 pow(_delta_phi(phi_up[index_up], phi_up_aver), 2));
554 index_up += 1;
555 } else {
556 cenbroad_down += object_pt_in[i]*sqrt(pow(object_eta_in[i]-eta_down, 2)+
557 pow(_delta_phi(phi_down[index_down], phi_down_aver), 2));
558 index_down += 1;
559 }
560 }
561
562 if (index_up == 0 || index_down ==0) _event_shapes[NEVTVAR] *= -1;
563
564 cenbroad_up=cenbroad_up/(2*pt_sum_c);
565 cenbroad_down=cenbroad_down/(2*pt_sum_c);
566
567 //central total jet broadening
568 double centotbroad = 0;
569 centotbroad = cenbroad_up + cenbroad_down;
570
571 _event_shapes[1] = centotbroad;
572
573 for (int ij = 0; ij < 5; ij++) {
574 if (_event_shapes[ij] < 1.e-20) _event_shapes[ij] = 1.e-20;
575 _event_shapes[ij] = log(_event_shapes[ij]);
576 }
577 }
578
579 return 1;
580 }
581
582
583 double EventShape::_three_jet_res(const vector<double>& in_object_px, const vector<double>& in_object_py, const vector<double>& in_object_pz, const vector<double>& in_object_e, int irap) {
584
585 size_t y3_length = (size_t)in_object_px.size();
586 if (((size_t) in_object_py.size()!=y3_length) ||
587 ((size_t) in_object_pz.size()!=y3_length) ||
588 (in_object_e.size()!=y3_length)) {
589 return 0.0;
590 }
591
592 vector<double> in_object_p, in_object_pt, in_object_eta, in_object_phi;
593 if (!in_object_p.empty()) {
594 in_object_p.clear();
595 in_object_pt.clear();
596 in_object_eta.clear();
597 in_object_phi.clear();
598 }
599 for (size_t j = 0; j < y3_length; j++) {
600 in_object_p.push_back(0.);
601 in_object_pt.push_back(0.);
602 in_object_eta.push_back(0.);
603 in_object_phi.push_back(0.);
604 }
605 double theta_y3_1st = 0;
606 for (size_t k =0; k<y3_length; k++) {
607 in_object_p[k] = sqrt(pow(in_object_px[k],2) + pow(in_object_py[k],2) + pow(in_object_pz[k],2));
608 in_object_pt[k] = sqrt(pow(in_object_px[k],2) + pow(in_object_py[k],2));
609
610 //calculates the pseudorapidity to prevent a division by zero
611 if (irap == 0) {
612 if (fabs(in_object_pz[k]) > 1E-5) {
613 theta_y3_1st = atan(in_object_pt[k]/(in_object_pz[k]));
614 } else {
615 theta_y3_1st = M_PI/2;
616 }
617 if (theta_y3_1st<0.) theta_y3_1st = theta_y3_1st + M_PI;
618 in_object_eta[k] = - log(tan(0.5*theta_y3_1st));
619 }
620 //calculates the real rapidity
621 if (irap == 1) {
622 in_object_eta[k]=0.5*log((in_object_e[k]+in_object_pz[k])/(in_object_e[k]-in_object_pz[k]));
623 }
624 in_object_phi[k] = atan2(in_object_py[k], in_object_px[k]);
625 }
626
627 //the three-jet resolution
628 //threshold y3
629 double y3 = 0;
630
631 //vector which will be filled with the
632 //minimum of the distances
633 double max_dmin_temp=0;
634
635 double max_dmin = 0;
636
637 //distance input object k, beam
638 double distance_jB = 0;
639 double distance_jB_min = 0;
640 //distance of input object k to l
641 double distance_jk = 0;
642 double distance_jk_min = 0;
643 //as we search the minimum of the distances
644 //give them values which are for sure higher
645 //than those we evaluate first in the for-loups
646
647 size_t index_jB = 0;
648 size_t index_j_jk = 0;
649 size_t index_k_jk = 0;
650
651 //to decide later if the minmum is a jB or jk
652 int decide_jB = -1;
653
654 vector<double> input_pt, input_px, input_py, input_pz;
655 vector<double> input_p, input_e, input_phi, input_eta;
656
657 if (!input_pt.empty()) {
658 input_pt.clear();
659 input_px.clear();
660 input_px.clear();
661 input_pz.clear();
662 input_p.clear();
663 input_e.clear();
664 input_phi.clear();
665 input_eta.clear();
666 }
667
668 for (size_t j = 0; j < y3_length; j++){
669 input_pt.push_back(in_object_pt[j]);
670 input_px.push_back(in_object_px[j]);
671 input_py.push_back(in_object_py[j]);
672 input_pz.push_back(in_object_pz[j]);
673 input_p.push_back(in_object_p[j]);
674 input_e.push_back(in_object_e[j]);
675 input_phi.push_back(in_object_phi[j]);
676 input_eta.push_back(in_object_eta[j]);
677 }
678 if (y3_length<3) {
679 return -1;
680 } else {
681 size_t rest = y3_length;
682 for (size_t i = 0; i<y3_length; i++) {
683 //make the minima at the initialization step
684 //of each looping bigger than the first values
685 distance_jB_min = 0.36*pow(input_pt[0],2) + 10;
686 //DELTA PHIs wanted not the pure difference
687 distance_jk_min = min(pow(input_pt[1], 2), pow(input_pt[0], 2)) *
688 (pow(input_eta[1]-input_eta[0], 2) +
689 pow(_delta_phi(input_phi[1], input_phi[0]), 2)) + 10;
690 //do the procedure only until we have only 2 objects left anymore
691 if (rest > 2) {
692 for (size_t j=0; j<rest;j++) {
693 //calculate the distance between object j and the beam
694 distance_jB = 0.36*pow(input_pt[j], 2);
695 if(distance_jB < distance_jB_min){
696 distance_jB_min = distance_jB;
697 index_jB = j;
698 }
699 if (j > 0) {
700 for(size_t k=0; k<j;k++){
701 //calculate the distance in delta eta and delta phi between object i and object j
702 distance_jk = min(pow(input_pt[j], 2),pow(input_pt[k], 2))*
703 (pow(input_eta[j]-input_eta[k], 2)+
704 pow(_delta_phi(input_phi[j],input_phi[k]), 2));
705 if (distance_jk<distance_jk_min) {
706 distance_jk_min = distance_jk;
707 index_j_jk = j;
708 index_k_jk =k;
709 }
710 }
711 }
712 }
713 //decide if the minimum is from a jB or jk combination
714 if (distance_jk_min<distance_jB_min) {
715 max_dmin_temp = max(distance_jk_min,max_dmin_temp);
716 decide_jB = 0;
717 } else {
718 max_dmin_temp = max(distance_jB_min,max_dmin_temp);
719 decide_jB=1;
720 }
721 //if we have only three jets left calculate
722 //the maxima of the dmin's
723 //if the minimum is a jB eliminate the input object
724 if (decide_jB == 1) {
725 //if index_jB is the last one nothing is to do
726 if (index_jB != rest-1) {
727 for (size_t i=index_jB; i<rest-1;i++) {
728 input_pt[i]=input_pt[i+1];
729 input_phi[i]=input_phi[i+1];
730 input_eta[i]=input_eta[i+1];
731 input_px[i]=input_px[i+1];
732 input_py[i]=input_py[i+1];
733 input_pz[i]=input_pz[i+1];
734 input_e[i]=input_e[i+1];
735 }
736 }
737 }
738 //if the minimum is a jk combine both input objects
739 if(decide_jB==0) {
740 input_px[index_k_jk] = input_px[index_k_jk]+input_px[index_j_jk];
741 input_py[index_k_jk] = input_py[index_k_jk]+input_py[index_j_jk];
742 input_pz[index_k_jk] = input_pz[index_k_jk]+input_pz[index_j_jk];
743 input_e[index_k_jk] = input_e[index_k_jk]+input_e[index_j_jk];
744 input_p[index_k_jk] = sqrt(pow(input_px[index_k_jk], 2)+
745 pow(input_py[index_k_jk], 2)+
746 pow(input_pz[index_k_jk], 2));
747 //calculate the pt, eta and phi of the new combined momenta k_jk
748 input_pt[index_k_jk] = sqrt(pow(input_px[index_k_jk], 2)+
749 pow(input_py[index_k_jk], 2));
750 //in the case of pseudorapidity
751 if (irap == 0) {
752 double theta_new =0;
753 if (fabs(input_pz[index_k_jk]) > 1E-5){
754 theta_new = atan(input_pt[index_k_jk]/(input_pz[index_k_jk]));
755 } else {
756 theta_new = M_PI/2;
757 }
758 if (theta_new < 0) {
759 theta_new = theta_new + M_PI;
760 }
761 input_eta[index_k_jk] = - log(tan(0.5*theta_new));
762 }
763 //in the real rapidity y is wanted
764 if (irap == 1) {
765 input_eta[index_k_jk] = 0.5 * log((input_e[index_k_jk]+
766 input_pz[index_k_jk]) /
767 (input_e[index_k_jk] -
768 input_pz[index_k_jk]));
769 }
770 input_phi[index_k_jk] = atan2(input_py[index_k_jk], input_px[index_k_jk]);
771 if (index_j_jk != rest-1) {
772 for (size_t i = index_j_jk; i<rest-1;i++) {
773 input_pt[i] = input_pt[i+1];
774 input_phi[i] = input_phi[i+1];
775 input_eta[i] = input_eta[i+1];
776 input_px[i] = input_px[i+1];
777 input_py[i] = input_py[i+1];
778 input_pz[i] = input_pz[i+1];
779 input_e[i] = input_e[i+1];
780 }
781 }
782 }
783 }
784 if (rest == 3) max_dmin = max_dmin_temp;
785 rest = rest-1;
786 }
787 }
788
789 double et2 = 0;
790 et2 = input_pt[0] + input_pt[1];
791 y3 = max_dmin/pow(et2,2);
792
793 return y3;
794 }
795
796
797 vector<double> EventShape::_thrust(const vector<double>& input_px, const vector<double>& input_py) {
798
799 double thrustmax_calc = 0;
800 double temp_calc = 0;
801 size_t length_thrust_calc = 0;
802 vector<double> thrust_values, thrust_axis_calc;
803 vector<double> p_thrust_max_calc, p_dec_1_calc, p_dec_2_calc, p_pt_beam_calc;
804
805 if (!thrust_values.empty()){
806 thrust_values.clear();
807 thrust_axis_calc.clear();
808 p_thrust_max_calc.clear();
809 p_dec_1_calc.clear();
810 p_dec_2_calc.clear();
811 p_pt_beam_calc.clear();
812 }
813
814 for (size_t j = 0; j < 3; j++){
815 p_pt_beam_calc.push_back(0.);
816 p_dec_1_calc.push_back(0.);
817 p_dec_2_calc.push_back(0.);
818 p_thrust_max_calc.push_back(0.);
819 thrust_axis_calc.push_back(0.);
820 }
821
822 for (size_t j = 0; j < 4; j++) {
823 thrust_values.push_back(0.);
824 }
825
826 length_thrust_calc = input_px.size();
827 if (input_py.size() != length_thrust_calc) {
828 /// @todo Change to exception or assert
829 cout << "ERROR in thrust calculation!!! Size of input vectors differs. Change that please!" << endl;
830 return thrust_values;
831 }
832
833 double pt_sum_calc =0;
834 for(size_t k=0;k<length_thrust_calc;k++){
835 pt_sum_calc+=sqrt(pow(input_px[k],2)+pow(input_py[k],2));
836 for(size_t j = 0; j < 3; j++){
837 p_thrust_max_calc[j]=0;
838 }
839 //get a vector perpendicular to the beam axis and
840 //perpendicular to the momentum of particle k
841 //per default beam axis b = (0,0,1)
842 p_pt_beam_calc[0] = input_py[k]*1;
843 p_pt_beam_calc[1] = - input_px[k]*1;
844 p_pt_beam_calc[2] = 0.; // GMA p_pt_beam_calc[3] = 0.;
845 for(size_t i=0;i<length_thrust_calc;i++){
846 if(i!=k){
847 if((input_px[i]*p_pt_beam_calc[0]+input_py[i]*p_pt_beam_calc[1])>=0){
848 p_thrust_max_calc[0]= p_thrust_max_calc[0] + input_px[i];
849 p_thrust_max_calc[1]= p_thrust_max_calc[1] + input_py[i];
850 }else{
851 p_thrust_max_calc[0]= p_thrust_max_calc[0] - input_px[i];
852 p_thrust_max_calc[1]= p_thrust_max_calc[1] - input_py[i];
853 }
854 }
855 }
856 p_dec_1_calc[0] = p_thrust_max_calc[0] + input_px[k];
857 p_dec_1_calc[1] = p_thrust_max_calc[1] + input_py[k];
858 p_dec_1_calc[2] = 0;
859 p_dec_2_calc[0] = p_thrust_max_calc[0] - input_px[k];
860 p_dec_2_calc[1] = p_thrust_max_calc[1] - input_py[k];
861 p_dec_2_calc[2] = 0;
862 temp_calc = pow(p_dec_1_calc[0], 2) + pow(p_dec_1_calc[1], 2);
863
864 if (temp_calc>thrustmax_calc) {
865 thrustmax_calc =temp_calc;
866 for (size_t i=0; i<3; i++) {
867 thrust_axis_calc[i] = p_dec_1_calc[i]/sqrt(thrustmax_calc);
868 }
869 }
870 temp_calc = pow(p_dec_2_calc[0], 2)+pow(p_dec_2_calc[1], 2);
871 if (temp_calc > thrustmax_calc) {
872 thrustmax_calc =temp_calc;
873 for (size_t i=0; i<3; i++) {
874 thrust_axis_calc[i] = p_dec_2_calc[i]/sqrt(thrustmax_calc);
875 }
876 }
877 }
878 for (size_t j = 0; j < 3; j++) thrust_values[j] = thrust_axis_calc[j];
879 const double thrust_calc = sqrt(thrustmax_calc)/pt_sum_calc;
880
881 // the variable which gets returned is not the thrust but tau=1-thrust
882 thrust_values[3] = 1 - thrust_calc;
883 if (thrust_values[3] < 1e-20) thrust_values[3] = 1e-20;
884
885 return thrust_values;
886 }
887
888
889 double EventShape::_delta_phi(double phi1, double phi2) {
890 double dphi = fabs(phi2 - phi1);
891 if (dphi > M_PI) dphi = 2*M_PI - dphi;
892 return dphi;
893 }
894
895
896 // Returns the scalar product between two 4 momenta
897 double EventShape::_lorentz_sp(const vector<double>& a, const vector<double>& b) {
898 size_t dim = (size_t) a.size();
899 if (a.size()!=b.size()) {
900 cout << "ERROR!!! Dimension of input vectors are different! Change that please!" << endl;
901 return 0;
902 } else {
903 double l_dot_product=a[dim-1]*b[dim-1];
904 for(size_t i=0; i<dim-1;i++){
905 l_dot_product-=a[i]*b[i];
906 }
907 return l_dot_product;
908 }
909 }
910
911
912 }
913
914
915}
|