1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6#include "Rivet/Projections/Sphericity.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// Rivet analysis class for ATLAS min bias event shapes
 12  class ATLAS_2012_I1124167 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    ATLAS_2012_I1124167()
 17      : Analysis("ATLAS_2012_I1124167") {  }
 18
 19
 20    /// Initialization, called once before running
 21    void init() {
 22      // Projections
 23      ChargedFinalState cfs(Cuts::abseta < 2.5 && Cuts::pT > 0.5*GeV);
 24      declare(cfs, "CFS");
 25
 26      // Book histograms
 27      book(_hist_T_05_25 ,1,1,1);
 28      book(_hist_T_05    ,2,1,1);
 29      book(_hist_T_25_50 ,1,1,2);
 30      book(_hist_T_25    ,2,1,2);
 31      book(_hist_T_50_75 ,1,1,3);
 32      book(_hist_T_50    ,2,1,3);
 33      book(_hist_T_75_100,1,1,4);
 34      book(_hist_T_75    ,2,1,4);
 35      book(_hist_T_100   ,2,1,5);
 36
 37      book(_hist_TM_05_25 ,3,1,1);
 38      book(_hist_TM_05    ,4,1,1);
 39      book(_hist_TM_25_50 ,3,1,2);
 40      book(_hist_TM_25    ,4,1,2);
 41      book(_hist_TM_50_75 ,3,1,3);
 42      book(_hist_TM_50    ,4,1,3);
 43      book(_hist_TM_75_100,3,1,4);
 44      book(_hist_TM_75    ,4,1,4);
 45      book(_hist_TM_100   ,4,1,5);
 46
 47      book(_hist_S_05_25 ,5,1,1);
 48      book(_hist_S_05    ,6,1,1);
 49      book(_hist_S_25_50 ,5,1,2);
 50      book(_hist_S_25    ,6,1,2);
 51      book(_hist_S_50_75 ,5,1,3);
 52      book(_hist_S_50    ,6,1,3);
 53      book(_hist_S_75_100,5,1,4);
 54      book(_hist_S_75    ,6,1,4);
 55      book(_hist_S_100   ,6,1,5);
 56
 57
 58      book(_hist_T_N  ,7,1,1);
 59      book(_hist_TM_N ,7,1,2);
 60      book(_hist_S_N  ,7,1,3);
 61
 62      book(_hist_T_S  ,8,1,1);
 63      book(_hist_TM_S ,8,1,2);
 64      book(_hist_S_S  ,8,1,3);
 65    }
 66
 67
 68    void analyze(const Event& event) {
 69
 70      // CFS projection and particles
 71      const Particles& particles500 = apply<ChargedFinalState>(event, "CFS").particlesByPt();
 72
 73      // Require at least 6 charged particles
 74      if (particles500.size() < 6) vetoEvent;
 75
 76      // Preparation for Thrust calculation
 77      vector<Vector3> momenta;
 78
 79      // Counters
 80      double num500 = 0;
 81      double ptSum500 = 0;
 82
 83      double pTlead = particles500[0].pT()/GeV;
 84
 85      // Loop over particles
 86      for (const Particle& p : particles500) {
 87        num500 += 1;
 88        ptSum500 += p.pT()/GeV;
 89
 90        // Transverse Thrust calculation requires p_z to be set to 0
 91        Vector3 mom = p.p3();
 92        mom.setZ(0.0);
 93        momenta.push_back(mom);
 94      }
 95
 96      // If only 2 particles, we need to use a ghost so that Thrust.calc() doesn't return 1.
 97      if (momenta.size() == 2) {
 98        momenta.push_back(Vector3(1e-10*MeV, 0., 0.));
 99      }
100
101      // Actual thrust calculation
102      Thrust thrust;
103      thrust.calc(momenta);
104
105      const double T  = 1.0 - thrust.thrust();
106      const double TM = thrust.thrustMajor();
107
108      Sphericity sphericity;
109      sphericity.calc(momenta);
110
111      double S = sphericity.transSphericity();
112      if ( std::isnan(S) )  S = -1.0; // put this in the underflow bin
113
114      // Fill histos, most inclusive first
115
116      // pTlead > 0.5
117      _hist_T_05->fill(T );
118      _hist_TM_05->fill(TM);
119      _hist_S_05->fill(S );
120
121      // pTlead 0.5 - 2.5
122      if (pTlead <= 2.5) {
123        _hist_T_05_25->fill(T );
124        _hist_TM_05_25->fill(TM);
125        _hist_S_05_25->fill(S );
126      }
127
128      // pTlead > 2.5
129      if (pTlead > 2.5) {
130        _hist_T_25->fill(T );
131        _hist_TM_25->fill(TM);
132        _hist_S_25->fill(S );
133      }
134
135      // pTlead 2.5 - .5
136      if (inRange(pTlead, 2.5, 5.0)) {
137        _hist_T_25_50->fill(T );
138        _hist_TM_25_50->fill(TM);
139        _hist_S_25_50->fill(S );
140      }
141
142      // pTlead > 5
143      if (pTlead > 5) {
144        _hist_T_50->fill(T );
145        _hist_TM_50->fill(TM);
146        _hist_S_50->fill(S );
147      }
148
149      // pTlead 5 - 7.5
150      if (inRange(pTlead, 5.0, 7.5)) {
151        _hist_T_50_75->fill(T );
152        _hist_TM_50_75->fill(TM);
153        _hist_S_50_75->fill(S );
154      }
155
156      // pTlead > 7.5
157      if (pTlead > 7.5) {
158        _hist_T_75->fill(T );
159        _hist_TM_75->fill(TM);
160        _hist_S_75->fill(S );
161      }
162
163      // pTlead 7.5 - 10
164      if (inRange(pTlead, 7.5, 10)) {
165        _hist_T_75_100->fill(T );
166        _hist_TM_75_100->fill(TM);
167        _hist_S_75_100->fill(S );
168      }
169
170      // pTlead > 10
171      if (pTlead > 10) {
172        _hist_T_100->fill(T );
173        _hist_TM_100->fill(TM);
174        _hist_S_100->fill(S );
175      }
176
177
178      // Profiles Nch vs. ES
179      _hist_T_N->fill(num500, T);
180      _hist_TM_N->fill(num500, TM);
181      _hist_S_N->fill(num500, S);
182
183      // Profiles pTsum vs. ES
184      _hist_T_S->fill(ptSum500, T);
185      _hist_TM_S->fill(ptSum500, TM);
186      _hist_S_S->fill(ptSum500, S);
187    }
188
189
190    void finalize() {
191      normalize(_hist_T_05_25);
192      normalize(_hist_T_05);
193      normalize(_hist_T_25_50);
194      normalize(_hist_T_25);
195      normalize(_hist_T_50_75);
196      normalize(_hist_T_50);
197      normalize(_hist_T_75_100);
198      normalize(_hist_T_75);
199      normalize(_hist_T_100);
200
201      normalize(_hist_TM_05_25);
202      normalize(_hist_TM_05);
203      normalize(_hist_TM_25_50);
204      normalize(_hist_TM_25);
205      normalize(_hist_TM_50_75);
206      normalize(_hist_TM_50);
207      normalize(_hist_TM_75_100);
208      normalize(_hist_TM_75);
209      normalize(_hist_TM_100);
210
211      normalize(_hist_S_05_25);
212      normalize(_hist_S_05);
213      normalize(_hist_S_25_50);
214      normalize(_hist_S_25);
215      normalize(_hist_S_50_75);
216      normalize(_hist_S_50);
217      normalize(_hist_S_75_100);
218      normalize(_hist_S_75);
219      normalize(_hist_S_100);
220    }
221
222
223  private:
224
225    Histo1DPtr _hist_T_05_25, _hist_T_05, _hist_T_25_50, _hist_T_25, _hist_T_50_75, _hist_T_50, _hist_T_75_100, _hist_T_75, _hist_T_100;
226    Histo1DPtr _hist_TM_05_25, _hist_TM_05, _hist_TM_25_50, _hist_TM_25, _hist_TM_50_75, _hist_TM_50, _hist_TM_75_100, _hist_TM_75, _hist_TM_100;
227    Histo1DPtr _hist_S_05_25, _hist_S_05, _hist_S_25_50, _hist_S_25, _hist_S_50_75, _hist_S_50, _hist_S_75_100, _hist_S_75, _hist_S_100;
228    Profile1DPtr _hist_T_N, _hist_TM_N, _hist_S_N;
229    Profile1DPtr _hist_T_S, _hist_TM_S, _hist_S_S;
230
231  };
232
233
234  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1124167);
235
236}