Index: contrib/contribs/ConstituentSubtractor/trunk/example.cc
===================================================================
--- contrib/contribs/ConstituentSubtractor/trunk/example.cc (revision 1224)
+++ contrib/contribs/ConstituentSubtractor/trunk/example.cc (revision 1225)
@@ -1,127 +0,0 @@
-// $Id$
-//
-//----------------------------------------------------------------------
-// Example on how to use this contribution
-//
-// run it with
-// ./example < ../data/Pythia-Zp2jets-lhc-pileup-1ev.dat
-//----------------------------------------------------------------------
-//
-//----------------------------------------------------------------------
-// This file is part of FastJet contrib.
-//
-// It is free software; you can redistribute it and/or modify it under
-// the terms of the GNU General Public License as published by the
-// Free Software Foundation; either version 2 of the License, or (at
-// your option) any later version.
-//
-// It is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
-// License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this code. If not, see .
-//----------------------------------------------------------------------
-
-
-#include "fastjet/ClusterSequenceArea.hh"
-#include "fastjet/Selector.hh"
-#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
-#include "ConstituentSubtractor.hh" // In external code, this should be fastjet/contrib/ConstituentSubtractor.hh
-
-#include "functions.hh"
-
-
-//----------------------------------------------------------------------
-int main(){
-
- //----------------------------------------------------------
- // read in input particles
- vector hard_event, full_event;
- read_event(hard_event, full_event);
-
- // keep the particles up to 4 units in rapidity
- hard_event = SelectorAbsRapMax(4.0)(hard_event);
- full_event = SelectorAbsRapMax(4.0)(full_event);
-
- cout << "# read an event with " << hard_event.size() << " signal particles and " << full_event.size() - hard_event.size() << " background particles with rapidity |y|<4" << endl;
-
-
- // do the clustering with ghosts and get the jets
- //----------------------------------------------------------
- JetDefinition jet_def(antikt_algorithm, 0.7);
- double ghost_area=0.01; // the density of ghosts can be changed through this variable
- AreaDefinition area_def(active_area_explicit_ghosts,GhostedAreaSpec(4.0,1,ghost_area)); // in this step, the ghosts are added among the constituents of the jets
-
- ClusterSequenceArea clust_seq_hard(hard_event, jet_def, area_def);
- ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
-
- Selector sel_jets = SelectorNHardest(2) * SelectorAbsRapMax(3.0);
-
- vector hard_jets = sel_jets(clust_seq_hard.inclusive_jets());
- vector full_jets = sel_jets(clust_seq_full.inclusive_jets());
-
- // create what we need for the background estimation
- //----------------------------------------------------------
- JetDefinition jet_def_for_rho(kt_algorithm, 0.4);
- Selector rho_range = SelectorAbsRapMax(3.0);
- ClusterSequenceArea clust_seq_rho(full_event, jet_def, area_def);
-
- JetMedianBackgroundEstimator bge_rho(rho_range, jet_def_for_rho, area_def);
- BackgroundJetScalarPtDensity *scalarPtDensity=new BackgroundJetScalarPtDensity();
- bge_rho.set_jet_density_class(scalarPtDensity); // this changes the computation of pt of patches from vector sum to scalar sum. The scalar sum seems more reasonable.
- bge_rho.set_particles(full_event);
-
- // subtractor:
- //----------------------------------------------------------
- contrib::ConstituentSubtractor subtractor(&bge_rho);
-
- // this sets the same background estimator to be used for deltaMass density, rho_m, as for pt density, rho:
- subtractor.set_common_bge_for_rho_and_rhom(); // for massless input particles it does not make any difference (rho_m is always zero)
- cout << subtractor.description() << endl;
-
- // shape variables:
- //----------------------------------------------------------
- JetWidth width;
-
- // subtract and print the result
- //----------------------------------------------------------
- cout << setprecision(4);
- cout << "# original hard jets" << endl;
- for (unsigned int i=0; i|max_eta| are removed and not used during the subtraction (they are not returned). The same parameter should be used for the GridMedianBackgroundEstimator as it is demonstrated in this example. If JetMedianBackgroundEstimator is used, then lower parameter should be used (to avoid including particles outside this range).
- double max_eta_jet=3; // the maximal pseudorapidity for selected jets. Not used for the subtraction - just for the final output jets.
-
- // background estimator
- GridMedianBackgroundEstimator bge_rho(max_eta,0.5);
-
- contrib::ConstituentSubtractor subtractor;
- subtractor.set_distance_type(contrib::ConstituentSubtractor::deltaR); // free parameter for the type of distance between particle i and ghost k. There are two options: "deltaR" or "angle" which are defined as deltaR=sqrt((y_i-y_k)^2+(phi_i-phi_k)^2) or Euclidean angle between the momenta
- subtractor.set_max_distance(0.3); // free parameter for the maximal allowed distance between particle i and ghost k
- subtractor.set_alpha(1); // free parameter for the distance measure (the exponent of particle pt). The larger the parameter alpha, the more are favoured the lower pt particles in the subtraction process
- subtractor.set_ghost_area(0.01); // free parameter for the density of ghosts. The smaller, the better - but also the computation is slower.
- subtractor.set_max_eta(max_eta); // parameter for the maximal eta cut
- subtractor.set_background_estimator(&bge_rho); // specify the background estimator to estimate rho. You can use also specify background estimator for the mass term as an additional parameter
-
- // For the treatment of masses, choose one of the following options:
- // 1. set all masses to zero. Nothing needs to be specified.
- // 2. keep the original masses. Use this function:
- // subtractor.set_keep_original_masses();
- // 3. do mass subtraction using the same background estimator as for rho. Use this function:
- subtractor.set_common_bge_for_rho_and_rhom(); // this must be used after set_background_estimator function.
- // 4. do mass subtraction using a separate background estimator - you need to specify it as an additional parameter in function set_background_estimator
-
- // example selector for ConstituentSubtractor:
- Selector sel_CS_correction = SelectorPhiRange(0,3) * SelectorEtaMin(-1.5) * SelectorEtaMax(0);
- // subtractor.set_selector(&sel_CS_correction); // uncomment in case you want to use selector. Only ghosts fulfilling this selector will be constructed. No selection on input particles is done. The selection on input particles is the responsibility of the user. However, the maximal eta cut (specified with "set_max_eta" function) is done always for both, ghosts and input particles.
-
- // If you want to use the information about hard proxies (typically charged tracks from primary vertex and/or high pt calorimeter clusters), then use this line:
- //subtractor.set_use_nearby_hard(0.2,2); // In this example, if there is a hard proxy within deltaR distance of 0.2, then the CS distance is multiplied by factor of 2, i.e. such particle is less favoured in the subtraction procedure. If you uncomment this line, then also uncomment line 106.
-
- subtractor.initialize();
- // print info (optional)
- cout << subtractor.description() << endl;
-
-
-
-
- // in event loop:
- // read in input particles
- vector hard_event, full_event;
- vector *hard_event_charged=new vector;
- vector *background_event_charged=new vector;
- read_event(hard_event, full_event, hard_event_charged, background_event_charged);
-
- hard_event = SelectorAbsEtaMax(max_eta)(hard_event);
- full_event = SelectorAbsEtaMax(max_eta)(full_event);
-
- cout << "# read an event with " << hard_event.size() << " signal particles and " << full_event.size() - hard_event.size() << " background particles with pseudo-rapidity |eta|<4" << endl;
-
-
- // jet definition and selector for jets
- JetDefinition jet_def(antikt_algorithm, 0.7);
- Selector sel_jets = SelectorNHardest(3) * SelectorAbsEtaMax(max_eta_jet);
-
- // clustering of the hard-scatter event and uncorrected event
- ClusterSequence clust_seq_hard(hard_event, jet_def);
- ClusterSequence clust_seq_full(full_event, jet_def);
- vector hard_jets = sel_jets(clust_seq_hard.inclusive_jets());
- vector full_jets = sel_jets(clust_seq_full.inclusive_jets());
-
- bge_rho.set_particles(full_event);
-
- // the correction of the whole event with ConstituentSubtractor
- vector corrected_event=subtractor.subtract_event(full_event);
- // if you want to use the information about hard proxies, use this version:
- // vector corrected_event=subtractor.subtract_event(full_event,hard_event_charged); // here all charged hard particles are used for hard proxies. In real experiment, this corresponds to charged tracks from primary vertex. Optionally, one can add among the hard proxies also high pt calorimeter clusters after some basic pileup correction.
-
- // clustering of the corrected event
- ClusterSequence clust_seq_corr(corrected_event, jet_def);
- vector corrected_jets = sel_jets(clust_seq_corr.inclusive_jets());
-
-
- ios::fmtflags f( cout.flags() );
- cout << setprecision(4) << fixed;
- cout << endl << "Corrected particles in the whole event:" << endl;
- for (unsigned int i=0; i|max_eta| are removed and not used during the subtraction (they are not returned). The same parameter should be used for the GridMedianBackgroundEstimator as it is demonstrated in this example. If JetMedianBackgroundEstimator is used, then lower parameter should be used (to avoid including particles outside this range).
- double max_eta_jet=3; // the maximal pseudorapidity for selected jets. Not important for the subtraction.
-
- // background estimator
- GridMedianBackgroundEstimator bge_rho(max_eta,0.5);
-
- contrib::IterativeConstituentSubtractor subtractor;
- subtractor.set_distance_type(contrib::ConstituentSubtractor::deltaR); // free parameter for the type of distance between particle i and ghost k. There are two options: "deltaR" or "angle" which are defined as deltaR=sqrt((y_i-y_k)^2+(phi_i-phi_k)^2) or Euclidean angle between the momenta
- vector max_distances;
- max_distances.push_back(0.1);
- max_distances.push_back(0.2);
- vector alphas;
- alphas.push_back(0);
- alphas.push_back(0);
- subtractor.set_parameters(max_distances,alphas); // in this example, 2 CS corrections will be performed: 1. correction with max_distance of 0.1 and alpha of 0, 2. correction with max_distance of 0.15 and alpha of 0. After the first correction, the scalar sum of pt from remaining ghosts is evaluated and redistributed uniformly accross the event.
- subtractor.set_remove_remaining_proxies(true); // set to true if the ghosts (proxies) which were not used in the previous CS procedure should be removed for the next CS procedure
- subtractor.set_ghost_area(0.004); // parameter for the density of ghosts. The smaller, the better - but also the computation is slower.
- subtractor.set_max_eta(max_eta); // parameter for maximal |eta| cut. It is pspecified above.
-
- subtractor.set_background_estimator(&bge_rho); // specify the background estimator to estimate rho. You can use also specify background estimator for the mass term as an additional parameter
- // For the treatment of masses, choose one of the following options:
- // 1. set all masses to zero. Nothing needs to be specified.
- // 2. keep the original masses. Use this function:
- // subtractor.set_keep_original_masses();
- // 3. do mass subtraction using the same background estimator as for rho. Use this function:
- subtractor.set_common_bge_for_rho_and_rhom();
- // 4. do mass subtraction using a separate background estimator - you need to specify it as an additional parameter in function set_background_estimator
-
-
- // example selector for ConstituentSubtractor:
- Selector sel_CS_correction = SelectorPhiRange(0,3) * SelectorEtaMin(-1.5) * SelectorEtaMax(0);
- // subtractor.set_selector(&sel_CS_correction); // uncomment in case you want to use selector. Only ghosts fulfilling this selector will be constructed. No selection on input particles is done. The selection on input particles is the responsibility of the user. However, the maximal eta cut (specified with "set_max_eta" function) is done always for both, ghosts and input particles.
-
- // If you want to use the information about hard proxies (typically charged tracks from primary vertex and/or high pt calorimeter clusters), then use the following lines:
- /*vector nearby_hard_radii;
- nearby_hard_radii.push_back(0.2);
- nearby_hard_radii.push_back(0.2);
- vector nearby_hard_factors;
- nearby_hard_factors.push_back(2);
- nearby_hard_factors.push_back(5);
- subtractor.set_nearby_hard_parameters(nearby_hard_radii,nearby_hard_factors); // In this example, during the first iteration, if there is a hard proxy within deltaR distance of 0.2, then the CS distance is multiplied by factor of 2, i.e. such particle is less favoured in the subtraction procedure. Similarly in the second iteration. If you uncomment these lines, then also uncomment line 120.
- */
-
- subtractor.initialize(); // this is new compared to previous usages of ConstituentSubtractor! It should be used after specifying all the parameters and before event loop.
- // print info (optional)
- cout << subtractor.description() << endl;
-
-
-
-
- // in event loop
- // read in input particles
- vector hard_event, full_event;
- vector *hard_event_charged=new vector;
- vector *background_event_charged=new vector;
- read_event(hard_event, full_event, hard_event_charged, background_event_charged);
-
- hard_event = SelectorAbsEtaMax(max_eta)(hard_event);
- full_event = SelectorAbsEtaMax(max_eta)(full_event);
-
- cout << "# read an event with " << hard_event.size() << " signal particles and " << full_event.size() - hard_event.size() << " background particles with pseudo-rapidity |eta|<4" << endl;
-
-
- // jet definition and selector for jets
- JetDefinition jet_def(antikt_algorithm, 0.7);
- Selector sel_jets = SelectorNHardest(3) * SelectorAbsEtaMax(max_eta_jet);
-
- // clustering of the hard-scatter event and uncorrected event
- ClusterSequence clust_seq_hard(hard_event, jet_def);
- ClusterSequence clust_seq_full(full_event, jet_def);
- vector hard_jets = sel_jets(clust_seq_hard.inclusive_jets());
- vector full_jets = sel_jets(clust_seq_full.inclusive_jets());
-
-
- bge_rho.set_particles(full_event);
-
- // the correction of the whole event with ConstituentSubtractor
- vector corrected_event=subtractor.subtract_event(full_event); // Do not use max_eta parameter here!!! It makes nothing and this parameter will be removed from here.
- // if you want to use the information about hard proxies, use this version:
- // vector corrected_event=subtractor.subtract_event(full_event,hard_event_charged); // here all charged hard particles are used for hard proxies. In real experiment, this corresponds to charged tracks from primary vertex. Optionally, one can add among the hard proxies also high pt calorimeter clusters after some basic pileup correction.
-
-
- // clustering of the corrected event
- ClusterSequence clust_seq_corr(corrected_event, jet_def);
- vector corrected_jets = sel_jets(clust_seq_corr.inclusive_jets());
-
-
- ios::fmtflags f( cout.flags() );
- cout << setprecision(4) << fixed;
- cout << endl << "Corrected particles in the whole event:" << endl;
- for (unsigned int i=0; i