#include <ClusterSequenceActiveArea.hh>
Inheritance diagram for fastjet::ClusterSequenceActiveArea:
Public Types | |
| enum | mean_pt_strategies { median = 0, non_ghost_median, pttot_over_areatot, pttot_over_areatot_cut, mean_ratio_cut, play } |
| enum providing a variety of tentative strategies for estimating the background (e.g. More... | |
Public Member Functions | |
| template<class L> | |
| ClusterSequenceActiveArea (const std::vector< L > &pseudojets, const JetDefinition &jet_def, const ActiveAreaSpec &area_spec, const bool &writeout_combinations=false) | |
| constructor based on JetDefinition and ActiveAreaSpec | |
| 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. | |
| void | parabolic_pt_per_unit_area (double &a, double &b, double raprange=-1.0, double exclude_above=-1.0) const |
| fits a form pt_per_unit_area(y) = a + b*y^2 in the range abs(y)<raprange (for negative raprange, it defaults to _etalim_for_area). | |
Private Member Functions | |
| void | _initialise_and_run_AA (const JetDefinition &jet_def, const ActiveAreaSpec &area_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 vector< int > &unique_hist_order, const ClusterSequenceActiveAreaExplicitGhosts &) |
| transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping. | |
| void | _extract_tree (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, valarray< bool > &, const valarray< int > &, 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, valarray< bool > &, const valarray< int > &, vector< int > &) const |
| do the part of the extraction associated with the parents of pos. | |
Private Attributes | |
| valarray< double > | _average_area |
| valarray< double > | _average_area2 |
| valarray< PseudoJet > | _average_area_4vector |
| double | _non_jet_area |
| double | _non_jet_area2 |
| double | _non_jet_number |
| double | _etamax_for_area |
| double | _etalim_for_area |
.. Figure out what to do about seeds later...)
Definition at line 49 of file ClusterSequenceActiveArea.hh.
|
|
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.
Definition at line 70 of file ClusterSequenceActiveArea.hh. 00070 {median=0, non_ghost_median, pttot_over_areatot,
00071 pttot_over_areatot_cut, mean_ratio_cut, play};
|
|
||||||||||||||||||||||||
|
constructor based on JetDefinition and ActiveAreaSpec
Definition at line 135 of file ClusterSequenceActiveArea.hh. 00138 {
00139
00140 // transfer the initial jets (type L) into our own array
00141 _transfer_input_jets(pseudojets);
00142
00143 // run the clustering for active areas
00144 _initialise_and_run_AA(jet_def, area_spec, writeout_combinations);
00145
00146 }
|
|
|
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
|
|
||||||||||||||||||||
|
do the part of the extraction associated with pos, working through its children and their parents
|
|
||||||||||||||||||||
|
do the part of the extraction associated with the parents of pos.
|
|
||||||||||||||||
|
does the initialisation and running specific to the active areas class
Definition at line 47 of file ClusterSequenceActiveArea.cc. References _average_area, _average_area2, _average_area_4vector, fastjet::ClusterSequence::_decant_options(), _etalim_for_area, _etamax_for_area, fastjet::ClusterSequence::_fill_initial_history(), fastjet::ClusterSequence::_initialise_and_run(), fastjet::ClusterSequence::_jets, _non_jet_area, _non_jet_area2, _non_jet_number, _transfer_areas(), _transfer_ghost_free_history(), fastjet::ActiveAreaSpec::ghost_etamax(), fastjet::JetDefinition::R(), fastjet::ActiveAreaSpec::repeat(), and fastjet::ClusterSequence::unique_history_order(). 00051 {
00052
00053 // initialize our local area information
00054 _average_area.resize(2*_jets.size()); _average_area = 0.0;
00055 _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
00056 _average_area_4vector.resize(2*_jets.size());
00057 _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
00058 _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
00059
00060 // for future reference...
00061 _etamax_for_area = area_spec.ghost_etamax();
00062 _etalim_for_area = _etamax_for_area - jet_def.R();
00063
00064 // Make sure we'll have at least one repetition -- then we can
00065 // deduce the unghosted clustering sequence from one of the ghosted
00066 // sequences. If we do not have any repetitions, then get the
00067 // unghosted sequence from the plain unghosted clustering.
00068 //
00069 // NB: all decanting and filling of initial history will then
00070 // be carried out by base-class routine
00071 if (area_spec.repeat() <= 0) {
00072 _initialise_and_run(jet_def, writeout_combinations);
00073 return;
00074 }
00075
00076 // transfer all relevant info into internal variables
00077 _decant_options(jet_def, writeout_combinations);
00078
00079 // set up the history entries for the initial particles (those
00080 // currently in _jets)
00081 _fill_initial_history();
00082
00083 // record the input jets as they are currently
00084 vector<PseudoJet> input_jets(_jets);
00085
00086 // code for testing the unique tree
00087 vector<int> unique_tree;
00088
00089
00090
00091
00092 // run the clustering multiple times so as to get areas of all the jets
00093 for (int irepeat = 0; irepeat < area_spec.repeat(); irepeat++) {
00094
00095 ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
00096 jet_def, area_spec);
00097
00098 if (irepeat == 0) {
00099 // take the non-ghost part of the history and put into our own
00100 // history.
00101 _transfer_ghost_free_history(clust_seq);
00102 // get the "unique" order that will be used for transferring all areas.
00103 unique_tree = unique_history_order();
00104 }
00105
00106 // transfer areas from clust_seq into our object
00107 _transfer_areas(unique_tree, clust_seq);
00108 }
00109
00110 _average_area /= area_spec.repeat();
00111 _average_area2 /= area_spec.repeat();
00112 if (area_spec.repeat() > 1) {
00113 _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/
00114 (area_spec.repeat()-1));
00115 } else {
00116 _average_area2 = 0.0;
00117 }
00118
00119 _non_jet_area /= area_spec.repeat();
00120 _non_jet_area2 /= area_spec.repeat();
00121 _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
00122 area_spec.repeat());
00123 _non_jet_number /= area_spec.repeat();
00124
00125 // following bizarre way of writing things is related to
00126 // poverty of operations on PseudoJet objects (as well as some confusion
00127 // in one or two places)
00128 for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
00129 _average_area_4vector[i] = (1.0/area_spec.repeat()) * _average_area_4vector[i];
00130 }
00131 //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
00132
00133
00134 }
|
|
||||||||||||
|
transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping. .. Definition at line 369 of file ClusterSequenceActiveArea.cc. References _average_area, _average_area2, _average_area_4vector, _etalim_for_area, fastjet::ClusterSequence::_history, fastjet::ClusterSequence::_initial_n, fastjet::ClusterSequence::_jets, _non_jet_area, _non_jet_area2, _non_jet_number, fastjet::ClusterSequenceActiveAreaExplicitGhosts::area(), area(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::area_4vector(), fastjet::ClusterSequence::BeamJet, fastjet::ClusterSequence::history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), fastjet::ClusterSequence::jets(), fastjet::ClusterSequence::n_particles(), fastjet::PseudoJet::perp(), fastjet::PseudoJet::perp2(), and fastjet::ClusterSequence::unique_history_order(). Referenced by _initialise_and_run_AA(). 00371 {
00372
00373 const vector<history_element> & gs_history = ghosted_seq.history();
00374 const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
00375 vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
00376
00377 const double tolerance = 1e-13; // to decide when two jets are the same
00378
00379 int j = -1;
00380 int hist_index = -1;
00381
00382 valarray<double> our_areas(_history.size());
00383 our_areas = 0.0;
00384
00385 valarray<PseudoJet> our_area_4vectors(_history.size());
00386 our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
00387
00388 for (unsigned i = 0; i < gs_history.size(); i++) {
00389 // only consider composite particles
00390 unsigned gs_hist_index = gs_unique_hist_order[i];
00391 if (gs_hist_index < ghosted_seq.n_particles()) continue;
00392 const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
00393 int parent1 = gs_hist.parent1;
00394 int parent2 = gs_hist.parent2;
00395
00396 if (parent2 == BeamJet) {
00397 // need to look at parent to get the actual jet
00398 const PseudoJet & jet =
00399 gs_jets[gs_history[parent1].jetp_index];
00400 double area = ghosted_seq.area(jet);
00401 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00402
00403 if (ghosted_seq.is_pure_ghost(parent1)) {
00404 if (abs(jet.rap()) < _etalim_for_area) {
00405 _non_jet_area += area;
00406 _non_jet_area2 += area*area;
00407 _non_jet_number += 1;
00408 }
00409 } else {
00410
00411 // get next "combined-particle" index in our own history
00412 // making sure we don't go beyond it's bounds (if we do
00413 // then we're in big trouble anyway...)
00414 while (++j < static_cast<int>(_history.size())) {
00415 hist_index = unique_hist_order[j];
00416 if (hist_index >= _initial_n) break;}
00417
00418 // sanity check
00419 const PseudoJet & refjet =
00420 _jets[_history[_history[hist_index].parent1].jetp_index];
00421 //if (jet.perp2() != refjet.perp2()) {
00422 if (abs(jet.perp2()-refjet.perp2()) >
00423 tolerance*max(jet.perp2(),refjet.perp2())) {
00424 cerr << jet.perp() << " " << refjet.perp() << " "<< jet.perp() - refjet.perp() << endl;
00425 throw Error("Could not match clustering sequence for an inclusive jet when reconstructing areas"); }
00426
00427 // set the area at this clustering stage
00428 our_areas[hist_index] = area;
00429 our_area_4vectors[hist_index] = ext_area;
00430
00431 // update the parent as well -- that way its area is the area
00432 // immediately before clustering (i.e. resolve an ambiguity in
00433 // the Cambridge case and ensure in the kt case that the original
00434 // particles get a correct area)
00435 our_areas[_history[hist_index].parent1] = area;
00436 our_area_4vectors[_history[hist_index].parent1] = ext_area;
00437
00438 }
00439 }
00440 else if (!ghosted_seq.is_pure_ghost(parent1) &&
00441 !ghosted_seq.is_pure_ghost(parent2)) {
00442
00443 // get next "combined-particle" index in our own history
00444 while (++j < static_cast<int>(_history.size())) {
00445 hist_index = unique_hist_order[j];
00446 if (hist_index >= _initial_n) break;}
00447
00448 const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
00449 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
00450
00451 // run sanity check
00452 if (abs(jet.perp2()-refjet.perp2()) >
00453 tolerance*max(jet.perp2(),refjet.perp2())) {
00454 cerr << jet.perp() << " " << refjet.perp() << " "<< jet.perp() - refjet.perp() << endl;
00455 throw Error("Could not match clustering sequence for an exclusive jet when reconstructing areas"); }
00456
00457 // update area and our local index (maybe redundant since later
00458 // the descendants will reupdate it?)
00459 double area = ghosted_seq.area(jet);
00460 our_areas[hist_index] += area;
00461
00462 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00463 our_area_4vectors[hist_index] = our_area_4vectors[hist_index] + ext_area;
00464
00465 // now update areas of parents (so that they becomes areas
00466 // immediately before clustering occurred). This is of use
00467 // because it allows us to set the areas of the original hard
00468 // particles in the kt algorithm; for the Cambridge case it
00469 // means a jet's area will be the area just before it clusters
00470 // with another hard jet.
00471 const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
00472 int our_parent1 = _history[hist_index].parent1;
00473 our_areas[our_parent1] = ghosted_seq.area(jet1);
00474 our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
00475
00476 const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
00477 int our_parent2 = _history[hist_index].parent2;
00478 our_areas[our_parent2] = ghosted_seq.area(jet2);
00479 our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
00480 }
00481
00482 }
00483
00484 _average_area += our_areas;
00485 _average_area2 += our_areas*our_areas;
00486
00487 // poverty of operators forces us to write things this way...
00488 //_average_area_4vector = _average_area_4vector + our_area_4vectors;
00489 _average_area_4vector += our_area_4vectors;
00490
00491 }
|
|
|
transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts
Definition at line 291 of file ClusterSequenceActiveArea.cc. References fastjet::ClusterSequence::_do_iB_recombination_step(), fastjet::ClusterSequence::_do_ij_recombination_step(), fastjet::ClusterSequence::_history, fastjet::ClusterSequence::_strategy, fastjet::ClusterSequence::BeamJet, fastjet::ClusterSequence::history(), fastjet::ClusterSequence::InexistentParent, fastjet::ClusterSequence::Invalid, fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), and fastjet::ClusterSequence::strategy_used(). Referenced by _initialise_and_run_AA(). 00292 {
00293
00294 const vector<history_element> & gs_history = ghosted_seq.history();
00295 vector<int> gs2self_hist_map(gs_history.size());
00296
00297 // work our way through to first non-trivial combination
00298 unsigned igs = 0;
00299 unsigned iself = 0;
00300 while (gs_history[igs].parent1 == InexistentParent) {
00301 // record correspondence
00302 if (!ghosted_seq.is_pure_ghost(igs)) {
00303 gs2self_hist_map[igs] = iself++;
00304 } else {
00305 gs2self_hist_map[igs] = Invalid;
00306 }
00307 igs++;
00308 };
00309
00310 // make sure the count of non-ghost initial jets is equal to
00311 // what we already have in terms of initial jets
00312 assert(iself == _history.size());
00313
00314 // now actually transfer things
00315 do {
00316 // if we are a pure ghost, then go on to next round
00317 if (ghosted_seq.is_pure_ghost(igs)) {
00318 gs2self_hist_map[igs] = Invalid;
00319 continue;
00320 }
00321
00322 const history_element & gs_hist_el = gs_history[igs];
00323
00324 bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
00325 bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);
00326
00327 // if exactly one parent is a ghost then maintain info about the
00328 // non-ghost correspondence for this jet, and then go on to next
00329 // recombination in the ghosted sequence
00330 if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
00331 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
00332 continue;
00333 }
00334 if (!parent1_is_ghost && parent2_is_ghost) {
00335 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
00336 continue;
00337 }
00338
00339 // no parents are ghosts...
00340 if (gs_hist_el.parent2 >= 0) {
00341 // recombination of two non-ghosts
00342 gs2self_hist_map[igs] = _history.size();
00343 // record the recombination in our own sequence
00344 int newjet_k; // dummy var -- not used
00345 //cerr << igs << " " << gs_hist_el.parent1 << " " << gs_hist_el.parent2 << endl;
00346 //cerr << gs2self_hist_map[gs_hist_el.parent1] << " " << gs2self_hist_map[gs_hist_el.parent2] << endl;
00347 int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
00348 int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
00349 //cerr << "recombining "<< jet_i << " and "<< jet_j << endl;
00350 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
00351 } else {
00352 // we have a non-ghost that has become a beam-jet
00353 assert(gs_history[igs].parent2 == BeamJet);
00354 // record position
00355 gs2self_hist_map[igs] = _history.size();
00356 // record the recombination in our own sequence
00357 _do_iB_recombination_step(
00358 _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
00359 gs_hist_el.dij);
00360 }
00361 } while (++igs < gs_history.size());
00362
00363 // finally transfer info about strategy used (which isn't necessarily
00364 // always the one that got asked for...)
00365 _strategy = ghosted_seq.strategy_used();
00366 }
|
|
|
return the area associated with the given jet; this base class returns 0.
Reimplemented from fastjet::ClusterSequenceWithArea. Definition at line 59 of file ClusterSequenceActiveArea.hh. References fastjet::PseudoJet::cluster_hist_index(). Referenced by _transfer_areas(), parabolic_pt_per_unit_area(), print_jets(), and pt_per_unit_area(). 00059 {
00060 return _average_area[jet.cluster_hist_index()];};
|
|
|
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 fastjet::ClusterSequenceWithArea. Definition at line 64 of file ClusterSequenceActiveArea.hh. References fastjet::PseudoJet::cluster_hist_index(). Referenced by print_jets(). 00064 {
00065 return _average_area_4vector[jet.cluster_hist_index()];};
|
|
|
return the error (uncertainty) associated with the determination of the area of this jet; this base class returns 0.
Reimplemented from fastjet::ClusterSequenceWithArea. Definition at line 61 of file ClusterSequenceActiveArea.hh. References fastjet::PseudoJet::cluster_hist_index(). 00061 {
00062 return _average_area2[jet.cluster_hist_index()];};
|
|
||||||||||||||||||||
|
fits a form pt_per_unit_area(y) = a + b*y^2 in the range abs(y)<raprange (for negative raprange, it defaults to _etalim_for_area).
Definition at line 239 of file ClusterSequenceActiveArea.cc. References _etalim_for_area, area(), and fastjet::ClusterSequence::inclusive_jets(). 00240 {
00241
00242 double this_raprange;
00243 if (raprange <= 0) {this_raprange = _etalim_for_area;}
00244 else {this_raprange = raprange;}
00245
00246 int n=0;
00247 int n_excluded = 0;
00248 double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
00249
00250 vector<PseudoJet> incl_jets = inclusive_jets();
00251
00252 for (unsigned i = 0; i < incl_jets.size(); i++) {
00253 if (abs(incl_jets[i].rap()) < this_raprange) {
00254 double this_area = area(incl_jets[i]);
00255 double f = incl_jets[i].perp()/this_area;
00256 if (exclude_above <= 0.0 || f < exclude_above) {
00257 double x = incl_jets[i].rap(); double x2 = x*x;
00258 mean_f += f;
00259 mean_x2 += x2;
00260 mean_x4 += x2*x2;
00261 mean_fx2 += f*x2;
00262 n++;
00263 } else {
00264 n_excluded++;
00265 }
00266 }
00267 }
00268
00269 if (n <= 1) {
00270 // meaningful results require at least two jets inside the
00271 // area -- mind you if there are empty jets we should be in
00272 // any case doing something special...
00273 a = 0.0;
00274 b = 0.0;
00275 } else {
00276 mean_f /= n;
00277 mean_x2 /= n;
00278 mean_x4 /= n;
00279 mean_fx2 /= n;
00280
00281 b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
00282 a = mean_f - b*mean_x2;
00283 }
00284 //cerr << "n_excluded = "<< n_excluded << endl;
00285 }
|
|
||||||||||||
|
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.
Definition at line 138 of file ClusterSequenceActiveArea.cc. References _etalim_for_area, _non_jet_area, _non_jet_number, area(), fastjet::ClusterSequence::inclusive_jets(), mean_ratio_cut, median, non_ghost_median, play, pttot_over_areatot, and pttot_over_areatot_cut. Referenced by print_jets(). 00139 {
00140
00141 vector<PseudoJet> incl_jets = inclusive_jets();
00142 vector<double> pt_over_areas;
00143
00144 for (unsigned i = 0; i < incl_jets.size(); i++) {
00145 if (abs(incl_jets[i].rap()) < _etalim_for_area) {
00146 double this_area = area(incl_jets[i]);
00147 pt_over_areas.push_back(incl_jets[i].perp()/this_area);
00148 }
00149 }
00150
00151 // there is nothing inside our region, so answer will always be zero
00152 if (pt_over_areas.size() == 0) {return 0.0;}
00153
00154 // get median (pt/area) [this is the "old" median definition. It considers
00155 // only the "real" jets in calculating the median, i.e. excluding the
00156 // only-ghost ones]
00157 sort(pt_over_areas.begin(), pt_over_areas.end());
00158 double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];
00159
00160 // new median definition that takes into account non-jet area (i.e.
00161 // jets composed only of ghosts), and for fractional median position
00162 // interpolates between the corresponding entries in the pt_over_areas array
00163 double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0;
00164 double nj_median_ratio;
00165 if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
00166 int int_nj_median = int(nj_median_pos);
00167 nj_median_ratio =
00168 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00169 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00170 } else {
00171 nj_median_ratio = 0.0;
00172 }
00173
00174
00175 // get various forms of mean (pt/area)
00176 double pt_sum = 0.0, pt_sum_with_cut = 0.0;
00177 double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area;
00178 double ratio_sum = 0.0;
00179 double ratio_n = _non_jet_number;
00180 for (unsigned i = 0; i < incl_jets.size(); i++) {
00181 if (abs(incl_jets[i].rap()) < _etalim_for_area) {
00182 double this_area = area(incl_jets[i]);
00183 pt_sum += incl_jets[i].perp();
00184 area_sum += this_area;
00185 double ratio = incl_jets[i].perp()/this_area;
00186 if (ratio < range*nj_median_ratio) {
00187 pt_sum_with_cut += incl_jets[i].perp();
00188 area_sum_with_cut += this_area;
00189 ratio_sum += ratio; ratio_n++;
00190 }
00191 }
00192 }
00193
00194 if (strat == play) {
00195 double trunc_sum = 0, trunc_sumsqr = 0;
00196 vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size());
00197 for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) {
00198 double ratio = pt_over_areas[i];
00199 trunc_sum += ratio;
00200 trunc_sumsqr += ratio*ratio;
00201 means[i] = trunc_sum / (i+1);
00202 sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1)));
00203 cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<<
00204 sd[i]/sqrt(i+1.0)<<endl;
00205 }
00206 cout << "-----------------------------------"<<endl;
00207 for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
00208 cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl;
00209 }
00210 cout << "Number of non-jets: "<<_non_jet_number<<endl;
00211 cout << "Area of non-jets: "<<_non_jet_area<<endl;
00212 cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl;
00213 cout << "NJ median position: " << nj_median_pos <<endl;
00214 cout << "NJ median value: " << nj_median_ratio <<endl;
00215 return 0.0;
00216 }
00217
00218 switch(strat) {
00219 case median:
00220 return nj_median_ratio;
00221 case non_ghost_median:
00222 return non_ghost_median_ratio;
00223 case pttot_over_areatot:
00224 return pt_sum / area_sum;
00225 case pttot_over_areatot_cut:
00226 return pt_sum_with_cut / area_sum_with_cut;
00227 case mean_ratio_cut:
00228 return ratio_sum/ratio_n;
00229 default:
00230 return nj_median_ratio;
00231 }
00232
00233 }
|
|
|
Definition at line 98 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas(). |
|
|
Definition at line 98 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas(). |
|
|
Definition at line 99 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas(). |
|
|
Definition at line 103 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), _transfer_areas(), parabolic_pt_per_unit_area(), and pt_per_unit_area(). |
|
|
Definition at line 102 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(). |
|
|
Definition at line 100 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), _transfer_areas(), and pt_per_unit_area(). |
|
|
Definition at line 100 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), and _transfer_areas(). |
|
|
Definition at line 100 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_and_run_AA(), _transfer_areas(), and pt_per_unit_area(). |
1.4.2