rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2014_I1305624

Study of hadronic event-shape variables in multijet final states in $pp$ collisions at $\sqrt{s} = 7$ TeV
Experiment: CMS (LHC)
Inspire ID: 1305624
Status: VALIDATED
Authors:
  • Sunanda Banerjee
  • Sandeep Bhowmik
  • Monoranjan Guchait
  • Gobinda Majumder
  • Manas Maity <
  • Debarati Roy
References:
  • 10.1007/JHEP10(2014)087
  • arXiv: 1407.2856
  • CMS-PAS-SMP-12-022
  • CERN-PH-EP-2014-146
Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • $pp$ QCD interactions at $\sqrt{s} = 7$ TeV. Data collected by CMS during the year 2011.

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), FastJets::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(30.0*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      const double weight = 1.0;
172      for (int ij = NJETPTMN-1; ij >= 0; --ij) {
173        if (leadingpt/GeV > LEADINGPTTHRESHOLD[ij]) {
174          if (inRange(eventvar[0], _alow1[ij], _ahgh1[ij])) _h_thrustc[ij]->fill(eventvar[0], weight);
175          if (inRange(eventvar[2], _alow3[ij], _ahgh3[ij])) _h_tot3dmass[ij]->fill(eventvar[2], weight);
176          if (inRange(eventvar[3], _alow4[ij], _ahgh4[ij])) _h_tottrnsmass[ij]->fill(eventvar[3], weight);
177          if (eventvar[NEVTVAR] >= 3) {
178            if (inRange(eventvar[1], _alow2[ij], _ahgh2[ij])) _h_broadt[ij]->fill(eventvar[1], weight);
179            if (inRange(eventvar[4], _alow5[ij], _ahgh5[ij])) _h_y23c[ij]->fill(eventvar[4], weight);
180          }
181          break;
182        }
183      }
184
185    }
186
187
188    /// Normalise histograms etc., after the run
189    void finalize() {
190      for (int ij = 0; ij < NJETPTMN; ij++) {
191        normalize(_h_thrustc[ij]);
192        normalize(_h_broadt[ij]);
193        normalize(_h_tot3dmass[ij]);
194        normalize(_h_tottrnsmass[ij]);
195        normalize(_h_y23c[ij]);
196      }
197    }
198
199
200  private:
201
202    /// @name Histograms
203    //@{
204    Histo1DPtr _h_thrustc[NJETPTMN];
205    Histo1DPtr _h_broadt[NJETPTMN];
206    Histo1DPtr _h_tot3dmass[NJETPTMN];
207    Histo1DPtr _h_tottrnsmass[NJETPTMN];
208    Histo1DPtr _h_y23c[NJETPTMN];
209    //@}
210
211    // Data members
212    bool _needBinInit;
213    double _alow1[NJETPTMN], _alow2[NJETPTMN], _alow3[NJETPTMN], _alow4[NJETPTMN], _alow5[NJETPTMN];
214    double _ahgh1[NJETPTMN], _ahgh2[NJETPTMN], _ahgh3[NJETPTMN], _ahgh4[NJETPTMN], _ahgh5[NJETPTMN];
215
216  };
217
218
219  RIVET_DECLARE_PLUGIN(CMS_2014_I1305624);
220
221
222
223  /////////////////////
224
225
226  namespace {
227
228    // EventShape helper class method implementations:
229
230    int EventShape::_calculate() {
231      if (!_event_shapes.empty() && !_thrust_axis.empty() && !_thrust_axis_c.empty())
232        return 1; //< return success if this appears to already have been run
233
234      const size_t length = (size_t) _object_px.size();
235
236      if (((size_t) _object_py.size() != length) ||
237          ((size_t) _object_pz.size() != length) ||
238          ((size_t) _object_e.size() != length)) {
239        /// @todo Change to exception or assert
240        // cout << "ERROR!!!! Input vectors differ in size! Change that please!" << '\n';
241        // cout<<"py_size: "<<_object_py.size()<<" ,pz_size: "<<_object_pz.size()
242        //     <<" ,px_size: "<<_object_px.size()<<" ,E_size: "<<_object_e.size()<<'\n';
243        return 0;
244      }
245
246      if (!_object_p.empty()) {
247        _object_p.clear();
248        _object_pt.clear();
249        _object_eta.clear();
250        _object_phi.clear();
251        _event_shapes.clear();
252        _thrust_axis.clear();
253        _thrust_axis_c.clear();
254      }
255
256      for (size_t j = 0; j < length; j++) {
257        _object_p.push_back(0.);
258        _object_pt.push_back(0.);
259        _object_eta.push_back(0.);
260        _object_phi.push_back(0.);
261      }
262
263      for (int j = 0; j < NEVTVAR; j++) {
264        _event_shapes.push_back(-50.);
265      }
266
267      _event_shapes.push_back(double(_object_px.size())); //< WTF?
268
269      for (size_t j = 0; j < 3; j++) {
270        _thrust_axis.push_back(0.);
271        _thrust_axis_c.push_back(0.);
272      }
273
274      double theta = 0;
275
276      for (size_t k = 0; k < length; k++) {
277        _object_p[k] = sqrt(pow(_object_px[k],2) + pow(_object_py[k],2) + pow(_object_pz[k],2));
278        _object_pt[k] = sqrt(pow(_object_px[k],2) + pow(_object_py[k],2));
279        if (_object_p[k] > _object_e[k] + 1e-4) {
280          /// @todo Change to exception or assert
281          // cout << "ERROR!!! object " << k <<" has P = " << _object_p[k]
282          //      << " which is bigger than E = " << _object_e[k] <<" "
283          //      << _object_px[k] <<" "<< _object_py[k] <<" "
284          //      << _object_pz[k] <<" of total length "<< length
285          //      << '\n';
286          return 0;
287        }
288
289        //to prevent a division by zero
290        if (_irap == 0) {
291          if (fabs(_object_pz[k]) > 1e-5) {
292            theta = atan(_object_pt[k]/(_object_pz[k]));
293          } else {
294            theta = M_PI/2;
295          }
296          if (theta < 0.) theta = theta + M_PI;
297          _object_eta[k] = -log(tan(0.5*theta));
298        }
299        if (_irap == 1) {
300          if (_object_pz[k] == _object_e[k]) {
301            /// @todo Change to exception or assert
302            // cout << "ERROR!!! object "<<k<<" has Pz "<< _object_pz[k] <<" which is equal to E = "<< _object_e[k] <<'\n';
303            return 0;
304          }
305          _object_eta[k]=0.5*log((_object_e[k]+_object_pz[k])/(_object_e[k]-_object_pz[k]));
306        }
307        if (_irap != 0 && _irap != 1) {
308          /// @todo Change to exception or assert
309          // cout << "ERROR!!!, The choice to use the rapidity y or the pseudorapidity eta is not set correctly! Change that please!" << '\n';
310          return 0;
311        }
312        _object_phi[k] = atan2(_object_py[k], _object_px[k]);
313      }
314
315      vector<double> object_px_in, object_py_in, object_pz_in, object_pt_in, object_e_in, object_et_in, object_eta_in;
316      vector<double> object_px_out, object_py_out, object_pz_out, object_e_out, object_pt_out, object_eta_out;
317      if (!object_px_in.empty()) { //< FFS, this is impossible: it's only just been created!
318        object_px_in.clear();
319        object_py_in.clear();
320        object_pz_in.clear();
321        object_pt_in.clear();
322        object_e_in.clear();
323        object_et_in.clear();
324        object_eta_in.clear();
325        object_px_out.clear();
326        object_py_out.clear();
327        object_pz_out.clear();
328        object_pt_out.clear();
329        object_e_out.clear();
330        object_eta_out.clear();
331      }
332
333      size_t nin = 0;
334
335      for (size_t j = 0; j < length; j++) {
336        if (fabs(_object_eta[j]) < _eta_c) {
337          object_px_in.push_back(_object_px[j]);
338          object_py_in.push_back(_object_py[j]);
339          object_pz_in.push_back(_object_pz[j]);
340          object_e_in.push_back(_object_e[j]);
341          object_pt_in.push_back(sqrt(pow(_object_px[j],2)+pow(_object_py[j],2)));
342          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))));
343          object_eta_in.push_back(_object_eta[j]);
344          nin += 1;
345      } else {
346          object_px_out.push_back(_object_px[j]);
347          object_py_out.push_back(_object_py[j]);
348          object_pz_out.push_back(_object_pz[j]);
349          object_e_out.push_back(_object_e[j]);
350          object_pt_out.push_back(sqrt(pow(_object_px[j],2)+pow(_object_py[j],2)));
351          object_eta_out.push_back(_object_eta[j]);
352        }
353      }
354
355      if (object_px_in.size() != nin) {
356        /// @todo Change to exception or assert
357        cout<<"ERROR!!! wrong dimension of 'in' momenta"<<'\n';
358        //return 0; ///< @todo Why not do this?
359      }
360      const size_t nout = length - nin;
361
362      if (nin < _nmnjet) {
363        for (int i = 0; i < NEVTVAR; i++) {
364          _event_shapes[i] = -50.0;
365        }
366      }
367
368      _event_shapes[NEVTVAR] = nin;
369
370      if (nin >= _nmnjet) {
371        double p_sum_c = 0; //GMA
372        double pt_sum_c = 0;
373        double eta_cw=0;
374        double px_sum_in = 0;
375        double py_sum_in = 0;
376        for (size_t j = 0; j < nin; j++) {
377          pt_sum_c += object_pt_in[j];
378          p_sum_c += sqrt(pow(object_pt_in[j],2.) + pow(object_pz_in[j], 2.0)); //GMA
379          eta_cw += object_pt_in[j]*object_eta_in[j];
380          px_sum_in += object_px_in[j];
381          py_sum_in += object_py_in[j];
382        }
383        eta_cw /= pt_sum_c;
384
385        double expTerm = 0;
386        for (size_t j = 0; j < nout; j++) {
387          expTerm += object_pt_out[j] * exp(-fabs(object_eta_out[j]-eta_cw));
388        }
389        expTerm /= pt_sum_c;
390
391        //the central global transverse thrust centrthr is calculated
392        double centrthr = 0;
393        vector<double> thrust_central = _thrust(object_px_in, object_py_in);
394
395        for (size_t l=0; l<3; l++) _thrust_axis_c[l] = thrust_central[l];
396        //the variable which gets resummed is not thrust
397        //but tau = 1 - thrust - see calculation
398        centrthr = thrust_central[3];
399        _event_shapes[0] = centrthr;
400
401        double alpha_c = atan2(_thrust_axis_c[1], _thrust_axis_c[0]);
402        //central jet masses
403        //define two jet masses in region U and D
404        double cenjm_up = 0;
405        double cenjm_down= 0;
406        double dot_product = 0;
407
408        vector<double> up_sum;
409        vector<double> down_sum;
410        for (size_t j=0; j<4;j++) {
411          up_sum.push_back(0.);
412          down_sum.push_back(0.);
413        }
414        for (size_t i=0;i<nin;i++) {
415          dot_product = object_px_in[i] * _thrust_axis_c[0] + object_py_in[i] * _thrust_axis_c[1];
416          if (dot_product >= 0) {
417            up_sum[0]+=object_px_in[i];
418            up_sum[1]+=object_py_in[i];
419            up_sum[2]+=object_pz_in[i];
420            up_sum[3]+=object_e_in[i];
421          } else {
422            down_sum[0]+=object_px_in[i];
423            down_sum[1]+=object_py_in[i];
424            down_sum[2]+=object_pz_in[i];
425            down_sum[3]+=object_e_in[i];
426          }
427        }
428        cenjm_up = _lorentz_sp(up_sum, up_sum) / pow(p_sum_c, 2.); //GMA pow(pt_sum_c,2);
429        cenjm_down = _lorentz_sp(down_sum, down_sum) / pow(p_sum_c, 2.); //GMA pow(pt_sum_c,2);
430
431        //central total jet mass centotjm
432        double centotjm=0;
433        centotjm = cenjm_up + cenjm_down;
434
435        _event_shapes[2]=centotjm;
436
437        double centrjm_up=0, centrjm_down=0;
438        vector<double> upsum;
439        vector<double> downsum;
440        for (size_t j = 0; j < 3; j++) {
441          upsum.push_back(0.);
442          downsum.push_back(0.);
443        }
444        for (size_t i = 0; i < nin; i++) {
445          dot_product = object_px_in[i]*_thrust_axis_c[0]+object_py_in[i]*_thrust_axis_c[1];
446          if (dot_product >= 0) {
447            upsum[0] += object_px_in[i];
448            upsum[1] += object_py_in[i];
449            upsum[2] += object_et_in[i];
450          } else {
451            downsum[0] += object_px_in[i];
452            downsum[1] += object_py_in[i];
453            downsum[2] += object_et_in[i];
454          }
455        }
456        centrjm_up = _lorentz_sp(upsum, upsum) / pow(pt_sum_c, 2);
457        centrjm_down = _lorentz_sp(downsum, downsum) / pow(pt_sum_c, 2);
458        double centottrjm = centrjm_up + centrjm_down;
459
460        _event_shapes[3] = centottrjm;
461
462        //central three-jet resolution threshold
463        double ceny3=0;
464        if (nin < 3) {
465          ceny3 = -1.0;
466        } else {
467          ceny3 = _three_jet_res(object_px_in, object_py_in, object_pz_in, object_e_in, _irap);
468        }
469
470        _event_shapes[4] = ceny3;
471
472        //the central jet broadenings in the up and down region
473        double cenbroad_up=0;
474        double cenbroad_down=0;
475
476        double eta_up=0;
477        size_t num_up=0;
478        double eta_down =0;
479        size_t num_down =0;
480        double phi_temp =0;
481        double phi_up_aver =0;
482        double phi_down_aver =0;
483        double pt_sum_up =0;
484        double pt_sum_down =0;
485        double dot_product_b =0;
486        vector<double> phi_up;
487        vector<double> phi_down;
488        double py_rot =0;
489        double px_rot =0;
490
491        for (size_t j = 0; j < 4; j++) {
492          up_sum.push_back(0.);
493          down_sum.push_back(0.);
494        }
495
496        for (size_t i=0;i<nin;i++) {
497          dot_product_b =sqrt(object_px_in[i]*_thrust_axis_c[0] + object_py_in[i]*_thrust_axis_c[1]);
498          if (dot_product_b>=0){
499            pt_sum_up += object_pt_in[i];
500            //rotate the coordinate system so that
501            //the central thrust axis is e_x
502            px_rot = cos(alpha_c)*object_px_in[i]+sin(alpha_c)*object_py_in[i];
503            py_rot = - sin(alpha_c)*object_px_in[i]+cos(alpha_c)*object_py_in[i];
504            //calculate the eta and phi in the rotated system
505            eta_up += object_pt_in[i]*object_eta_in[i];
506            phi_temp = atan2(py_rot,px_rot);
507
508            if(phi_temp > M_PI/2){
509              phi_temp = phi_temp - M_PI/2;
510            }
511            if (phi_temp < -M_PI/2){
512              phi_temp = phi_temp + M_PI/2;
513            }
514            phi_up.push_back(phi_temp);
515            phi_up_aver += object_pt_in[i]*phi_temp;
516            num_up += 1;
517          } else {
518            eta_down += object_pt_in[i]*object_eta_in[i];
519            pt_sum_down += object_pt_in[i];
520            px_rot = cos(alpha_c)*object_px_in[i]+sin(alpha_c)*object_py_in[i];
521            py_rot = - sin(alpha_c)*object_px_in[i]+cos(alpha_c)*object_py_in[i];
522            phi_temp = atan2(py_rot,px_rot);
523            if (phi_temp > M_PI/2) {
524              //if phi is bigger than pi/2 in the new system calculate
525              //the difference to the thrust axis
526              phi_temp = M_PI -phi_temp;
527            }
528            if (phi_temp<-M_PI/2) {
529              //if phi is smaller than
530              phi_temp = -M_PI-phi_temp;
531            }
532            phi_down.push_back(phi_temp);
533            //calculate the pt-weighted phi
534            phi_down_aver += object_pt_in[i]*phi_temp;
535            num_down += 1;
536          }
537        }
538        if (num_up!=0){
539          eta_up = eta_up/pt_sum_up;
540          phi_up_aver = phi_up_aver/pt_sum_up;
541        }
542        if (num_down!=0) {
543          eta_down = eta_down/pt_sum_down;
544          phi_down_aver = phi_down_aver/pt_sum_down;
545        }
546
547        size_t index_up=0, index_down=0;
548        for (size_t i = 0; i < nin; i++) {
549          dot_product_b = object_px_in[i]*_thrust_axis_c[0] + object_py_in[i]*_thrust_axis_c[1];
550          if (dot_product_b >= 0) {
551            //calculate the broadenings of the regions with the rotated system
552            //and the pt-weighted average of phi in the rotated system
553            cenbroad_up += object_pt_in[i]*sqrt(pow(object_eta_in[i]-eta_up, 2) +
554                                                pow(_delta_phi(phi_up[index_up], phi_up_aver), 2));
555            index_up += 1;
556          } else {
557            cenbroad_down += object_pt_in[i]*sqrt(pow(object_eta_in[i]-eta_down, 2)+
558                                                  pow(_delta_phi(phi_down[index_down], phi_down_aver), 2));
559            index_down += 1;
560          }
561        }
562
563        if (index_up == 0 || index_down ==0) _event_shapes[NEVTVAR] *= -1;
564
565        cenbroad_up=cenbroad_up/(2*pt_sum_c);
566        cenbroad_down=cenbroad_down/(2*pt_sum_c);
567
568        //central total jet broadening
569        double centotbroad = 0;
570        centotbroad = cenbroad_up + cenbroad_down;
571
572        _event_shapes[1] = centotbroad;
573
574        for (int ij = 0; ij < 5; ij++) {
575          if (_event_shapes[ij] < 1.e-20) _event_shapes[ij] = 1.e-20;
576          _event_shapes[ij] = log(_event_shapes[ij]);
577        }
578      }
579
580      return 1;
581    }
582
583
584    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) {
585
586      size_t y3_length = (size_t)in_object_px.size();
587      if (((size_t) in_object_py.size()!=y3_length) ||
588          ((size_t) in_object_pz.size()!=y3_length) ||
589          (in_object_e.size()!=y3_length)) {
590        // cout << "ERROR!!!! Input vectors differ in size! Change that please!" << '\n';
591        // cout<<"py_size: "<<in_object_py.size()<<" ,pz_size: "<<in_object_pz.size()
592        //     <<" ,px_size: "<<in_object_px.size()<<" , E_size: "<<in_object_e.size() <<'\n';
593        return 0.0;
594      }
595
596      vector<double> in_object_p, in_object_pt, in_object_eta, in_object_phi;
597      if (!in_object_p.empty()) {
598        in_object_p.clear();
599        in_object_pt.clear();
600        in_object_eta.clear();
601        in_object_phi.clear();
602      }
603      for (size_t j = 0; j < y3_length; j++) {
604        in_object_p.push_back(0.);
605        in_object_pt.push_back(0.);
606        in_object_eta.push_back(0.);
607        in_object_phi.push_back(0.);
608      }
609      double theta_y3_1st = 0;
610      for (size_t k =0; k<y3_length; k++) {
611        in_object_p[k] = sqrt(pow(in_object_px[k],2) + pow(in_object_py[k],2) + pow(in_object_pz[k],2));
612        in_object_pt[k] = sqrt(pow(in_object_px[k],2) + pow(in_object_py[k],2));
613
614        //calculates the pseudorapidity to prevent a division by zero
615        if (irap == 0) {
616          if (fabs(in_object_pz[k]) > 1E-5) {
617            theta_y3_1st = atan(in_object_pt[k]/(in_object_pz[k]));
618          } else {
619            theta_y3_1st = M_PI/2;
620          }
621          if (theta_y3_1st<0.) theta_y3_1st = theta_y3_1st + M_PI;
622          in_object_eta[k] = - log(tan(0.5*theta_y3_1st));
623        }
624        //calculates the real rapidity
625        if (irap == 1) {
626          in_object_eta[k]=0.5*log((in_object_e[k]+in_object_pz[k])/(in_object_e[k]-in_object_pz[k]));
627        }
628        in_object_phi[k] = atan2(in_object_py[k], in_object_px[k]);
629      }
630
631      //the three-jet resolution
632      //threshold y3
633      double y3 = 0;
634
635      //vector which will be filled with the
636      //minimum of the distances
637      double max_dmin_temp=0;
638
639      double max_dmin = 0;
640
641      //distance input object k, beam
642      double distance_jB = 0;
643      double distance_jB_min = 0;
644      //distance of input object k to l
645      double distance_jk = 0;
646      double distance_jk_min = 0;
647      //as we search the minimum of the distances
648      //give them values which are for sure higher
649      //than those we evaluate first in the for-loups
650
651      size_t index_jB = 0;
652      size_t index_j_jk = 0;
653      size_t index_k_jk = 0;
654
655      //to decide later if the minmum is a jB or jk
656      int decide_jB = -1;
657
658      vector<double> input_pt, input_px, input_py, input_pz;
659      vector<double> input_p, input_e, input_phi, input_eta;
660
661      if (!input_pt.empty()) {
662        input_pt.clear();
663        input_px.clear();
664        input_px.clear();
665        input_pz.clear();
666        input_p.clear();
667        input_e.clear();
668        input_phi.clear();
669        input_eta.clear();
670      }
671
672      for (size_t j = 0; j < y3_length; j++){
673        input_pt.push_back(in_object_pt[j]);
674        input_px.push_back(in_object_px[j]);
675        input_py.push_back(in_object_py[j]);
676        input_pz.push_back(in_object_pz[j]);
677        input_p.push_back(in_object_p[j]);
678        input_e.push_back(in_object_e[j]);
679        input_phi.push_back(in_object_phi[j]);
680        input_eta.push_back(in_object_eta[j]);
681      }
682      if (y3_length<3) {
683        return -1;
684      } else {
685        size_t rest = y3_length;
686        for (size_t i = 0; i<y3_length; i++) {
687          //make the minima at the initialization step
688          //of each looping bigger than the first values
689          distance_jB_min = 0.36*pow(input_pt[0],2) + 10;
690          //DELTA PHIs wanted not the pure difference
691          distance_jk_min = min(pow(input_pt[1], 2), pow(input_pt[0], 2)) *
692            (pow(input_eta[1]-input_eta[0], 2) +
693             pow(_delta_phi(input_phi[1], input_phi[0]), 2)) + 10;
694          //do the procedure only until we have only 2 objects left anymore
695          if (rest > 2) {
696            for (size_t j=0; j<rest;j++) {
697              //calculate the distance between object j and the beam
698              distance_jB = 0.36*pow(input_pt[j], 2);
699              if(distance_jB < distance_jB_min){
700                distance_jB_min = distance_jB;
701                index_jB = j;
702              }
703              if (j > 0) {
704                for(size_t k=0; k<j;k++){
705                  //calculate the distance in delta eta and delta phi between object i and object j
706                  distance_jk = min(pow(input_pt[j], 2),pow(input_pt[k], 2))*
707                    (pow(input_eta[j]-input_eta[k], 2)+
708                     pow(_delta_phi(input_phi[j],input_phi[k]), 2));
709                  if (distance_jk<distance_jk_min) {
710                    distance_jk_min = distance_jk;
711                    index_j_jk = j;
712                    index_k_jk =k;
713                  }
714                }
715              }
716            }
717            //decide if the minimum is from a jB or jk combination
718            if (distance_jk_min<distance_jB_min) {
719              max_dmin_temp = max(distance_jk_min,max_dmin_temp);
720              decide_jB = 0;
721            } else {
722              max_dmin_temp = max(distance_jB_min,max_dmin_temp);
723              decide_jB=1;
724            }
725            //if we have only three jets left calculate
726            //the maxima of the dmin's
727            //if the minimum is a jB eliminate the input object
728            if (decide_jB == 1) {
729              //if index_jB is the last one nothing is to do
730              if (index_jB != rest-1) {
731                for (size_t i=index_jB; i<rest-1;i++) {
732                  input_pt[i]=input_pt[i+1];
733                  input_phi[i]=input_phi[i+1];
734                  input_eta[i]=input_eta[i+1];
735                  input_px[i]=input_px[i+1];
736                  input_py[i]=input_py[i+1];
737                  input_pz[i]=input_pz[i+1];
738                  input_e[i]=input_e[i+1];
739                }
740              }
741            }
742            //if the minimum is a jk combine both input objects
743            if(decide_jB==0) {
744              input_px[index_k_jk] = input_px[index_k_jk]+input_px[index_j_jk];
745              input_py[index_k_jk] = input_py[index_k_jk]+input_py[index_j_jk];
746              input_pz[index_k_jk] = input_pz[index_k_jk]+input_pz[index_j_jk];
747              input_e[index_k_jk] = input_e[index_k_jk]+input_e[index_j_jk];
748              input_p[index_k_jk] = sqrt(pow(input_px[index_k_jk], 2)+
749                                         pow(input_py[index_k_jk], 2)+
750                                         pow(input_pz[index_k_jk], 2));
751              //calculate the pt, eta and phi of the new combined momenta k_jk
752              input_pt[index_k_jk] = sqrt(pow(input_px[index_k_jk], 2)+
753                                          pow(input_py[index_k_jk], 2));
754              //in the case of pseudorapidity
755              if (irap == 0) {
756                double theta_new =0;
757                if (fabs(input_pz[index_k_jk]) > 1E-5){
758                  theta_new = atan(input_pt[index_k_jk]/(input_pz[index_k_jk]));
759                } else {
760                  theta_new = M_PI/2;
761                }
762                if (theta_new < 0) {
763                  theta_new = theta_new + M_PI;
764                }
765                input_eta[index_k_jk] = - log(tan(0.5*theta_new));
766              }
767              //in the real rapidity y is wanted
768              if (irap == 1) {
769                input_eta[index_k_jk] = 0.5 * log((input_e[index_k_jk]+
770                                                   input_pz[index_k_jk]) /
771                                                  (input_e[index_k_jk] -
772                                                   input_pz[index_k_jk]));
773              }
774              input_phi[index_k_jk] = atan2(input_py[index_k_jk], input_px[index_k_jk]);
775              if (index_j_jk != rest-1) {
776                for (size_t i = index_j_jk; i<rest-1;i++) {
777                  input_pt[i] = input_pt[i+1];
778                  input_phi[i] = input_phi[i+1];
779                  input_eta[i] = input_eta[i+1];
780                  input_px[i] = input_px[i+1];
781                  input_py[i] = input_py[i+1];
782                  input_pz[i] = input_pz[i+1];
783                  input_e[i] = input_e[i+1];
784                }
785              }
786            }
787          }
788          if (rest == 3) max_dmin = max_dmin_temp;
789          rest = rest-1;
790        }
791      }
792
793      double et2 = 0;
794      et2 = input_pt[0] + input_pt[1];
795      y3 = max_dmin/pow(et2,2);
796
797      return y3;
798    }
799
800
801    vector<double> EventShape::_thrust(const vector<double>& input_px, const vector<double>& input_py) {
802
803      double thrustmax_calc = 0;
804      double temp_calc = 0;
805      size_t length_thrust_calc = 0;
806      vector<double> thrust_values, thrust_axis_calc;
807      vector<double> p_thrust_max_calc, p_dec_1_calc,  p_dec_2_calc, p_pt_beam_calc;
808
809      if (!thrust_values.empty()){
810        thrust_values.clear();
811        thrust_axis_calc.clear();
812        p_thrust_max_calc.clear();
813        p_dec_1_calc.clear();
814        p_dec_2_calc.clear();
815        p_pt_beam_calc.clear();
816      }
817
818      for (size_t j = 0; j < 3; j++){
819        p_pt_beam_calc.push_back(0.);
820        p_dec_1_calc.push_back(0.);
821        p_dec_2_calc.push_back(0.);
822        p_thrust_max_calc.push_back(0.);
823        thrust_axis_calc.push_back(0.);
824      }
825
826      for (size_t j = 0; j < 4; j++) {
827        thrust_values.push_back(0.);
828      }
829
830      length_thrust_calc = input_px.size();
831      if (input_py.size() != length_thrust_calc) {
832        /// @todo Change to exception or assert
833        cout<<"ERROR in thrust calculation!!! Size of input vectors differs. Change that please!"<<'\n';
834        return thrust_values;
835      }
836
837      double pt_sum_calc =0;
838      for(size_t k=0;k<length_thrust_calc;k++){
839        pt_sum_calc+=sqrt(pow(input_px[k],2)+pow(input_py[k],2));
840        for(size_t j = 0; j < 3; j++){
841          p_thrust_max_calc[j]=0;
842        }
843        //get a vector perpendicular to the beam axis and
844        //perpendicular to the momentum of particle k
845        //per default beam axis b = (0,0,1)
846        p_pt_beam_calc[0] = input_py[k]*1;
847        p_pt_beam_calc[1] = - input_px[k]*1;
848        p_pt_beam_calc[2] = 0.; // GMA p_pt_beam_calc[3] = 0.;
849        for(size_t i=0;i<length_thrust_calc;i++){
850          if(i!=k){
851            if((input_px[i]*p_pt_beam_calc[0]+input_py[i]*p_pt_beam_calc[1])>=0){
852              p_thrust_max_calc[0]= p_thrust_max_calc[0] + input_px[i];
853              p_thrust_max_calc[1]= p_thrust_max_calc[1] + input_py[i];
854            }else{
855              p_thrust_max_calc[0]= p_thrust_max_calc[0] - input_px[i];
856              p_thrust_max_calc[1]= p_thrust_max_calc[1] - input_py[i];
857            }
858          }
859        }
860        p_dec_1_calc[0] = p_thrust_max_calc[0] + input_px[k];
861        p_dec_1_calc[1] = p_thrust_max_calc[1] + input_py[k];
862        p_dec_1_calc[2] = 0;
863        p_dec_2_calc[0] = p_thrust_max_calc[0] - input_px[k];
864        p_dec_2_calc[1] = p_thrust_max_calc[1] - input_py[k];
865        p_dec_2_calc[2] = 0;
866        temp_calc = pow(p_dec_1_calc[0], 2) + pow(p_dec_1_calc[1], 2);
867
868        if (temp_calc>thrustmax_calc) {
869          thrustmax_calc =temp_calc;
870          for (size_t i=0; i<3; i++) {
871            thrust_axis_calc[i] = p_dec_1_calc[i]/sqrt(thrustmax_calc);
872          }
873        }
874        temp_calc = pow(p_dec_2_calc[0], 2)+pow(p_dec_2_calc[1], 2);
875        if (temp_calc > thrustmax_calc) {
876          thrustmax_calc =temp_calc;
877          for (size_t i=0; i<3; i++) {
878            thrust_axis_calc[i] = p_dec_2_calc[i]/sqrt(thrustmax_calc);
879          }
880        }
881      }
882      for (size_t j = 0; j < 3; j++) thrust_values[j] = thrust_axis_calc[j];
883      const double thrust_calc = sqrt(thrustmax_calc)/pt_sum_calc;
884
885      // the variable which gets returned is not the thrust but tau=1-thrust
886      thrust_values[3] = 1 - thrust_calc;
887      if (thrust_values[3] < 1e-20) thrust_values[3] = 1e-20;
888
889      return thrust_values;
890    }
891
892
893    double EventShape::_delta_phi(double phi1, double phi2) {
894      double dphi = fabs(phi2 - phi1);
895      if (dphi > M_PI) dphi = 2*M_PI - dphi;
896      return dphi;
897    }
898
899
900    // Returns the scalar product between two 4 momenta
901    double EventShape::_lorentz_sp(const vector<double>& a, const vector<double>& b) {
902      size_t dim = (size_t) a.size();
903      if (a.size()!=b.size()) {
904        cout<<"ERROR!!! Dimension of input vectors are different! Change that please!"<<'\n';
905        return 0;
906      } else {
907        double l_dot_product=a[dim-1]*b[dim-1];
908        for(size_t i=0; i<dim-1;i++){
909          l_dot_product-=a[i]*b[i];
910        }
911        return l_dot_product;
912      }
913    }
914
915
916  }
917
918
919}