fastjet 2.4.5
fastjet_subtraction.cc
Go to the documentation of this file.
00001 
00002 //STARTHEADER
00003 // $Id: fastjet_subtraction.cc 958 2007-11-28 17:43:48Z cacciari $
00004 //
00005 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00027 //  Foundation, Inc.:
00028 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00029 //----------------------------------------------------------------------
00030 //ENDHEADER
00031 
00032 
00033 
00034 
00035 
00036 
00037 
00038 //----------------------------------------------------------------------
00039 // fastjet subtraction example program. 
00040 //
00041 // Compile it with: make fastjet_subtraction
00042 // run it with    : ./fastjet_subtraction < data/Pythia-dijet-ptmin100-lhc-pileup-1ev.dat
00043 //
00044 //----------------------------------------------------------------------
00045 #include "fastjet/PseudoJet.hh"
00046 //#include "fastjet/ClusterSequenceActiveArea.hh"
00047 #include "fastjet/ClusterSequenceArea.hh"
00048 //#include "../areas/ClusterSequencePassiveArea.hh"
00049 #include<iostream> // needed for io
00050 #include<sstream>  // needed for internal io
00051 #include<vector> 
00052 
00053 using namespace std;
00054 
00055 // A declaration of a function that pretty prints a list of jets
00056 // The subtraction is also performed inside this function
00057 void print_jets (const fastjet::ClusterSequenceAreaBase &, 
00058                  const vector<fastjet::PseudoJet> &);
00059 
00061 int main (int argc, char ** argv) {
00062   
00063   // since we use here simulated data we can split the hard event
00064   // from the full (i.e. with pileup added) one
00065   vector<fastjet::PseudoJet> hard_event, full_event;
00066   
00067   // maximum rapidity that we accept for the input particles
00068   double etamax = 6.0;
00069   
00070   // read in input particles. Keep the hard event generated by PYTHIA
00071   // separated from the full event, so as to be able to gauge the
00072   // "goodness" of the subtraction from the full event which also
00073   // includes pileup
00074   string line;
00075   int  nsub  = 0;
00076   while (getline(cin, line)) {
00077     istringstream linestream(line);
00078     // take substrings to avoid problems when there are extra "pollution"
00079     // characters (e.g. line-feed).
00080     if (line.substr(0,4) == "#END") {break;}
00081     if (line.substr(0,9) == "#SUBSTART") {
00082       // if more sub events follow, make copy of hard one here
00083       if (nsub == 1) hard_event = full_event;
00084       nsub += 1;
00085     }
00086     if (line.substr(0,1) == "#") {continue;}
00087     valarray<double> fourvec(4);
00088     linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
00089     fastjet::PseudoJet psjet(fourvec);
00090     psjet.set_user_index(0);
00091     // push event onto back of full_event vector
00092     if (abs(psjet.rap()) < etamax) {full_event.push_back(psjet);}
00093   }
00094 
00095   // if we have read in only one event, copy it across here...
00096   if (nsub == 1) hard_event = full_event;
00097 
00098   // if there was nothing in the event 
00099   if (nsub == 0) {
00100     cerr << "Error: read empty event\n";
00101     exit(-1);
00102   }
00103   
00104   
00105   // create an object that represents your choice of jet algorithm and 
00106   // the associated parameters
00107   double R = 0.7;
00108   fastjet::Strategy strategy = fastjet::Best;
00109   fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R, strategy);
00110   //fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, R, strategy);
00111   //fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R, strategy);
00112 
00113   // create an object that specifies how we to define the (active) area
00114   double ghost_etamax = 6.0;
00115   int    active_area_repeats = 1;
00116   double ghost_area    = 0.01;
00117   fastjet::ActiveAreaSpec area_spec(ghost_etamax, active_area_repeats, 
00118                                     ghost_area);
00119   fastjet::AreaDefinition area_def(area_spec);
00120   //fastjet::AreaDefinition area_def(fastjet::VoronoiAreaSpec(1.0));
00121 
00122   // run the jet clustering with the above jet definition. hard event first
00123   fastjet::ClusterSequenceArea clust_seq(hard_event, jet_def, area_def);
00124   //fastjet::ClusterSequencePassiveArea clust_seq(hard_event, jet_def);
00125 
00126 
00127   // extract the inclusive jets with pt > 5 GeV, sorted by pt
00128   double ptmin = 5.0;
00129   vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
00130 
00131   // print them out
00132   cout << "" <<endl;
00133   cout << "Hard event only"<<endl;
00134   cout << "Number of input particles: "<<hard_event.size()<<endl;
00135   cout << "Strategy used: "<<clust_seq.strategy_string()<<endl;
00136   cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n";
00137   cout << "---------------------------------------\n";
00138   print_jets(clust_seq, inclusive_jets);
00139   cout << endl;
00140 
00141 
00142   // repeat everything on the full event
00143 
00144   // run the jet clustering with the above jet definition
00145   //fastjet::ClusterSequenceActiveArea clust_seq_full(full_event, jet_def, area_spec);
00146   fastjet::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
00147   //fastjet::ClusterSequencePassiveArea clust_seq_full(full_event, jet_def,0.9);
00148 
00149   // extract the inclusive jets with pt > 20 GeV, sorted by pt
00150   ptmin = 20.0;
00151   inclusive_jets = clust_seq_full.inclusive_jets(ptmin);
00152 
00153   // print them out
00154   cout << "Full event, with pileup, and its subtraction"<<endl;
00155   cout << "Number of input particles: "<<full_event.size()<<endl;
00156   cout << "Strategy used: "<<clust_seq_full.strategy_string()<<endl;
00157   cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV (before subtraction)\n";
00158   cout << "---------------------------------------\n";
00159   print_jets(clust_seq_full, inclusive_jets);
00160   cout << endl;
00161 
00162 
00163 }
00164 
00165 
00166 //----------------------------------------------------------------------
00170 void print_jets (const fastjet::ClusterSequenceAreaBase & clust_seq, 
00171                  const vector<fastjet::PseudoJet> & unsorted_jets ) {
00172 
00173   // sort jets into increasing pt
00174   vector<fastjet::PseudoJet> jets = sorted_by_pt(unsorted_jets);  
00175 
00176   // the corrected jets will go in here
00177   vector<fastjet::PseudoJet> corrected_jets(jets.size());
00178   
00179   // get median pt per unit area -- you must decide the rapidity range
00180   // in which you measure the activity. (If doing a ghost-based active
00181   // area, make sure that your ghosts go up at least to \sim range+R).
00182   double range = 5.0;
00183   double median_pt_per_area = clust_seq.median_pt_per_unit_area(range);
00184   double median_pt_per_area4vector = clust_seq.median_pt_per_unit_area_4vector(range);
00185 
00186   printf(" ijet     rap     phi        Pt    area  Pt corr  (rap corr phi corr Pt corr)ext\n");
00187   for (size_t j = 0; j < jets.size(); j++) {
00188 
00189     // get area of each jet
00190     double area     = clust_seq.area(jets[j]);
00191 
00192     // "standard" correction. Subtract only the Pt
00193     double pt_corr  = jets[j].perp() - area*median_pt_per_area;
00194 
00195     // "extended" correction
00196     fastjet::PseudoJet sub_4vect = 
00197                        median_pt_per_area4vector*clust_seq.area_4vector(jets[j]);
00198     if (sub_4vect.perp2() >= jets[j].perp2() || 
00199         sub_4vect.E()     >= jets[j].E()) {
00200       // if the correction is too large, set the jet to zero
00201       corrected_jets[j] =  0.0 * jets[j];
00202     } else {
00203       // otherwise do an E-scheme subtraction
00204       corrected_jets[j] = jets[j] - sub_4vect;
00205     }
00206     // NB We could also leave out the above "if": 
00207     // corrected_jets[j] = jets[j] - sub_4vect;
00208     // but the result would be different, since we would not avoid
00209     // jets with negative Pt or energy
00210     
00211     printf("%5u %7.3f %7.3f %9.3f %7.3f %9.3f %7.3f %7.3f %9.3f\n",
00212      j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, pt_corr,
00213      corrected_jets[j].rap(),corrected_jets[j].phi(), corrected_jets[j].perp());
00214   }
00215 
00216   cout << endl;
00217   cout << "median pt_over_area = " << median_pt_per_area << endl;
00218   cout << "median pt_over_area4vector = " << median_pt_per_area4vector << endl << endl;
00219 
00220 
00221 }
00222 
00223 
00224 
00225 
00226 
00227 
00228 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines