ClusterSequenceActiveArea Class Reference

Class that behaves essentially like ClusterSequence except that it also provides access to the area of a jet (which will be a random quantity. More...

#include <ClusterSequenceActiveArea.hh>

Inheritance diagram for ClusterSequenceActiveArea:
Inheritance graph
[legend]
Collaboration diagram for ClusterSequenceActiveArea:
Collaboration graph
[legend]

List of all members.

Classes

class  GhostJet
 a class for our internal storage of ghost jets More...

Public Types

enum  mean_pt_strategies {
  median = 0, non_ghost_median, pttot_over_areatot, pttot_over_areatot_cut,
  mean_ratio_cut, play, median_4vector
}
 

enum providing a variety of tentative strategies for estimating the background (e.g.

More...

Public Member Functions

 ClusterSequenceActiveArea ()
 default constructor
template<class L >
 ClusterSequenceActiveArea (const std::vector< L > &pseudojets, const JetDefinition &jet_def, const GhostedAreaSpec &ghost_spec, const bool &writeout_combinations=false)
 constructor based on JetDefinition and GhostedAreaSpec
virtual double area (const PseudoJet &jet) const
 return the area associated with the given jet; this base class returns 0.
virtual double area_error (const PseudoJet &jet) const
 return the error (uncertainty) associated with the determination of the area of this jet; this base class returns 0.
virtual PseudoJet area_4vector (const PseudoJet &jet) const
 return a PseudoJet whose 4-vector is defined by the following integral
double pt_per_unit_area (mean_pt_strategies strat=median, double range=2.0) const
 return the transverse momentum per unit area according to one of the above strategies; for some strategies (those with "cut" in their name) the parameter "range" allows one to exclude a subset of the jets for the background estimation, those that have pt/area > median(pt/area)*range.
virtual double empty_area (const RangeDefinition &range) const
 rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, in the given y-phi range, that consists of ghost jets or unclustered ghosts
virtual double n_empty_jets (const RangeDefinition &range) const
 return the true number of empty jets (replaces ClusterSequenceAreaBase::n_empty_jets(.

Protected Member Functions

void _resize_and_zero_AA ()
void _initialise_AA (const JetDefinition &jet_def, const GhostedAreaSpec &ghost_spec, const bool &writeout_combinations, bool &continue_running)
void _run_AA (const GhostedAreaSpec &ghost_spec)
void _postprocess_AA (const GhostedAreaSpec &ghost_spec)
 run the postprocessing for the active area (and derived classes)
void _initialise_and_run_AA (const JetDefinition &jet_def, const GhostedAreaSpec &ghost_spec, const bool &writeout_combinations=false)
 does the initialisation and running specific to the active areas class
void _transfer_ghost_free_history (const ClusterSequenceActiveAreaExplicitGhosts &clust_seq)
 transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts
void _transfer_areas (const std::vector< int > &unique_hist_order, const ClusterSequenceActiveAreaExplicitGhosts &)
 transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping.
bool has_dangerous_particles () const
 returns true if there are any particles whose transverse momentum if so low that there's a risk of the ghosts having modified the clustering sequence

Protected Attributes

std::valarray< double > _average_area
 child classes benefit from having these at their disposal
std::valarray< double > _average_area2
std::valarray< PseudoJet > _average_area_4vector

Private Member Functions

void _extract_tree (std::vector< int > &) const
 routine for extracting the tree in an order that will be independent of any degeneracies in the recombination sequence that don't affect the composition of the final jets
void _extract_tree_children (int pos, std::valarray< bool > &, const std::valarray< int > &, std::vector< int > &) const
 do the part of the extraction associated with pos, working through its children and their parents
void _extract_tree_parents (int pos, std::valarray< bool > &, const std::valarray< int > &, std::vector< int > &) const
 do the part of the extraction associated with the parents of pos.
void _throw_unless_jets_have_same_perp_or_E (const PseudoJet &jet, const PseudoJet &refjet, double tolerance, const ClusterSequenceActiveAreaExplicitGhosts &jets_ghosted_seq) const
 check if two jets have the same momentum to within the tolerance (and if pt's are not the same we're forgiving and look to see if the energy is the same)

Private Attributes

double _non_jet_area
double _non_jet_area2
double _non_jet_number
double _maxrap_for_area
double _safe_rap_for_area
bool _has_dangerous_particles
int _ghost_spec_repeat
 since we are playing nasty games with seeds, we should warn the user a few times
std::vector< GhostJet > _ghost_jets
std::vector< GhostJet > _unclustered_ghosts

Detailed Description

Class that behaves essentially like ClusterSequence except that it also provides access to the area of a jet (which will be a random quantity.

.. Figure out what to do about seeds later...)

Definition at line 56 of file ClusterSequenceActiveArea.hh.

Member Enumeration Documentation

enum ClusterSequenceActiveArea::mean_pt_strategies

enum providing a variety of tentative strategies for estimating the background (e.g.

non-jet) activity in a highly populated event; the one that has been most extensively tested is median.

Enumerator:
median 
non_ghost_median 
pttot_over_areatot 
pttot_over_areatot_cut 
mean_ratio_cut 
play 
median_4vector 

Definition at line 80 of file ClusterSequenceActiveArea.hh.

00080                          {median=0, non_ghost_median, pttot_over_areatot, 
00081                           pttot_over_areatot_cut, mean_ratio_cut, play,
00082                           median_4vector};

Constructor & Destructor Documentation

ClusterSequenceActiveArea::ClusterSequenceActiveArea (  )  [inline]

default constructor

Definition at line 60 of file ClusterSequenceActiveArea.hh.

00060 {}

template<class L >
ClusterSequenceActiveArea::ClusterSequenceActiveArea ( const std::vector< L > &  pseudojets,
const JetDefinition &  jet_def,
const GhostedAreaSpec &  ghost_spec,
const bool &  writeout_combinations = false 
) [inline]

constructor based on JetDefinition and GhostedAreaSpec

Definition at line 204 of file ClusterSequenceActiveArea.hh.

00207                                      {
00208 
00209   // transfer the initial jets (type L) into our own array
00210   _transfer_input_jets(pseudojets);
00211 
00212   // run the clustering for active areas
00213   _initialise_and_run_AA(jet_def, ghost_spec, writeout_combinations);
00214 
00215 }

Member Function Documentation

void ClusterSequenceActiveArea::_extract_tree ( std::vector< int > &   )  const [private]

routine for extracting the tree in an order that will be independent of any degeneracies in the recombination sequence that don't affect the composition of the final jets

void ClusterSequenceActiveArea::_extract_tree_children ( int  pos,
std::valarray< bool > &  ,
const std::valarray< int > &  ,
std::vector< int > &   
) const [private]

do the part of the extraction associated with pos, working through its children and their parents

Reimplemented from ClusterSequence.

void ClusterSequenceActiveArea::_extract_tree_parents ( int  pos,
std::valarray< bool > &  ,
const std::valarray< int > &  ,
std::vector< int > &   
) const [private]

do the part of the extraction associated with the parents of pos.

Reimplemented from ClusterSequence.

void ClusterSequenceActiveArea::_initialise_AA ( const JetDefinition &  jet_def,
const GhostedAreaSpec &  ghost_spec,
const bool &  writeout_combinations,
bool &  continue_running 
) [protected]

Definition at line 77 of file ClusterSequenceActiveArea.cc.

References ClusterSequence::_decant_options(), ClusterSequence::_fill_initial_history(), _ghost_spec_repeat, _has_dangerous_particles, ClusterSequence::_initialise_and_run(), _maxrap_for_area, _resize_and_zero_AA(), _safe_rap_for_area, GhostedAreaSpec::ghost_maxrap(), JetDefinition::R(), and GhostedAreaSpec::repeat().

Referenced by ClusterSequence1GhostPassiveArea::_initialise_and_run_1GPA(), and _initialise_and_run_AA().

00082 {
00083 
00084   // store this for future use
00085   _ghost_spec_repeat = ghost_spec.repeat();
00086 
00087   // make sure placeholders are there & zeroed
00088   _resize_and_zero_AA();
00089      
00090   // for future reference...
00091   _maxrap_for_area = ghost_spec.ghost_maxrap();
00092   _safe_rap_for_area = _maxrap_for_area - jet_def.R();
00093 
00094   // Make sure we'll have at least one repetition -- then we can
00095   // deduce the unghosted clustering sequence from one of the ghosted
00096   // sequences. If we do not have any repetitions, then get the
00097   // unghosted sequence from the plain unghosted clustering.
00098   //
00099   // NB: all decanting and filling of initial history will then
00100   // be carried out by base-class routine
00101   if (ghost_spec.repeat() <= 0) {
00102     _initialise_and_run(jet_def, writeout_combinations);
00103     continue_running = false;
00104     return;
00105   }
00106 
00107   // transfer all relevant info into internal variables
00108   _decant_options(jet_def, writeout_combinations);
00109 
00110   // set up the history entries for the initial particles (those
00111   // currently in _jets)
00112   _fill_initial_history();
00113 
00114   // by default it does not...
00115   _has_dangerous_particles = false;
00116   
00117   continue_running = true;
00118 }

void ClusterSequenceActiveArea::_initialise_and_run_AA ( const JetDefinition &  jet_def,
const GhostedAreaSpec &  ghost_spec,
const bool &  writeout_combinations = false 
) [protected]

does the initialisation and running specific to the active areas class

global routine for running active area

Definition at line 53 of file ClusterSequenceActiveArea.cc.

References _initialise_AA(), _postprocess_AA(), and _run_AA().

Referenced by ClusterSequencePassiveArea::_initialise_and_run_PA().

00056                                                     {
00057 
00058   bool continue_running;
00059   _initialise_AA(jet_def,  ghost_spec, writeout_combinations, continue_running);
00060   if (continue_running) {
00061     _run_AA(ghost_spec);
00062     _postprocess_AA(ghost_spec);
00063   }
00064 }

void ClusterSequenceActiveArea::_postprocess_AA ( const GhostedAreaSpec &  ghost_spec  )  [protected]

run the postprocessing for the active area (and derived classes)

Definition at line 152 of file ClusterSequenceActiveArea.cc.

References _average_area, _average_area2, _average_area_4vector, _non_jet_area, _non_jet_area2, _non_jet_number, and GhostedAreaSpec::repeat().

Referenced by ClusterSequence1GhostPassiveArea::_initialise_and_run_1GPA(), and _initialise_and_run_AA().

00152                                                                                    {
00153   _average_area  /= ghost_spec.repeat();
00154   _average_area2 /= ghost_spec.repeat();
00155   if (ghost_spec.repeat() > 1) {
00156     // the VC compiler complains if one puts everything on a single line.
00157     // An alternative solution would be to use -1.0 (+single line)
00158     const double tmp = ghost_spec.repeat()-1;
00159     _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/tmp);
00160   } else {
00161     _average_area2 = 0.0;
00162   }
00163 
00164   _non_jet_area  /= ghost_spec.repeat();
00165   _non_jet_area2 /= ghost_spec.repeat();
00166   _non_jet_area2  = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
00167                          ghost_spec.repeat());
00168   _non_jet_number /= ghost_spec.repeat();
00169 
00170   // following bizarre way of writing things is related to 
00171   // poverty of operations on PseudoJet objects (as well as some confusion
00172   // in one or two places)
00173   for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
00174     _average_area_4vector[i] = (1.0/ghost_spec.repeat()) * _average_area_4vector[i];
00175   }
00176   //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
00177 }

void ClusterSequenceActiveArea::_resize_and_zero_AA (  )  [protected]

Definition at line 67 of file ClusterSequenceActiveArea.cc.

References _average_area, _average_area2, _average_area_4vector, ClusterSequence::_jets, _non_jet_area, _non_jet_area2, and _non_jet_number.

Referenced by _initialise_AA(), and ClusterSequencePassiveArea::_initialise_and_run_PA().

00067                                                      {
00068   // initialize our local area information
00069   _average_area.resize(2*_jets.size());  _average_area  = 0.0;
00070   _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
00071   _average_area_4vector.resize(2*_jets.size()); 
00072   _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
00073   _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
00074 }

void ClusterSequenceActiveArea::_run_AA ( const GhostedAreaSpec &  ghost_spec  )  [protected]

Definition at line 122 of file ClusterSequenceActiveArea.cc.

References _has_dangerous_particles, ClusterSequence::_jets, _transfer_areas(), _transfer_ghost_free_history(), ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles(), ClusterSequence::jet_def(), GhostedAreaSpec::repeat(), and ClusterSequence::unique_history_order().

Referenced by _initialise_and_run_AA().

00122                                                                            {
00123   // record the input jets as they are currently
00124   vector<PseudoJet> input_jets(_jets);
00125 
00126   // code for testing the unique tree
00127   vector<int> unique_tree;
00128 
00129   // run the clustering multiple times so as to get areas of all the jets
00130   for (int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) {
00131 
00132     ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, 
00133                                                       jet_def(), ghost_spec);
00134 
00135     _has_dangerous_particles |= clust_seq.has_dangerous_particles();
00136     if (irepeat == 0) {
00137       // take the non-ghost part of the history and put into our own
00138       // history.
00139       _transfer_ghost_free_history(clust_seq);
00140       // get the "unique" order that will be used for transferring all areas. 
00141       unique_tree = unique_history_order();
00142     }
00143 
00144     // transfer areas from clust_seq into our object
00145     _transfer_areas(unique_tree, clust_seq);
00146   }
00147 }

void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E ( const PseudoJet &  jet,
const PseudoJet &  refjet,
double  tolerance,
const ClusterSequenceActiveAreaExplicitGhosts &  jets_ghosted_seq 
) const [private]

check if two jets have the same momentum to within the tolerance (and if pt's are not the same we're forgiving and look to see if the energy is the same)

Definition at line 749 of file ClusterSequenceActiveArea.cc.

References PseudoJet::E(), ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles(), PseudoJet::perp2(), PseudoJet::px(), PseudoJet::py(), and PseudoJet::pz().

00754         {
00755 
00756   if (abs(jet.perp2()-refjet.perp2()) > 
00757       tolerance*max(jet.perp2(),refjet.perp2())
00758       && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
00759     ostringstream ostr;
00760     ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas. See FAQ for possible explanations." << endl;
00761     ostr << "  Ref-Jet: "
00762          << refjet.px() << " " 
00763          << refjet.py() << " " 
00764          << refjet.pz() << " " 
00765          << refjet.E() << endl;
00766     ostr << "  New-Jet: "
00767          << jet.px() << " " 
00768          << jet.py() << " " 
00769          << jet.pz() << " " 
00770          << jet.E() << endl;
00771     if (jets_ghosted_seq.has_dangerous_particles()) {
00772       ostr << "  NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
00773     //ostr << jet.perp() << " " << refjet.perp() << " "
00774     //     << jet.perp() - refjet.perp() << endl;
00775     throw Error(ostr.str());
00776   }
00777 }

void ClusterSequenceActiveArea::_transfer_areas ( const std::vector< int > &  unique_hist_order,
const ClusterSequenceActiveAreaExplicitGhosts &   
) [protected]

transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping.

..

Referenced by ClusterSequence1GhostPassiveArea::_run_1GPA(), and _run_AA().

void ClusterSequenceActiveArea::_transfer_ghost_free_history ( const ClusterSequenceActiveAreaExplicitGhosts &  clust_seq  )  [protected]

transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts

Definition at line 476 of file ClusterSequenceActiveArea.cc.

References ClusterSequence::_do_iB_recombination_step(), ClusterSequence::_do_ij_recombination_step(), ClusterSequence::_history, ClusterSequence::_strategy, ClusterSequence::BeamJet, ClusterSequence::history_element::dij, ClusterSequence::history(), ClusterSequence::InexistentParent, ClusterSequence::Invalid, ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), ClusterSequence::history_element::parent1, ClusterSequence::history_element::parent2, and ClusterSequence::strategy_used().

Referenced by ClusterSequence1GhostPassiveArea::_run_1GPA(), and _run_AA().

00477                                                                           {
00478   
00479   const vector<history_element> & gs_history  = ghosted_seq.history();
00480   vector<int> gs2self_hist_map(gs_history.size());
00481 
00482   // first transfer info about strategy used (which isn't necessarily
00483   // always the one that got asked for...)
00484   _strategy = ghosted_seq.strategy_used();
00485 
00486   // work our way through to first non-trivial combination
00487   unsigned igs = 0;
00488   unsigned iself = 0;
00489   while (igs < gs_history.size() && gs_history[igs].parent1 == InexistentParent) {
00490     // record correspondence 
00491     if (!ghosted_seq.is_pure_ghost(igs)) {
00492       gs2self_hist_map[igs] = iself++; 
00493     } else {
00494       gs2self_hist_map[igs] = Invalid; 
00495     }
00496     igs++;
00497   };
00498 
00499   // make sure the count of non-ghost initial jets is equal to
00500   // what we already have in terms of initial jets
00501   assert(iself == _history.size());
00502 
00503   // if there was no clustering in this event (e.g. SISCone passive
00504   // area with zero input particles, or with a pt cut on stable cones
00505   // that kills all jets), then don't bother with the rest (which
00506   // would crash!)
00507   if (igs == gs_history.size()) return;
00508   
00509   // now actually transfer things
00510   do  {
00511     // if we are a pure ghost, then go on to next round
00512     if (ghosted_seq.is_pure_ghost(igs)) {
00513       gs2self_hist_map[igs] = Invalid;
00514       continue;
00515     }
00516 
00517     const history_element & gs_hist_el = gs_history[igs];
00518 
00519     bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
00520     bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);
00521 
00522     // if exactly one parent is a ghost then maintain info about the
00523     // non-ghost correspondence for this jet, and then go on to next
00524     // recombination in the ghosted sequence
00525     if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
00526       gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
00527       continue;
00528     }
00529     if (!parent1_is_ghost && parent2_is_ghost) {
00530       gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
00531       continue;
00532     }
00533 
00534     // no parents are ghosts...
00535     if (gs_hist_el.parent2 >= 0) {
00536       // recombination of two non-ghosts
00537       gs2self_hist_map[igs] = _history.size();
00538       // record the recombination in our own sequence
00539       int newjet_k; // dummy var -- not used
00540       int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
00541       int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
00542       //cerr << "recombining "<< jet_i << " and "<< jet_j << endl;
00543       _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
00544     } else {
00545       // we have a non-ghost that has become a beam-jet
00546       assert(gs_history[igs].parent2 == BeamJet);
00547       // record position
00548       gs2self_hist_map[igs] = _history.size();
00549       // record the recombination in our own sequence
00550       _do_iB_recombination_step(
00551              _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
00552              gs_hist_el.dij);
00553     }
00554   } while (++igs < gs_history.size());
00555 
00556 }

virtual double ClusterSequenceActiveArea::area ( const PseudoJet &   )  const [inline, virtual]

return the area associated with the given jet; this base class returns 0.

Reimplemented from ClusterSequenceAreaBase.

Definition at line 69 of file ClusterSequenceActiveArea.hh.

References _average_area, and PseudoJet::cluster_hist_index().

Referenced by pt_per_unit_area().

00069                                                     {
00070                              return _average_area[jet.cluster_hist_index()];};

virtual PseudoJet ClusterSequenceActiveArea::area_4vector ( const PseudoJet &   )  const [inline, virtual]

return a PseudoJet whose 4-vector is defined by the following integral

drap d PseudoJet("rap,phi,pt=one") * Theta("rap,phi inside jet boundary")

where PseudoJet("rap,phi,pt=one") is a 4-vector with the given rapidity (rap), azimuth (phi) and pt=1, while Theta("rap,phi inside jet boundary") is a function that is 1 when rap,phi define a direction inside the jet boundary and 0 otherwise.

This base class returns a null 4-vector.

Reimplemented from ClusterSequenceAreaBase.

Definition at line 74 of file ClusterSequenceActiveArea.hh.

References _average_area_4vector, and PseudoJet::cluster_hist_index().

Referenced by pt_per_unit_area().

00074                                                                {
00075                     return _average_area_4vector[jet.cluster_hist_index()];};

virtual double ClusterSequenceActiveArea::area_error ( const PseudoJet &   )  const [inline, virtual]

return the error (uncertainty) associated with the determination of the area of this jet; this base class returns 0.

Reimplemented from ClusterSequenceAreaBase.

Definition at line 71 of file ClusterSequenceActiveArea.hh.

References _average_area2, and PseudoJet::cluster_hist_index().

00071                                                           {
00072                              return _average_area2[jet.cluster_hist_index()];};

double ClusterSequenceActiveArea::empty_area ( const RangeDefinition &  range  )  const [virtual]

rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, in the given y-phi range, that consists of ghost jets or unclustered ghosts

Reimplemented from ClusterSequenceAreaBase.

Reimplemented in ClusterSequencePassiveArea.

Definition at line 445 of file ClusterSequenceActiveArea.cc.

References _ghost_jets, _ghost_spec_repeat, _unclustered_ghosts, and RangeDefinition::is_in_range().

00445                                                                                 {
00446   double empty = 0.0;
00447   // first deal with ghost jets
00448   for (unsigned  i = 0; i < _ghost_jets.size(); i++) {
00449     if (range.is_in_range(_ghost_jets[i])) {
00450       empty += _ghost_jets[i].area;
00451     }
00452   }
00453   // then deal with unclustered ghosts
00454   for (unsigned  i = 0; i < _unclustered_ghosts.size(); i++) {
00455     if (range.is_in_range(_unclustered_ghosts[i])) {
00456       empty += _unclustered_ghosts[i].area;
00457     }
00458   }
00459   empty /= _ghost_spec_repeat;
00460   return empty;
00461 }

bool ClusterSequenceActiveArea::has_dangerous_particles (  )  const [inline, protected]

returns true if there are any particles whose transverse momentum if so low that there's a risk of the ghosts having modified the clustering sequence

Definition at line 149 of file ClusterSequenceActiveArea.hh.

References _has_dangerous_particles.

00149 {return _has_dangerous_particles;}

double ClusterSequenceActiveArea::n_empty_jets ( const RangeDefinition &  range  )  const [virtual]

return the true number of empty jets (replaces ClusterSequenceAreaBase::n_empty_jets(.

..))

Reimplemented from ClusterSequenceAreaBase.

Reimplemented in ClusterSequence1GhostPassiveArea.

Definition at line 464 of file ClusterSequenceActiveArea.cc.

References _ghost_jets, _ghost_spec_repeat, and RangeDefinition::is_in_range().

00464                                                                                   {
00465   double inrange = 0;
00466   for (unsigned  i = 0; i < _ghost_jets.size(); i++) {
00467     if (range.is_in_range(_ghost_jets[i])) inrange++;
00468   }
00469   inrange /= _ghost_spec_repeat;
00470   return inrange;
00471 }

double ClusterSequenceActiveArea::pt_per_unit_area ( mean_pt_strategies  strat = median,
double  range = 2.0 
) const

return the transverse momentum per unit area according to one of the above strategies; for some strategies (those with "cut" in their name) the parameter "range" allows one to exclude a subset of the jets for the background estimation, those that have pt/area > median(pt/area)*range.

NB: This call is OBSOLETE; use media_pt_per_unit_area from the

Definition at line 276 of file ClusterSequenceActiveArea.cc.

References _non_jet_area, _non_jet_number, _safe_rap_for_area, area(), area_4vector(), ClusterSequence::inclusive_jets(), mean_ratio_cut, median, median_4vector, non_ghost_median, PseudoJet::perp(), play, pttot_over_areatot, and pttot_over_areatot_cut.

00277                                                                      {
00278   
00279   vector<PseudoJet> incl_jets = inclusive_jets();
00280   vector<double> pt_over_areas;
00281 
00282   for (unsigned i = 0; i < incl_jets.size(); i++) {
00283     if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
00284       double this_area;
00285       if ( strat == median_4vector ) {
00286           this_area = area_4vector(incl_jets[i]).perp();
00287       } else {
00288           this_area = area(incl_jets[i]);
00289       }
00290       pt_over_areas.push_back(incl_jets[i].perp()/this_area);
00291     }
00292   }
00293   
00294   // there is nothing inside our region, so answer will always be zero
00295   if (pt_over_areas.size() == 0) {return 0.0;}
00296   
00297   // get median (pt/area) [this is the "old" median definition. It considers
00298   // only the "real" jets in calculating the median, i.e. excluding the
00299   // only-ghost ones]
00300   sort(pt_over_areas.begin(), pt_over_areas.end());
00301   double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];
00302 
00303   // new median definition that takes into account non-jet area (i.e.
00304   // jets composed only of ghosts), and for fractional median position 
00305   // interpolates between the corresponding entries in the pt_over_areas array
00306   double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0;
00307   double nj_median_ratio;
00308   if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
00309     int int_nj_median = int(nj_median_pos);
00310     nj_median_ratio = 
00311       pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00312       + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00313   } else {
00314     nj_median_ratio = 0.0;
00315   }
00316 
00317 
00318   // get various forms of mean (pt/area)
00319   double pt_sum = 0.0, pt_sum_with_cut = 0.0;
00320   double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area;
00321   double ratio_sum = 0.0; 
00322   double ratio_n = _non_jet_number;
00323   for (unsigned i = 0; i < incl_jets.size(); i++) {
00324     if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
00325       double this_area;
00326       if ( strat == median_4vector ) {
00327           this_area = area_4vector(incl_jets[i]).perp();
00328       } else {
00329           this_area = area(incl_jets[i]);
00330       }
00331       pt_sum   += incl_jets[i].perp();
00332       area_sum += this_area;
00333       double ratio = incl_jets[i].perp()/this_area;
00334       if (ratio < range*nj_median_ratio) {
00335         pt_sum_with_cut   += incl_jets[i].perp();
00336         area_sum_with_cut += this_area;
00337         ratio_sum += ratio; ratio_n++;
00338       }
00339     }
00340   }
00341   
00342   if (strat == play) {
00343     double trunc_sum = 0, trunc_sumsqr = 0;
00344     vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size());
00345     for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) {
00346       double ratio = pt_over_areas[i];
00347       trunc_sum += ratio;
00348       trunc_sumsqr += ratio*ratio;
00349       means[i] = trunc_sum / (i+1);
00350       sd[i]    = sqrt(abs(means[i]*means[i]  - trunc_sumsqr/(i+1)));
00351       cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<<
00352         sd[i]/sqrt(i+1.0)<<endl;
00353     }
00354     cout << "-----------------------------------"<<endl;
00355     for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
00356       cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl;
00357     }
00358     cout << "Number of non-jets: "<<_non_jet_number<<endl;
00359     cout << "Area of non-jets: "<<_non_jet_area<<endl;
00360     cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl;
00361     cout << "NJ median position: " << nj_median_pos <<endl;
00362     cout << "NJ median value: " << nj_median_ratio <<endl;
00363     return 0.0;
00364   }
00365 
00366   switch(strat) {
00367   case median:
00368   case median_4vector:
00369     return nj_median_ratio;
00370   case non_ghost_median:
00371     return non_ghost_median_ratio; 
00372   case pttot_over_areatot:
00373     return pt_sum / area_sum;
00374   case pttot_over_areatot_cut:
00375     return pt_sum_with_cut / area_sum_with_cut;
00376   case mean_ratio_cut:
00377     return ratio_sum/ratio_n;
00378   default:
00379     return nj_median_ratio;
00380   }
00381 
00382 }

Member Data Documentation

std::valarray<double> ClusterSequenceActiveArea::_average_area [protected]

child classes benefit from having these at their disposal

Definition at line 143 of file ClusterSequenceActiveArea.hh.

Referenced by ClusterSequencePassiveArea::_initialise_and_run_PA(), _postprocess_AA(), _resize_and_zero_AA(), ClusterSequence1GhostPassiveArea::_run_1GPA(), and area().

std::valarray<double> ClusterSequenceActiveArea::_average_area2 [protected]

Definition at line 143 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), _resize_and_zero_AA(), ClusterSequence1GhostPassiveArea::_run_1GPA(), and area_error().

std::valarray<PseudoJet> ClusterSequenceActiveArea::_average_area_4vector [protected]

Definition at line 144 of file ClusterSequenceActiveArea.hh.

Referenced by ClusterSequencePassiveArea::_initialise_and_run_PA(), _postprocess_AA(), _resize_and_zero_AA(), and area_4vector().

std::vector<GhostJet> ClusterSequenceActiveArea::_ghost_jets [private]

Definition at line 196 of file ClusterSequenceActiveArea.hh.

Referenced by empty_area(), and n_empty_jets().

int ClusterSequenceActiveArea::_ghost_spec_repeat [private]

since we are playing nasty games with seeds, we should warn the user a few times

Definition at line 187 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA(), empty_area(), and n_empty_jets().

bool ClusterSequenceActiveArea::_has_dangerous_particles [private]

Definition at line 159 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA(), _run_AA(), and has_dangerous_particles().

double ClusterSequenceActiveArea::_maxrap_for_area [private]

Definition at line 156 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA().

double ClusterSequenceActiveArea::_non_jet_area [private]

Definition at line 154 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), _resize_and_zero_AA(), and pt_per_unit_area().

double ClusterSequenceActiveArea::_non_jet_area2 [private]

Definition at line 154 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), and _resize_and_zero_AA().

double ClusterSequenceActiveArea::_non_jet_number [private]

Definition at line 154 of file ClusterSequenceActiveArea.hh.

Referenced by _postprocess_AA(), _resize_and_zero_AA(), and pt_per_unit_area().

double ClusterSequenceActiveArea::_safe_rap_for_area [private]

Definition at line 157 of file ClusterSequenceActiveArea.hh.

Referenced by _initialise_AA(), and pt_per_unit_area().

std::vector<GhostJet> ClusterSequenceActiveArea::_unclustered_ghosts [private]

Definition at line 197 of file ClusterSequenceActiveArea.hh.

Referenced by empty_area().

The documentation for this class was generated from the following files:
Generated on 26 Feb 2010 for fastjet by  doxygen 1.6.1