FastJet 3.0.5
11-filter.cc
Go to the documentation of this file.
00001 //----------------------------------------------------------------------
00002 /// \file
00003 /// \page Example11 11 - use of filtering
00004 ///
00005 /// fastjet example program to illustrate the use of the
00006 /// fastjet::Filter class
00007 ///
00008 /// We apply different filter examples to either the hardest jet of
00009 /// the given event, or to the composition of the two hardest jets:
00010 ///
00011 ///   - two examples of a filter keeping a fixed number of subjets (as 
00012 ///     in arXiv:0802.2470)
00013 ///   - a "trimmer" i.e. a filter keeping subjets carrying at least a given 
00014 ///     fraction of the pt of the jet (arXiv:0912.1342).
00015 ///   - two examples of filter in combination with background subtraction
00016 ///
00017 /// run it with    : ./11-filter < data/single-event.dat
00018 ///
00019 /// Source code: 11-filter.cc
00020 //----------------------------------------------------------------------
00021 
00022 //STARTHEADER
00023 // $Id: 11-filter.cc 2684 2011-11-14 07:41:44Z soyez $
00024 //
00025 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00026 //
00027 //----------------------------------------------------------------------
00028 // This file is part of FastJet.
00029 //
00030 //  FastJet is free software; you can redistribute it and/or modify
00031 //  it under the terms of the GNU General Public License as published by
00032 //  the Free Software Foundation; either version 2 of the License, or
00033 //  (at your option) any later version.
00034 //
00035 //  The algorithms that underlie FastJet have required considerable
00036 //  development and are described in hep-ph/0512210. If you use
00037 //  FastJet as part of work towards a scientific publication, please
00038 //  include a citation to the FastJet paper.
00039 //
00040 //  FastJet is distributed in the hope that it will be useful,
00041 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00042 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00043 //  GNU General Public License for more details.
00044 //
00045 //  You should have received a copy of the GNU General Public License
00046 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00047 //----------------------------------------------------------------------
00048 //ENDHEADER
00049 
00050 #include <fastjet/PseudoJet.hh>
00051 #include <fastjet/ClusterSequence.hh>
00052 #include <fastjet/Selector.hh>
00053 #include <iostream>
00054 #include "fastjet/tools/Filter.hh"
00055 // the following includes are only needed when combining filtering with subtraction
00056 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
00057 #include "fastjet/ClusterSequenceArea.hh"
00058 #include "fastjet/tools/Subtractor.hh"
00059 
00060 #include <cstdio>   // needed for io
00061 
00062 using namespace fastjet;
00063 using namespace std;
00064 
00065 // a function returning
00066 //   min(Rmax, deltaR_factor * deltaR(j1,j2))
00067 // where j1 and j2 are the 2 subjets of j
00068 // if the jet does not have 2 exactly pieces, Rmax is used.
00069 class DynamicRfilt : public FunctionOfPseudoJet<double>{
00070 public:
00071   // default ctor 
00072   DynamicRfilt(double Rmax, double deltaR_factor) : _Rmax(Rmax), _deltaR_factor(deltaR_factor){}
00073 
00074   // action of the function
00075   double result(const PseudoJet &j) const{
00076     if (! j.has_pieces()) return _Rmax;
00077 
00078     vector<PseudoJet> pieces = j.pieces();
00079     if (! pieces.size()==2) return _Rmax;
00080 
00081     double deltaR = pieces[0].delta_R(pieces[1]);
00082     return min(_Rmax, _deltaR_factor * deltaR);
00083   }
00084 
00085 private:
00086   double _Rmax, _deltaR_factor;
00087 };
00088 
00089 /// an example program showing how to use Filter in FastJet
00090 int main(){
00091   // read in input particles
00092   //----------------------------------------------------------
00093   vector<PseudoJet> input_particles;
00094   
00095   double px, py , pz, E;
00096   while (cin >> px >> py >> pz >> E) {
00097     // create a PseudoJet with these components and put it onto
00098     // back of the input_particles vector
00099     input_particles.push_back(PseudoJet(px,py,pz,E)); 
00100   }
00101  
00102   // get the resulting jets ordered in pt
00103   //----------------------------------------------------------
00104   JetDefinition jet_def(cambridge_algorithm, 1.2);
00105   // the use of a ClusterSequenceArea (instead of a plain ClusterSequence)
00106   // is only needed because we will later combine filtering with area-based
00107   // subtraction
00108   ClusterSequenceArea clust_seq(input_particles, jet_def, 
00109                                 AreaDefinition(active_area_explicit_ghosts));
00110   vector<PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(5.0));
00111 
00112   // label the columns
00113   printf("%5s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt");
00114  
00115   // print out the details for each jet
00116   for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
00117     printf("%5u %15.8f %15.8f %15.8f\n",
00118            i, inclusive_jets[i].rap(), inclusive_jets[i].phi(),
00119            inclusive_jets[i].perp());
00120   }
00121 
00122   // simple test to avoid that the example below crashes:
00123   // make sure there is at least 3 jets above our 5 GeV
00124   if (inclusive_jets.size()<3){
00125     cout << "Please provide an event with at least 3 jets above 5 GeV" << endl;
00126     return 1;
00127   }
00128 
00129   // the sample PseudoJet that we will filter
00130   //  - the hardest jet of the event
00131   //  - the composition of the second and third hardest jets 
00132   ///   (this shows that the Filter can also be applied to a composite jet)
00133   //----------------------------------------------------------
00134   vector<PseudoJet> candidates;
00135   candidates.push_back(inclusive_jets[0]);
00136   candidates.push_back(join(inclusive_jets[1],inclusive_jets[2]));
00137 
00138 
00139   // create 5 filters
00140   //----------------------------------------------------------
00141   vector<Filter> filters;
00142   
00143   // 1.
00144   // the Cambridge/Aachen filter with Rfilt=0.3 (simpliefied version of arXiv:0802.2470)
00145   filters.push_back(Filter(JetDefinition(cambridge_algorithm, 0.3), SelectorNHardest(3)));
00146 
00147   // 2.
00148   // the Cambridge/Aachen filter with Rfilt=min(0.3, 0.5*Rbb) as in arXiv:0802.2470
00149   SharedPtr<DynamicRfilt> dynamic_Rfilt(new DynamicRfilt(0.3, 0.5));
00150   filters.push_back(Filter(dynamic_Rfilt.get(), SelectorNHardest(3)));
00151 
00152   // 3.
00153   // Filtering with a pt cut as for trimming (arXiv:0912.1342)
00154   filters.push_back(Filter(JetDefinition(kt_algorithm, 0.2), SelectorPtFractionMin(0.03)));
00155 
00156   // 4.
00157   // First example of filtering with subtraction of the background: provide rho
00158   // First, estimate the background for the given event
00159   GridMedianBackgroundEstimator bkgd(4.5, 0.55); // uses particles up to |y|=4.5
00160   bkgd.set_particles(input_particles);
00161   double rho = bkgd.rho();
00162   // Then, define the filter
00163   filters.push_back(Filter(JetDefinition(cambridge_algorithm, 0.3), SelectorNHardest(3), rho));
00164 
00165   // 5.
00166   // Second example of filtering with subtraction of the background: set a subtractor
00167   // First, define a subtractor from a background estimator
00168   Subtractor subtractor(&bkgd);
00169   // Then, define the filter
00170   Filter filt(JetDefinition(cambridge_algorithm, 0.3), SelectorNHardest(3));
00171   // Finally, tell the filter about the subtractor
00172   filt.set_subtractor(&subtractor);
00173   filters.push_back(filt);
00174 
00175 
00176   // apply the various filters to the test PseudoJet
00177   // and show the result
00178   //----------------------------------------------------------
00179 
00180   // print out original jet candidates
00181   cout << "\nOriginal jets that will be filtered: " << endl;
00182   for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
00183     const PseudoJet & c = *jit;
00184     cout << "  rap = " << c.rap() << ", phi = " << c.phi() << ", pt = " << c.perp() 
00185          << "  [" << c.description() << "]" <<  endl;
00186   }
00187 
00188   // loop on filters
00189   for (vector<Filter>::iterator it=filters.begin(); it!=filters.end(); it++){
00190     const Filter & f = *it;
00191     cout << "\nUsing filter: " << f.description() << endl;
00192     
00193     // loop on jet candidates
00194     for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
00195       const PseudoJet & c = *jit;
00196       
00197       // apply filter f to jet c      
00198       PseudoJet j = f(c);
00199       
00200       // access properties specific to the Filter
00201       //
00202       // We first make sure that the jet indeed has a structure
00203       // compatible with the result of a Filter (using
00204       // has_structure_of()), and then retrieve the pieces rejected by the
00205       // filter (using structure_of())
00206       assert(j.has_structure_of<Filter>());
00207       const Filter::StructureType & fj_struct = j.structure_of<Filter>();
00208       
00209       // write out result
00210       cout << "  rap = " << j.rap() << ", phi = " << j.phi() << ", pt = " << j.perp() 
00211            << "  [kept: " << j.pieces().size() << ", rejected: "
00212            << fj_struct.rejected().size() << "]" << endl;
00213     }
00214   }
00215 
00216   return 0;
00217 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends