JETAN/fastjet/siscone/split_merge.h

   1 // -*- C++ -*-
   2 ///////////////////////////////////////////////////////////////////////////////
   3 // File: split_merge.h                                                       //
   4 // Description: header file for splitting/merging (contains the CJet class)  //
   5 // This file is part of the SISCone project.                                 //
   6 // For more details, see http://projects.hepforge.org/siscone                //
   7 //                                                                           //
   8 // Copyright (c) 2006 Gavin Salam and Gregory Soyez                          //
   9 //                                                                           //
  10 // This program is free software; you can redistribute it and/or modify      //
  11 // it under the terms of the GNU General Public License as published by      //
  12 // the Free Software Foundation; either version 2 of the License, or         //
  13 // (at your option) any later version.                                       //
  14 //                                                                           //
  15 // This program is distributed in the hope that it will be useful,           //
  16 // but WITHOUT ANY WARRANTY; without even the implied warranty of            //
  17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
  18 // GNU General Public License for more details.                              //
  19 //                                                                           //
  20 // You should have received a copy of the GNU General Public License         //
  21 // along with this program; if not, write to the Free Software               //
  22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
  23 //                                                                           //
  24 // $Revision:: 268                                                          $//
  25 // $Date:: 2009-03-12 21:24:16 +0100 (Thu, 12 Mar 2009)                     $//
  26 ///////////////////////////////////////////////////////////////////////////////
  27
  28 #ifndef __SPLIT_MERGE_H__
  29 #define __SPLIT_MERGE_H__
  30
  31 #include "defines.h"
  32 #include "geom_2d.h"
  33 #include "momentum.h"
  34 #include <stdio.h>
  35 #include <vector>
  36 #include <set>
  37 #include <memory>
  38 #include <string>
  39
  40 namespace siscone{
  41
  42 /**
  43  * \class Cjet
  44  * real Jet information.
  45  *
  46  * This class contains information for one single jet.
  47  * That is, first, its momentum carrying information
  48  * about its centre and pT, and second, its particle
  49  * contents
  50  */
  51 class Cjet{
  52  public:
  53   /// default ctor
  54   Cjet();
  55
  56   /// default dtor
  57   ~Cjet();
  58
  59   Cmomentum v;               ///< jet momentum
  60   double pt_tilde;           ///< p-scheme pt
  61   int n;                     ///< number of particles inside
  62   std::vector<int> contents; ///< particle contents (list of indices)
  63
  64   /// ordering variable used for ordering and overlap in the
  65   /// split--merge. This variable is automatically set either to
  66   /// pt_tilde, or to mt or to pt, depending on the siscone
  67   /// parameter. Note that the default behaviour is pt_tilde and that
  68   /// other chices may lead to infrared unsafe situations.
  69   /// Note: we use the square of the varible rather than the variable itself
  70   double sm_var2;
  71
  72   /// covered range in eta-phi
  73   Ceta_phi_range range;
  74
  75   /// pass at which the jet has been found
  76   /// It starts at 0 (first pass), -1 means infinite rapidity
  77   int pass;
  78 };
  79
  80 /// ordering of jets in pt (e.g. used in final jets ordering)
  81 bool jets_pt_less(const Cjet &j1, const Cjet &j2);
  82
  83
  84 /// the choices of scale variable that can be used in the split-merge
  85 /// step, both for ordering the protojets and for measuing their
  86 /// overlap; pt, Et and mt=sqrt(pt^2+m^2) are all defined in E-scheme
  87 /// (4-momentum) recombination; pttilde = \sum_{i\in jet} |p_{t,i}|
  88 ///
  89 /// NB: if one changes the order here, one _MUST_ also change the order
  90 ///     in the SISCone plugin
  91 enum Esplit_merge_scale {
  92            SM_pt,     ///< transverse momentum (E-scheme), IR unsafe
  93            SM_Et,     ///< transverse energy (E-scheme), not long. boost inv.
  94                       ///< original run-II choice [may not be implemented]
  95            SM_mt,     ///< transverse mass (E-scheme), IR safe except
  96                       ///< in decays of two identical narrow heavy particles
  97            SM_pttilde ///< pt-scheme pt = \sum_{i in jet} |p_{ti}|, should
  98                       ///< be IR safe in all cases
  99 };
 100
 101 /// return the name of the split-merge scale choice
 102 std::string split_merge_scale_name(Esplit_merge_scale sms);
 103
 104 /**
 105  * \class Csplit_merge_ptcomparison
 106  * comparison of jets for split--merge ordering
 107  *
 108  * a class that allows us to carry out comparisons of pt of jets, using
 109  * information from exact particle contents where necessary.
 110  */
 111 class Csplit_merge_ptcomparison{
 112 public:
 113   /// default ctor
 114   Csplit_merge_ptcomparison() :
 115     particles(0), split_merge_scale(SM_pttilde){};
 116
 117   /// return the name corresponding to the SM scale variable
 118   std::string SM_scale_name() const {
 119     return split_merge_scale_name(split_merge_scale);}
 120
 121   std::vector<Cmomentum> * particles; ///< pointer to the list of particles
 122   std::vector<double> * pt;           ///< pointer to the pt of the particles
 123
 124   /// comparison between 2 jets
 125   bool operator()(const Cjet &jet1, const Cjet &jet2) const;
 126
 127   /**
 128    * get the difference between 2 jets, calculated such that rounding
 129    * errors will not affect the result even if the two jets have
 130    * almost the same content (so that the difference is below the
 131    * rounding errors)
 132    *
 133    * \param j1        first jet
 134    * \param j2        second jet
 135    * \param v         jet1-jet2
 136    * \param pt_tilde  jet1-jet2 pt_tilde
 137    */
 138   void get_difference(const Cjet &j1, const Cjet &j2, Cmomentum *v, double *pt_tilde) const;
 139
 140   /// the following parameter controls the variable we're using for
 141   /// the split-merge process i.e. the variable we use for
 142   ///  1. ordering jet candidates;
 143   ///  2. computing te overlap fraction of two candidates.
 144   /// The default value uses pttile (p-scheme pt). Other alternatives are
 145   /// pt, mt=sqrt(pt^2+m^2)=sqrt(E^2-pz^2) or Et.
 146   /// NOTE: Modifying the default choice can have nasty effects:
 147   /// - using pt leads to some IR unsafety when we have two jets,
 148   ///   e.g. back-to-back, with the same pt. In that case, their ordering
 149   ///   in pt is random and can be affected by the addition of a
 150   ///   soft particle.  Hence, we highly recommand to keep this to
 151   ///   the default value i.e.  to use pt only for the purpose of
 152   ///   investigating the IR issue
 153   /// - using Et is safe but do not respect boost invariance
 154   /// - using mt solves the IR unsafety issues with the pt variable
 155   ///   for QCD jets but the IR unsafety remains for nack-to-back
 156   ///   jets of unstable narrow-width particles (e.g. Higgs).
 157   /// Therefore, keeping the default value is strongly advised.
 158   Esplit_merge_scale split_merge_scale;
 159 };
 160
 161
 162 // iterator types
 163 /// iterator definition for the jet candidates structure
 164 typedef std::multiset<siscone::Cjet,Csplit_merge_ptcomparison>::iterator cjet_iterator;
 165
 166 /// iterator definition for the jet structure
 167 typedef std::vector<siscone::Cjet>::iterator jet_iterator;
 168
 169
 170
 171 /**
 172  * \class Csplit_merge
 173  * Class used to split and merge jets.
 174  */
 175 class Csplit_merge{
 176  public:
 177   /// default ctor
 178   Csplit_merge();
 179
 180   /// default dtor
 181   ~Csplit_merge();
 182
 183
 184   //////////////////////////////
 185   // initialisation functions //
 186   //////////////////////////////
 187
 188   /**
 189    * initialisation function
 190    * \param _particles  list of particles
 191    * \param protocones  list of protocones (initial jet candidates)
 192    * \param R2          cone radius (squared)
 193    * \param ptmin       minimal pT allowed for jets
 194    * \return 0 on success, 1 on error
 195    */
 196   int init(std::vector<Cmomentum> &_particles, std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
 197
 198   /**
 199    * initialisation function for particle list
 200    * \param _particles  list of particles
 201    * \return 0 on success, 1 on error
 202    */
 203   int init_particles(std::vector<Cmomentum> &_particles);
 204
 205   /**
 206    * build initial list of left particles
 207    */
 208   int init_pleft();
 209
 210   /**
 211    * use a pt-dependent boundary for splitting
 212    * When called with true, the criterium for splitting two protojets
 213    * will be to compare D1^2/kt1^2 vs. D2^2/kt2^2, the (anti-)kt-weighted
 214    * distance instead of the plain distance D1^2 vs. D2^2.
 215    * This can be set in order to produce more circular hard jets,
 216    * with the same underlying philosophy as for the anti-kt algorithm.
 217    * We thus expect a behaviour closer to the IterativeCone one.
 218    * By default, we use the standard D1^2 vs. D2^2 comparison and this
 219    * function is not called.
 220    */
 221   inline int set_pt_weighted_splitting(bool _use_pt_weighted_splitting){
 222     use_pt_weighted_splitting = _use_pt_weighted_splitting;
 223     return 0;
 224   }
 225
 226   ////////////////////////
 227   // cleaning functions //
 228   ////////////////////////
 229
 230   /// partial clearance
 231   int partial_clear();
 232
 233   /// full clearance
 234   int full_clear();
 235
 236
 237   /////////////////////////////////
 238   // main parts of the algorithm //
 239   /////////////////////////////////
 240
 241   /**
 242    * build the list 'p_uncol_hard' from p_remain by clustering
 243    * collinear particles and removing particles softer than
 244    * stable_cone_soft_pt2_cutoff
 245    * note that thins in only used for stable-cone detection
 246    * so the parent_index field is unnecessary
 247    */
 248   int merge_collinear_and_remove_soft();
 249
 250   /**
 251    * add a list of protocones
 252    * \param protocones  list of protocones (initial jet candidates)
 253    * \param R2          cone radius (squared)
 254    * \param ptmin       minimal pT allowed for jets
 255    * \return 0 on success, 1 on error
 256    */
 257   int add_protocones(std::vector<Cmomentum> *protocones, double R2, double ptmin=0.0);
 258
 259   /**
 260    * really do the splitting and merging
 261    * At the end, the vector jets is filled with the jets found.
 262    * the 'contents' field of each jets contains the indices
 263    * of the particles included in that jet.
 264    * \param overlap_tshold  threshold for splitting/merging transition
 265    * \param ptmin           minimal pT allowed for jets
 266    * \return the number of jets is returned
 267    */
 268   int perform(double overlap_tshold, double ptmin=0.0);
 269
 270
 271   //////////////////////////////
 272   // save and debug functions //
 273   //////////////////////////////
 274
 275   /// save final jets
 276   /// \param flux   stream to save the jet contentss
 277   int save_contents(FILE *flux);
 278
 279   /// show jets/candidates status
 280   int show();
 281
 282   // particle information
 283   int n;                               ///< number of particles
 284   std::vector<Cmomentum> particles;    ///< list of particles
 285   std::vector<double> pt;              ///< list of particles' pt
 286   int n_left;                          ///< numer of particles that does not belong to any jet
 287   std::vector<Cmomentum> p_remain;     ///< list of particles remaining to deal with
 288   std::vector<Cmomentum> p_uncol_hard; ///< list of particles remaining with collinear clustering
 289   int n_pass;                          ///< index of the run
 290
 291   /// minimal difference in squared distance between a particle and
 292   /// two overlapping protojets when doing a split (useful when
 293   /// testing approx. collinear safety)
 294   double most_ambiguous_split;
 295
 296   // jets information
 297   std::vector<Cjet> jets;            ///< list of jets
 298
 299   // working entries
 300   int *indices;                      ///< maximal size array for indices works
 301   int idx_size;                      ///< number of elements in indices1
 302
 303   /// The following flag indicates that identical protocones
 304   /// are to be merged automatically each time around the split-merge
 305   /// loop and before anything else happens.
 306   ///
 307   /// This flag is only effective if ALLOW_MERGE_IDENTICAL_PROTOCONES
 308   /// is set in 'defines.h'
 309   /// Note that this lead to infrared-unsafety so it is disabled
 310   /// by default
 311   bool merge_identical_protocones;
 312
 313   /// member used for detailed comparisons of pt's
 314   Csplit_merge_ptcomparison ptcomparison;
 315
 316   /// stop split--merge when the SM_var of the hardest protojet
 317   /// is below this cut-off.
 318   /// This is not collinear-safe so you should not use this
 319   /// variable unless you really know what you are doing
 320   /// Note that the cut-off is set on the variable squared.
 321   double SM_var2_hardest_cut_off;
 322
 323   /// pt cutoff for the particles to put in p_uncol_hard
 324   /// this is meant to allow removing soft particles in the
 325   /// stable-cone search.
 326   double stable_cone_soft_pt2_cutoff;
 327
 328  private:
 329   /**
 330    * get the overlap between 2 jets
 331    * \param j1   first jet
 332    * \param j2   second jet
 333    * \param v    returned overlap^2 (determined by the choice of SM variable)
 334    * \return true if overlapping, false if disjoint
 335    */
 336   bool get_overlap(const Cjet &j1, const Cjet &j2, double *v);
 337
 338
 339   /**
 340    * split the two given jets.
 341    * during this procedure, the jets j1 & j2 are replaced
 342    * by 2 new jets. Common particles are associted to the
 343    * closest initial jet.
 344    * \param it_j1  iterator of the first jet in 'candidates'
 345    * \param it_j2  iterator of the second jet in 'candidates'
 346    * \param j1     first jet (Cjet instance)
 347    * \param j2     second jet (Cjet instance)
 348    * \return true on success, false on error
 349    */
 350   bool split(cjet_iterator &it_j1, cjet_iterator &it_j2);
 351
 352   /**
 353    * merge the two given jet.
 354    * during this procedure, the jets j1 & j2 are replaced
 355    * by 1 single jets containing both of them.
 356    * \param it_j1  iterator of the first jet in 'candidates'
 357    * \param it_j2  iterator of the second jet in 'candidates'
 358    * \return true on success, false on error
 359    */
 360   bool merge(cjet_iterator &it_j1, cjet_iterator &it_j2);
 361
 362   /**
 363    * Check whether or not a jet has to be inserted in the
 364    * list of protojets. If it has, set its sm_variable and
 365    * insert it to the list of protojets.
 366    * \param jet    jet to insert
 367    */
 368   bool insert(Cjet &jet);
 369
 370   /**
 371    * given a 4-momentum and its associated pT, return the
 372    * variable tht has to be used for SM
 373    * \param v          4 momentum of the protojet
 374    * \param pt_tilde   pt_tilde of the protojet
 375    */
 376   double get_sm_var2(Cmomentum &v, double &pt_tilde);
 377
 378   // jet information
 379   /// list of jet candidates
 380   std::auto_ptr<std::multiset<Cjet,Csplit_merge_ptcomparison> > candidates;
 381
 382   /// minimal pt2
 383   double pt_min2;
 384
 385   /**
 386    * do we have or not to use the pt-weighted splitting
 387    * (see description for set_pt_weighted_splitting)
 388    * This will be false by default
 389    */
 390   bool use_pt_weighted_splitting;
 391
 392 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
 393   /// checkxor for the candidates (to avoid having twice the same contents)
 394   std::set<Creference> cand_refs;
 395 #endif
 396 };
 397
 398 }
 399
 400
 401 #endif