FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: ClusterSequenceActiveArea.cc 2687 2011-11-14 11:17:51Z soyez $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 #include "fastjet/PseudoJet.hh" 00030 #include "fastjet/ClusterSequence.hh" 00031 #include "fastjet/ClusterSequenceActiveArea.hh" 00032 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh" 00033 #include<iostream> 00034 #include<vector> 00035 #include<sstream> 00036 #include<algorithm> 00037 #include<cmath> 00038 #include<valarray> 00039 00040 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00041 00042 00043 using namespace std; 00044 00045 00046 //int ClusterSequenceActiveArea::_n_seed_warnings = 0; 00047 //const int _max_seed_warnings = 10; 00048 00049 //---------------------------------------------------------------------- 00050 /// global routine for running active area 00051 void ClusterSequenceActiveArea::_initialise_and_run_AA ( 00052 const JetDefinition & jet_def_in, 00053 const GhostedAreaSpec & ghost_spec, 00054 const bool & writeout_combinations) { 00055 00056 bool continue_running; 00057 _initialise_AA(jet_def_in, ghost_spec, writeout_combinations, continue_running); 00058 if (continue_running) { 00059 _run_AA(ghost_spec); 00060 _postprocess_AA(ghost_spec); 00061 } 00062 } 00063 00064 //---------------------------------------------------------------------- 00065 void ClusterSequenceActiveArea::_resize_and_zero_AA () { 00066 // initialize our local area information 00067 _average_area.resize(2*_jets.size()); _average_area = 0.0; 00068 _average_area2.resize(2*_jets.size()); _average_area2 = 0.0; 00069 _average_area_4vector.resize(2*_jets.size()); 00070 _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0); 00071 _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0; 00072 } 00073 00074 //---------------------------------a------------------------------------- 00075 void ClusterSequenceActiveArea::_initialise_AA ( 00076 const JetDefinition & jet_def_in, 00077 const GhostedAreaSpec & ghost_spec, 00078 const bool & writeout_combinations, 00079 bool & continue_running) 00080 { 00081 00082 // store this for future use 00083 _ghost_spec_repeat = ghost_spec.repeat(); 00084 00085 // make sure placeholders are there & zeroed 00086 _resize_and_zero_AA(); 00087 00088 // for future reference... 00089 _maxrap_for_area = ghost_spec.ghost_maxrap(); 00090 _safe_rap_for_area = _maxrap_for_area - jet_def_in.R(); 00091 00092 // Make sure we'll have at least one repetition -- then we can 00093 // deduce the unghosted clustering sequence from one of the ghosted 00094 // sequences. If we do not have any repetitions, then get the 00095 // unghosted sequence from the plain unghosted clustering. 00096 // 00097 // NB: all decanting and filling of initial history will then 00098 // be carried out by base-class routine 00099 if (ghost_spec.repeat() <= 0) { 00100 _initialise_and_run(jet_def_in, writeout_combinations); 00101 continue_running = false; 00102 return; 00103 } 00104 00105 // transfer all relevant info into internal variables 00106 _decant_options(jet_def_in, writeout_combinations); 00107 00108 // set up the history entries for the initial particles (those 00109 // currently in _jets) 00110 _fill_initial_history(); 00111 00112 // by default it does not... 00113 _has_dangerous_particles = false; 00114 00115 continue_running = true; 00116 } 00117 00118 00119 //---------------------------------------------------------------------- 00120 void ClusterSequenceActiveArea::_run_AA (const GhostedAreaSpec & ghost_spec) { 00121 // record the input jets as they are currently 00122 vector<PseudoJet> input_jets(_jets); 00123 00124 // code for testing the unique tree 00125 vector<int> unique_tree; 00126 00127 // run the clustering multiple times so as to get areas of all the jets 00128 for (int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) { 00129 00130 ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, 00131 jet_def(), ghost_spec); 00132 00133 _has_dangerous_particles |= clust_seq.has_dangerous_particles(); 00134 if (irepeat == 0) { 00135 // take the non-ghost part of the history and put into our own 00136 // history. 00137 _transfer_ghost_free_history(clust_seq); 00138 // get the "unique" order that will be used for transferring all areas. 00139 unique_tree = unique_history_order(); 00140 } 00141 00142 // transfer areas from clust_seq into our object 00143 _transfer_areas(unique_tree, clust_seq); 00144 } 00145 } 00146 00147 00148 //---------------------------------------------------------------------- 00149 /// run the postprocessing for the active area (and derived classes) 00150 void ClusterSequenceActiveArea::_postprocess_AA (const GhostedAreaSpec & ghost_spec) { 00151 _average_area /= ghost_spec.repeat(); 00152 _average_area2 /= ghost_spec.repeat(); 00153 if (ghost_spec.repeat() > 1) { 00154 // the VC compiler complains if one puts everything on a single line. 00155 // An alternative solution would be to use -1.0 (+single line) 00156 const double tmp = ghost_spec.repeat()-1; 00157 _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/tmp); 00158 } else { 00159 _average_area2 = 0.0; 00160 } 00161 00162 _non_jet_area /= ghost_spec.repeat(); 00163 _non_jet_area2 /= ghost_spec.repeat(); 00164 _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/ 00165 ghost_spec.repeat()); 00166 _non_jet_number /= ghost_spec.repeat(); 00167 00168 // following bizarre way of writing things is related to 00169 // poverty of operations on PseudoJet objects (as well as some confusion 00170 // in one or two places) 00171 for (unsigned i = 0; i < _average_area_4vector.size(); i++) { 00172 _average_area_4vector[i] = (1.0/ghost_spec.repeat()) * _average_area_4vector[i]; 00173 } 00174 //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl; 00175 } 00176 00177 00178 // //---------------------------------------------------------------------- 00179 // void ClusterSequenceActiveArea::_initialise_and_run_AA ( 00180 // const JetDefinition & jet_def, 00181 // const GhostedAreaSpec & ghost_spec, 00182 // const bool & writeout_combinations) 00183 // { 00184 // 00185 // // store this for future use 00186 // _ghost_spec_repeat = ghost_spec.repeat(); 00187 // 00188 // // initialize our local area information 00189 // _average_area.resize(2*_jets.size()); _average_area = 0.0; 00190 // _average_area2.resize(2*_jets.size()); _average_area2 = 0.0; 00191 // _average_area_4vector.resize(2*_jets.size()); 00192 // _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0); 00193 // _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0; 00194 // 00195 // // for future reference... 00196 // _maxrap_for_area = ghost_spec.ghost_maxrap(); 00197 // _safe_rap_for_area = _maxrap_for_area - jet_def.R(); 00198 // 00199 // // Make sure we'll have at least one repetition -- then we can 00200 // // deduce the unghosted clustering sequence from one of the ghosted 00201 // // sequences. If we do not have any repetitions, then get the 00202 // // unghosted sequence from the plain unghosted clustering. 00203 // // 00204 // // NB: all decanting and filling of initial history will then 00205 // // be carried out by base-class routine 00206 // if (ghost_spec.repeat() <= 0) { 00207 // _initialise_and_run(jet_def, writeout_combinations); 00208 // return; 00209 // } 00210 // 00211 // // transfer all relevant info into internal variables 00212 // _decant_options(jet_def, writeout_combinations); 00213 // 00214 // // set up the history entries for the initial particles (those 00215 // // currently in _jets) 00216 // _fill_initial_history(); 00217 // 00218 // // record the input jets as they are currently 00219 // vector<PseudoJet> input_jets(_jets); 00220 // 00221 // // code for testing the unique tree 00222 // vector<int> unique_tree; 00223 // 00224 // 00225 // 00226 // 00227 // // run the clustering multiple times so as to get areas of all the jets 00228 // for (int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) { 00229 // 00230 // ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, 00231 // jet_def, ghost_spec); 00232 // 00233 // if (irepeat == 0) { 00234 // // take the non-ghost part of the history and put into our own 00235 // // history. 00236 // _transfer_ghost_free_history(clust_seq); 00237 // // get the "unique" order that will be used for transferring all areas. 00238 // unique_tree = unique_history_order(); 00239 // } 00240 // 00241 // // transfer areas from clust_seq into our object 00242 // _transfer_areas(unique_tree, clust_seq); 00243 // } 00244 // 00245 // _average_area /= ghost_spec.repeat(); 00246 // _average_area2 /= ghost_spec.repeat(); 00247 // if (ghost_spec.repeat() > 1) { 00248 // _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/ 00249 // (ghost_spec.repeat()-1)); 00250 // } else { 00251 // _average_area2 = 0.0; 00252 // } 00253 // 00254 // _non_jet_area /= ghost_spec.repeat(); 00255 // _non_jet_area2 /= ghost_spec.repeat(); 00256 // _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/ 00257 // ghost_spec.repeat()); 00258 // _non_jet_number /= ghost_spec.repeat(); 00259 // 00260 // // following bizarre way of writing things is related to 00261 // // poverty of operations on PseudoJet objects (as well as some confusion 00262 // // in one or two places) 00263 // for (unsigned i = 0; i < _average_area_4vector.size(); i++) { 00264 // _average_area_4vector[i] = (1.0/ghost_spec.repeat()) * _average_area_4vector[i]; 00265 // } 00266 // //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl; 00267 // 00268 // 00269 // } 00270 // 00271 00272 00273 //---------------------------------------------------------------------- 00274 double ClusterSequenceActiveArea::pt_per_unit_area( 00275 mean_pt_strategies strat, double range) const { 00276 00277 vector<PseudoJet> incl_jets = inclusive_jets(); 00278 vector<double> pt_over_areas; 00279 00280 for (unsigned i = 0; i < incl_jets.size(); i++) { 00281 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { 00282 double this_area; 00283 if ( strat == median_4vector ) { 00284 this_area = area_4vector(incl_jets[i]).perp(); 00285 } else { 00286 this_area = area(incl_jets[i]); 00287 } 00288 pt_over_areas.push_back(incl_jets[i].perp()/this_area); 00289 } 00290 } 00291 00292 // there is nothing inside our region, so answer will always be zero 00293 if (pt_over_areas.size() == 0) {return 0.0;} 00294 00295 // get median (pt/area) [this is the "old" median definition. It considers 00296 // only the "real" jets in calculating the median, i.e. excluding the 00297 // only-ghost ones] 00298 sort(pt_over_areas.begin(), pt_over_areas.end()); 00299 double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2]; 00300 00301 // new median definition that takes into account non-jet area (i.e. 00302 // jets composed only of ghosts), and for fractional median position 00303 // interpolates between the corresponding entries in the pt_over_areas array 00304 double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0; 00305 double nj_median_ratio; 00306 if (nj_median_pos >= 0 && pt_over_areas.size() > 1) { 00307 int int_nj_median = int(nj_median_pos); 00308 nj_median_ratio = 00309 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos) 00310 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median); 00311 } else { 00312 nj_median_ratio = 0.0; 00313 } 00314 00315 00316 // get various forms of mean (pt/area) 00317 double pt_sum = 0.0, pt_sum_with_cut = 0.0; 00318 double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area; 00319 double ratio_sum = 0.0; 00320 double ratio_n = _non_jet_number; 00321 for (unsigned i = 0; i < incl_jets.size(); i++) { 00322 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { 00323 double this_area; 00324 if ( strat == median_4vector ) { 00325 this_area = area_4vector(incl_jets[i]).perp(); 00326 } else { 00327 this_area = area(incl_jets[i]); 00328 } 00329 pt_sum += incl_jets[i].perp(); 00330 area_sum += this_area; 00331 double ratio = incl_jets[i].perp()/this_area; 00332 if (ratio < range*nj_median_ratio) { 00333 pt_sum_with_cut += incl_jets[i].perp(); 00334 area_sum_with_cut += this_area; 00335 ratio_sum += ratio; ratio_n++; 00336 } 00337 } 00338 } 00339 00340 if (strat == play) { 00341 double trunc_sum = 0, trunc_sumsqr = 0; 00342 vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size()); 00343 for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) { 00344 double ratio = pt_over_areas[i]; 00345 trunc_sum += ratio; 00346 trunc_sumsqr += ratio*ratio; 00347 means[i] = trunc_sum / (i+1); 00348 sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1))); 00349 cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<< 00350 sd[i]/sqrt(i+1.0)<<endl; 00351 } 00352 cout << "-----------------------------------"<<endl; 00353 for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) { 00354 cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl; 00355 } 00356 cout << "Number of non-jets: "<<_non_jet_number<<endl; 00357 cout << "Area of non-jets: "<<_non_jet_area<<endl; 00358 cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl; 00359 cout << "NJ median position: " << nj_median_pos <<endl; 00360 cout << "NJ median value: " << nj_median_ratio <<endl; 00361 return 0.0; 00362 } 00363 00364 switch(strat) { 00365 case median: 00366 case median_4vector: 00367 return nj_median_ratio; 00368 case non_ghost_median: 00369 return non_ghost_median_ratio; 00370 case pttot_over_areatot: 00371 return pt_sum / area_sum; 00372 case pttot_over_areatot_cut: 00373 return pt_sum_with_cut / area_sum_with_cut; 00374 case mean_ratio_cut: 00375 return ratio_sum/ratio_n; 00376 default: 00377 return nj_median_ratio; 00378 } 00379 00380 } 00381 00382 00383 // The following functionality is now provided by the base class 00384 // //---------------------------------------------------------------------- 00385 // // fit a parabola to pt/area as a function of rapidity, using the 00386 // // formulae of CCN28-36 (which actually fits f = a+b*x^2) 00387 // void ClusterSequenceActiveArea::parabolic_pt_per_unit_area( 00388 // double & a, double & b, double raprange, double exclude_above, 00389 // bool use_area_4vector) const { 00390 // 00391 // double this_raprange; 00392 // if (raprange <= 0) {this_raprange = _safe_rap_for_area;} 00393 // else {this_raprange = raprange;} 00394 // 00395 // int n=0; 00396 // int n_excluded = 0; 00397 // double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0; 00398 // 00399 // vector<PseudoJet> incl_jets = inclusive_jets(); 00400 // 00401 // for (unsigned i = 0; i < incl_jets.size(); i++) { 00402 // if (abs(incl_jets[i].rap()) < this_raprange) { 00403 // double this_area; 00404 // if ( use_area_4vector ) { 00405 // this_area = area_4vector(incl_jets[i]).perp(); 00406 // } else { 00407 // this_area = area(incl_jets[i]); 00408 // } 00409 // double f = incl_jets[i].perp()/this_area; 00410 // if (exclude_above <= 0.0 || f < exclude_above) { 00411 // double x = incl_jets[i].rap(); double x2 = x*x; 00412 // mean_f += f; 00413 // mean_x2 += x2; 00414 // mean_x4 += x2*x2; 00415 // mean_fx2 += f*x2; 00416 // n++; 00417 // } else { 00418 // n_excluded++; 00419 // } 00420 // } 00421 // } 00422 // 00423 // if (n <= 1) { 00424 // // meaningful results require at least two jets inside the 00425 // // area -- mind you if there are empty jets we should be in 00426 // // any case doing something special... 00427 // a = 0.0; 00428 // b = 0.0; 00429 // } else { 00430 // mean_f /= n; 00431 // mean_x2 /= n; 00432 // mean_x4 /= n; 00433 // mean_fx2 /= n; 00434 // 00435 // b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4); 00436 // a = mean_f - b*mean_x2; 00437 // } 00438 // //cerr << "n_excluded = "<< n_excluded << endl; 00439 // } 00440 00441 00442 //---------------------------------------------------------------------- 00443 double ClusterSequenceActiveArea::empty_area(const Selector & selector) const { 00444 // make sure that the selector applies jet by jet 00445 if (! selector.applies_jet_by_jet()){ 00446 throw Error("ClusterSequenceActiveArea: empty area can only be computed from selectors applying jet by jet"); 00447 } 00448 00449 double empty = 0.0; 00450 // first deal with ghost jets 00451 for (unsigned i = 0; i < _ghost_jets.size(); i++) { 00452 if (selector.pass(_ghost_jets[i])) { 00453 empty += _ghost_jets[i].area; 00454 } 00455 } 00456 // then deal with unclustered ghosts 00457 for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) { 00458 if (selector.pass(_unclustered_ghosts[i])) { 00459 empty += _unclustered_ghosts[i].area; 00460 } 00461 } 00462 empty /= _ghost_spec_repeat; 00463 return empty; 00464 } 00465 00466 //---------------------------------------------------------------------- 00467 double ClusterSequenceActiveArea::n_empty_jets(const Selector & selector) const { 00468 _check_selector_good_for_median(selector); 00469 00470 double inrange = 0; 00471 for (unsigned i = 0; i < _ghost_jets.size(); i++) { 00472 if (selector.pass(_ghost_jets[i])) inrange++; 00473 } 00474 inrange /= _ghost_spec_repeat; 00475 return inrange; 00476 } 00477 00478 //---------------------------------------------------------------------- 00479 /// transfer the history (and jet-momenta) from clust_seq to our 00480 /// own internal structure while removing ghosts 00481 void ClusterSequenceActiveArea::_transfer_ghost_free_history( 00482 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) { 00483 00484 const vector<history_element> & gs_history = ghosted_seq.history(); 00485 vector<int> gs2self_hist_map(gs_history.size()); 00486 00487 // first transfer info about strategy used (which isn't necessarily 00488 // always the one that got asked for...) 00489 _strategy = ghosted_seq.strategy_used(); 00490 00491 // work our way through to first non-trivial combination 00492 unsigned igs = 0; 00493 unsigned iself = 0; 00494 while (igs < gs_history.size() && gs_history[igs].parent1 == InexistentParent) { 00495 // record correspondence 00496 if (!ghosted_seq.is_pure_ghost(igs)) { 00497 gs2self_hist_map[igs] = iself++; 00498 } else { 00499 gs2self_hist_map[igs] = Invalid; 00500 } 00501 igs++; 00502 }; 00503 00504 // make sure the count of non-ghost initial jets is equal to 00505 // what we already have in terms of initial jets 00506 assert(iself == _history.size()); 00507 00508 // if there was no clustering in this event (e.g. SISCone passive 00509 // area with zero input particles, or with a pt cut on stable cones 00510 // that kills all jets), then don't bother with the rest (which 00511 // would crash!) 00512 if (igs == gs_history.size()) return; 00513 00514 // now actually transfer things 00515 do { 00516 // if we are a pure ghost, then go on to next round 00517 if (ghosted_seq.is_pure_ghost(igs)) { 00518 gs2self_hist_map[igs] = Invalid; 00519 continue; 00520 } 00521 00522 const history_element & gs_hist_el = gs_history[igs]; 00523 00524 bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1); 00525 bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2); 00526 00527 // if exactly one parent is a ghost then maintain info about the 00528 // non-ghost correspondence for this jet, and then go on to next 00529 // recombination in the ghosted sequence 00530 if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) { 00531 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2]; 00532 continue; 00533 } 00534 if (!parent1_is_ghost && parent2_is_ghost) { 00535 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1]; 00536 continue; 00537 } 00538 00539 // no parents are ghosts... 00540 if (gs_hist_el.parent2 >= 0) { 00541 // recombination of two non-ghosts 00542 gs2self_hist_map[igs] = _history.size(); 00543 // record the recombination in our own sequence 00544 int newjet_k; // dummy var -- not used 00545 int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index; 00546 int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index; 00547 //cerr << "recombining "<< jet_i << " and "<< jet_j << endl; 00548 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k); 00549 } else { 00550 // we have a non-ghost that has become a beam-jet 00551 assert(gs_history[igs].parent2 == BeamJet); 00552 // record position 00553 gs2self_hist_map[igs] = _history.size(); 00554 // record the recombination in our own sequence 00555 _do_iB_recombination_step( 00556 _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index, 00557 gs_hist_el.dij); 00558 } 00559 } while (++igs < gs_history.size()); 00560 00561 } 00562 00563 //---------------------------------------------------------------------- 00564 void ClusterSequenceActiveArea::_transfer_areas( 00565 const vector<int> & unique_hist_order, 00566 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq ) { 00567 00568 const vector<history_element> & gs_history = ghosted_seq.history(); 00569 const vector<PseudoJet> & gs_jets = ghosted_seq.jets(); 00570 vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order(); 00571 00572 const double tolerance = 1e-11; // to decide when two jets are the same 00573 00574 int j = -1; 00575 int hist_index = -1; 00576 00577 valarray<double> our_areas(_history.size()); 00578 our_areas = 0.0; 00579 00580 valarray<PseudoJet> our_area_4vectors(_history.size()); 00581 our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0); 00582 00583 for (unsigned i = 0; i < gs_history.size(); i++) { 00584 // only consider composite particles 00585 unsigned gs_hist_index = gs_unique_hist_order[i]; 00586 if (gs_hist_index < ghosted_seq.n_particles()) continue; 00587 const history_element & gs_hist = gs_history[gs_unique_hist_order[i]]; 00588 int parent1 = gs_hist.parent1; 00589 int parent2 = gs_hist.parent2; 00590 00591 if (parent2 == BeamJet) { 00592 // need to look at parent to get the actual jet 00593 const PseudoJet & jet = 00594 gs_jets[gs_history[parent1].jetp_index]; 00595 double area_local = ghosted_seq.area(jet); 00596 PseudoJet ext_area = ghosted_seq.area_4vector(jet); 00597 00598 if (ghosted_seq.is_pure_ghost(parent1)) { 00599 // record the existence of the pure ghost jet for future use 00600 _ghost_jets.push_back(GhostJet(jet,area_local)); 00601 if (abs(jet.rap()) < _safe_rap_for_area) { 00602 _non_jet_area += area_local; 00603 _non_jet_area2 += area_local*area_local; 00604 _non_jet_number += 1; 00605 } 00606 } else { 00607 00608 // get next "combined-particle" index in our own history 00609 // making sure we don't go beyond its bounds (if we do 00610 // then we're in big trouble anyway...) 00611 while (++j < int(_history.size())) { 00612 hist_index = unique_hist_order[j]; 00613 if (hist_index >= _initial_n) break;} 00614 00615 // sanity checking -- do not overrun 00616 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching"); 00617 00618 // sanity check -- make sure we are taking about the same 00619 // jet in reference and new sequences 00620 int refjet_index = _history[_history[hist_index].parent1].jetp_index; 00621 assert(refjet_index >= 0 && refjet_index < int(_jets.size())); 00622 const PseudoJet & refjet = _jets[refjet_index]; 00623 00624 //cerr << "Inclusive" << endl; 00625 //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl; 00626 //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl; 00627 00628 // If pt disagrees check E; if they both disagree there's a 00629 // problem here... NB: a massive particle with zero pt may 00630 // have its pt changed when a ghost is added -- this is why we 00631 // also require the energy to be wrong before complaining 00632 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance, 00633 ghosted_seq); 00634 00635 // set the area at this clustering stage 00636 our_areas[hist_index] = area_local; 00637 our_area_4vectors[hist_index] = ext_area; 00638 00639 // update the parent as well -- that way its area is the area 00640 // immediately before clustering (i.e. resolve an ambiguity in 00641 // the Cambridge case and ensure in the kt case that the original 00642 // particles get a correct area) 00643 our_areas[_history[hist_index].parent1] = area_local; 00644 our_area_4vectors[_history[hist_index].parent1] = ext_area; 00645 00646 } 00647 } 00648 else if (!ghosted_seq.is_pure_ghost(parent1) && 00649 !ghosted_seq.is_pure_ghost(parent2)) { 00650 00651 // get next "combined-particle" index in our own history 00652 while (++j < int(_history.size())) { 00653 hist_index = unique_hist_order[j]; 00654 if (hist_index >= _initial_n) break;} 00655 00656 // sanity checking -- do not overrun 00657 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching"); 00658 00659 // make sure that our reference history entry is also for 00660 // an exclusive (dij) clustering (otherwise the comparison jet 00661 // will not exist) 00662 if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)"); 00663 00664 //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl; 00665 //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl; 00666 //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl; 00667 00668 const PseudoJet & jet = gs_jets[gs_hist.jetp_index]; 00669 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index]; 00670 00671 // run sanity check 00672 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance, 00673 ghosted_seq); 00674 00675 // update area and our local index (maybe redundant since later 00676 // the descendants will reupdate it?) 00677 double area_local = ghosted_seq.area(jet); 00678 our_areas[hist_index] += area_local; 00679 00680 PseudoJet ext_area = ghosted_seq.area_4vector(jet); 00681 00682 // GPS TMP debugging (jetclu) ----------------------- 00683 //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100); 00684 //our_area_4vectors[hist_index] = ext_area; 00685 //cout << "aa " 00686 // << our_area_4vectors[hist_index].px() << " " 00687 // << our_area_4vectors[hist_index].py() << " " 00688 // << our_area_4vectors[hist_index].pz() << " " 00689 // << our_area_4vectors[hist_index].E() << endl; 00690 //cout << "bb " 00691 // << ext_area.px() << " " 00692 // << ext_area.py() << " " 00693 // << ext_area.pz() << " " 00694 // << ext_area.E() << endl; 00695 //--------------------------------------------------- 00696 00697 _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area); 00698 00699 // now update areas of parents (so that they becomes areas 00700 // immediately before clustering occurred). This is of use 00701 // because it allows us to set the areas of the original hard 00702 // particles in the kt algorithm; for the Cambridge case it 00703 // means a jet's area will be the area just before it clusters 00704 // with another hard jet. 00705 const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index]; 00706 int our_parent1 = _history[hist_index].parent1; 00707 our_areas[our_parent1] = ghosted_seq.area(jet1); 00708 our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1); 00709 00710 const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index]; 00711 int our_parent2 = _history[hist_index].parent2; 00712 our_areas[our_parent2] = ghosted_seq.area(jet2); 00713 our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2); 00714 } 00715 00716 } 00717 00718 // now add unclustered ghosts to the relevant list so that we can 00719 // calculate empty area later. 00720 vector<PseudoJet> unclust = ghosted_seq.unclustered_particles(); 00721 for (unsigned iu = 0; iu < unclust.size(); iu++) { 00722 if (ghosted_seq.is_pure_ghost(unclust[iu])) { 00723 double area_local = ghosted_seq.area(unclust[iu]); 00724 _unclustered_ghosts.push_back(GhostJet(unclust[iu],area_local)); 00725 } 00726 } 00727 00728 /* 00729 * WARNING: 00730 * _average_area has explicitly been sized initially to 2*jets().size() 00731 * which can be bigger than our_areas (of size _history.size() 00732 * if there are some unclustered particles. 00733 * So we must take care about boundaries 00734 */ 00735 00736 for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){ 00737 _average_area[area_index] += our_areas[area_index]; 00738 _average_area2[area_index] += our_areas[area_index]*our_areas[area_index]; 00739 } 00740 00741 //_average_area_4vector += our_area_4vectors; 00742 // Use the proper recombination scheme when averaging the area_4vectors 00743 // over multiple ghost runs (i.e. the repeat stage); 00744 for (unsigned i = 0; i < our_area_4vectors.size(); i++) { 00745 _jet_def.recombiner()->plus_equal(_average_area_4vector[i], 00746 our_area_4vectors[i]); 00747 } 00748 } 00749 00750 00751 /// check if two jets have the same momentum to within the 00752 /// tolerance (and if pt's are not the same we're forgiving and 00753 /// look to see if the energy is the same) 00754 void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E( 00755 const PseudoJet & jet, 00756 const PseudoJet & refjet, 00757 double tolerance, 00758 const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq 00759 ) const { 00760 00761 if (abs(jet.perp2()-refjet.perp2()) > 00762 tolerance*max(jet.perp2(),refjet.perp2()) 00763 && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) { 00764 ostringstream ostr; 00765 ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas. See FAQ for possible explanations." << endl; 00766 ostr << " Ref-Jet: " 00767 << refjet.px() << " " 00768 << refjet.py() << " " 00769 << refjet.pz() << " " 00770 << refjet.E() << endl; 00771 ostr << " New-Jet: " 00772 << jet.px() << " " 00773 << jet.py() << " " 00774 << jet.pz() << " " 00775 << jet.E() << endl; 00776 if (jets_ghosted_seq.has_dangerous_particles()) { 00777 ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;} 00778 //ostr << jet.perp() << " " << refjet.perp() << " " 00779 // << jet.perp() - refjet.perp() << endl; 00780 throw Error(ostr.str()); 00781 } 00782 } 00783 00784 FASTJET_END_NAMESPACE 00785