FastJet 3.0.2
|
00001 //---------------------------------------------------------------------- 00002 /// \file 00003 /// \page Example07 07 - subtracting jet background contamination 00004 /// 00005 /// fastjet subtraction example program. 00006 /// 00007 /// run it with : ./07-subtraction < data/Pythia-Zp2jets-lhc-pileup-1ev.dat 00008 /// 00009 /// Source code: 07-subtraction.cc 00010 //---------------------------------------------------------------------- 00011 00012 //STARTHEADER 00013 // $Id: 07-subtraction.cc 2684 2011-11-14 07:41:44Z soyez $ 00014 // 00015 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00016 // 00017 //---------------------------------------------------------------------- 00018 // This file is part of FastJet. 00019 // 00020 // FastJet is free software; you can redistribute it and/or modify 00021 // it under the terms of the GNU General Public License as published by 00022 // the Free Software Foundation; either version 2 of the License, or 00023 // (at your option) any later version. 00024 // 00025 // The algorithms that underlie FastJet have required considerable 00026 // development and are described in hep-ph/0512210. If you use 00027 // FastJet as part of work towards a scientific publication, please 00028 // include a citation to the FastJet paper. 00029 // 00030 // FastJet is distributed in the hope that it will be useful, 00031 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00032 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00033 // GNU General Public License for more details. 00034 // 00035 // You should have received a copy of the GNU General Public License 00036 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00037 //---------------------------------------------------------------------- 00038 //ENDHEADER 00039 00040 #include "fastjet/PseudoJet.hh" 00041 #include "fastjet/ClusterSequenceArea.hh" 00042 #include "fastjet/Selector.hh" 00043 #include "fastjet/tools/JetMedianBackgroundEstimator.hh" 00044 #include "fastjet/tools/Subtractor.hh" 00045 #include <iostream> // needed for io 00046 00047 using namespace std; 00048 using namespace fastjet; 00049 00050 int main(){ 00051 00052 // read in input particles 00053 // 00054 // since we use here simulated data we can split the hard event 00055 // from the full (i.e. with pileup added) one 00056 //---------------------------------------------------------- 00057 00058 vector<PseudoJet> hard_event, full_event; 00059 00060 // read in input particles. Keep the hard event generated by PYTHIA 00061 // separated from the full event, so as to be able to gauge the 00062 // "goodness" of the subtraction from the full event, which also 00063 // includes pileup 00064 double particle_maxrap = 5.0; 00065 00066 string line; 00067 int nsub = 0; // counter to keep track of which sub-event we're reading 00068 while (getline(cin, line)) { 00069 istringstream linestream(line); 00070 // take substrings to avoid problems when there are extra "pollution" 00071 // characters (e.g. line-feed). 00072 if (line.substr(0,4) == "#END") {break;} 00073 if (line.substr(0,9) == "#SUBSTART") { 00074 // if more sub events follow, make copy of first one (the hard one) here 00075 if (nsub == 1) hard_event = full_event; 00076 nsub += 1; 00077 } 00078 if (line.substr(0,1) == "#") {continue;} 00079 double px,py,pz,E; 00080 linestream >> px >> py >> pz >> E; 00081 // you can construct 00082 PseudoJet particle(px,py,pz,E); 00083 00084 // push event onto back of full_event vector 00085 if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle); 00086 } 00087 00088 // if we have read in only one event, copy it across here... 00089 if (nsub == 1) hard_event = full_event; 00090 00091 // if there was nothing in the event 00092 if (nsub == 0) { 00093 cerr << "Error: read empty event\n"; 00094 exit(-1); 00095 } 00096 00097 00098 // create a jet definition for the clustering 00099 // We use the anti-kt algorithm with a radius of 0.5 00100 //---------------------------------------------------------- 00101 double R = 0.5; 00102 JetDefinition jet_def(antikt_algorithm, R); 00103 00104 // create an area definition for the clustering 00105 //---------------------------------------------------------- 00106 // ghosts should go up to the acceptance of the detector or 00107 // (with infinite acceptance) at least 2R beyond the region 00108 // where you plan to investigate jets. 00109 double ghost_maxrap = 6.0; 00110 GhostedAreaSpec area_spec(ghost_maxrap); 00111 AreaDefinition area_def(active_area, area_spec); 00112 00113 // run the jet clustering with the above jet and area definitions 00114 // for both the hard and full event 00115 // 00116 // We retrieve the jets above 7 GeV in both case (note that the 00117 // 7-GeV cut we be applied again later on after we subtract the jets 00118 // from the full event) 00119 // ---------------------------------------------------------- 00120 ClusterSequenceArea clust_seq_hard(hard_event, jet_def, area_def); 00121 ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def); 00122 00123 double ptmin = 7.0; 00124 vector<PseudoJet> hard_jets = sorted_by_pt(clust_seq_hard.inclusive_jets(ptmin)); 00125 vector<PseudoJet> full_jets = sorted_by_pt(clust_seq_full.inclusive_jets(ptmin)); 00126 00127 // Now turn to the estimation of the background (for the full event) 00128 // 00129 // There are different ways to do that. In general, this also 00130 // requires clustering the particles that will be handled internally 00131 // in FastJet. 00132 // 00133 // The suggested way to proceed is to use a BackgroundEstimator 00134 // constructed from the following 3 arguments: 00135 // - a jet definition used to cluster the particles. 00136 // . We strongly recommend using the kt or Cambridge/Aachen 00137 // algorithm (a warning will be issued otherwise) 00138 // . The choice of the radius is a bit more subtle. R=0.4 has 00139 // been chosen to limit the impact of hard jets; in samples of 00140 // dominantly sparse events it may cause the UE/pileup to be 00141 // underestimated a little, a slightly larger value (0.5 or 00142 // 0.6) may be better. 00143 // - An area definition for which we recommend the use of explicit 00144 // ghosts (i.e. active_area_explicit_ghosts) 00145 // As mentionned in the area example (06-area.cc), ghosts should 00146 // extend sufficiently far in rapidity to cover the jets used in 00147 // the computation of the background (see also the comment below) 00148 // - A Selector specifying the range over which we will keep the 00149 // jets entering the estimation of the background (you should 00150 // thus make sure the ghosts extend far enough in rapidity to 00151 // cover the range, a warning will be issued otherwise). 00152 // In this particular example, the two hardest jets in the event 00153 // are removed from the background estimation 00154 // ---------------------------------------------------------- 00155 JetDefinition jet_def_bkgd(kt_algorithm, 0.4); 00156 AreaDefinition area_def_bkgd(active_area_explicit_ghosts, 00157 GhostedAreaSpec(ghost_maxrap)); 00158 Selector selector = SelectorAbsRapMax(4.5) * (!SelectorNHardest(2)); 00159 JetMedianBackgroundEstimator bkgd_estimator(selector, jet_def_bkgd, area_def_bkgd); 00160 00161 // To help manipulate the background estimator, we also provide a 00162 // transformer that allows to apply directly the background 00163 // subtraction on the jets. This will use the background estimator 00164 // to compute rho for the jets to be subtracted. 00165 // ---------------------------------------------------------- 00166 Subtractor subtractor(&bkgd_estimator); 00167 00168 // Finally, once we have an event, we can just tell the background 00169 // estimator to use that list of particles 00170 // This could be done directly when declaring the background 00171 // estimator but the usage below can more easily be accomodated to a 00172 // loop over a set of events. 00173 // ---------------------------------------------------------- 00174 bkgd_estimator.set_particles(full_event); 00175 00176 // show a summary of what was done so far 00177 // - the description of the algorithms, areas and ranges used 00178 // - the background properties 00179 // - the jets in the hard event 00180 //---------------------------------------------------------- 00181 cout << "Main clustering:" << endl; 00182 cout << " Ran: " << jet_def.description() << endl; 00183 cout << " Area: " << area_def.description() << endl; 00184 cout << " Particles up to |y|=" << particle_maxrap << endl; 00185 cout << endl; 00186 00187 cout << "Background estimation:" << endl; 00188 cout << " " << bkgd_estimator.description() << endl << endl;; 00189 cout << " Giving, for the full event" << endl; 00190 cout << " rho = " << bkgd_estimator.rho() << endl; 00191 cout << " sigma = " << bkgd_estimator.sigma() << endl; 00192 cout << endl; 00193 00194 cout << "Jets above " << ptmin << " GeV in the hard event (" << hard_event.size() << " particles)" << endl; 00195 cout << "---------------------------------------\n"; 00196 printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area"); 00197 for (unsigned int i = 0; i < hard_jets.size(); i++) { 00198 printf("%5u %15.8f %15.8f %15.8f %15.8f\n", i, 00199 hard_jets[i].rap(), hard_jets[i].phi(), hard_jets[i].perp(), 00200 hard_jets[i].area()); 00201 } 00202 cout << endl; 00203 00204 // Once the background properties have been computed, subtraction 00205 // can be applied on the jets. Subtraction is performed on the 00206 // full 4-vector 00207 // 00208 // We output the jets before and after subtraction 00209 // ---------------------------------------------------------- 00210 cout << "Jets above " << ptmin << " GeV in the full event (" << full_event.size() << " particles)" << endl; 00211 cout << "---------------------------------------\n"; 00212 printf("%5s %15s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area", "rap_sub", "phi_sub", "pt_sub"); 00213 unsigned int idx=0; 00214 00215 // get the subtracted jets 00216 vector<PseudoJet> subtracted_jets = subtractor(full_jets); 00217 00218 for (unsigned int i=0; i<full_jets.size(); i++){ 00219 // re-apply the pt cut 00220 if (subtracted_jets[i].perp2() >= ptmin*ptmin){ 00221 printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", idx, 00222 full_jets[i].rap(), full_jets[i].phi(), full_jets[i].perp(), 00223 full_jets[i].area(), 00224 subtracted_jets[i].rap(), subtracted_jets[i].phi(), 00225 subtracted_jets[i].perp()); 00226 idx++; 00227 } 00228 } 00229 00230 return 0; 00231 }