FastJet 3.0.2
|
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 }