1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  class ATLAS_2012_I1091481 : public Analysis {
  9  public:
 10
 11    /// Constructor
 12    ATLAS_2012_I1091481()
 13      : Analysis("ATLAS_2012_I1091481")
 14    {   }
 15
 16
 17    /// Book histograms and initialise projections before the run
 18    void init() {
 19
 20      ChargedFinalState cfs100(Cuts::abseta < 2.5 && Cuts::pT > 0.1*GeV);
 21      declare(cfs100,"CFS100");
 22      ChargedFinalState cfs500(Cuts::abseta < 2.5 && Cuts::pT > 0.5*GeV);
 23      declare(cfs500,"CFS500");
 24
 25      // collision energy
 26      int isqrts = -1;
 27      if (isCompatibleWithSqrtS(900*GeV))  isqrts = 2;
 28      if (isCompatibleWithSqrtS(7000*GeV)) isqrts = 1;
 29      assert(isqrts > 0);
 30
 31      book(_sE_10_100   ,isqrts, 1, 1);
 32      book(_sE_1_100    ,isqrts, 1, 2);
 33      book(_sE_10_500   ,isqrts, 1, 3);
 34
 35      book(_sEta_10_100 ,isqrts, 2, 1);
 36      book(_sEta_1_100  ,isqrts, 2, 2);
 37      book(_sEta_10_500 ,isqrts, 2, 3);
 38
 39      book(norm_inclusive, "norm_inclusive");
 40      book(norm_lowPt, "norm_lowPt");
 41      book(norm_pt500, "norm_pt500");
 42    }
 43
 44
 45    // Recalculate particle energy assuming pion mass
 46    double getPionEnergy(const Particle& p) {
 47      double m_pi = 0.1396*GeV;
 48      double p2 = p.p3().mod2()/(GeV*GeV);
 49      return sqrt(sqr(m_pi) + p2);
 50    }
 51
 52
 53    // S_eta core for one event
 54    //
 55    //  -1 + 1/Nch * |sum_j^Nch exp[i*(xi eta_j - Phi_j)]|^2
 56    //
 57    double getSeta(const Particles& part, double xi) {
 58      std::complex<double> c_eta (0.0, 0.0);
 59      for (const Particle& p : part) {
 60        double eta = p.eta();
 61        double phi = p.phi();
 62        double arg = xi*eta-phi;
 63         std::complex<double> temp(cos(arg), sin(arg));
 64         c_eta += temp;
 65      }
 66      return std::norm(c_eta)/part.size() - 1.0;
 67    }
 68
 69
 70    // S_E core for one event
 71    //
 72    //  -1 + 1/Nch * |sum_j^Nch exp[i*(omega X_j - Phi_j)]|^2
 73    //
 74    double getSE(const Particles& part, double omega) {
 75      double Xj = 0.0;
 76      std::complex<double> c_E (0.0, 0.0);
 77      for (unsigned int i=0; i < part.size(); ++i) {
 78        Xj += 0.5*getPionEnergy(part[i]);
 79        double phi = part[i].phi();
 80        double arg = omega*Xj - phi;
 81        std::complex<double> temp(cos(arg), sin(arg));
 82        c_E += temp;
 83        Xj += 0.5*getPionEnergy(part[i]);
 84      }
 85      return std::norm(c_E)/part.size() - 1.0;
 86    }
 87
 88
 89    // Convenient fill function
 90    void fillS(Histo1DPtr h, const Particles& part, bool SE=true) {
 91      // Loop over bins, take bin centers as parameter values
 92      for(size_t i=0; i < h->numBins(); ++i) {
 93        double x = h->bin(i).xMid();
 94        double width = h->bin(i).xMax() - h->bin(i).xMin();
 95        double y;
 96        if(SE)  y = getSE(part,   x);
 97        else    y = getSeta(part, x);
 98        h->fill(x, y * width);
 99        // Histo1D objects will be converted to Estimate1D objects for plotting
100        // As part of this conversion, Rivet will divide by bin width
101        // However, we want the (x,y) of the Estimate1D to be the (binCenter, sumW) of
102        // the current Histo1D. This is why in the above line we multiply by bin width,
103        // so as to undo later division by bin width.
104        //
105        // Could have used Estimate1D objects in the first place, but they cannot be merged
106        // as easily as Histo1Ds can using yodamerge (missing ScaledBy attribute)
107        }
108    }
109
110
111    /// Perform the per-event analysis
112    void analyze(const Event& event) {
113      // Charged fs
114      const ChargedFinalState& cfs100  = apply<ChargedFinalState>(event, "CFS100");
115      const Particles          part100 = cfs100.particles(cmpMomByEta);
116      const ChargedFinalState& cfs500  = apply<ChargedFinalState>(event, "CFS500");
117      const Particles&         part500 = cfs500.particles(cmpMomByEta);
118
119      // Veto event if the most inclusive phase space has less than 10 particles and the max pT is > 10 GeV
120      if (part100.size() < 11) vetoEvent;
121      double ptmax = cfs100.particlesByPt()[0].pT()/GeV;
122      if (ptmax > 10.0) vetoEvent;
123
124      // Fill the pt>100, pTmax<10 GeV histos
125      fillS(_sE_10_100, part100, true);
126      fillS(_sEta_10_100, part100, false);
127      norm_inclusive->fill();
128
129      // Fill the pt>100, pTmax<1 GeV histos
130      if (ptmax < 1.0) {
131        fillS(_sE_1_100,   part100, true);
132        fillS(_sEta_1_100, part100, false);
133        norm_lowPt->fill();
134      }
135
136      // Fill the pt>500, pTmax<10 GeV histos
137      if (part500.size() > 10) {
138        fillS(_sE_10_500,   part500, true );
139        fillS(_sEta_10_500, part500, false);
140        norm_pt500->fill();
141      }
142    }
143
144    /// Normalise histograms etc., after the run
145    void finalize() {
146      // The scaling takes the multiple fills per event into account
147      scale(_sE_10_100, 1.0/ *norm_inclusive);
148      scale(_sE_1_100 , 1.0/ *norm_lowPt);
149      scale(_sE_10_500, 1.0/ *norm_pt500);
150
151      scale(_sEta_10_100, 1.0/ *norm_inclusive);
152      scale(_sEta_1_100 , 1.0/ *norm_lowPt);
153      scale(_sEta_10_500, 1.0/ *norm_pt500);
154    }
155
156
157  private:
158
159    Histo1DPtr _sE_10_100;
160    Histo1DPtr _sE_1_100;
161    Histo1DPtr _sE_10_500;
162
163    Histo1DPtr _sEta_10_100;
164    Histo1DPtr _sEta_1_100;
165    Histo1DPtr _sEta_10_500;
166
167    CounterPtr norm_inclusive;
168    CounterPtr norm_lowPt;
169    CounterPtr norm_pt500;
170  };
171
172
173  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1091481);
174
175}