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