FastJet 3.0.0
RestFrameNSubjettinessTagger.cc
00001 //STARTHEADER
00002 // $Id: RestFrameNSubjettinessTagger.cc 2616 2011-09-30 18:03:40Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 #include <fastjet/tools/RestFrameNSubjettinessTagger.hh>
00030 #include <fastjet/tools/Boost.hh>
00031 #include <fastjet/ClusterSequence.hh>
00032 #include <sstream>
00033 
00034 using namespace std;
00035 
00036 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00037 
00038 //------------------------------------------------------------------------
00039 // RestFrameNSubjettinessTagger class implementation
00040 //------------------------------------------------------------------------
00041 
00042 //------------------------------------------------------------------------
00043 // tagger description
00044 string RestFrameNSubjettinessTagger::description() const{ 
00045   ostringstream oss;
00046   oss << "RestFrameNSubjettiness tagger that performs clustering in the jet rest frame with " 
00047       << _subjet_def.description() 
00048       << ", supplemented with cuts tau_2 < " << _t2cut 
00049       << " and cos(theta_s) < " << _costscut;
00050   return oss.str();
00051 }
00052 
00053 
00054 //------------------------------------------------------------------------
00055 // action on a single jet
00056 // returns the tagged PseudoJet if successful, 0 otherwise
00057 PseudoJet RestFrameNSubjettinessTagger::result(const PseudoJet & jet) const{
00058   // make sure that the jet has constituents
00059   if (!jet.has_constituents())
00060     throw("The jet you try to tag needs to have accessible constituents");
00061    
00062   // get the constituents and boost them into the rest frame of the jet
00063   vector<PseudoJet> rest_input = jet.constituents();
00064   for (unsigned int i=0; i<rest_input.size(); i++)
00065     rest_input[i].unboost(jet);
00066 
00067   ClusterSequence cs_rest(rest_input, _subjet_def);
00068   vector<PseudoJet> subjets = (_use_exclusive)
00069     ? cs_rest.exclusive_jets(2)
00070     : sorted_by_E(cs_rest.inclusive_jets());
00071 
00072   // impose the cuts in the rest-frame
00073   if (subjets.size()<2) return PseudoJet();
00074 
00075   const PseudoJet &j0 = subjets[0];
00076   const PseudoJet &j1 = subjets[1];
00077 
00078   /// impose the cut on cos(theta_s)
00079   double ct0 = (j0.px()*jet.px() + j0.py()*jet.py() + j0.pz()*jet.pz())
00080     /sqrt(j0.modp2()*jet.modp2());
00081   double ct1 = (j1.px()*jet.px() + j1.py()*jet.py() + j1.pz()*jet.pz())
00082     /sqrt(j1.modp2()*jet.modp2());
00083   if ((ct0 > _costscut) || (ct1 > _costscut)) return PseudoJet();
00084   
00085   // ccompute the 2-subjettiness and impose the coresponding cut
00086   double tau2 = 0.0;
00087   for (unsigned int i=0; i<rest_input.size(); i++)
00088     tau2 += min(dot_product(rest_input[i], j0), 
00089                 dot_product(rest_input[i], j1));
00090 
00091   tau2 *= (2.0/jet.m2());
00092 
00093   if (tau2 > _t2cut) return PseudoJet();
00094 
00095   // We have a positive tag, 
00096   //  - boost everything back into the lab frame
00097   //  - record the info in the interface
00098   // Note that in order to point to the correct Clustersequence, the
00099   // subjets must be taken from the boosted one. We extract that
00100   // through the history index of the rest-frame subjets
00101   ClusterSequence * cs_structure = new ClusterSequence();
00102   Boost boost(jet);
00103   cs_structure->transfer_from_sequence(cs_rest, &boost);
00104   PseudoJet subjet_lab1 = cs_structure->jets()[cs_rest.history()[subjets[0].cluster_hist_index()].jetp_index];
00105   PseudoJet subjet_lab2 = cs_structure->jets()[cs_rest.history()[subjets[0].cluster_hist_index()].jetp_index];
00106     
00107   PseudoJet result = join<StructureType>(subjet_lab1,subjet_lab2);
00108   StructureType * s = (StructureType *) result.structure_non_const_ptr();
00109 //  s->_original_jet = jet;
00110   s->_tau2 = tau2;
00111   s->_costhetas = max(ct0, ct1);
00112 
00113   // keep the rest-frame CS alive
00114   cs_structure->delete_self_when_unused();
00115 
00116   return result;
00117 }
00118 
00119 
00120 FASTJET_END_NAMESPACE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends