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
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 double ClusterSequenceActiveArea::empty_area(const RangeDefinition & range) const {
00441 double empty = 0.0;
00442
00443 for (unsigned i = 0; i < _ghost_jets.size(); i++) {
00444 if (range.is_in_range(_ghost_jets[i])) {
00445 empty += _ghost_jets[i].area;
00446 }
00447 }
00448
00449 for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) {
00450 if (range.is_in_range(_unclustered_ghosts[i])) {
00451 empty += _unclustered_ghosts[i].area;
00452 }
00453 }
00454 empty /= _area_spec_repeat;
00455 return empty;
00456 }
00457
00458
00459 double ClusterSequenceActiveArea::n_empty_jets(const RangeDefinition & range) const {
00460 double inrange = 0;
00461 for (unsigned i = 0; i < _ghost_jets.size(); i++) {
00462 if (range.is_in_range(_ghost_jets[i])) inrange++;
00463 }
00464 inrange /= _area_spec_repeat;
00465 return inrange;
00466 }
00467
00468
00471 void ClusterSequenceActiveArea::_transfer_ghost_free_history(
00472 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) {
00473
00474 const vector<history_element> & gs_history = ghosted_seq.history();
00475 vector<int> gs2self_hist_map(gs_history.size());
00476
00477
00478 unsigned igs = 0;
00479 unsigned iself = 0;
00480 while (gs_history[igs].parent1 == InexistentParent) {
00481
00482 if (!ghosted_seq.is_pure_ghost(igs)) {
00483 gs2self_hist_map[igs] = iself++;
00484 } else {
00485 gs2self_hist_map[igs] = Invalid;
00486 }
00487 igs++;
00488 };
00489
00490
00491
00492 assert(iself == _history.size());
00493
00494
00495 do {
00496
00497 if (ghosted_seq.is_pure_ghost(igs)) {
00498 gs2self_hist_map[igs] = Invalid;
00499 continue;
00500 }
00501
00502 const history_element & gs_hist_el = gs_history[igs];
00503
00504 bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
00505 bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);
00506
00507
00508
00509
00510 if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
00511 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
00512 continue;
00513 }
00514 if (!parent1_is_ghost && parent2_is_ghost) {
00515 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
00516 continue;
00517 }
00518
00519
00520 if (gs_hist_el.parent2 >= 0) {
00521
00522 gs2self_hist_map[igs] = _history.size();
00523
00524 int newjet_k;
00525 int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
00526 int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
00527
00528 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
00529 } else {
00530
00531 assert(gs_history[igs].parent2 == BeamJet);
00532
00533 gs2self_hist_map[igs] = _history.size();
00534
00535 _do_iB_recombination_step(
00536 _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
00537 gs_hist_el.dij);
00538 }
00539 } while (++igs < gs_history.size());
00540
00541
00542
00543 _strategy = ghosted_seq.strategy_used();
00544 }
00545
00546
00547 void ClusterSequenceActiveArea::_transfer_areas(
00548 const vector<int> & unique_hist_order,
00549 const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq ) {
00550
00551 const vector<history_element> & gs_history = ghosted_seq.history();
00552 const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
00553 vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
00554
00555 const double tolerance = 1e-11;
00556
00557 int j = -1;
00558 int hist_index = -1;
00559
00560 valarray<double> our_areas(_history.size());
00561 our_areas = 0.0;
00562
00563 valarray<PseudoJet> our_area_4vectors(_history.size());
00564 our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
00565
00566 for (unsigned i = 0; i < gs_history.size(); i++) {
00567
00568 unsigned gs_hist_index = gs_unique_hist_order[i];
00569 if (gs_hist_index < ghosted_seq.n_particles()) continue;
00570 const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
00571 int parent1 = gs_hist.parent1;
00572 int parent2 = gs_hist.parent2;
00573
00574 if (parent2 == BeamJet) {
00575
00576 const PseudoJet & jet =
00577 gs_jets[gs_history[parent1].jetp_index];
00578 double area = ghosted_seq.area(jet);
00579 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00580
00581 if (ghosted_seq.is_pure_ghost(parent1)) {
00582
00583 _ghost_jets.push_back(GhostJet(jet,area));
00584 if (abs(jet.rap()) < _safe_rap_for_area) {
00585 _non_jet_area += area;
00586 _non_jet_area2 += area*area;
00587 _non_jet_number += 1;
00588 }
00589 } else {
00590
00591
00592
00593
00594 while (++j < int(_history.size())) {
00595 hist_index = unique_hist_order[j];
00596 if (hist_index >= _initial_n) break;}
00597
00598
00599 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
00600
00601
00602
00603 int refjet_index = _history[_history[hist_index].parent1].jetp_index;
00604 assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
00605 const PseudoJet & refjet = _jets[refjet_index];
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00616 ghosted_seq);
00617
00618
00619 our_areas[hist_index] = area;
00620 our_area_4vectors[hist_index] = ext_area;
00621
00622
00623
00624
00625
00626 our_areas[_history[hist_index].parent1] = area;
00627 our_area_4vectors[_history[hist_index].parent1] = ext_area;
00628
00629 }
00630 }
00631 else if (!ghosted_seq.is_pure_ghost(parent1) &&
00632 !ghosted_seq.is_pure_ghost(parent2)) {
00633
00634
00635 while (++j < int(_history.size())) {
00636 hist_index = unique_hist_order[j];
00637 if (hist_index >= _initial_n) break;}
00638
00639
00640 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
00641
00642
00643
00644
00645 if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
00646
00647
00648
00649
00650
00651 const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
00652 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
00653
00654
00655 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00656 ghosted_seq);
00657
00658
00659
00660 double area = ghosted_seq.area(jet);
00661 our_areas[hist_index] += area;
00662
00663 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680 _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
00681
00682
00683
00684
00685
00686
00687
00688 const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
00689 int our_parent1 = _history[hist_index].parent1;
00690 our_areas[our_parent1] = ghosted_seq.area(jet1);
00691 our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
00692
00693 const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
00694 int our_parent2 = _history[hist_index].parent2;
00695 our_areas[our_parent2] = ghosted_seq.area(jet2);
00696 our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
00697 }
00698
00699 }
00700
00701
00702
00703 vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
00704 for (unsigned iu = 0; iu < unclust.size(); iu++) {
00705 if (ghosted_seq.is_pure_ghost(unclust[iu])) {
00706 double area = ghosted_seq.area(unclust[iu]);
00707 _unclustered_ghosts.push_back(GhostJet(unclust[iu],area));
00708 }
00709 }
00710
00711 _average_area += our_areas;
00712 _average_area2 += our_areas*our_areas;
00713
00714
00715
00716
00717 for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
00718 _jet_def.recombiner()->plus_equal(_average_area_4vector[i],
00719 our_area_4vectors[i]);
00720 }
00721 }
00722
00723
00727 void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E(
00728 const PseudoJet & jet,
00729 const PseudoJet & refjet,
00730 double tolerance,
00731 const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
00732 ) const {
00733
00734 if (abs(jet.perp2()-refjet.perp2()) >
00735 tolerance*max(jet.perp2(),refjet.perp2())
00736 && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
00737 ostringstream ostr;
00738 ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas" << endl;
00739 ostr << " Ref-Jet: "
00740 << refjet.px() << " "
00741 << refjet.py() << " "
00742 << refjet.pz() << " "
00743 << refjet.E() << endl;
00744 ostr << " New-Jet: "
00745 << jet.px() << " "
00746 << jet.py() << " "
00747 << jet.pz() << " "
00748 << jet.E() << endl;
00749 if (jets_ghosted_seq.has_dangerous_particles()) {
00750 ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
00751
00752
00753 throw Error(ostr.str());
00754 }
00755 }
00756
00757 FASTJET_END_NAMESPACE
00758