FastJet  3.3.3
RestFrameNSubjettinessTagger.cc
1 //FJSTARTHEADER
2 // $Id: RestFrameNSubjettinessTagger.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2005-2019, 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. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 #include <fastjet/tools/RestFrameNSubjettinessTagger.hh>
32 #include <fastjet/tools/Boost.hh>
33 #include <fastjet/ClusterSequence.hh>
34 #include <sstream>
35 
36 using namespace std;
37 
38 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39 
40 //------------------------------------------------------------------------
41 // RestFrameNSubjettinessTagger class implementation
42 //------------------------------------------------------------------------
43 
44 //------------------------------------------------------------------------
45 // tagger description
46 string RestFrameNSubjettinessTagger::description() const{
47  ostringstream oss;
48  oss << "RestFrameNSubjettiness tagger that performs clustering in the jet rest frame with "
49  << _subjet_def.description()
50  << ", supplemented with cuts tau_2 < " << _t2cut
51  << " and cos(theta_s) < " << _costscut;
52  return oss.str();
53 }
54 
55 
56 //------------------------------------------------------------------------
57 // action on a single jet
58 // returns the tagged PseudoJet if successful, 0 otherwise
59 PseudoJet RestFrameNSubjettinessTagger::result(const PseudoJet & jet) const{
60  // make sure that the jet has constituents
61  if (!jet.has_constituents())
62  throw("The jet you try to tag needs to have accessible constituents");
63 
64  // get the constituents and boost them into the rest frame of the jet
65  vector<PseudoJet> rest_input = jet.constituents();
66  for (unsigned int i=0; i<rest_input.size(); i++)
67  rest_input[i].unboost(jet);
68 
69  ClusterSequence cs_rest(rest_input, _subjet_def);
70  vector<PseudoJet> subjets = (_use_exclusive)
71  ? cs_rest.exclusive_jets(2)
72  : sorted_by_E(cs_rest.inclusive_jets());
73 
74  // impose the cuts in the rest-frame
75  if (subjets.size()<2) return PseudoJet();
76 
77  const PseudoJet &j0 = subjets[0];
78  const PseudoJet &j1 = subjets[1];
79 
80  /// impose the cut on cos(theta_s)
81  double ct0 = (j0.px()*jet.px() + j0.py()*jet.py() + j0.pz()*jet.pz())
82  /sqrt(j0.modp2()*jet.modp2());
83  double ct1 = (j1.px()*jet.px() + j1.py()*jet.py() + j1.pz()*jet.pz())
84  /sqrt(j1.modp2()*jet.modp2());
85  if ((ct0 > _costscut) || (ct1 > _costscut)) return PseudoJet();
86 
87  // ccompute the 2-subjettiness and impose the coresponding cut
88  double tau2 = 0.0;
89  for (unsigned int i=0; i<rest_input.size(); i++)
90  tau2 += min(dot_product(rest_input[i], j0),
91  dot_product(rest_input[i], j1));
92 
93  tau2 *= (2.0/jet.m2());
94 
95  if (tau2 > _t2cut) return PseudoJet();
96 
97  // We have a positive tag,
98  // - boost everything back into the lab frame
99  // - record the info in the interface
100  // Note that in order to point to the correct Clustersequence, the
101  // subjets must be taken from the boosted one. We extract that
102  // through the history index of the rest-frame subjets
103  ClusterSequence * cs_structure = new ClusterSequence();
104  Boost boost(jet);
105  cs_structure->transfer_from_sequence(cs_rest, &boost);
106  PseudoJet subjet_lab1 = cs_structure->jets()[cs_rest.history()[subjets[0].cluster_hist_index()].jetp_index];
107  PseudoJet subjet_lab2 = cs_structure->jets()[cs_rest.history()[subjets[0].cluster_hist_index()].jetp_index];
108 
109  PseudoJet result_local = join<StructureType>(subjet_lab1,subjet_lab2);
110  StructureType * s = (StructureType *) result_local.structure_non_const_ptr();
111 // s->_original_jet = jet;
112  s->_tau2 = tau2;
113  s->_costhetas = max(ct0, ct1);
114 
115  // keep the rest-frame CS alive
116  cs_structure->delete_self_when_unused();
117 
118  return result_local;
119 }
120 
121 
122 FASTJET_END_NAMESPACE
double m2() const
returns the squared invariant mass // like CLHEP
Definition: PseudoJet.hh:150
deals with clustering
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
the structure returned by the RestFrameNSubjettinessTagger transformer.
PseudoJetStructureBase * structure_non_const_ptr()
return a non-const pointer to the structure (of type PseudoJetStructureBase*) associated with this Ps...
Definition: PseudoJet.cc:498
std::vector< PseudoJet > inclusive_jets(const double ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
double modp2() const
return the squared 3-vector modulus = px^2+py^2+pz^2
Definition: PseudoJet.hh:165
double dot_product(const PseudoJet &a, const PseudoJet &b)
returns the 4-vector dot product of a and b
Definition: PseudoJet.hh:858
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:584
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition: PseudoJet.cc:790
virtual bool has_constituents() const
returns true if the PseudoJet has constituents
Definition: PseudoJet.cc:578
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
Class to boost a PseudoJet.
Definition: Boost.hh:47
std::vector< PseudoJet > exclusive_jets(const double dcut) const
return a vector of all jets (in the sense of the exclusive algorithm) that would be obtained when run...
void transfer_from_sequence(const ClusterSequence &from_seq, const FunctionOfPseudoJet< PseudoJet > *action_on_jets=0)
transfer the sequence contained in other_seq into our own; any plugin "extras" contained in the from_...