FastJet 3.0.5
14-groomers.cc
Go to the documentation of this file.
00001 //----------------------------------------------------------------------
00002 /// \file
00003 /// \page Example14 14 - unified use of transformers
00004 ///
00005 /// fastjet example program to illustrate the use of the
00006 /// fastjet::Filter and fastjet::Pruner classes in a unified way
00007 /// through their derivation from fastjet::Transformer
00008 ///
00009 /// The two hardest jets in a boosted top event, clustered with an
00010 /// (abnormally) large R, are then groomed using different tools. One
00011 /// notes the reduction in the mass of the jets after grooming.
00012 ///
00013 /// run it with    : ./14-groomers < data/boosted_top_event.dat
00014 ///
00015 /// Source code: 14-groomers.cc
00016 //----------------------------------------------------------------------
00017 
00018 //STARTHEADER
00019 // $Id: 14-groomers.cc 2684 2011-11-14 07:41:44Z soyez $
00020 //
00021 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00022 //
00023 //----------------------------------------------------------------------
00024 // This file is part of FastJet.
00025 //
00026 //  FastJet is free software; you can redistribute it and/or modify
00027 //  it under the terms of the GNU General Public License as published by
00028 //  the Free Software Foundation; either version 2 of the License, or
00029 //  (at your option) any later version.
00030 //
00031 //  The algorithms that underlie FastJet have required considerable
00032 //  development and are described in hep-ph/0512210. If you use
00033 //  FastJet as part of work towards a scientific publication, please
00034 //  include a citation to the FastJet paper.
00035 //
00036 //  FastJet is distributed in the hope that it will be useful,
00037 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00038 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00039 //  GNU General Public License for more details.
00040 //
00041 //  You should have received a copy of the GNU General Public License
00042 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00043 //----------------------------------------------------------------------
00044 //ENDHEADER
00045 
00046 #include <fastjet/PseudoJet.hh>
00047 #include <fastjet/ClusterSequence.hh>
00048 #include <fastjet/Selector.hh>
00049 #include <iostream>
00050 #include "fastjet/tools/Filter.hh"
00051 #include "fastjet/tools/Pruner.hh"
00052 
00053 #include <cstdio>   // needed for io
00054 
00055 using namespace fastjet;
00056 using namespace std;
00057 
00058 /// an example program showing how to use Filter and Pruner in FastJet
00059 int main(){
00060   // read in input particles
00061   //----------------------------------------------------------
00062   vector<PseudoJet> input_particles;
00063   
00064   double px, py , pz, E;
00065   while (cin >> px >> py >> pz >> E) {
00066     // create a PseudoJet with these components and put it onto
00067     // back of the input_particles vector
00068     input_particles.push_back(PseudoJet(px,py,pz,E)); 
00069   }
00070  
00071   // get the resulting jets ordered in pt
00072   //----------------------------------------------------------
00073   JetDefinition jet_def(cambridge_algorithm, 1.5);
00074   ClusterSequence clust_seq(input_particles, jet_def);
00075   vector<PseudoJet> inclusive_jets = 
00076                              sorted_by_pt(clust_seq.inclusive_jets(5.0));
00077 
00078   // label the columns
00079   printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "mass");
00080  
00081   // print out the details for each jet
00082   for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
00083     printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
00084            i, inclusive_jets[i].rap(), inclusive_jets[i].phi(),
00085            inclusive_jets[i].perp(),inclusive_jets[i].m());
00086   }
00087 
00088   // simple test to avoid that the example below crashes:
00089   // make sure there is at least 3 jets above our 5 GeV
00090   if (inclusive_jets.size()<2){
00091     cout << "Please provide an event with at least 2 jets above 5 GeV" << endl;
00092     return 1;
00093   }
00094 
00095   // We will groom the two hardest jets of the event
00096   //----------------------------------------------------------
00097   vector<PseudoJet> candidates = SelectorNHardest(2)(inclusive_jets);
00098 
00099   // create 3 groomers
00100   //----------------------------------------------------------
00101   vector<Transformer *> groomers;
00102   
00103   // 1.
00104   // the Cambridge/Aachen filter with Rfilt=0.3 
00105   // (simplified version of arXiv:0802.2470)
00106   double Rfilt = 0.3;
00107   unsigned int nfilt = 3;
00108   groomers.push_back(new Filter(JetDefinition(cambridge_algorithm, Rfilt), 
00109                                 SelectorNHardest(nfilt) ) );
00110 
00111   // 2.
00112   // Filtering with a pt cut as for trimming (arXiv:0912.1342)
00113   double Rtrim = 0.2;
00114   double ptfrac = 0.03;
00115   groomers.push_back(new Filter(JetDefinition(kt_algorithm, Rtrim), 
00116                                 SelectorPtFractionMin(ptfrac) ) );
00117 
00118   // 3.
00119   // Pruning (arXiv:0903.5081)
00120   double zcut = 0.1;
00121   double rcut_factor = 0.5;
00122   groomers.push_back(new Pruner(cambridge_algorithm, zcut, rcut_factor));
00123 
00124   // apply the various groomers to the test PseudoJet's
00125   // and show the result
00126   //----------------------------------------------------------
00127 
00128   // print out original jet candidates
00129   cout << "\nOriginal jets that will be grooomed: " << endl;
00130   for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
00131     const PseudoJet & c = *jit;
00132     cout << "  rap = " << c.rap() << ", phi = " << c.phi() << ", pt = " << c.perp() 
00133          << ", mass = " << c.m() 
00134          << "  [" << c.description() << "]" <<  endl;
00135   }
00136 
00137   // loop on groomers
00138   for (unsigned int i=0; i < groomers.size(); i++){
00139     const Transformer & f = *groomers[i];
00140     cout << "\nUsing groomer: " << f.description() << endl;
00141     
00142     // loop on jet candidates
00143     for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
00144       const PseudoJet & c = *jit;
00145       
00146       // apply groomer f to jet c      
00147       PseudoJet j = f(c);
00148       
00149       // access properties specific to the given transformer
00150       //
00151       // We first make sure that the jet indeed has a structure
00152       // compatible with the result of a Filter or Pruner (using
00153       // has_structure_of()), and then retrieve the pieces rejected by the
00154       // groomer (using structure_of())
00155       int n_rejected;         
00156       if (j.has_structure_of<Filter>()) {
00157          const Filter::StructureType & fj_struct = j.structure_of<Filter>(); 
00158          n_rejected = fj_struct.rejected().size();
00159       }else { 
00160          assert(j.has_structure_of<Pruner>());  // make sure
00161          const Pruner::StructureType & fj_struct = j.structure_of<Pruner>();
00162          n_rejected = fj_struct.rejected().size();
00163       }
00164       
00165       // write out result
00166       cout << "  rap = " << j.rap() << ", phi = " << j.phi() << ", pt = " << j.perp()
00167            << " mass = " << j.m() 
00168            << "  [kept: " << j.pieces().size() << ", rejected: "
00169            << n_rejected << "]" << endl;
00170     }
00171   }
00172 
00173   // a bit of memory cleaning
00174   for (unsigned int i=0; i < groomers.size(); i++) delete groomers[i];
00175 
00176   return 0;
00177 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends