FastJet 3.0.2
ClusterSequenceActiveArea.cc
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends