ClusterSequenceActiveArea.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: ClusterSequenceActiveArea.cc 957 2007-11-22 18:58:45Z salam $
00003 //
00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
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      // defined in fastjet/internal/base.hh
00040 
00041 
00042 using namespace std;
00043 
00044 
00045 //int ClusterSequenceActiveArea::_n_seed_warnings = 0;
00046 //const int _max_seed_warnings = 10;
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   // initialize our local area information
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 //---------------------------------a-------------------------------------
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   // store this for future use
00082   _area_spec_repeat = area_spec.repeat();
00083 
00084   // make sure placeholders are there & zeroed
00085   _resize_and_zero_AA();
00086      
00087   // for future reference...
00088   _maxrap_for_area = area_spec.ghost_maxrap();
00089   _safe_rap_for_area = _maxrap_for_area - jet_def.R();
00090 
00091   // Make sure we'll have at least one repetition -- then we can
00092   // deduce the unghosted clustering sequence from one of the ghosted
00093   // sequences. If we do not have any repetitions, then get the
00094   // unghosted sequence from the plain unghosted clustering.
00095   //
00096   // NB: all decanting and filling of initial history will then
00097   // be carried out by base-class routine
00098   if (area_spec.repeat() <= 0) {
00099     _initialise_and_run(jet_def, writeout_combinations);
00100     continue_running = false;
00101     return;
00102   }
00103 
00104   // transfer all relevant info into internal variables
00105   _decant_options(jet_def, writeout_combinations);
00106 
00107   // set up the history entries for the initial particles (those
00108   // currently in _jets)
00109   _fill_initial_history();
00110 
00111   // by default it does not...
00112   _has_dangerous_particles = false;
00113   
00114   continue_running = true;
00115 }
00116 
00117 
00118 //----------------------------------------------------------------------
00119 void ClusterSequenceActiveArea::_run_AA (const GhostedAreaSpec & area_spec) {
00120   // record the input jets as they are currently
00121   vector<PseudoJet> input_jets(_jets);
00122 
00123   // code for testing the unique tree
00124   vector<int> unique_tree;
00125 
00126   // run the clustering multiple times so as to get areas of all the jets
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       // take the non-ghost part of the history and put into our own
00135       // history.
00136       _transfer_ghost_free_history(clust_seq);
00137       // get the "unique" order that will be used for transferring all areas. 
00138       unique_tree = unique_history_order();
00139     }
00140 
00141     // transfer areas from clust_seq into our object
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   // following bizarre way of writing things is related to 
00166   // poverty of operations on PseudoJet objects (as well as some confusion
00167   // in one or two places)
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   //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
00172 }
00173 
00174 
00175 // //----------------------------------------------------------------------
00176 // void ClusterSequenceActiveArea::_initialise_and_run_AA (
00177 //              const JetDefinition & jet_def,
00178 //              const GhostedAreaSpec & area_spec,
00179 //              const bool & writeout_combinations) 
00180 // {
00181 // 
00182 //   // store this for future use
00183 //   _area_spec_repeat = area_spec.repeat();
00184 // 
00185 //   // initialize our local area information
00186 //   _average_area.resize(2*_jets.size());  _average_area  = 0.0;
00187 //   _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
00188 //   _average_area_4vector.resize(2*_jets.size()); 
00189 //   _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
00190 //   _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
00191 //      
00192 //   // for future reference...
00193 //   _maxrap_for_area = area_spec.ghost_maxrap();
00194 //   _safe_rap_for_area = _maxrap_for_area - jet_def.R();
00195 // 
00196 //   // Make sure we'll have at least one repetition -- then we can
00197 //   // deduce the unghosted clustering sequence from one of the ghosted
00198 //   // sequences. If we do not have any repetitions, then get the
00199 //   // unghosted sequence from the plain unghosted clustering.
00200 //   //
00201 //   // NB: all decanting and filling of initial history will then
00202 //   // be carried out by base-class routine
00203 //   if (area_spec.repeat() <= 0) {
00204 //     _initialise_and_run(jet_def, writeout_combinations);
00205 //     return;
00206 //   }
00207 // 
00208 //   // transfer all relevant info into internal variables
00209 //   _decant_options(jet_def, writeout_combinations);
00210 // 
00211 //   // set up the history entries for the initial particles (those
00212 //   // currently in _jets)
00213 //   _fill_initial_history();
00214 // 
00215 //   // record the input jets as they are currently
00216 //   vector<PseudoJet> input_jets(_jets);
00217 // 
00218 //   // code for testing the unique tree
00219 //   vector<int> unique_tree;
00220 // 
00221 //   
00222 //   
00223 // 
00224 //   // run the clustering multiple times so as to get areas of all the jets
00225 //   for (int irepeat = 0; irepeat < area_spec.repeat(); irepeat++) {
00226 // 
00227 //     ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, 
00228 //                                                       jet_def, area_spec);
00229 // 
00230 //     if (irepeat == 0) {
00231 //       // take the non-ghost part of the history and put into our own
00232 //       // history.
00233 //       _transfer_ghost_free_history(clust_seq);
00234 //       // get the "unique" order that will be used for transferring all areas. 
00235 //       unique_tree = unique_history_order();
00236 //     }
00237 // 
00238 //     // transfer areas from clust_seq into our object
00239 //     _transfer_areas(unique_tree, clust_seq);
00240 //   }
00241 //   
00242 //   _average_area  /= area_spec.repeat();
00243 //   _average_area2 /= area_spec.repeat();
00244 //   if (area_spec.repeat() > 1) {
00245 //     _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/
00246 //                           (area_spec.repeat()-1));
00247 //   } else {
00248 //     _average_area2 = 0.0;
00249 //   }
00250 // 
00251 //   _non_jet_area  /= area_spec.repeat();
00252 //   _non_jet_area2 /= area_spec.repeat();
00253 //   _non_jet_area2  = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
00254 //                       area_spec.repeat());
00255 //   _non_jet_number /= area_spec.repeat();
00256 // 
00257 //   // following bizarre way of writing things is related to 
00258 //   // poverty of operations on PseudoJet objects (as well as some confusion
00259 //   // in one or two places)
00260 //   for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
00261 //     _average_area_4vector[i] = (1.0/area_spec.repeat()) * _average_area_4vector[i];
00262 //   }
00263 //   //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
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   // there is nothing inside our region, so answer will always be zero
00290   if (pt_over_areas.size() == 0) {return 0.0;}
00291   
00292   // get median (pt/area) [this is the "old" median definition. It considers
00293   // only the "real" jets in calculating the median, i.e. excluding the
00294   // only-ghost ones]
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   // new median definition that takes into account non-jet area (i.e.
00299   // jets composed only of ghosts), and for fractional median position 
00300   // interpolates between the corresponding entries in the pt_over_areas array
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   // get various forms of mean (pt/area)
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 // fit a parabola to pt/area as a function of rapidity, using the
00382 // formulae of CCN28-36 (which actually fits f = a+b*x^2)
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     // meaningful results require at least two jets inside the
00421     // area -- mind you if there are empty jets we should be in 
00422     // any case doing something special...
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   //cerr << "n_excluded = "<< n_excluded << endl;
00435 }
00436 
00437 
00438 //----------------------------------------------------------------------
00439 double ClusterSequenceActiveArea::empty_area(const RangeDefinition & range) const {
00440   double empty = 0.0;
00441   // first deal with ghost jets
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   // then deal with unclustered ghosts
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   // work our way through to first non-trivial combination
00477   unsigned igs = 0;
00478   unsigned iself = 0;
00479   while (gs_history[igs].parent1 == InexistentParent) {
00480     // record correspondence 
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   // make sure the count of non-ghost initial jets is equal to
00490   // what we already have in terms of initial jets
00491   assert(iself == _history.size());
00492 
00493   // now actually transfer things
00494   do  {
00495     // if we are a pure ghost, then go on to next round
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     // if exactly one parent is a ghost then maintain info about the
00507     // non-ghost correspondence for this jet, and then go on to next
00508     // recombination in the ghosted sequence
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     // no parents are ghosts...
00519     if (gs_hist_el.parent2 >= 0) {
00520       // recombination of two non-ghosts
00521       gs2self_hist_map[igs] = _history.size();
00522       // record the recombination in our own sequence
00523       int newjet_k; // dummy var -- not used
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       //cerr << "recombining "<< jet_i << " and "<< jet_j << endl;
00527       _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
00528     } else {
00529       // we have a non-ghost that has become a beam-jet
00530       assert(gs_history[igs].parent2 == BeamJet);
00531       // record position
00532       gs2self_hist_map[igs] = _history.size();
00533       // record the recombination in our own sequence
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   // finally transfer info about strategy used (which isn't necessarily
00541   // always the one that got asked for...)
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; // to decide when two jets are the same
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     // only consider composite particles
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       // need to look at parent to get the actual jet
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         // record the existence of the pure ghost jet for future use
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         // get next "combined-particle" index in our own history
00591         // making sure we don't go beyond its bounds (if we do
00592         // then we're in big trouble anyway...)
00593         while (++j < int(_history.size())) {
00594           hist_index = unique_hist_order[j];
00595           if (hist_index >= _initial_n) break;}
00596 
00597         // sanity checking -- do not overrun
00598         if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
00599 
00600         // sanity check -- make sure we are taking about the same 
00601         // jet in reference and new sequences
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       //cerr << "Inclusive" << endl;
00607       //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl;
00608       //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl;
00609 
00610         // If pt disagrees check E; if they both disagree there's a
00611         // problem here... NB: a massive particle with zero pt may
00612         // have its pt changed when a ghost is added -- this is why we
00613         // also require the energy to be wrong before complaining
00614         _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00615                                                ghosted_seq);
00616 
00617         // set the area at this clustering stage
00618         our_areas[hist_index]  = area; 
00619         our_area_4vectors[hist_index]  = ext_area; 
00620 
00621         // update the parent as well -- that way its area is the area
00622         // immediately before clustering (i.e. resolve an ambiguity in
00623         // the Cambridge case and ensure in the kt case that the original
00624         // particles get a correct area)
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       // get next "combined-particle" index in our own history
00634       while (++j < int(_history.size())) {
00635         hist_index = unique_hist_order[j];
00636         if (hist_index >= _initial_n) break;}
00637       
00638       // sanity checking -- do not overrun
00639       if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
00640 
00641       // make sure that our reference history entry is also for
00642       // an exclusive (dij) clustering (otherwise the comparison jet
00643       // will not exist)
00644       if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
00645 
00646       //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl;
00647       //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl;
00648       //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl;
00649 
00650       const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
00651       const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
00652 
00653       // run sanity check 
00654       _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00655                                              ghosted_seq);
00656 
00657       // update area and our local index (maybe redundant since later
00658       // the descendants will reupdate it?)
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       // GPS TMP debugging (jetclu) -----------------------
00665       //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100);
00666       //our_area_4vectors[hist_index] = ext_area;
00667       //cout << "aa " 
00668       //     << our_area_4vectors[hist_index].px() << " "
00669       //     << our_area_4vectors[hist_index].py() << " "
00670       //     << our_area_4vectors[hist_index].pz() << " "
00671       //     << our_area_4vectors[hist_index].E() << endl;
00672       //cout << "bb " 
00673       //     << ext_area.px() << " "
00674       //     << ext_area.py() << " "
00675       //     << ext_area.pz() << " "
00676       //     << ext_area.E() << endl;
00677       //---------------------------------------------------
00678 
00679       _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
00680 
00681       // now update areas of parents (so that they becomes areas
00682       // immediately before clustering occurred). This is of use
00683       // because it allows us to set the areas of the original hard
00684       // particles in the kt algorithm; for the Cambridge case it
00685       // means a jet's area will be the area just before it clusters
00686       // with another hard jet.
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   // now add unclustered ghosts to the relevant list so that we can
00701   // calculate empty area later.
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   //_average_area_4vector += our_area_4vectors;
00714   // Use the proper recombination scheme when averaging the area_4vectors
00715   // over multiple ghost runs (i.e. the repeat stage);
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     //ostr << jet.perp() << " " << refjet.perp() << " "
00751     //     << jet.perp() - refjet.perp() << endl;
00752     throw Error(ostr.str());
00753   }
00754 }
00755 
00756 FASTJET_END_NAMESPACE
00757 

Generated on Tue Dec 18 17:05:03 2007 for fastjet by  doxygen 1.5.2