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  /// ATLAS jet substructure measurement
 10  class ATLAS_2012_I1094564 : public Analysis {
 11  public:
 12
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1094564);
 14
 15
 16    // Returns constituents to make it easier to do the filtering
 17    PseudoJets splitjet(fastjet::PseudoJet jet, double& last_R, const FastJets& fj, bool& unclustered) const {
 18
 19      // Build a new cluster sequence just using the constituents of this jet.
 20      fastjet::ClusterSequence cs(fj.clusterSeq()->constituents(jet), fastjet::JetDefinition(fastjet::cambridge_algorithm, M_PI/2.));
 21
 22      // Get the jet back again
 23      vector<fastjet::PseudoJet> remadeJets = cs.inclusive_jets(0.);
 24
 25      if ( remadeJets.size() != 1 ) return remadeJets;
 26
 27      fastjet::PseudoJet remadeJet = remadeJets[0];
 28      fastjet::PseudoJet parent1, parent2;
 29      unclustered = false;
 30
 31      while ( cs.has_parents(remadeJet, parent1, parent2) ) {
 32        if (parent1.squared_distance(parent2) < 0.09) break;
 33        if (parent1.m2() < parent2.m2()) {
 34          fastjet::PseudoJet tmp;
 35          tmp = parent1; parent1 = parent2; parent2 = tmp;
 36        }
 37
 38        double ktdist = parent1.kt_distance(parent2);
 39        double rtycut2 = 0.3*0.3;
 40        if (parent1.m() < 0.67*remadeJet.m() && ktdist > rtycut2*remadeJet.m2()) {
 41          unclustered = true;
 42          break;
 43        } else {
 44          remadeJet = parent1;
 45        }
 46      }
 47
 48      last_R = 0.5 * sqrt(parent1.squared_distance(parent2));
 49      return cs.constituents(remadeJet);
 50    }
 51
 52
 53    fastjet::PseudoJet filterjet(PseudoJets jets, double& stingy_R, const double def_R) const {
 54      if (stingy_R == 0.0) stingy_R = def_R;
 55      stingy_R = def_R < stingy_R ? def_R : stingy_R;
 56      fastjet::JetDefinition stingy_jet_def(fastjet::cambridge_algorithm, stingy_R);
 57      fastjet::ClusterSequence scs(jets, stingy_jet_def);
 58      vector<fastjet::PseudoJet> stingy_jets = sorted_by_pt(scs.inclusive_jets(0));
 59      fastjet::PseudoJet reconst_jet(0, 0, 0, 0);
 60      for (size_t isj = 0; isj < std::min((size_t) 3, stingy_jets.size()); ++isj) {
 61        reconst_jet += stingy_jets[isj];
 62      }
 63      return reconst_jet;
 64    }
 65
 66
 67    // These are custom functions for n-subjettiness.
 68    PseudoJets jetGetAxes(int n_jets, const PseudoJets& inputJets, double subR) const {
 69      // Sanity check
 70      if (inputJets.size() < (size_t) n_jets) {
 71        MSG_ERROR("Not enough input particles.");
 72        return inputJets;
 73      }
 74
 75      // Get subjets, return
 76      fastjet::ClusterSequence sub_clust_seq(inputJets, fastjet::JetDefinition(fastjet::kt_algorithm, subR, fastjet::E_scheme, fastjet::Best));
 77      return sub_clust_seq.exclusive_jets(n_jets);
 78    }
 79
 80
 81    double jetTauValue(double beta, double jet_rad, const PseudoJets& particles, const PseudoJets& axes, double Rcut) const {
 82      double tauNum = 0.0;
 83      double tauDen = 0.0;
 84
 85      if (particles.size() == 0) return 0.0;
 86
 87      for (size_t i = 0; i < particles.size(); i++) {
 88        // find minimum distance (set R large to begin)
 89        double minR = 10000.0;
 90        for (size_t j = 0; j < axes.size(); j++) {
 91          double tempR = sqrt(particles[i].squared_distance(axes[j]));
 92          if (tempR < minR) minR = tempR;
 93        }
 94        if (minR > Rcut) minR = Rcut;
 95        // calculate nominator and denominator
 96        tauNum += particles[i].perp() * pow(minR,beta);
 97        tauDen += particles[i].perp() * pow(jet_rad,beta);
 98      }
 99
100      // return N-subjettiness (or 0 if denominator is 0)
101      return safediv(tauNum, tauDen, 0);
102    }
103
104
105    void jetUpdateAxes(double beta, const PseudoJets& particles, PseudoJets& axes) const {
106      vector<int> belongsto;
107      for (size_t i = 0; i < particles.size(); i++) {
108        // find minimum distance axis
109        int assign = 0;
110        double minR = 10000.0;
111        for (size_t j = 0; j < axes.size(); j++) {
112          double tempR = sqrt(particles[i].squared_distance(axes[j]));
113          if (tempR < minR) {
114            minR = tempR;
115            assign = j;
116          }
117        }
118        belongsto.push_back(assign);
119      }
120
121      // iterative step
122      double deltaR2, distphi;
123      vector<double> ynom, phinom, den;
124
125      ynom.resize(axes.size());
126      phinom.resize(axes.size());
127      den.resize(axes.size());
128
129      for (size_t i = 0; i < particles.size(); i++) {
130        distphi = particles[i].phi() - axes[belongsto[i]].phi();
131        deltaR2 = particles[i].squared_distance(axes[belongsto[i]]);
132        if (deltaR2 == 0.) continue;
133        if (abs(distphi) <= M_PI) phinom.at(belongsto[i]) += particles[i].perp() * particles[i].phi() * pow(deltaR2, (beta-2)/2);
134        else if ( distphi > M_PI) phinom.at(belongsto[i]) += particles[i].perp() * (-2 * M_PI + particles[i].phi()) * pow(deltaR2, (beta-2)/2);
135        else if ( distphi < M_PI) phinom.at(belongsto[i]) += particles[i].perp() * (+2 * M_PI + particles[i].phi()) * pow(deltaR2, (beta-2)/2);
136        ynom.at(belongsto[i]) += particles[i].perp() * particles[i].rap() * pow(deltaR2, (beta-2)/2);
137        den.at(belongsto[i])  += particles[i].perp() * pow(deltaR2, (beta-2)/2);
138      }
139
140      // reset to new axes
141      for (size_t j = 0; j < axes.size(); j++) {
142        if (den[j] == 0.) axes.at(j) = axes[j];
143        else {
144          double phi_new = fmod( 2*M_PI + (phinom[j] / den[j]), 2*M_PI );
145          double pt_new  = axes[j].perp();
146          double y_new   = ynom[j] / den[j];
147          double px = pt_new * cos(phi_new);
148          double py = pt_new * sin(phi_new);
149          double pz = pt_new * sinh(y_new);
150          axes.at(j).reset(px, py, pz, axes[j].perp()/2);
151        }
152      }
153    }
154
155
156    void init() {
157
158      /// Projections:
159      FinalState fs((Cuts::etaIn(-4.5, 4.5) && Cuts::pT >=  0.*GeV));
160      declare(fs, "FS");
161      declare(FastJets(fs, JetAlg::ANTIKT, 1.0), "AKT");
162      declare(FastJets(fs, JetAlg::CAM, 1.2)   , "CA" );
163
164      /// Histograms:
165      book(_h_camass,   {200., 300., 400., 500., 600.},
166                        {"d01-x01-y01", "d02-x01-y01", "d03-x01-y01", "d04-x01-y01"});
167      book(_h_filtmass, {200., 300., 400., 500., 600.},
168                        {"d05-x01-y01", "d06-x01-y01", "d07-x01-y01", "d08-x01-y01"});
169      book(_h_ktmass,   {200., 300., 400., 500., 600.},
170                        {"d09-x01-y01", "d10-x01-y01", "d11-x01-y01", "d12-x01-y01"});
171      book(_h_ktd12,    {200., 300., 400., 500., 600.},
172                        {"d13-x01-y01", "d14-x01-y01", "d15-x01-y01", "d16-x01-y01"});
173      book(_h_ktd23,    {200., 300., 400., 500., 600.},
174                        {"d17-x01-y01", "d18-x01-y01", "d19-x01-y01", "d20-x01-y01"});
175      book(_h_cat21,    {200., 300., 400., 500., 600.},
176                        {"d21-x01-y01", "d22-x01-y01", "d23-x01-y01", "d24-x01-y01"});
177      book(_h_cat32,    {200., 300., 400., 500., 600.},
178                        {"d25-x01-y01", "d26-x01-y01", "d27-x01-y01", "d28-x01-y01"});
179      book(_h_ktt21,    {200., 300., 400., 500., 600.},
180                        {"d29-x01-y01", "d30-x01-y01", "d31-x01-y01", "d32-x01-y01"});
181      book(_h_ktt32,    {200., 300., 400., 500., 600.},
182                        {"d33-x01-y01", "d34-x01-y01", "d35-x01-y01", "d36-x01-y01"});
183
184    }
185
186
187    /// Perform the per-event analysis
188    void analyze(const Event& event) {
189
190      using namespace fastjet;
191
192      // Get anti-kt jets with p_T > 200 GeV, check abs(y) < 2, and fill mass histograms
193      const FastJets& ktfj = apply<FastJets>(event, "AKT");
194      PseudoJets ktjets = ktfj.pseudojetsByPt(200*GeV);
195      for (const PseudoJet & ajet : ktjets) {
196        if (abs(ajet.rap()) < 2) {
197          _h_ktmass->fill(ajet.perp(), ajet.m());
198        }
199      }
200
201      // Same as above but C/A jets
202      const FastJets& cafj = apply<FastJets>(event, "CA");
203      PseudoJets cajets = cafj.pseudojetsByPt(200*GeV);
204      for (const PseudoJet & ajet : cajets) {
205        if (abs(ajet.rap()) < 2) {
206          _h_camass->fill(ajet.perp(), ajet.m());
207        }
208      }
209
210      // Split and filter.
211      // Only do this to C/A jets in this analysis.
212      for (const PseudoJet & pjet : cajets) {
213        if ( pjet.perp() > 600 || abs(pjet.rap()) > 2) continue;
214        double dR = 0;
215        bool unclustered = false;
216        PseudoJets split_jets = splitjet(pjet, dR, cafj, unclustered);
217        if ( (dR < 0.15) || (unclustered == false) ) continue;
218        PseudoJet filt_jet = filterjet(split_jets, dR, 0.3);
219        _h_filtmass->fill(filt_jet.perp(), filt_jet.m());
220      }
221
222      // Use the two last stages of clustering to get sqrt(d_12) and sqrt(d_23).
223      // Only use anti-kt jets in this analysis.
224      for (const PseudoJet & pjet : ktjets) {
225        if (pjet.perp() > 600 || abs(pjet.rap()) > 2) continue;
226        ClusterSequence subjet_cseq(ktfj.clusterSeq()->constituents(pjet), JetDefinition(kt_algorithm, M_PI/2.));
227        double d_12 = subjet_cseq.exclusive_dmerge(1) * M_PI*M_PI/4.;
228        double d_23 = subjet_cseq.exclusive_dmerge(2) * M_PI*M_PI/4.;
229        _h_ktd12->fill(pjet.perp(), sqrt(d_12));
230        _h_ktd23->fill(pjet.perp(), sqrt(d_23));
231      }
232
233      // N-subjettiness, use beta = 1 (no rationale given).
234      // Uses the functions defined above.
235      // C/A jets first, anti-kt after.
236      double beta = 1.;
237      //Rcut is used for particles that are very far from the closest axis. At 10
238      //is has no impact on the outcome of the calculation
239      double Rcut = 10.;
240      for (const PseudoJet & pjet : cajets) {
241        if (pjet.perp() > 600*GeV || fabs(pjet.rap()) > 2) continue;
242
243        const PseudoJets constituents = cafj.clusterSeq()->constituents(pjet);
244        if (constituents.size() < 3) continue;
245
246        const PseudoJets axis1 = jetGetAxes(1, constituents, M_PI/2.0);
247        const PseudoJets axis2 = jetGetAxes(2, constituents, M_PI/2.0);
248        const PseudoJets axis3 = jetGetAxes(3, constituents, M_PI/2.0);
249
250        const double radius = 1.2;
251        const double tau1 = jetTauValue(beta, radius, constituents, axis1, Rcut);
252        const double tau2 = jetTauValue(beta, radius, constituents, axis2, Rcut);
253        const double tau3 = jetTauValue(beta, radius, constituents, axis3, Rcut);
254
255        if (tau1 == 0 || tau2 == 0) continue;
256        _h_cat21->fill(pjet.perp(), tau2/tau1);
257        _h_cat32->fill(pjet.perp(), tau3/tau2);
258      }
259
260      for (const PseudoJet & pjet : ktjets) {
261        if (pjet.perp() > 600*GeV || fabs(pjet.rap()) > 2) continue;
262
263        const PseudoJets constituents = ktfj.clusterSeq()->constituents(pjet);
264        if (constituents.size() < 3) continue;
265
266        const PseudoJets axis1 = jetGetAxes(1, constituents, M_PI/2.0);
267        const PseudoJets axis2 = jetGetAxes(2, constituents, M_PI/2.0);
268        const PseudoJets axis3 = jetGetAxes(3, constituents, M_PI/2.0);
269
270        const double radius = 1.0;
271        const double tau1 = jetTauValue(beta, radius, constituents, axis1, Rcut);
272        const double tau2 = jetTauValue(beta, radius, constituents, axis2, Rcut);
273        const double tau3 = jetTauValue(beta, radius, constituents, axis3, Rcut);
274        if (tau1 == 0 || tau2 == 0) continue;
275
276        _h_ktt21->fill(pjet.perp(), tau2/tau1);
277        _h_ktt32->fill(pjet.perp(), tau3/tau2);
278      }
279    }
280
281
282    /// Normalise histograms etc., after the run
283    void finalize() {
284      normalize(_h_camass);
285      normalize(_h_filtmass);
286      normalize(_h_ktmass);
287      normalize(_h_ktd12);
288      normalize(_h_ktd23);
289      normalize(_h_cat21);
290      normalize(_h_cat32);
291      normalize(_h_ktt21);
292      normalize(_h_ktt32);
293    }
294
295  private:
296
297    Histo1DGroupPtr _h_camass;
298    Histo1DGroupPtr _h_filtmass;
299    Histo1DGroupPtr _h_ktmass;
300    Histo1DGroupPtr _h_ktd12;
301    Histo1DGroupPtr _h_ktd23;
302    Histo1DGroupPtr _h_cat21;
303    Histo1DGroupPtr _h_cat32;
304    Histo1DGroupPtr _h_ktt21;
305    Histo1DGroupPtr _h_ktt32;
306
307  };
308
309
310  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1094564);
311
312}