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: 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), 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}