FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
RestFrameNSubjettinessTagger.cc
1 //STARTHEADER
2 // $Id: RestFrameNSubjettinessTagger.cc 2689 2011-11-14 14:51:06Z soyez $
3 //
4 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #include <fastjet/tools/RestFrameNSubjettinessTagger.hh>
30 #include <fastjet/tools/Boost.hh>
31 #include <fastjet/ClusterSequence.hh>
32 #include <sstream>
33 
34 using namespace std;
35 
36 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37 
38 //------------------------------------------------------------------------
39 // RestFrameNSubjettinessTagger class implementation
40 //------------------------------------------------------------------------
41 
42 //------------------------------------------------------------------------
43 // tagger description
44 string RestFrameNSubjettinessTagger::description() const{
45  ostringstream oss;
46  oss << "RestFrameNSubjettiness tagger that performs clustering in the jet rest frame with "
47  << _subjet_def.description()
48  << ", supplemented with cuts tau_2 < " << _t2cut
49  << " and cos(theta_s) < " << _costscut;
50  return oss.str();
51 }
52 
53 
54 //------------------------------------------------------------------------
55 // action on a single jet
56 // returns the tagged PseudoJet if successful, 0 otherwise
57 PseudoJet RestFrameNSubjettinessTagger::result(const PseudoJet & jet) const{
58  // make sure that the jet has constituents
59  if (!jet.has_constituents())
60  throw("The jet you try to tag needs to have accessible constituents");
61 
62  // get the constituents and boost them into the rest frame of the jet
63  vector<PseudoJet> rest_input = jet.constituents();
64  for (unsigned int i=0; i<rest_input.size(); i++)
65  rest_input[i].unboost(jet);
66 
67  ClusterSequence cs_rest(rest_input, _subjet_def);
68  vector<PseudoJet> subjets = (_use_exclusive)
69  ? cs_rest.exclusive_jets(2)
70  : sorted_by_E(cs_rest.inclusive_jets());
71 
72  // impose the cuts in the rest-frame
73  if (subjets.size()<2) return PseudoJet();
74 
75  const PseudoJet &j0 = subjets[0];
76  const PseudoJet &j1 = subjets[1];
77 
78  /// impose the cut on cos(theta_s)
79  double ct0 = (j0.px()*jet.px() + j0.py()*jet.py() + j0.pz()*jet.pz())
80  /sqrt(j0.modp2()*jet.modp2());
81  double ct1 = (j1.px()*jet.px() + j1.py()*jet.py() + j1.pz()*jet.pz())
82  /sqrt(j1.modp2()*jet.modp2());
83  if ((ct0 > _costscut) || (ct1 > _costscut)) return PseudoJet();
84 
85  // ccompute the 2-subjettiness and impose the coresponding cut
86  double tau2 = 0.0;
87  for (unsigned int i=0; i<rest_input.size(); i++)
88  tau2 += min(dot_product(rest_input[i], j0),
89  dot_product(rest_input[i], j1));
90 
91  tau2 *= (2.0/jet.m2());
92 
93  if (tau2 > _t2cut) return PseudoJet();
94 
95  // We have a positive tag,
96  // - boost everything back into the lab frame
97  // - record the info in the interface
98  // Note that in order to point to the correct Clustersequence, the
99  // subjets must be taken from the boosted one. We extract that
100  // through the history index of the rest-frame subjets
101  ClusterSequence * cs_structure = new ClusterSequence();
102  Boost boost(jet);
103  cs_structure->transfer_from_sequence(cs_rest, &boost);
104  PseudoJet subjet_lab1 = cs_structure->jets()[cs_rest.history()[subjets[0].cluster_hist_index()].jetp_index];
105  PseudoJet subjet_lab2 = cs_structure->jets()[cs_rest.history()[subjets[0].cluster_hist_index()].jetp_index];
106 
107  PseudoJet result_local = join<StructureType>(subjet_lab1,subjet_lab2);
108  StructureType * s = (StructureType *) result_local.structure_non_const_ptr();
109 // s->_original_jet = jet;
110  s->_tau2 = tau2;
111  s->_costhetas = max(ct0, ct1);
112 
113  // keep the rest-frame CS alive
114  cs_structure->delete_self_when_unused();
115 
116  return result_local;
117 }
118 
119 
120 FASTJET_END_NAMESPACE