FastJet 3.0.2
fastjet_subtraction.cc
00001 
00002 //STARTHEADER
00003 // $Id: fastjet_subtraction.cc 2577 2011-09-13 15:11:38Z salam $
00004 //
00005 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00006 //
00007 //----------------------------------------------------------------------
00008 // This file is part of FastJet.
00009 //
00010 //  FastJet is free software; you can redistribute it and/or modify
00011 //  it under the terms of the GNU General Public License as published by
00012 //  the Free Software Foundation; either version 2 of the License, or
00013 //  (at your option) any later version.
00014 //
00015 //  The algorithms that underlie FastJet have required considerable
00016 //  development and are described in hep-ph/0512210. If you use
00017 //  FastJet as part of work towards a scientific publication, please
00018 //  include a citation to the FastJet paper.
00019 //
00020 //  FastJet is distributed in the hope that it will be useful,
00021 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00022 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00023 //  GNU General Public License for more details.
00024 //
00025 //  You should have received a copy of the GNU General Public License
00026 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00027 //----------------------------------------------------------------------
00028 //ENDHEADER
00029 
00030 
00031 
00032 
00033 
00034 
00035 
00036 //----------------------------------------------------------------------
00037 // fastjet subtraction example program. 
00038 //
00039 // Compile it with: make fastjet_subtraction
00040 // run it with    : ./fastjet_subtraction < data/Pythia-dijet-ptmin100-lhc-pileup-1ev.dat
00041 //
00042 //----------------------------------------------------------------------
00043 #include "fastjet/PseudoJet.hh"
00044 #include "fastjet/ClusterSequenceArea.hh"
00045 #include<iostream> // needed for io
00046 #include<sstream>  // needed for internal io
00047 #include<vector> 
00048 
00049 using namespace std;
00050 
00051 // A declaration of a function that pretty prints a list of jets
00052 // The subtraction is also performed inside this function
00053 void print_jets (const fastjet::ClusterSequenceAreaBase &, 
00054                  const vector<fastjet::PseudoJet> &);
00055 
00056 /// an example program showing how to use fastjet
00057 int main (int argc, char ** argv) {
00058   
00059   // since we use here simulated data we can split the hard event
00060   // from the full (i.e. with pileup added) one
00061   vector<fastjet::PseudoJet> hard_event, full_event;
00062   
00063   // maximum rapidity that we accept for the input particles
00064   double etamax = 6.0;
00065   
00066   // read in input particles. Keep the hard event generated by PYTHIA
00067   // separated from the full event, so as to be able to gauge the
00068   // "goodness" of the subtraction from the full event which also
00069   // includes pileup
00070   string line;
00071   int  nsub  = 0;
00072   while (getline(cin, line)) {
00073     istringstream linestream(line);
00074     // take substrings to avoid problems when there are extra "pollution"
00075     // characters (e.g. line-feed).
00076     if (line.substr(0,4) == "#END") {break;}
00077     if (line.substr(0,9) == "#SUBSTART") {
00078       // if more sub events follow, make copy of hard one here
00079       if (nsub == 1) hard_event = full_event;
00080       nsub += 1;
00081     }
00082     if (line.substr(0,1) == "#") {continue;}
00083     valarray<double> fourvec(4);
00084     linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00085     fastjet::PseudoJet psjet(fourvec);
00086     psjet.set_user_index(0);
00087     // push event onto back of full_event vector
00088     if (abs(psjet.rap()) < etamax) {full_event.push_back(psjet);}
00089   }
00090 
00091   // if we have read in only one event, copy it across here...
00092   if (nsub == 1) hard_event = full_event;
00093 
00094   // if there was nothing in the event 
00095   if (nsub == 0) {
00096     cerr << "Error: read empty event\n";
00097     exit(-1);
00098   }
00099   
00100   
00101   // create an object that represents your choice of jet algorithm and 
00102   // the associated parameters
00103   double R = 0.7;
00104   fastjet::Strategy strategy = fastjet::Best;
00105   fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R, strategy);
00106   //fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, R, strategy);
00107   //fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R, strategy);
00108 
00109   // create an object that specifies how we to define the (active) area
00110   double ghost_etamax = 6.0;
00111   int    active_area_repeats = 1;
00112   double ghost_area    = 0.01;
00113   fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, 
00114                                     ghost_area);
00115   fastjet::AreaDefinition area_def(ghost_spec);
00116   //fastjet::AreaDefinition area_def(fastjet::VoronoiAreaSpec(1.0));
00117 
00118   // run the jet clustering with the above jet definition. hard event first
00119   fastjet::ClusterSequenceArea clust_seq(hard_event, jet_def, area_def);
00120   //fastjet::ClusterSequencePassiveArea clust_seq(hard_event, jet_def);
00121 
00122 
00123   // extract the inclusive jets with pt > 5 GeV, sorted by pt
00124   double ptmin = 5.0;
00125   vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
00126 
00127   // print them out
00128   cout << "" <<endl;
00129   cout << "Hard event only"<<endl;
00130   cout << "Number of input particles: "<<hard_event.size()<<endl;
00131   cout << "Strategy used: "<<clust_seq.strategy_string()<<endl;
00132   cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n";
00133   cout << "---------------------------------------\n";
00134   print_jets(clust_seq, inclusive_jets);
00135   cout << endl;
00136 
00137 
00138   // repeat everything on the full event
00139 
00140   // run the jet clustering with the above jet definition
00141   fastjet::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
00142 
00143   // extract the inclusive jets with pt > 20 GeV, sorted by pt
00144   ptmin = 20.0;
00145   inclusive_jets = clust_seq_full.inclusive_jets(ptmin);
00146 
00147   // print them out
00148   cout << "Full event, with pileup, and its subtraction"<<endl;
00149   cout << "Number of input particles: "<<full_event.size()<<endl;
00150   cout << "Strategy used: "<<clust_seq_full.strategy_string()<<endl;
00151   cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV (before subtraction)\n";
00152   cout << "---------------------------------------\n";
00153   print_jets(clust_seq_full, inclusive_jets);
00154   cout << endl;
00155 
00156 
00157 }
00158 
00159 
00160 //----------------------------------------------------------------------
00161 /// a function that pretty prints a list of jets, and performs the subtraction
00162 /// in two different ways, using a generic ClusterSequenceAreaBase
00163 /// type object.
00164 void print_jets (const fastjet::ClusterSequenceAreaBase & clust_seq, 
00165                  const vector<fastjet::PseudoJet> & unsorted_jets ) {
00166 
00167   // sort jets into increasing pt
00168   vector<fastjet::PseudoJet> jets = sorted_by_pt(unsorted_jets);  
00169 
00170   // the corrected jets will go in here
00171   vector<fastjet::PseudoJet> corrected_jets(jets.size());
00172   
00173   // get median pt per unit area -- you must decide the rapidity range
00174   // in which you measure the activity. (If doing a ghost-based active
00175   // area, make sure that your ghosts go up at least to \sim range+R).
00176   double range = 5.0;
00177   double median_pt_per_area = clust_seq.median_pt_per_unit_area(range);
00178   double median_pt_per_area4vector = clust_seq.median_pt_per_unit_area_4vector(range);
00179 
00180   printf(" ijet     rap     phi        Pt    area  Pt corr  (rap corr phi corr Pt corr)ext\n");
00181   for (size_t j = 0; j < jets.size(); j++) {
00182 
00183     // get area of each jet
00184     double area     = jets[j].area();
00185 
00186     // "standard" correction. Subtract only the Pt
00187     double pt_corr  = jets[j].perp() - area*median_pt_per_area;
00188 
00189     // "extended" correction
00190     fastjet::PseudoJet sub_4vect = 
00191                        median_pt_per_area4vector*jets[j].area_4vector();
00192     if (sub_4vect.perp2() >= jets[j].perp2() || 
00193         sub_4vect.E()     >= jets[j].E()) {
00194       // if the correction is too large, set the jet to zero
00195       corrected_jets[j] =  0.0 * jets[j];
00196     } else {
00197       // otherwise do an E-scheme subtraction
00198       corrected_jets[j] = jets[j] - sub_4vect;
00199     }
00200     // NB We could also leave out the above "if": 
00201     // corrected_jets[j] = jets[j] - sub_4vect;
00202     // but the result would be different, since we would not avoid
00203     // jets with negative Pt or energy
00204     
00205     printf("%5u %7.3f %7.3f %9.3f %7.3f %9.3f %7.3f %7.3f %9.3f\n",
00206      j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, pt_corr,
00207      corrected_jets[j].rap(),corrected_jets[j].phi(), corrected_jets[j].perp());
00208   }
00209 
00210   cout << endl;
00211   cout << "median pt_over_area = " << median_pt_per_area << endl;
00212   cout << "median pt_over_area4vector = " << median_pt_per_area4vector << endl << endl;
00213 
00214 
00215 }
00216 
00217 
00218 
00219 
00220 
00221 
00222 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends