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