FastJet 3.0.5
07-subtraction.cc
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends