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
00039 FASTJET_BEGIN_NAMESPACE
00040
00041
00042 using namespace std;
00043
00044
00045
00046
00047
00048
00050 void ClusterSequenceActiveArea::_initialise_and_run_AA (
00051 const JetDefinition & jet_def,
00052 const GhostedAreaSpec & area_spec,
00053 const bool & writeout_combinations) {
00054
00055 bool continue_running;
00056 _initialise_AA(jet_def, area_spec, writeout_combinations, continue_running);
00057 if (continue_running) {
00058 _run_AA(area_spec);
00059 _postprocess_AA(area_spec);
00060 }
00061 }
00062
00063
00064 void ClusterSequenceActiveArea::_resize_and_zero_AA () {
00065
00066 _average_area.resize(2*_jets.size()); _average_area = 0.0;
00067 _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
00068 _average_area_4vector.resize(2*_jets.size());
00069 _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
00070 _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
00071 }
00072
00073
00074 void ClusterSequenceActiveArea::_initialise_AA (
00075 const JetDefinition & jet_def,
00076 const GhostedAreaSpec & area_spec,
00077 const bool & writeout_combinations,
00078 bool & continue_running)
00079 {
00080
00081
00082 _area_spec_repeat = area_spec.repeat();
00083
00084
00085 _resize_and_zero_AA();
00086
00087
00088 _maxrap_for_area = area_spec.ghost_maxrap();
00089 _safe_rap_for_area = _maxrap_for_area - jet_def.R();
00090
00091
00092
00093
00094
00095
00096
00097
00098 if (area_spec.repeat() <= 0) {
00099 _initialise_and_run(jet_def, writeout_combinations);
00100 continue_running = false;
00101 return;
00102 }
00103
00104
00105 _decant_options(jet_def, writeout_combinations);
00106
00107
00108
00109 _fill_initial_history();
00110
00111
00112 _has_dangerous_particles = false;
00113
00114 continue_running = true;
00115 }
00116
00117
00118
00119 void ClusterSequenceActiveArea::_run_AA (const GhostedAreaSpec & area_spec) {
00120
00121 vector<PseudoJet> input_jets(_jets);
00122
00123
00124 vector<int> unique_tree;
00125
00126
00127 for (int irepeat = 0; irepeat < area_spec.repeat(); irepeat++) {
00128
00129 ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
00130 jet_def(), area_spec);
00131
00132 _has_dangerous_particles |= clust_seq.has_dangerous_particles();
00133 if (irepeat == 0) {
00134
00135
00136 _transfer_ghost_free_history(clust_seq);
00137
00138 unique_tree = unique_history_order();
00139 }
00140
00141
00142 _transfer_areas(unique_tree, clust_seq);
00143 }
00144 }
00145
00146
00147
00149 void ClusterSequenceActiveArea::_postprocess_AA (const GhostedAreaSpec & area_spec) {
00150 _average_area /= area_spec.repeat();
00151 _average_area2 /= area_spec.repeat();
00152 if (area_spec.repeat() > 1) {
00153 _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/
00154 (area_spec.repeat()-1));
00155 } else {
00156 _average_area2 = 0.0;
00157 }
00158
00159 _non_jet_area /= area_spec.repeat();
00160 _non_jet_area2 /= area_spec.repeat();
00161 _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
00162 area_spec.repeat());
00163 _non_jet_number /= area_spec.repeat();
00164
00165
00166
00167
00168 for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
00169 _average_area_4vector[i] = (1.0/area_spec.repeat()) * _average_area_4vector[i];
00170 }
00171
00172 }
00173
00174
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 double ClusterSequenceActiveArea::pt_per_unit_area(
00272 mean_pt_strategies strat, double range) const {
00273
00274 vector<PseudoJet> incl_jets = inclusive_jets();
00275 vector<double> pt_over_areas;
00276
00277 for (unsigned i = 0; i < incl_jets.size(); i++) {
00278 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
00279 double this_area;
00280 if ( strat == median_4vector ) {
00281 this_area = area_4vector(incl_jets[i]).perp();
00282 } else {
00283 this_area = area(incl_jets[i]);
00284 }
00285 pt_over_areas.push_back(incl_jets[i].perp()/this_area);
00286 }
00287 }
00288
00289
00290 if (pt_over_areas.size() == 0) {return 0.0;}
00291
00292
00293
00294
00295 sort(pt_over_areas.begin(), pt_over_areas.end());
00296 double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];
00297
00298
00299
00300
00301 double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0;
00302 double nj_median_ratio;
00303 if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
00304 int int_nj_median = int(nj_median_pos);
00305 nj_median_ratio =
00306 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00307 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00308 } else {
00309 nj_median_ratio = 0.0;
00310 }
00311
00312
00313
00314 double pt_sum = 0.0, pt_sum_with_cut = 0.0;
00315 double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area;
00316 double ratio_sum = 0.0;
00317 double ratio_n = _non_jet_number;
00318 for (unsigned i = 0; i < incl_jets.size(); i++) {
00319 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
00320 double this_area;
00321 if ( strat == median_4vector ) {
00322 this_area = area_4vector(incl_jets[i]).perp();
00323 } else {
00324 this_area = area(incl_jets[i]);
00325 }
00326 pt_sum += incl_jets[i].perp();
00327 area_sum += this_area;
00328 double ratio = incl_jets[i].perp()/this_area;
00329 if (ratio < range*nj_median_ratio) {
00330 pt_sum_with_cut += incl_jets[i].perp();
00331 area_sum_with_cut += this_area;
00332 ratio_sum += ratio; ratio_n++;
00333 }
00334 }
00335 }
00336
00337 if (strat == play) {
00338 double trunc_sum = 0, trunc_sumsqr = 0;
00339 vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size());
00340 for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) {
00341 double ratio = pt_over_areas[i];
00342 trunc_sum += ratio;
00343 trunc_sumsqr += ratio*ratio;
00344 means[i] = trunc_sum / (i+1);
00345 sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1)));
00346 cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<<
00347 sd[i]/sqrt(i+1.0)<<endl;
00348 }
00349 cout << "-----------------------------------"<<endl;
00350 for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
00351 cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl;
00352 }
00353 cout << "Number of non-jets: "<<_non_jet_number<<endl;
00354 cout << "Area of non-jets: "<<_non_jet_area<<endl;
00355 cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl;
00356 cout << "NJ median position: " << nj_median_pos <<endl;
00357 cout << "NJ median value: " << nj_median_ratio <<endl;
00358 return 0.0;
00359 }
00360
00361 switch(strat) {
00362 case median:
00363 case median_4vector:
00364 return nj_median_ratio;
00365 case non_ghost_median:
00366 return non_ghost_median_ratio;
00367 case pttot_over_areatot:
00368 return pt_sum / area_sum;
00369 case pttot_over_areatot_cut:
00370 return pt_sum_with_cut / area_sum_with_cut;
00371 case mean_ratio_cut:
00372 return ratio_sum/ratio_n;
00373 default:
00374 return nj_median_ratio;
00375 }
00376
00377 }
00378
00379
00380
00381
00382
00383 void ClusterSequenceActiveArea::parabolic_pt_per_unit_area(
00384 double & a, double & b, double raprange, double exclude_above,
00385 bool use_area_4vector) const {
00386
00387 double this_raprange;
00388 if (raprange <= 0) {this_raprange = _safe_rap_for_area;}
00389 else {this_raprange = raprange;}
00390
00391 int n=0;
00392 int n_excluded = 0;
00393 double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
00394
00395 vector<PseudoJet> incl_jets = inclusive_jets();
00396
00397 for (unsigned i = 0; i < incl_jets.size(); i++) {
00398 if (abs(incl_jets[i].rap()) < this_raprange) {
00399 double this_area;
00400 if ( use_area_4vector ) {
00401 this_area = area_4vector(incl_jets[i]).perp();
00402 } else {
00403 this_area = area(incl_jets[i]);
00404 }
00405 double f = incl_jets[i].perp()/this_area;
00406 if (exclude_above <= 0.0 || f < exclude_above) {
00407 double x = incl_jets[i].rap(); double x2 = x*x;
00408 mean_f += f;
00409 mean_x2 += x2;
00410 mean_x4 += x2*x2;
00411 mean_fx2 += f*x2;
00412 n++;
00413 } else {
00414 n_excluded++;
00415 }
00416 }
00417 }
00418
00419 if (n <= 1) {
00420
00421
00422
00423 a = 0.0;
00424 b = 0.0;
00425 } else {
00426 mean_f /= n;
00427 mean_x2 /= n;
00428 mean_x4 /= n;
00429 mean_fx2 /= n;
00430
00431 b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
00432 a = mean_f - b*mean_x2;
00433 }
00434
00435 }
00436
00437
00438
00439 double ClusterSequenceActiveArea::empty_area(const RangeDefinition & range) const {
00440 double empty = 0.0;
00441
00442 for (unsigned i = 0; i < _ghost_jets.size(); i++) {
00443 if (range.is_in_range(_ghost_jets[i])) {
00444 empty += _ghost_jets[i].area;
00445 }
00446 }
00447
00448 for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) {
00449 if (range.is_in_range(_unclustered_ghosts[i])) {
00450 empty += _unclustered_ghosts[i].area;
00451 }
00452 }
00453 empty /= _area_spec_repeat;
00454 return empty;
00455 }
00456
00457
00458 double ClusterSequenceActiveArea::n_empty_jets(const RangeDefinition & range) const {
00459 double inrange = 0;
00460 for (unsigned i = 0; i < _ghost_jets.size(); i++) {
00461 if (range.is_in_range(_ghost_jets[i])) inrange++;
00462 }
00463 inrange /= _area_spec_repeat;
00464 return inrange;
00465 }
00466
00467
00470 void ClusterSequenceActiveArea::_transfer_ghost_free_history(
00471 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) {
00472
00473 const vector<history_element> & gs_history = ghosted_seq.history();
00474 vector<int> gs2self_hist_map(gs_history.size());
00475
00476
00477 unsigned igs = 0;
00478 unsigned iself = 0;
00479 while (gs_history[igs].parent1 == InexistentParent) {
00480
00481 if (!ghosted_seq.is_pure_ghost(igs)) {
00482 gs2self_hist_map[igs] = iself++;
00483 } else {
00484 gs2self_hist_map[igs] = Invalid;
00485 }
00486 igs++;
00487 };
00488
00489
00490
00491 assert(iself == _history.size());
00492
00493
00494 do {
00495
00496 if (ghosted_seq.is_pure_ghost(igs)) {
00497 gs2self_hist_map[igs] = Invalid;
00498 continue;
00499 }
00500
00501 const history_element & gs_hist_el = gs_history[igs];
00502
00503 bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
00504 bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);
00505
00506
00507
00508
00509 if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
00510 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
00511 continue;
00512 }
00513 if (!parent1_is_ghost && parent2_is_ghost) {
00514 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
00515 continue;
00516 }
00517
00518
00519 if (gs_hist_el.parent2 >= 0) {
00520
00521 gs2self_hist_map[igs] = _history.size();
00522
00523 int newjet_k;
00524 int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
00525 int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
00526
00527 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
00528 } else {
00529
00530 assert(gs_history[igs].parent2 == BeamJet);
00531
00532 gs2self_hist_map[igs] = _history.size();
00533
00534 _do_iB_recombination_step(
00535 _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
00536 gs_hist_el.dij);
00537 }
00538 } while (++igs < gs_history.size());
00539
00540
00541
00542 _strategy = ghosted_seq.strategy_used();
00543 }
00544
00545
00546 void ClusterSequenceActiveArea::_transfer_areas(
00547 const vector<int> & unique_hist_order,
00548 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq ) {
00549
00550 const vector<history_element> & gs_history = ghosted_seq.history();
00551 const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
00552 vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
00553
00554 const double tolerance = 1e-11;
00555
00556 int j = -1;
00557 int hist_index = -1;
00558
00559 valarray<double> our_areas(_history.size());
00560 our_areas = 0.0;
00561
00562 valarray<PseudoJet> our_area_4vectors(_history.size());
00563 our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
00564
00565 for (unsigned i = 0; i < gs_history.size(); i++) {
00566
00567 unsigned gs_hist_index = gs_unique_hist_order[i];
00568 if (gs_hist_index < ghosted_seq.n_particles()) continue;
00569 const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
00570 int parent1 = gs_hist.parent1;
00571 int parent2 = gs_hist.parent2;
00572
00573 if (parent2 == BeamJet) {
00574
00575 const PseudoJet & jet =
00576 gs_jets[gs_history[parent1].jetp_index];
00577 double area = ghosted_seq.area(jet);
00578 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00579
00580 if (ghosted_seq.is_pure_ghost(parent1)) {
00581
00582 _ghost_jets.push_back(GhostJet(jet,area));
00583 if (abs(jet.rap()) < _safe_rap_for_area) {
00584 _non_jet_area += area;
00585 _non_jet_area2 += area*area;
00586 _non_jet_number += 1;
00587 }
00588 } else {
00589
00590
00591
00592
00593 while (++j < int(_history.size())) {
00594 hist_index = unique_hist_order[j];
00595 if (hist_index >= _initial_n) break;}
00596
00597
00598 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
00599
00600
00601
00602 int refjet_index = _history[_history[hist_index].parent1].jetp_index;
00603 assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
00604 const PseudoJet & refjet = _jets[refjet_index];
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00615 ghosted_seq);
00616
00617
00618 our_areas[hist_index] = area;
00619 our_area_4vectors[hist_index] = ext_area;
00620
00621
00622
00623
00624
00625 our_areas[_history[hist_index].parent1] = area;
00626 our_area_4vectors[_history[hist_index].parent1] = ext_area;
00627
00628 }
00629 }
00630 else if (!ghosted_seq.is_pure_ghost(parent1) &&
00631 !ghosted_seq.is_pure_ghost(parent2)) {
00632
00633
00634 while (++j < int(_history.size())) {
00635 hist_index = unique_hist_order[j];
00636 if (hist_index >= _initial_n) break;}
00637
00638
00639 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
00640
00641
00642
00643
00644 if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
00645
00646
00647
00648
00649
00650 const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
00651 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
00652
00653
00654 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00655 ghosted_seq);
00656
00657
00658
00659 double area = ghosted_seq.area(jet);
00660 our_areas[hist_index] += area;
00661
00662 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679 _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
00680
00681
00682
00683
00684
00685
00686
00687 const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
00688 int our_parent1 = _history[hist_index].parent1;
00689 our_areas[our_parent1] = ghosted_seq.area(jet1);
00690 our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
00691
00692 const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
00693 int our_parent2 = _history[hist_index].parent2;
00694 our_areas[our_parent2] = ghosted_seq.area(jet2);
00695 our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
00696 }
00697
00698 }
00699
00700
00701
00702 vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
00703 for (unsigned iu = 0; iu < unclust.size(); iu++) {
00704 if (ghosted_seq.is_pure_ghost(unclust[iu])) {
00705 double area = ghosted_seq.area(unclust[iu]);
00706 _unclustered_ghosts.push_back(GhostJet(unclust[iu],area));
00707 }
00708 }
00709
00710 _average_area += our_areas;
00711 _average_area2 += our_areas*our_areas;
00712
00713
00714
00715
00716 for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
00717 _jet_def.recombiner()->plus_equal(_average_area_4vector[i],
00718 our_area_4vectors[i]);
00719 }
00720 }
00721
00722
00726 void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E(
00727 const PseudoJet & jet,
00728 const PseudoJet & refjet,
00729 double tolerance,
00730 const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
00731 ) const {
00732
00733 if (abs(jet.perp2()-refjet.perp2()) >
00734 tolerance*max(jet.perp2(),refjet.perp2())
00735 && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
00736 ostringstream ostr;
00737 ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas" << endl;
00738 ostr << " Ref-Jet: "
00739 << refjet.px() << " "
00740 << refjet.py() << " "
00741 << refjet.pz() << " "
00742 << refjet.E() << endl;
00743 ostr << " New-Jet: "
00744 << jet.px() << " "
00745 << jet.py() << " "
00746 << jet.pz() << " "
00747 << jet.E() << endl;
00748 if (jets_ghosted_seq.has_dangerous_particles()) {
00749 ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
00750
00751
00752 throw Error(ostr.str());
00753 }
00754 }
00755
00756 FASTJET_END_NAMESPACE
00757