00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "fastjet/PseudoJet.hh"
00032 #include "fastjet/ClusterSequence.hh"
00033 #include "fastjet/ClusterSequenceActiveArea.hh"
00034 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
00035 #include<iostream>
00036 #include<vector>
00037 #include<sstream>
00038 #include<algorithm>
00039 #include<cmath>
00040 #include<valarray>
00041
00042 FASTJET_BEGIN_NAMESPACE
00043
00044
00045 using namespace std;
00046
00047
00048
00049
00050
00051
00053 void ClusterSequenceActiveArea::_initialise_and_run_AA (
00054 const JetDefinition & jet_def,
00055 const GhostedAreaSpec & ghost_spec,
00056 const bool & writeout_combinations) {
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 }
00065
00066
00067 void ClusterSequenceActiveArea::_resize_and_zero_AA () {
00068
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 }
00075
00076
00077 void ClusterSequenceActiveArea::_initialise_AA (
00078 const JetDefinition & jet_def,
00079 const GhostedAreaSpec & ghost_spec,
00080 const bool & writeout_combinations,
00081 bool & continue_running)
00082 {
00083
00084
00085 _ghost_spec_repeat = ghost_spec.repeat();
00086
00087
00088 _resize_and_zero_AA();
00089
00090
00091 _maxrap_for_area = ghost_spec.ghost_maxrap();
00092 _safe_rap_for_area = _maxrap_for_area - jet_def.R();
00093
00094
00095
00096
00097
00098
00099
00100
00101 if (ghost_spec.repeat() <= 0) {
00102 _initialise_and_run(jet_def, writeout_combinations);
00103 continue_running = false;
00104 return;
00105 }
00106
00107
00108 _decant_options(jet_def, writeout_combinations);
00109
00110
00111
00112 _fill_initial_history();
00113
00114
00115 _has_dangerous_particles = false;
00116
00117 continue_running = true;
00118 }
00119
00120
00121
00122 void ClusterSequenceActiveArea::_run_AA (const GhostedAreaSpec & ghost_spec) {
00123
00124 vector<PseudoJet> input_jets(_jets);
00125
00126
00127 vector<int> unique_tree;
00128
00129
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
00138
00139 _transfer_ghost_free_history(clust_seq);
00140
00141 unique_tree = unique_history_order();
00142 }
00143
00144
00145 _transfer_areas(unique_tree, clust_seq);
00146 }
00147 }
00148
00149
00150
00152 void ClusterSequenceActiveArea::_postprocess_AA (const GhostedAreaSpec & ghost_spec) {
00153 _average_area /= ghost_spec.repeat();
00154 _average_area2 /= ghost_spec.repeat();
00155 if (ghost_spec.repeat() > 1) {
00156
00157
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
00171
00172
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
00177 }
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276 double ClusterSequenceActiveArea::pt_per_unit_area(
00277 mean_pt_strategies strat, double range) const {
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
00295 if (pt_over_areas.size() == 0) {return 0.0;}
00296
00297
00298
00299
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
00304
00305
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
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 }
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 double ClusterSequenceActiveArea::empty_area(const RangeDefinition & range) const {
00446 double empty = 0.0;
00447
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
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 }
00462
00463
00464 double ClusterSequenceActiveArea::n_empty_jets(const RangeDefinition & range) const {
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 }
00472
00473
00476 void ClusterSequenceActiveArea::_transfer_ghost_free_history(
00477 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) {
00478
00479 const vector<history_element> & gs_history = ghosted_seq.history();
00480 vector<int> gs2self_hist_map(gs_history.size());
00481
00482
00483
00484 _strategy = ghosted_seq.strategy_used();
00485
00486
00487 unsigned igs = 0;
00488 unsigned iself = 0;
00489 while (igs < gs_history.size() && gs_history[igs].parent1 == InexistentParent) {
00490
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
00500
00501 assert(iself == _history.size());
00502
00503
00504
00505
00506
00507 if (igs == gs_history.size()) return;
00508
00509
00510 do {
00511
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
00523
00524
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
00535 if (gs_hist_el.parent2 >= 0) {
00536
00537 gs2self_hist_map[igs] = _history.size();
00538
00539 int newjet_k;
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
00543 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
00544 } else {
00545
00546 assert(gs_history[igs].parent2 == BeamJet);
00547
00548 gs2self_hist_map[igs] = _history.size();
00549
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 }
00557
00558
00559 void ClusterSequenceActiveArea::_transfer_areas(
00560 const vector<int> & unique_hist_order,
00561 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq ) {
00562
00563 const vector<history_element> & gs_history = ghosted_seq.history();
00564 const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
00565 vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
00566
00567 const double tolerance = 1e-11;
00568
00569 int j = -1;
00570 int hist_index = -1;
00571
00572 valarray<double> our_areas(_history.size());
00573 our_areas = 0.0;
00574
00575 valarray<PseudoJet> our_area_4vectors(_history.size());
00576 our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
00577
00578 for (unsigned i = 0; i < gs_history.size(); i++) {
00579
00580 unsigned gs_hist_index = gs_unique_hist_order[i];
00581 if (gs_hist_index < ghosted_seq.n_particles()) continue;
00582 const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
00583 int parent1 = gs_hist.parent1;
00584 int parent2 = gs_hist.parent2;
00585
00586 if (parent2 == BeamJet) {
00587
00588 const PseudoJet & jet =
00589 gs_jets[gs_history[parent1].jetp_index];
00590 double area = ghosted_seq.area(jet);
00591 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00592
00593 if (ghosted_seq.is_pure_ghost(parent1)) {
00594
00595 _ghost_jets.push_back(GhostJet(jet,area));
00596 if (abs(jet.rap()) < _safe_rap_for_area) {
00597 _non_jet_area += area;
00598 _non_jet_area2 += area*area;
00599 _non_jet_number += 1;
00600 }
00601 } else {
00602
00603
00604
00605
00606 while (++j < int(_history.size())) {
00607 hist_index = unique_hist_order[j];
00608 if (hist_index >= _initial_n) break;}
00609
00610
00611 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
00612
00613
00614
00615 int refjet_index = _history[_history[hist_index].parent1].jetp_index;
00616 assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
00617 const PseudoJet & refjet = _jets[refjet_index];
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00628 ghosted_seq);
00629
00630
00631 our_areas[hist_index] = area;
00632 our_area_4vectors[hist_index] = ext_area;
00633
00634
00635
00636
00637
00638 our_areas[_history[hist_index].parent1] = area;
00639 our_area_4vectors[_history[hist_index].parent1] = ext_area;
00640
00641 }
00642 }
00643 else if (!ghosted_seq.is_pure_ghost(parent1) &&
00644 !ghosted_seq.is_pure_ghost(parent2)) {
00645
00646
00647 while (++j < int(_history.size())) {
00648 hist_index = unique_hist_order[j];
00649 if (hist_index >= _initial_n) break;}
00650
00651
00652 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
00653
00654
00655
00656
00657 if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
00658
00659
00660
00661
00662
00663 const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
00664 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
00665
00666
00667 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00668 ghosted_seq);
00669
00670
00671
00672 double area = ghosted_seq.area(jet);
00673 our_areas[hist_index] += area;
00674
00675 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692 _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
00693
00694
00695
00696
00697
00698
00699
00700 const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
00701 int our_parent1 = _history[hist_index].parent1;
00702 our_areas[our_parent1] = ghosted_seq.area(jet1);
00703 our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
00704
00705 const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
00706 int our_parent2 = _history[hist_index].parent2;
00707 our_areas[our_parent2] = ghosted_seq.area(jet2);
00708 our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
00709 }
00710
00711 }
00712
00713
00714
00715 vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
00716 for (unsigned iu = 0; iu < unclust.size(); iu++) {
00717 if (ghosted_seq.is_pure_ghost(unclust[iu])) {
00718 double area = ghosted_seq.area(unclust[iu]);
00719 _unclustered_ghosts.push_back(GhostJet(unclust[iu],area));
00720 }
00721 }
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731 for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){
00732 _average_area[area_index] += our_areas[area_index];
00733 _average_area2[area_index] += our_areas[area_index]*our_areas[area_index];
00734 }
00735
00736
00737
00738
00739 for (unsigned i = 0; i < our_area_4vectors.size(); i++) {
00740 _jet_def.recombiner()->plus_equal(_average_area_4vector[i],
00741 our_area_4vectors[i]);
00742 }
00743 }
00744
00745
00749 void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E(
00750 const PseudoJet & jet,
00751 const PseudoJet & refjet,
00752 double tolerance,
00753 const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
00754 ) const {
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
00774
00775 throw Error(ostr.str());
00776 }
00777 }
00778
00779 FASTJET_END_NAMESPACE
00780