← Back to branch summary

~maddevelopers/mg5amcnlo/2.5.4_run_py8_at_evtgen

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/fjcore.cc

  • Committer: olivier Mattelaer
  • Date: 2016-05-12 11:00:18 UTC
  • mfrom: (262.1.150 2.3.4)
  • Revision ID: olivier.mattelaer@uclouvain.be-20160512110018-sevb79f0wm4g8mpp
pass to 2.4.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/fjcore.cc
1
 
// fjcore -- extracted from FastJet v3.0.5 (http://fastjet.fr)
 
1
// fjcore -- extracted from FastJet v3.1.3 (http://fastjet.fr)
2
2
//
3
3
// fjcore constitutes a digest of the main FastJet functionality.
4
4
// The files fjcore.hh and fjcore.cc are meant to provide easy access to these 
19
19
//
20
20
//   - access to all native pp and ee algorithms, kt, anti-kt, C/A.
21
21
//     For C/A, the NlnN method is available, while anti-kt and kt
22
 
//     are limited to the N^2 one (still the fastest for N < 20k particles)
 
22
//     are limited to the N^2 one (still the fastest for N < 100k particles)
23
23
//   - access to selectors, for implementing cuts and selections
24
24
//   - access to all functionalities related to pseudojets (e.g. a jet's
25
25
//     structure or user-defined information)
42
42
// Like FastJet, fjcore is released under the terms of the GNU General Public
43
43
// License version 2 (GPLv2). If you use this code as part of work towards a
44
44
// scientific publication, whether directly or contained within another program
45
 
// (e.g. Delphes, SpartyJet, Rivet, LHC collaboration software frameworks, 
 
45
// (e.g. Delphes, MadGraph, SpartyJet, Rivet, LHC collaboration software frameworks, 
46
46
// etc.), you should include a citation to
47
47
// 
48
48
//   EPJC72(2012)1896 [arXiv:1111.6097] (FastJet User Manual)
49
49
//   and, optionally, Phys.Lett.B641 (2006) 57 [arXiv:hep-ph/0512210]
50
50
//
51
 
// Copyright (c) 2005-2013, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
 
51
// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
52
52
//
53
53
//----------------------------------------------------------------------
54
54
// This file is part of FastJet.
171
171
    if (second >= twopi) second -= twopi;
172
172
  }
173
173
};
174
 
class DnnError {
 
174
class DnnError : public Error {
175
175
public:
176
 
  DnnError() {;};
177
 
  DnnError(const std::string & message_in) {
178
 
    _message = message_in; std::cerr << message_in << std::endl;};
179
 
  std::string message() const {return _message;};
180
 
private:
181
 
  std::string _message;
 
176
  DnnError(const std::string & message_in) : Error(message_in) {}
182
177
};
183
178
class DynamicNearestNeighbours {
184
179
public:
185
 
  virtual int NearestNeighbourIndex(const int & ii) const = 0;
186
 
  virtual double NearestNeighbourDistance(const int & ii) const = 0;
187
 
  virtual bool Valid(const int & index) const = 0;
 
180
  virtual int NearestNeighbourIndex(const int ii) const = 0;
 
181
  virtual double NearestNeighbourDistance(const int ii) const = 0;
 
182
  virtual bool Valid(const int index) const = 0;
188
183
  virtual void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
189
184
                          const std::vector<EtaPhi> & points_to_add,
190
185
                          std::vector<int> & indices_added,
191
186
                          std::vector<int> & indices_of_updated_neighbours) = 0;
192
 
  inline void RemovePoint (const int & index,
 
187
  inline void RemovePoint (const int index,
193
188
                           std::vector<int> & indices_of_updated_neighbours) {
194
189
    std::vector<int> indices_added;
195
190
    std::vector<EtaPhi> points_to_add;
199
194
                       indices_of_updated_neighbours
200
195
                       );};
201
196
  inline void RemoveCombinedAddCombination(
202
 
                        const int & index1, const int & index2,
 
197
                        const int index1, const int index2,
203
198
                        const EtaPhi & newpoint,
204
199
                        int & index3,
205
200
                        std::vector<int> & indices_of_updated_neighbours) {
290
285
}
291
286
template<class T> class SearchTree<T>::circulator{
292
287
public:
293
 
  template<class U> friend class SearchTree<U>::const_circulator;
 
288
  friend class SearchTree<T>::const_circulator;
294
289
  friend class SearchTree<T>;
295
290
  circulator() : _node(NULL) {}
296
291
  circulator(Node * node) : _node(node) {}
644
639
class MinHeap {
645
640
public:
646
641
  MinHeap (const std::vector<double> & values, unsigned int max_size) :
647
 
    _heap(max_size) {_initialise(values);};
 
642
    _heap(max_size) {initialise(values);}
 
643
  MinHeap (unsigned int max_size) : _heap(max_size) {}
648
644
  MinHeap (const std::vector<double> & values) :
649
 
    _heap(values.size()) {_initialise(values);};
 
645
    _heap(values.size()) {initialise(values);}
 
646
  void initialise(const std::vector<double> & values);
650
647
  inline unsigned int minloc() const {
651
 
    return (_heap[0].minloc) - &(_heap[0]);};
652
 
  inline double       minval() const {return _heap[0].minloc->value;};
653
 
  inline double operator[](int i) const {return _heap[i].value;};
 
648
    return (_heap[0].minloc) - &(_heap[0]);}
 
649
  inline double       minval() const {return _heap[0].minloc->value;}
 
650
  inline double operator[](int i) const {return _heap[i].value;}
654
651
  void remove(unsigned int loc) {
655
652
    update(loc,std::numeric_limits<double>::max());};
656
653
  void update(unsigned int, double);
660
657
    ValueLoc * minloc;
661
658
  };
662
659
  std::vector<ValueLoc> _heap;
663
 
  void _initialise(const std::vector<double> & values);
664
660
};
665
661
FJCORE_END_NAMESPACE
666
662
#endif // __FJCORE_MINHEAP__HH__
818
814
}
819
815
FJCORE_END_NAMESPACE
820
816
#endif // __FJCORE_CLOSESTPAIR2D__HH__
 
817
#ifndef __FJCORE_LAZYTILING9ALT_HH__
 
818
#define __FJCORE_LAZYTILING9ALT_HH__
 
819
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
820
const double tile_edge_security_margin=1.0e-7;
 
821
class TiledJet {
 
822
public:
 
823
  double     eta, phi, kt2, NN_dist;
 
824
  TiledJet * NN, *previous, * next; 
 
825
  int        _jets_index, tile_index;
 
826
  bool _minheap_update_needed;
 
827
  inline void label_minheap_update_needed() {_minheap_update_needed = true;}
 
828
  inline void label_minheap_update_done()   {_minheap_update_needed = false;}
 
829
  inline bool minheap_update_needed() const {return _minheap_update_needed;}
 
830
};
 
831
const int n_tile_neighbours = 9;
 
832
class Tile {
 
833
public:
 
834
  typedef double (Tile::*DistToTileFn)(const TiledJet*) const;
 
835
  typedef std::pair<Tile *, DistToTileFn> TileFnPair;
 
836
  TileFnPair begin_tiles[n_tile_neighbours]; 
 
837
  TileFnPair *  surrounding_tiles; 
 
838
  TileFnPair *  RH_tiles;  
 
839
  TileFnPair *  end_tiles; 
 
840
  TiledJet * head;    
 
841
  bool     tagged;    
 
842
  bool     use_periodic_delta_phi;
 
843
  double max_NN_dist;
 
844
  double eta_min, eta_max, phi_min, phi_max;
 
845
  double distance_to_centre(const TiledJet *) const {return 0;}
 
846
  double distance_to_left(const TiledJet * jet) const {
 
847
    double deta = jet->eta - eta_min;
 
848
    return deta*deta;
 
849
  }
 
850
  double distance_to_right(const TiledJet * jet) const {
 
851
    double deta = jet->eta - eta_max;
 
852
    return deta*deta;
 
853
  }
 
854
  double distance_to_bottom(const TiledJet * jet) const {
 
855
    double dphi = jet->phi - phi_min;
 
856
    return dphi*dphi;
 
857
  }
 
858
  double distance_to_top(const TiledJet * jet) const {
 
859
    double dphi = jet->phi - phi_max;
 
860
    return dphi*dphi;
 
861
  }
 
862
  double distance_to_left_top(const TiledJet * jet) const {
 
863
    double deta = jet->eta - eta_min;
 
864
    double dphi = jet->phi - phi_max;
 
865
    return deta*deta + dphi*dphi;
 
866
  }
 
867
  double distance_to_left_bottom(const TiledJet * jet) const {
 
868
    double deta = jet->eta - eta_min;
 
869
    double dphi = jet->phi - phi_min;
 
870
    return deta*deta + dphi*dphi;
 
871
  }
 
872
  double distance_to_right_top(const TiledJet * jet) const {
 
873
    double deta = jet->eta - eta_max;
 
874
    double dphi = jet->phi - phi_max;
 
875
    return deta*deta + dphi*dphi;
 
876
  }
 
877
  double distance_to_right_bottom(const TiledJet * jet) const {
 
878
    double deta = jet->eta - eta_max;
 
879
    double dphi = jet->phi - phi_min;
 
880
    return deta*deta + dphi*dphi;
 
881
  }
 
882
};
 
883
class LazyTiling9Alt {
 
884
public:
 
885
  LazyTiling9Alt(ClusterSequence & cs);
 
886
  void run();
 
887
protected:
 
888
  ClusterSequence & _cs;
 
889
  const std::vector<PseudoJet> & _jets;
 
890
  std::vector<Tile> _tiles;
 
891
  double _Rparam, _R2, _invR2;
 
892
  double _tiles_eta_min, _tiles_eta_max;
 
893
  double _tile_size_eta, _tile_size_phi;
 
894
  double _tile_half_size_eta, _tile_half_size_phi;
 
895
  int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
 
896
  std::vector<TiledJet *> _jets_for_minheap;
 
897
  void _initialise_tiles();
 
898
  inline int _tile_index (int ieta, int iphi) const {
 
899
    return (ieta-_tiles_ieta_min)*_n_tiles_phi
 
900
                  + (iphi+_n_tiles_phi) % _n_tiles_phi;
 
901
  }
 
902
  void  _bj_remove_from_tiles(TiledJet * const jet);
 
903
  int _tile_index(const double eta, const double phi) const;
 
904
  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
 
905
  void _print_tiles(TiledJet * briefjets ) const;
 
906
  void _add_neighbours_to_tile_union(const int tile_index, 
 
907
                 std::vector<int> & tile_union, int & n_near_tiles) const;
 
908
  void _add_untagged_neighbours_to_tile_union(const int tile_index, 
 
909
                 std::vector<int> & tile_union, int & n_near_tiles);
 
910
  void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet, 
 
911
                 std::vector<int> & tile_union, int & n_near_tiles);
 
912
  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
 
913
  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
 
914
  template <class J> double _bj_diJ(const J * const jet) const {
 
915
    double kt2 = jet->kt2;
 
916
    if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
 
917
    return jet->NN_dist * kt2;
 
918
  }
 
919
  template <class J> inline void _bj_set_jetinfo(
 
920
                            J * const jetA, const int _jets_index) const {
 
921
    jetA->eta  = _jets[_jets_index].rap();
 
922
    jetA->phi  = _jets[_jets_index].phi_02pi();
 
923
    jetA->kt2  = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
 
924
    jetA->_jets_index = _jets_index;
 
925
    jetA->NN_dist = _R2;
 
926
    jetA->NN      = NULL;
 
927
  }
 
928
  template <class J> inline double _bj_dist(
 
929
                const J * const jetA, const J * const jetB) const {
 
930
    double dphi = std::abs(jetA->phi - jetB->phi);
 
931
    double deta = (jetA->eta - jetB->eta);
 
932
    if (dphi > pi) {dphi = twopi - dphi;}
 
933
    return dphi*dphi + deta*deta;
 
934
  }
 
935
  template <class J> inline double _bj_dist_not_periodic(
 
936
                const J * const jetA, const J * const jetB) const {
 
937
    double dphi = jetA->phi - jetB->phi;
 
938
    double deta = (jetA->eta - jetB->eta);
 
939
    return dphi*dphi + deta*deta;
 
940
  }
 
941
};
 
942
FJCORE_END_NAMESPACE
 
943
#endif // __FJCORE_LAZYTILING9ALT_HH__
 
944
#ifndef __FJCORE_LAZYTILING9_HH__
 
945
#define __FJCORE_LAZYTILING9_HH__
 
946
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
947
template<int NN>
 
948
class Tile2Base {
 
949
public:
 
950
  Tile2Base *   begin_tiles[NN]; 
 
951
  Tile2Base **  surrounding_tiles; 
 
952
  Tile2Base **  RH_tiles;  
 
953
  Tile2Base **  end_tiles; 
 
954
  TiledJet * head;    
 
955
  bool     tagged;    
 
956
  bool     use_periodic_delta_phi;
 
957
  double max_NN_dist;
 
958
  double eta_centre, phi_centre;
 
959
  int jet_count() const {
 
960
    int count = 0;
 
961
    const TiledJet * jet = head;
 
962
    while (jet != 0) {
 
963
      count++;
 
964
      jet = jet->next;
 
965
    }
 
966
    return count;
 
967
  }
 
968
};
 
969
typedef Tile2Base<9> Tile2;
 
970
class LazyTiling9 {
 
971
public:
 
972
  LazyTiling9(ClusterSequence & cs);
 
973
  void run();
 
974
protected:
 
975
  ClusterSequence & _cs;
 
976
  const std::vector<PseudoJet> & _jets;
 
977
  std::vector<Tile2> _tiles;
 
978
#ifdef INSTRUMENT2
 
979
  int _ncall; // GPS tmp
 
980
  int _ncall_dtt; // GPS tmp
 
981
#endif // INSTRUMENT2
 
982
  double _Rparam, _R2, _invR2;
 
983
  double _tiles_eta_min, _tiles_eta_max;
 
984
  double _tile_size_eta, _tile_size_phi;
 
985
  double _tile_half_size_eta, _tile_half_size_phi;
 
986
  int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
 
987
  std::vector<TiledJet *> _jets_for_minheap;
 
988
  void _initialise_tiles();
 
989
  inline int _tile_index (int ieta, int iphi) const {
 
990
    return (ieta-_tiles_ieta_min)*_n_tiles_phi
 
991
                  + (iphi+_n_tiles_phi) % _n_tiles_phi;
 
992
  }
 
993
  void  _bj_remove_from_tiles(TiledJet * const jet);
 
994
  int _tile_index(const double eta, const double phi) const;
 
995
  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
 
996
  void _print_tiles(TiledJet * briefjets ) const;
 
997
  void _add_neighbours_to_tile_union(const int tile_index, 
 
998
                 std::vector<int> & tile_union, int & n_near_tiles) const;
 
999
  void _add_untagged_neighbours_to_tile_union(const int tile_index, 
 
1000
                 std::vector<int> & tile_union, int & n_near_tiles);
 
1001
  void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet, 
 
1002
                 std::vector<int> & tile_union, int & n_near_tiles);
 
1003
  double _distance_to_tile(const TiledJet * bj, const Tile2 *) 
 
1004
#ifdef INSTRUMENT2
 
1005
    ;
 
1006
#else
 
1007
    const;
 
1008
#endif 
 
1009
  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
 
1010
  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
 
1011
  template <class J> double _bj_diJ(const J * const jet) const {
 
1012
    double kt2 = jet->kt2;
 
1013
    if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
 
1014
    return jet->NN_dist * kt2;
 
1015
  }
 
1016
  template <class J> inline void _bj_set_jetinfo(
 
1017
                            J * const jetA, const int _jets_index) const {
 
1018
    jetA->eta  = _jets[_jets_index].rap();
 
1019
    jetA->phi  = _jets[_jets_index].phi_02pi();
 
1020
    jetA->kt2  = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
 
1021
    jetA->_jets_index = _jets_index;
 
1022
    jetA->NN_dist = _R2;
 
1023
    jetA->NN      = NULL;
 
1024
  }
 
1025
  template <class J> inline double _bj_dist(
 
1026
                const J * const jetA, const J * const jetB) 
 
1027
#ifdef INSTRUMENT2
 
1028
    {
 
1029
    _ncall++; // GPS tmp
 
1030
#else
 
1031
    const {
 
1032
#endif 
 
1033
    double dphi = std::abs(jetA->phi - jetB->phi);
 
1034
    double deta = (jetA->eta - jetB->eta);
 
1035
    if (dphi > pi) {dphi = twopi - dphi;}
 
1036
    return dphi*dphi + deta*deta;
 
1037
  }
 
1038
  template <class J> inline double _bj_dist_not_periodic(
 
1039
                const J * const jetA, const J * const jetB)
 
1040
#ifdef INSTRUMENT2
 
1041
    {
 
1042
    _ncall++; // GPS tmp
 
1043
#else
 
1044
    const {
 
1045
#endif 
 
1046
    double dphi = jetA->phi - jetB->phi;
 
1047
    double deta = (jetA->eta - jetB->eta);
 
1048
    return dphi*dphi + deta*deta;
 
1049
  }
 
1050
};
 
1051
FJCORE_END_NAMESPACE
 
1052
#endif // __FJCORE_LAZYTILING9_HH__
 
1053
#ifndef __FJCORE_LAZYTILING25_HH__
 
1054
#define __FJCORE_LAZYTILING25_HH__
 
1055
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
1056
typedef Tile2Base<25> Tile25;
 
1057
class LazyTiling25 {
 
1058
public:
 
1059
  LazyTiling25(ClusterSequence & cs);
 
1060
  void run();
 
1061
protected:
 
1062
  ClusterSequence & _cs;
 
1063
  const std::vector<PseudoJet> & _jets;
 
1064
  std::vector<Tile25> _tiles;
 
1065
#ifdef INSTRUMENT2
 
1066
  int _ncall; // GPS tmp
 
1067
  int _ncall_dtt; // GPS tmp
 
1068
#endif // INSTRUMENT2
 
1069
  double _Rparam, _R2, _invR2;
 
1070
  double _tiles_eta_min, _tiles_eta_max;
 
1071
  double _tile_size_eta, _tile_size_phi;
 
1072
  double _tile_half_size_eta, _tile_half_size_phi;
 
1073
  int    _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
 
1074
  std::vector<TiledJet *> _jets_for_minheap;
 
1075
  void _initialise_tiles();
 
1076
  inline int _tile_index (int ieta, int iphi) const {
 
1077
    return (ieta-_tiles_ieta_min)*_n_tiles_phi
 
1078
                  + (iphi+_n_tiles_phi) % _n_tiles_phi;
 
1079
  }
 
1080
  void  _bj_remove_from_tiles(TiledJet * const jet);
 
1081
  int _tile_index(const double eta, const double phi) const;
 
1082
  void _tj_set_jetinfo(TiledJet * const jet, const int _jets_index);
 
1083
  void _print_tiles(TiledJet * briefjets ) const;
 
1084
  void _add_neighbours_to_tile_union(const int tile_index, 
 
1085
                 std::vector<int> & tile_union, int & n_near_tiles) const;
 
1086
  void _add_untagged_neighbours_to_tile_union(const int tile_index, 
 
1087
                 std::vector<int> & tile_union, int & n_near_tiles);
 
1088
  void _add_untagged_neighbours_to_tile_union_using_max_info(const TiledJet * const jet, 
 
1089
                 std::vector<int> & tile_union, int & n_near_tiles);
 
1090
  double _distance_to_tile(const TiledJet * bj, const Tile25 *) 
 
1091
#ifdef INSTRUMENT2
 
1092
    ;
 
1093
#else
 
1094
    const;
 
1095
#endif 
 
1096
  void _update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
 
1097
  void _set_NN(TiledJet * jetI, std::vector<TiledJet *> & jets_for_minheap);
 
1098
  template <class J> double _bj_diJ(const J * const jet) const {
 
1099
    double kt2 = jet->kt2;
 
1100
    if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
 
1101
    return jet->NN_dist * kt2;
 
1102
  }
 
1103
  template <class J> inline void _bj_set_jetinfo(
 
1104
                            J * const jetA, const int _jets_index) const {
 
1105
    jetA->eta  = _jets[_jets_index].rap();
 
1106
    jetA->phi  = _jets[_jets_index].phi_02pi();
 
1107
    jetA->kt2  = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
 
1108
    jetA->_jets_index = _jets_index;
 
1109
    jetA->NN_dist = _R2;
 
1110
    jetA->NN      = NULL;
 
1111
  }
 
1112
  template <class J> inline double _bj_dist(
 
1113
                const J * const jetA, const J * const jetB) 
 
1114
#ifdef INSTRUMENT2
 
1115
    {
 
1116
    _ncall++; // GPS tmp
 
1117
#else
 
1118
    const {
 
1119
#endif 
 
1120
    double dphi = std::abs(jetA->phi - jetB->phi);
 
1121
    double deta = (jetA->eta - jetB->eta);
 
1122
    if (dphi > pi) {dphi = twopi - dphi;}
 
1123
    return dphi*dphi + deta*deta;
 
1124
  }
 
1125
  template <class J> inline double _bj_dist_not_periodic(
 
1126
                const J * const jetA, const J * const jetB)
 
1127
#ifdef INSTRUMENT2
 
1128
    {
 
1129
    _ncall++; // GPS tmp
 
1130
#else
 
1131
    const {
 
1132
#endif 
 
1133
    double dphi = jetA->phi - jetB->phi;
 
1134
    double deta = (jetA->eta - jetB->eta);
 
1135
    return dphi*dphi + deta*deta;
 
1136
  }
 
1137
};
 
1138
FJCORE_END_NAMESPACE
 
1139
#endif // __FJCORE_LAZYTILING25_HH__
 
1140
#ifndef __FJCORE_TILINGEXTENT_HH__
 
1141
#define __FJCORE_TILINGEXTENT_HH__
 
1142
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
1143
class TilingExtent {
 
1144
public:
 
1145
  TilingExtent(ClusterSequence & cs);
 
1146
  double minrap() const {return _minrap;}
 
1147
  double maxrap() const {return _maxrap;}
 
1148
  double sum_of_binned_squared_multiplicity() const {return _cumul2;}
 
1149
private:
 
1150
  double _minrap, _maxrap, _cumul2;
 
1151
  void _determine_rapidity_extent(const std::vector<PseudoJet> & particles);
 
1152
};
 
1153
FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
 
1154
#endif // __FJCORE_TILINGEXTENT_HH__
821
1155
#include<limits>
822
1156
#include<iostream>
823
1157
#include<iomanip>
824
1158
#include<algorithm>
825
1159
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
826
 
const unsigned int huge_unsigned = 4294967295U;
827
1160
const unsigned int twopow31      = 2147483648U;
828
1161
using namespace std;
829
1162
void ClosestPair2D::_point2shuffle(Point & point, Shuffle & shuffle, 
1120
1453
    throw Error("A ClusterSequence cannot be created with an uninitialised JetDefinition");
1121
1454
  }
1122
1455
  if (_strategy == Best) {
 
1456
    _strategy = _best_strategy();
 
1457
#ifdef __FJCORE_DROP_CGAL
 
1458
    if (_strategy == NlnN) _strategy = N2MHTLazy25;
 
1459
#endif  // __FJCORE_DROP_CGAL
 
1460
  } else if (_strategy == BestFJ30) {
1123
1461
    int N = _jets.size();
1124
1462
    if (min(1.0,max(0.1,_Rparam)*3.3)*N <= 30) {
1125
1463
      _strategy = N2Plain;
1163
1501
    this->_faster_tiled_N2_cluster();
1164
1502
  } else if (_strategy == N2MinHeapTiled) {
1165
1503
    this->_minheap_faster_tiled_N2_cluster();
 
1504
  } else if (_strategy == N2MHTLazy9Alt) {
 
1505
    _plugin_activated = true;
 
1506
    LazyTiling9Alt tiling(*this);
 
1507
    tiling.run();
 
1508
    _plugin_activated = false;
 
1509
  } else if (_strategy == N2MHTLazy25) {
 
1510
    _plugin_activated = true;
 
1511
    LazyTiling25 tiling(*this);
 
1512
    tiling.run();
 
1513
    _plugin_activated = false;
 
1514
  } else if (_strategy == N2MHTLazy9) {
 
1515
    _plugin_activated = true;
 
1516
    LazyTiling9 tiling(*this);
 
1517
    tiling.run();
 
1518
    _plugin_activated = false;
 
1519
  } else if (_strategy == N2MHTLazy9AntiKtSeparateGhosts) {
 
1520
    throw Error("N2MHTLazy9AntiKtSeparateGhosts strategy not supported with FJCORE");
1166
1521
  } else if (_strategy == NlnN) {
1167
1522
    this->_delaunay_cluster();
1168
1523
  } else if (_strategy == NlnNCam) {
1184
1539
  }
1185
1540
}
1186
1541
bool ClusterSequence::_first_time = true;
1187
 
int ClusterSequence::_n_exclusive_warnings = 0;
 
1542
LimitedWarning ClusterSequence::_exclusive_warnings;
1188
1543
string fastjet_version_string() {
1189
1544
  return "FastJet version "+string(fastjet_version)+" [fjcore]";
1190
1545
}
1266
1621
    strategy = "N2MinHeapTiled"; break;
1267
1622
  case N2PoorTiled:
1268
1623
    strategy = "N2PoorTiled"; break;
 
1624
  case N2MHTLazy9:
 
1625
    strategy = "N2MHTLazy9"; break;
 
1626
  case N2MHTLazy9Alt:
 
1627
    strategy = "N2MHTLazy9Alt"; break;
 
1628
  case N2MHTLazy25:
 
1629
    strategy = "N2MHTLazy25"; break;
 
1630
  case N2MHTLazy9AntiKtSeparateGhosts:
 
1631
    strategy = "N2MHTLazy9AntiKtSeparateGhosts"; break;
1269
1632
  case N3Dumb:
1270
1633
    strategy = "N3Dumb"; break;
1271
1634
  case NlnNCam4pi:
1301
1664
    } else {return 1.0;}
1302
1665
  } else {throw Error("Unrecognised jet algorithm");}
1303
1666
}
 
1667
Strategy ClusterSequence::_best_strategy() const {
 
1668
  int N = _jets.size();
 
1669
  double bounded_R = max(_Rparam, 0.1);
 
1670
  if (N <= 30 || N <= 39.0/(bounded_R + 0.6)) {
 
1671
    return N2Plain;
 
1672
  } 
 
1673
  const static _Parabola N_Tiled_to_MHT_lowR             (-45.4947,54.3528,44.6283);
 
1674
  const static _Parabola L_MHT_to_MHTLazy9_lowR          (0.677807,-1.05006,10.6994);
 
1675
  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_lowR(0.169967,-0.512589,12.1572);
 
1676
  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_lowR (0.16237,-0.484612,12.3373);
 
1677
  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_lowR = L_MHTLazy9_to_MHTLazy25_kt_lowR;
 
1678
  const static _Parabola L_MHTLazy25_to_NlnN_akt_lowR    (0.0472051,-0.22043,15.9196);
 
1679
  const static _Parabola L_MHTLazy25_to_NlnN_kt_lowR     (0.118609,-0.326811,14.8287);
 
1680
  const static _Parabola L_MHTLazy25_to_NlnN_cam_lowR    (0.10119,-0.295748,14.3924);
 
1681
  const static _Line     L_Tiled_to_MHTLazy9_medR         (-1.31304,7.29621);
 
1682
  const static _Parabola L_MHTLazy9_to_MHTLazy25_akt_medR = L_MHTLazy9_to_MHTLazy25_akt_lowR;
 
1683
  const static _Parabola L_MHTLazy9_to_MHTLazy25_kt_medR  = L_MHTLazy9_to_MHTLazy25_kt_lowR;
 
1684
  const static _Parabola L_MHTLazy9_to_MHTLazy25_cam_medR = L_MHTLazy9_to_MHTLazy25_cam_lowR;
 
1685
  const static _Parabola L_MHTLazy25_to_NlnN_akt_medR     = L_MHTLazy25_to_NlnN_akt_lowR;
 
1686
  const static _Parabola L_MHTLazy25_to_NlnN_kt_medR      = L_MHTLazy25_to_NlnN_kt_lowR;
 
1687
  const static _Parabola L_MHTLazy25_to_NlnN_cam_medR     = L_MHTLazy25_to_NlnN_cam_lowR;
 
1688
  const static double    N_Plain_to_MHTLazy9_largeR         = 75;
 
1689
  const static double    N_MHTLazy9_to_MHTLazy25_akt_largeR = 700;
 
1690
  const static double    N_MHTLazy9_to_MHTLazy25_kt_largeR  = 1000;
 
1691
  const static double    N_MHTLazy9_to_MHTLazy25_cam_largeR = 1000;
 
1692
  const static double    N_MHTLazy25_to_NlnN_akt_largeR     = 100000;
 
1693
  const static double    N_MHTLazy25_to_NlnN_kt_largeR      = 40000;
 
1694
  const static double    N_MHTLazy25_to_NlnN_cam_largeR     = 15000;
 
1695
  JetAlgorithm jet_algorithm;
 
1696
  if (_jet_algorithm == genkt_algorithm) {
 
1697
    double p   = jet_def().extra_param();
 
1698
    if (p < 0.0) jet_algorithm = antikt_algorithm;
 
1699
    else         jet_algorithm =     kt_algorithm;
 
1700
  } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
 
1701
    jet_algorithm = kt_algorithm;
 
1702
  } else {
 
1703
    jet_algorithm = _jet_algorithm;
 
1704
  }
 
1705
  if (bounded_R < 0.65) {
 
1706
    if          (N    < N_Tiled_to_MHT_lowR(bounded_R))              return N2Tiled;
 
1707
    double logN = log(double(N));
 
1708
    if          (logN < L_MHT_to_MHTLazy9_lowR(bounded_R))           return N2MinHeapTiled;
 
1709
    else {
 
1710
      if (jet_algorithm == antikt_algorithm){
 
1711
        if      (logN < L_MHTLazy9_to_MHTLazy25_akt_lowR(bounded_R)) return N2MHTLazy9;
 
1712
        else if (logN < L_MHTLazy25_to_NlnN_akt_lowR(bounded_R))     return N2MHTLazy25;
 
1713
        else                                                         return NlnN;
 
1714
      } else if (jet_algorithm == kt_algorithm){
 
1715
        if      (logN < L_MHTLazy9_to_MHTLazy25_kt_lowR(bounded_R))  return N2MHTLazy9;
 
1716
        else if (logN < L_MHTLazy25_to_NlnN_kt_lowR(bounded_R))      return N2MHTLazy25;
 
1717
        else                                                         return NlnN;
 
1718
      } else if (jet_algorithm == cambridge_algorithm)  {
 
1719
        if      (logN < L_MHTLazy9_to_MHTLazy25_cam_lowR(bounded_R)) return N2MHTLazy9;
 
1720
        else if (logN < L_MHTLazy25_to_NlnN_cam_lowR(bounded_R))     return N2MHTLazy25;
 
1721
        else                                                         return NlnNCam;
 
1722
      }
 
1723
    }
 
1724
  } else if (bounded_R < 0.5*pi) {
 
1725
    double logN = log(double(N));
 
1726
    if      (logN < L_Tiled_to_MHTLazy9_medR(bounded_R))             return N2Tiled;
 
1727
    else {
 
1728
      if (jet_algorithm == antikt_algorithm){
 
1729
        if      (logN < L_MHTLazy9_to_MHTLazy25_akt_medR(bounded_R)) return N2MHTLazy9;
 
1730
        else if (logN < L_MHTLazy25_to_NlnN_akt_medR(bounded_R))     return N2MHTLazy25;
 
1731
        else                                                         return NlnN;
 
1732
      } else if (jet_algorithm == kt_algorithm){
 
1733
        if      (logN < L_MHTLazy9_to_MHTLazy25_kt_medR(bounded_R))  return N2MHTLazy9;
 
1734
        else if (logN < L_MHTLazy25_to_NlnN_kt_medR(bounded_R))      return N2MHTLazy25;
 
1735
        else                                                         return NlnN;
 
1736
      } else if (jet_algorithm == cambridge_algorithm)  {
 
1737
        if      (logN < L_MHTLazy9_to_MHTLazy25_cam_medR(bounded_R)) return N2MHTLazy9;
 
1738
        else if (logN < L_MHTLazy25_to_NlnN_cam_medR(bounded_R))     return N2MHTLazy25;
 
1739
        else                                                         return NlnNCam;
 
1740
      }
 
1741
    }
 
1742
  } else {
 
1743
    if      (N    < N_Plain_to_MHTLazy9_largeR)                      return N2Plain;
 
1744
    else {
 
1745
      if (jet_algorithm == antikt_algorithm){
 
1746
        if      (N < N_MHTLazy9_to_MHTLazy25_akt_largeR)             return N2MHTLazy9;
 
1747
        else if (N < N_MHTLazy25_to_NlnN_akt_largeR)                 return N2MHTLazy25;
 
1748
        else                                                         return NlnN;
 
1749
      } else if (jet_algorithm == kt_algorithm){
 
1750
        if      (N < N_MHTLazy9_to_MHTLazy25_kt_largeR)              return N2MHTLazy9;
 
1751
        else if (N < N_MHTLazy25_to_NlnN_kt_largeR)                  return N2MHTLazy25;
 
1752
        else                                                         return NlnN;
 
1753
      } else if (jet_algorithm == cambridge_algorithm)  {
 
1754
        if      (N < N_MHTLazy9_to_MHTLazy25_cam_largeR)             return N2MHTLazy9;
 
1755
        else if (N < N_MHTLazy25_to_NlnN_cam_largeR)                 return N2MHTLazy25;
 
1756
        else                                                         return NlnNCam;
 
1757
      }
 
1758
    }
 
1759
  }
 
1760
  bool code_should_never_reach_here = false;
 
1761
  assert(code_should_never_reach_here); 
 
1762
  return N2MHTLazy9;
 
1763
}
1304
1764
void ClusterSequence::transfer_from_sequence(const ClusterSequence & from_seq,
1305
1765
                                             const FunctionOfPseudoJet<PseudoJet> * action_on_jets){
1306
1766
  if (will_delete_self_when_unused()) 
1342
1802
  _jets[newjet_k].set_cluster_hist_index(tmp_index);
1343
1803
  _set_structure_shared_ptr(_jets[newjet_k]);
1344
1804
}
1345
 
vector<PseudoJet> ClusterSequence::inclusive_jets (const double & ptmin) const{
 
1805
vector<PseudoJet> ClusterSequence::inclusive_jets (const double ptmin) const{
1346
1806
  double dcut = ptmin*ptmin;
1347
1807
  int i = _history.size() - 1; // last jet
1348
1808
  vector<PseudoJet> jets_local;
1380
1840
  } else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
1381
1841
  return jets_local;
1382
1842
}
1383
 
int ClusterSequence::n_exclusive_jets (const double & dcut) const {
 
1843
int ClusterSequence::n_exclusive_jets (const double dcut) const {
1384
1844
  int i = _history.size() - 1; // last jet
1385
1845
  while (i >= 0) {
1386
1846
    if (_history[i].max_dij_so_far <= dcut) {break;}
1390
1850
  int njets = 2*_initial_n - stop_point;
1391
1851
  return njets;
1392
1852
}
1393
 
vector<PseudoJet> ClusterSequence::exclusive_jets (const double & dcut) const {
 
1853
vector<PseudoJet> ClusterSequence::exclusive_jets (const double dcut) const {
1394
1854
  int njets = n_exclusive_jets(dcut);
1395
1855
  return exclusive_jets(njets);
1396
1856
}
1397
 
vector<PseudoJet> ClusterSequence::exclusive_jets (const int & njets) const {
 
1857
vector<PseudoJet> ClusterSequence::exclusive_jets (const int njets) const {
1398
1858
  if (njets > _initial_n) {
1399
1859
    ostringstream err;
1400
1860
    err << "Requested " << njets << " exclusive jets, but there were only " 
1403
1863
  }
1404
1864
  return exclusive_jets_up_to(njets);
1405
1865
}
1406
 
vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int & njets) const {
 
1866
vector<PseudoJet> ClusterSequence::exclusive_jets_up_to (const int njets) const {
1407
1867
  if (( _jet_def.jet_algorithm() != kt_algorithm) &&
1408
1868
      ( _jet_def.jet_algorithm() != cambridge_algorithm) &&
1409
1869
      ( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
1411
1871
        (_jet_def.jet_algorithm() != ee_genkt_algorithm)) || 
1412
1872
       (_jet_def.extra_param() <0)) &&
1413
1873
      ((_jet_def.jet_algorithm() != plugin_algorithm) ||
1414
 
       (!_jet_def.plugin()->exclusive_sequence_meaningful())) &&
1415
 
      (_n_exclusive_warnings < 5)) {
1416
 
    _n_exclusive_warnings++;
1417
 
    cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
 
1874
       (!_jet_def.plugin()->exclusive_sequence_meaningful()))) {
 
1875
    _exclusive_warnings.warn("dcut and exclusive jets for jet-finders other than kt, C/A or genkt with p>=0 should be interpreted with care.");
1418
1876
  }
1419
1877
  int stop_point = 2*_initial_n - njets;
1420
1878
  if (stop_point < _initial_n) stop_point = _initial_n;
1443
1901
  }
1444
1902
  return jets_local;
1445
1903
}
1446
 
double ClusterSequence::exclusive_dmerge (const int & njets) const {
 
1904
double ClusterSequence::exclusive_dmerge (const int njets) const {
1447
1905
  assert(njets >= 0);
1448
1906
  if (njets >= _initial_n) {return 0.0;}
1449
1907
  return _history[2*_initial_n-njets-1].dij;
1450
1908
}
1451
 
double ClusterSequence::exclusive_dmerge_max (const int & njets) const {
 
1909
double ClusterSequence::exclusive_dmerge_max (const int njets) const {
1452
1910
  assert(njets >= 0);
1453
1911
  if (njets >= _initial_n) {return 0.0;}
1454
1912
  return _history[2*_initial_n-njets-1].max_dij_so_far;
1455
1913
}
1456
1914
std::vector<PseudoJet> ClusterSequence::exclusive_subjets 
1457
 
   (const PseudoJet & jet, const double & dcut) const {
 
1915
   (const PseudoJet & jet, const double dcut) const {
1458
1916
  set<const history_element*> subhist;
1459
1917
  get_subhist_set(subhist, jet, dcut, 0);
1460
1918
  vector<PseudoJet> subjets;
1466
1924
  return subjets;
1467
1925
}
1468
1926
int ClusterSequence::n_exclusive_subjets(const PseudoJet & jet, 
1469
 
                        const double & dcut) const {
 
1927
                        const double dcut) const {
1470
1928
  set<const history_element*> subhist;
1471
1929
  get_subhist_set(subhist, jet, dcut, 0);
1472
1930
  return subhist.size();
1659
2117
  }
1660
2118
}
1661
2119
void ClusterSequence::_add_step_to_history (
1662
 
               const int & step_number, const int & parent1, 
1663
 
               const int & parent2, const int & jetp_index,
1664
 
               const double & dij) {
 
2120
               const int step_number, const int parent1, 
 
2121
               const int parent2, const int jetp_index,
 
2122
               const double dij) {
1665
2123
  history_element element;
1666
2124
  element.parent1 = parent1;
1667
2125
  element.parent2 = parent2;
1673
2131
  int local_step = _history.size()-1;
1674
2132
  assert(local_step == step_number);
1675
2133
  assert(parent1 >= 0);
 
2134
  if (_history[parent1].child != Invalid){
 
2135
    throw InternalError("trying to recomine an object that has previsously been recombined");
 
2136
  }
1676
2137
  _history[parent1].child = local_step;
1677
 
  if (parent2 >= 0) {_history[parent2].child = local_step;}
 
2138
  if (parent2 >= 0) {
 
2139
    if (_history[parent2].child != Invalid){
 
2140
      throw InternalError("trying to recomine an object that has previsously been recombined");
 
2141
    }
 
2142
    _history[parent2].child = local_step;
 
2143
  }
1678
2144
  if (jetp_index != Invalid) {
1679
2145
    assert(jetp_index >= 0);
1680
2146
    _jets[jetp_index].set_cluster_hist_index(local_step);
1761
2227
  }
1762
2228
}
1763
2229
void ClusterSequence::_do_ij_recombination_step(
1764
 
                               const int & jet_i, const int & jet_j, 
1765
 
                               const double & dij, 
 
2230
                               const int jet_i, const int jet_j, 
 
2231
                               const double dij, 
1766
2232
                               int & newjet_k) {
1767
2233
  PseudoJet newjet(false); 
1768
2234
  _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
1776
2242
                       newjet_k, dij);
1777
2243
}
1778
2244
void ClusterSequence::_do_iB_recombination_step(
1779
 
                                  const int & jet_i, const double & diB) {
 
2245
                                  const int jet_i, const double diB) {
1780
2246
  int newstep_k = _history.size();
1781
2247
  _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
1782
2248
                       Invalid, diB);
1989
2455
    points[i].sanitize(); // make sure things are in the right range
1990
2456
  }
1991
2457
  auto_ptr<DynamicNearestNeighbours> DNN;
 
2458
  bool verbose = false;
1992
2459
#ifndef __FJCORE_DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
1993
 
  bool verbose = false;
1994
2460
  bool ignore_nearest_is_mirror = (_Rparam < twopi);
1995
2461
  if (_strategy == NlnN4pi) {
1996
2462
    DNN.reset(new Dnn4piCylinder(points,verbose));
2005
2471
    err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
2006
2472
    err << "       supported because FastJet was compiled without CGAL"<<endl;
2007
2473
    throw Error(err.str());
2008
 
  }
 
2474
  } else
2009
2475
#endif // __FJCORE_DROP_CGAL
2010
2476
  {
2011
 
    ostringstream err;
2012
 
    err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
2013
2477
    assert(false);
2014
 
    throw Error(err.str());
2015
2478
  }
2016
2479
  DistMap DijMap;
2017
2480
  for (int ii = 0; ii < n; ii++) {
2028
2491
      SmallestDijPair = DijMap.begin()->second;
2029
2492
      jet_i = SmallestDijPair.first;
2030
2493
      jet_j = SmallestDijPair.second;
 
2494
      if (verbose) cout << "CS_Delaunay found recombination candidate: " << jet_i << " " << jet_j << " " << SmallestDij << endl; // GPS debugging
2031
2495
      DijMap.erase(DijMap.begin());
2032
2496
      recombine_with_beam = (jet_j == BeamJet);
2033
2497
      if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);} 
2034
2498
      else {Valid2 = true;}
 
2499
      if (verbose) cout << "CS_Delaunay validities i & j: " << DNN->Valid(jet_i) << " " << Valid2 << endl;
2035
2500
    } while ( !DNN->Valid(jet_i) || !Valid2);
2036
2501
    if (! recombine_with_beam) {
2037
2502
      int nn; // will be index of new jet
 
2503
      if (verbose) cout << "CS_Delaunay call _do_ij_recomb: " << jet_i << " " << jet_j << " " << SmallestDij << endl; // GPS debug
2038
2504
      _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
2039
2505
      EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
2040
2506
      newpoint.sanitize(); // make sure it is in correct range
2041
2507
      points.push_back(newpoint);
2042
2508
    } else {
 
2509
      if (verbose) cout << "CS_Delaunay call _do_iB_recomb: " << jet_i << " " << SmallestDij << endl; // GPS debug
2043
2510
      _do_iB_recombination_step(jet_i, SmallestDij);
2044
2511
    }
2045
2512
    if (i == n-1) {break;}
2062
2529
  } // end clustering loop 
2063
2530
}
2064
2531
void ClusterSequence::_add_ktdistance_to_map(
2065
 
                          const int & ii, 
 
2532
                          const int ii, 
2066
2533
                          DistMap & DijMap,
2067
2534
                          const DynamicNearestNeighbours * DNN) {
2068
2535
  double yiB = jet_scale_for_algorithm(_jets[ii]);
2283
2750
  _tile_size_eta = default_size;
2284
2751
  _n_tiles_phi   = max(3,int(floor(twopi/default_size)));
2285
2752
  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
2286
 
  _tiles_eta_min = 0.0;
2287
 
  _tiles_eta_max = 0.0;
2288
 
  const double maxrap = 7.0;
2289
 
  for(unsigned int i = 0; i < _jets.size(); i++) {
2290
 
    double eta = _jets[i].rap();
2291
 
    if (abs(eta) < maxrap) {
2292
 
      if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
2293
 
      if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
2294
 
    }
2295
 
  }
 
2753
  TilingExtent tiling_analysis(*this);
 
2754
  _tiles_eta_min = tiling_analysis.minrap();
 
2755
  _tiles_eta_max = tiling_analysis.maxrap();
2296
2756
  _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
2297
2757
  _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
2298
2758
  _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
2331
2791
    }
2332
2792
  }
2333
2793
}
2334
 
int ClusterSequence::_tile_index(const double & eta, const double & phi) const {
 
2794
int ClusterSequence::_tile_index(const double eta, const double phi) const {
2335
2795
  int ieta, iphi;
2336
2796
  if      (eta <= _tiles_eta_min) {ieta = 0;}
2337
2797
  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
2392
2852
  TiledJet * briefjets = new TiledJet[n];
2393
2853
  TiledJet * jetA = briefjets, * jetB;
2394
2854
  TiledJet oldB;
2395
 
  oldB.tile_index=0; // prevents a gcc warning  
 
2855
  oldB.tile_index=0; // prevents a gcc warning
2396
2856
  vector<int> tile_union(3*n_tile_neighbours);
2397
2857
  for (int i = 0; i< n; i++) {
2398
2858
    _tj_set_jetinfo(jetA, i);
2546
3006
  TiledJet * briefjets = new TiledJet[n];
2547
3007
  TiledJet * jetA = briefjets, * jetB;
2548
3008
  TiledJet oldB;
2549
 
  oldB.tile_index=0; // prevents a gcc warning  
 
3009
  oldB.tile_index=0; // prevents a gcc warning
2550
3010
  vector<int> tile_union(3*n_tile_neighbours);
2551
3011
  for (int i = 0; i< n; i++) {
2552
3012
    _tj_set_jetinfo(jetA, i);
2846
3306
}
2847
3307
FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
2848
3308
#include <sstream>
2849
 
#ifdef FJCORE_HAVE_EXECINFO_H
2850
 
#include <execinfo.h>
2851
 
#include <cstdlib>
2852
 
#endif
2853
3309
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
2854
3310
using namespace std;
2855
3311
bool Error::_print_errors = true;
2856
3312
bool Error::_print_backtrace = false;
2857
3313
ostream * Error::_default_ostr = & cerr;
 
3314
#if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
 
3315
  LimitedWarning Error::_execinfo_undefined;
 
3316
#endif
2858
3317
Error::Error(const std::string & message_in) {
2859
3318
  _message = message_in; 
2860
3319
  if (_print_errors && _default_ostr){
2861
3320
    ostringstream oss;
2862
3321
    oss << "fjcore::Error:  "<< message_in << endl;
2863
 
#ifdef FJCORE_HAVE_EXECINFO_H
2864
 
    if (_print_backtrace){
2865
 
      void * array[10];
2866
 
      char ** messages;
2867
 
      int size = backtrace(array, 10);
2868
 
      messages = backtrace_symbols(array, size);
2869
 
      oss << "stack:" << endl;
2870
 
      for (int i = 1; i < size && messages != NULL; ++i){
2871
 
        oss << "  #" << i << ": " << messages[i] << endl;
2872
 
      }
2873
 
      free(messages);
2874
 
    }
2875
 
#endif
2876
3322
    *_default_ostr << oss.str();
2877
3323
    _default_ostr->flush(); 
2878
3324
  }
2879
3325
}
 
3326
void Error::set_print_backtrace(bool enabled) {
 
3327
#if (!defined(FJCORE_HAVE_EXECINFO_H)) || defined(__FJCORE__)
 
3328
  if (enabled) {
 
3329
    _execinfo_undefined.warn("Error::set_print_backtrace(true) will not work with this build of FastJet");
 
3330
  }
 
3331
#endif    
 
3332
  _print_backtrace = enabled;
 
3333
}
2880
3334
FJCORE_END_NAMESPACE
2881
3335
#include <string>
2882
3336
#include <sstream>
2902
3356
      throw Error(oss.str());
2903
3357
    }
2904
3358
  }
2905
 
  switch (jet_algorithm_in) {
2906
 
  case ee_kt_algorithm:
2907
 
    if (nparameters != 0) {
2908
 
      ostringstream oss;
2909
 
      oss << "ee_kt_algorithm should be constructed with 0 parameters but was called with " 
2910
 
          << nparameters << " parameter(s)\n";
2911
 
      throw Error(oss.str()); 
2912
 
    }
2913
 
    break;
2914
 
  case genkt_algorithm: 
2915
 
  case ee_genkt_algorithm: 
2916
 
    if (nparameters != 2) {
2917
 
      ostringstream oss;
2918
 
      oss << "(ee_)genkt_algorithm should be constructed with 2 parameters but was called with " 
2919
 
          << nparameters << " parameter(s)\n";
2920
 
      throw Error(oss.str()); 
2921
 
    }
2922
 
    break;
2923
 
  default:
2924
 
    if (nparameters != 1) {
2925
 
      ostringstream oss;
2926
 
      oss << "The jet algorithm you requested ("
2927
 
          << jet_algorithm_in << ") should be constructed with 1 parameter but was called with " 
2928
 
          << nparameters << " parameter(s)\n";
2929
 
      throw Error(oss.str()); 
2930
 
    }
 
3359
  unsigned int nparameters_expected = n_parameters_for_algorithm(jet_algorithm_in);
 
3360
  if (nparameters != (int) nparameters_expected){
 
3361
    ostringstream oss;
 
3362
    oss << "The jet algorithm you requested ("
 
3363
        << jet_algorithm_in << ") should be constructed with " << nparameters_expected 
 
3364
        << " parameter(s) but was called with " << nparameters << " parameter(s)\n";
 
3365
    throw Error(oss.str()); 
2931
3366
  }
2932
3367
  assert (_strategy  != plugin_strategy);
2933
3368
  _plugin = NULL;
2934
3369
  set_recombination_scheme(recomb_scheme_in);
2935
3370
  set_extra_param(0.0); // make sure it's defined
2936
3371
}
 
3372
bool JetDefinition::is_spherical() const {
 
3373
  if (jet_algorithm() == plugin_algorithm) {
 
3374
    return plugin()->is_spherical();
 
3375
  } else {
 
3376
    return (jet_algorithm() == ee_kt_algorithm ||  // as of 2013-02-14, the two
 
3377
            jet_algorithm() == ee_genkt_algorithm  // native spherical algorithms
 
3378
            );
 
3379
  }
 
3380
}
2937
3381
string JetDefinition::description() const {
2938
3382
  ostringstream name;
 
3383
  name << description_no_recombiner();
 
3384
  if ((jet_algorithm() == plugin_algorithm) || (jet_algorithm() == undefined_jet_algorithm)){
 
3385
    return name.str();
 
3386
  }
 
3387
  if (n_parameters_for_algorithm(jet_algorithm()) == 0)
 
3388
    name << " with ";
 
3389
  else 
 
3390
    name << " and ";
 
3391
  name << recombiner()->description();
 
3392
  return name.str();
 
3393
}
 
3394
string JetDefinition::description_no_recombiner() const {
 
3395
  ostringstream name;
2939
3396
  if (jet_algorithm() == plugin_algorithm) {
2940
3397
    return plugin()->description();
2941
 
  } else if (jet_algorithm() == kt_algorithm) {
2942
 
    name << "Longitudinally invariant kt algorithm with R = " << R();
2943
 
    name << " and " << recombiner()->description();
2944
 
  } else if (jet_algorithm() == cambridge_algorithm) {
2945
 
    name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 
2946
 
         << R() ;
2947
 
    name << " and " << recombiner()->description();
2948
 
  } else if (jet_algorithm() == antikt_algorithm) {
2949
 
    name << "Longitudinally invariant anti-kt algorithm with R = " 
2950
 
         << R() ;
2951
 
    name << " and " << recombiner()->description();
2952
 
  } else if (jet_algorithm() == genkt_algorithm) {
2953
 
    name << "Longitudinally invariant generalised kt algorithm with R = " 
2954
 
         << R() << ", p = " << extra_param();
2955
 
    name << " and " << recombiner()->description();
2956
 
  } else if (jet_algorithm() == cambridge_for_passive_algorithm) {
2957
 
    name << "Longitudinally invariant Cambridge/Aachen algorithm with R = " 
2958
 
         << R() << "and a special hack whereby particles with kt < " 
2959
 
         << extra_param() << "are treated as passive ghosts";
2960
 
  } else if (jet_algorithm() == ee_kt_algorithm) {
2961
 
    name << "e+e- kt (Durham) algorithm (NB: no R)";
2962
 
    name << " with " << recombiner()->description();
2963
 
  } else if (jet_algorithm() == ee_genkt_algorithm) {
2964
 
    name << "e+e- generalised kt algorithm with R = " 
2965
 
         << R() << ", p = " << extra_param();
2966
 
    name << " and " << recombiner()->description();
2967
3398
  } else if (jet_algorithm() == undefined_jet_algorithm) {
2968
 
    name << "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
2969
 
  } else {
2970
 
    throw Error("JetDefinition::description(): unrecognized jet_algorithm");
 
3399
    return "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
2971
3400
  }
 
3401
  name << algorithm_description(jet_algorithm());
 
3402
  switch (n_parameters_for_algorithm(jet_algorithm())){
 
3403
  case 0: name << " (NB: no R)"; break;
 
3404
  case 1: name << " with R = " << R(); break; // the parameter is always R
 
3405
  case 2: 
 
3406
    name << " with R = " << R();
 
3407
    if (jet_algorithm() == cambridge_for_passive_algorithm){
 
3408
      name << "and a special hack whereby particles with kt < " 
 
3409
           << extra_param() << "are treated as passive ghosts";
 
3410
    } else {
 
3411
      name << ", p = " << extra_param();
 
3412
    }
 
3413
  };
2972
3414
  return name.str();
2973
3415
}
 
3416
string JetDefinition::algorithm_description(const JetAlgorithm jet_alg){
 
3417
  ostringstream name;
 
3418
  switch (jet_alg){
 
3419
  case plugin_algorithm:                return "plugin algorithm";
 
3420
  case kt_algorithm:                    return "Longitudinally invariant kt algorithm";
 
3421
  case cambridge_algorithm:             return "Longitudinally invariant Cambridge/Aachen algorithm";
 
3422
  case antikt_algorithm:                return "Longitudinally invariant anti-kt algorithm";
 
3423
  case genkt_algorithm:                 return "Longitudinally invariant generalised kt algorithm";
 
3424
  case cambridge_for_passive_algorithm: return "Longitudinally invariant Cambridge/Aachen algorithm";
 
3425
  case ee_kt_algorithm:                 return "e+e- kt (Durham) algorithm (NB: no R)";
 
3426
  case ee_genkt_algorithm:              return "e+e- generalised kt algorithm";
 
3427
  case undefined_jet_algorithm:         return "undefined jet algorithm";
 
3428
  default:
 
3429
    throw Error("JetDefinition::algorithm_description(): unrecognized jet_algorithm");
 
3430
  };
 
3431
}
 
3432
unsigned int JetDefinition::n_parameters_for_algorithm(const JetAlgorithm jet_alg){
 
3433
  switch (jet_alg) {
 
3434
  case ee_kt_algorithm:    return 0;
 
3435
  case genkt_algorithm:
 
3436
  case ee_genkt_algorithm: return 2;
 
3437
  default:                 return 1;
 
3438
  };
 
3439
}
2974
3440
void JetDefinition::set_recombination_scheme(
2975
3441
                               RecombinationScheme recomb_scheme) {
2976
3442
  _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
2977
 
  if (_recombiner_shared()) _recombiner_shared.reset();
 
3443
  if (_shared_recombiner()) _shared_recombiner.reset();
2978
3444
  _recombiner = 0;
2979
3445
}
 
3446
void JetDefinition::set_recombiner(const JetDefinition &other_jet_def){
 
3447
  assert(other_jet_def._recombiner || 
 
3448
         other_jet_def.recombination_scheme() != external_scheme);
 
3449
  if (other_jet_def._recombiner == 0){
 
3450
    set_recombination_scheme(other_jet_def.recombination_scheme());
 
3451
    return;
 
3452
  }
 
3453
  _recombiner = other_jet_def._recombiner;
 
3454
  _default_recombiner = DefaultRecombiner(external_scheme);
 
3455
  _shared_recombiner.reset(other_jet_def._shared_recombiner);
 
3456
}
2980
3457
bool JetDefinition::has_same_recombiner(const JetDefinition &other_jd) const{
2981
3458
  const RecombinationScheme & scheme = recombination_scheme();
2982
3459
  if (other_jd.recombination_scheme() != scheme) return false;
2986
3463
void JetDefinition::delete_recombiner_when_unused(){
2987
3464
  if (_recombiner == 0){
2988
3465
    throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
 
3466
  } else if (_shared_recombiner.get()) {
 
3467
    throw Error("Error in JetDefinition::delete_recombiner_when_unused: the recombiner is already scheduled for deletion when unused (or was already set as shared)");
2989
3468
  }
2990
 
  _recombiner_shared.reset(_recombiner);
 
3469
  _shared_recombiner.reset(_recombiner);
2991
3470
}
2992
3471
void JetDefinition::delete_plugin_when_unused(){
2993
3472
  if (_plugin == 0){
3011
3490
    return "boost-invariant pt scheme recombination";
3012
3491
  case BIpt2_scheme:
3013
3492
    return "boost-invariant pt2 scheme recombination";
 
3493
  case WTA_pt_scheme:
 
3494
    return "pt-ordered Winner-Takes-All recombination";
 
3495
  case WTA_modp_scheme:
 
3496
    return "|3-momentum|-ordered Winner-Takes-All recombination";
3014
3497
  default:
3015
3498
    ostringstream err;
3016
3499
    err << "DefaultRecombiner: unrecognized recombination scheme " 
3041
3524
    weighta = pa.perp2(); 
3042
3525
    weightb = pb.perp2();
3043
3526
    break;
 
3527
  case WTA_pt_scheme:{
 
3528
    const PseudoJet & phard = (pa.pt2() >= pb.pt2()) ? pa : pb;
 
3529
    pab.reset_PtYPhiM(pa.pt()+pb.pt(), 
 
3530
                      phard.rap(), phard.phi(), phard.m());
 
3531
    return;}
 
3532
  case WTA_modp_scheme:{
 
3533
    bool a_hardest = (pa.modp2() >= pb.modp2());
 
3534
    const PseudoJet & phard = a_hardest ? pa : pb;
 
3535
    const PseudoJet & psoft = a_hardest ? pb : pa;
 
3536
    double modp_hard = phard.modp();
 
3537
    double modp_ab = modp_hard + psoft.modp();
 
3538
    if (phard.modp2()==0.0){
 
3539
      pab.reset(0.0, 0.0, 0.0, phard.m());
 
3540
    } else {
 
3541
      double scale = modp_ab/modp_hard;
 
3542
      pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale,
 
3543
                sqrt(modp_ab*modp_ab + phard.m2()));
 
3544
    }
 
3545
    return;}
3044
3546
  default:
3045
3547
    ostringstream err;
3046
3548
    err << "DefaultRecombiner: unrecognized recombination scheme " 
3064
3566
  case E_scheme:
3065
3567
  case BIpt_scheme:
3066
3568
  case BIpt2_scheme:
 
3569
  case WTA_pt_scheme:
 
3570
  case WTA_modp_scheme:
3067
3571
    break;
3068
3572
  case pt_scheme:
3069
3573
  case pt2_scheme:
3136
3640
ostream * LimitedWarning::_default_ostr = &cerr;
3137
3641
std::list< LimitedWarning::Summary > LimitedWarning::_global_warnings_summary;
3138
3642
int LimitedWarning::_max_warn_default = 5;
3139
 
void LimitedWarning::warn(const std::string & warning) {
3140
 
  warn(warning, _default_ostr);
3141
 
}
3142
 
void LimitedWarning::warn(const std::string & warning, std::ostream * ostr) {
 
3643
void LimitedWarning::warn(const char * warning, std::ostream * ostr) {
3143
3644
  if (_this_warning_summary == 0) {
3144
3645
    _global_warnings_summary.push_back(Summary(warning, 0));
3145
3646
    _this_warning_summary = & (_global_warnings_summary.back());
3146
3647
  }
3147
3648
  if (_n_warn_so_far < _max_warn) {
3148
3649
    ostringstream warnstr;
3149
 
    warnstr << "WARNING: ";
 
3650
    warnstr << "WARNING from FastJet: ";
3150
3651
    warnstr << warning;
3151
3652
    _n_warn_so_far++;
3152
3653
    if (_n_warn_so_far == _max_warn) warnstr << " (LAST SUCH WARNING)";
3174
3675
#include<limits>
3175
3676
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
3176
3677
using namespace std;
3177
 
void MinHeap::_initialise(const std::vector<double> & values){
 
3678
void MinHeap::initialise(const std::vector<double> & values){
3178
3679
  for (unsigned i = values.size(); i < _heap.size(); i++) {
3179
3680
    _heap[i].value = std::numeric_limits<double>::max();
3180
3681
    _heap[i].minloc = &(_heap[i]);
3305
3806
                   jet1.E() -jet2.E()  );
3306
3807
} 
3307
3808
PseudoJet operator* (double coeff, const PseudoJet & jet) {
 
3809
  jet._ensure_valid_rap_phi(); 
3308
3810
  PseudoJet coeff_times_jet(jet);
3309
3811
  coeff_times_jet *= coeff;
3310
3812
  return coeff_times_jet;
3316
3818
  return (1.0/coeff)*jet;
3317
3819
} 
3318
3820
void PseudoJet::operator*=(double coeff) {
 
3821
  _ensure_valid_rap_phi(); 
3319
3822
  _px *= coeff;
3320
3823
  _py *= coeff;
3321
3824
  _pz *= coeff;
3504
4007
bool PseudoJet::has_exclusive_subjets() const{
3505
4008
  return (_structure()) && (_structure->has_exclusive_subjets());
3506
4009
}
3507
 
std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double & dcut) const {
 
4010
std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double dcut) const {
3508
4011
  return validated_structure_ptr()->exclusive_subjets(*this, dcut);
3509
4012
}
3510
 
int PseudoJet::n_exclusive_subjets(const double & dcut) const {
 
4013
int PseudoJet::n_exclusive_subjets(const double dcut) const {
3511
4014
  return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
3512
4015
}
3513
4016
std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const {
3588
4091
}
3589
4092
PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
3590
4093
  vector<PseudoJet> pieces;
 
4094
  pieces.reserve(2);
3591
4095
  pieces.push_back(j1);
3592
4096
  pieces.push_back(j2);
3593
4097
  return join(pieces);
3594
4098
}
3595
4099
PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
3596
4100
  vector<PseudoJet> pieces;
 
4101
  pieces.reserve(3);
3597
4102
  pieces.push_back(j1);
3598
4103
  pieces.push_back(j2);
3599
4104
  pieces.push_back(j3);
3601
4106
}
3602
4107
PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
3603
4108
  vector<PseudoJet> pieces;
 
4109
  pieces.reserve(4);
3604
4110
  pieces.push_back(j1);
3605
4111
  pieces.push_back(j2);
3606
4112
  pieces.push_back(j3);
3693
4199
  }
3694
4200
  return n;
3695
4201
}
 
4202
PseudoJet Selector::sum(const std::vector<PseudoJet> & jets) const {
 
4203
  PseudoJet this_sum(0,0,0,0);
 
4204
  const SelectorWorker * worker_local = validated_worker();
 
4205
  if (worker_local->applies_jet_by_jet()) {
 
4206
    for (unsigned i = 0; i < jets.size(); i++) {
 
4207
      if (worker_local->pass(jets[i])) this_sum += jets[i];
 
4208
    }
 
4209
  } else {
 
4210
    std::vector<const PseudoJet *> jetptrs(jets.size());
 
4211
    for (unsigned i = 0; i < jets.size(); i++) {
 
4212
      jetptrs[i] = & jets[i];
 
4213
    }
 
4214
    worker_local->terminator(jetptrs);
 
4215
    for (unsigned i = 0; i < jetptrs.size(); i++) {
 
4216
      if (jetptrs[i]) this_sum += jets[i];
 
4217
    }
 
4218
  }
 
4219
  return this_sum;
 
4220
}
 
4221
double Selector::scalar_pt_sum(const std::vector<PseudoJet> & jets) const {
 
4222
  double this_sum = 0.0;
 
4223
  const SelectorWorker * worker_local = validated_worker();
 
4224
  if (worker_local->applies_jet_by_jet()) {
 
4225
    for (unsigned i = 0; i < jets.size(); i++) {
 
4226
      if (worker_local->pass(jets[i])) this_sum += jets[i].pt();
 
4227
    }
 
4228
  } else {
 
4229
    std::vector<const PseudoJet *> jetptrs(jets.size());
 
4230
    for (unsigned i = 0; i < jets.size(); i++) {
 
4231
      jetptrs[i] = & jets[i];
 
4232
    }
 
4233
    worker_local->terminator(jetptrs);
 
4234
    for (unsigned i = 0; i < jetptrs.size(); i++) {
 
4235
      if (jetptrs[i]) this_sum += jets[i].pt();
 
4236
    }
 
4237
  }
 
4238
  return this_sum;
 
4239
}
3696
4240
void Selector::sift(const std::vector<PseudoJet> & jets,
3697
4241
                    std::vector<PseudoJet> & jets_that_pass,
3698
4242
                    std::vector<PseudoJet> & jets_that_fail
4232
4776
};
4233
4777
class SW_Circle : public SW_WithReference {
4234
4778
public:
4235
 
  SW_Circle(const double &radius) : _radius2(radius*radius) {}
 
4779
  SW_Circle(const double radius) : _radius2(radius*radius) {}
4236
4780
  virtual SelectorWorker* copy(){ return new SW_Circle(*this);}
4237
4781
  virtual bool pass(const PseudoJet & jet) const {
4238
4782
    if (! _is_initialised)
4259
4803
protected:
4260
4804
  double _radius2;
4261
4805
};
4262
 
Selector SelectorCircle(const double & radius) {
 
4806
Selector SelectorCircle(const double radius) {
4263
4807
  return Selector(new SW_Circle(radius));
4264
4808
}
4265
4809
class SW_Doughnut : public SW_WithReference {
4266
4810
public:
4267
 
  SW_Doughnut(const double &radius_in, const double &radius_out)
 
4811
  SW_Doughnut(const double radius_in, const double radius_out)
4268
4812
    : _radius_in2(radius_in*radius_in), _radius_out2(radius_out*radius_out) {}
4269
4813
  virtual SelectorWorker* copy(){ return new SW_Doughnut(*this);}
4270
4814
  virtual bool pass(const PseudoJet & jet) const {
4293
4837
protected:
4294
4838
  double _radius_in2, _radius_out2;
4295
4839
};
4296
 
Selector SelectorDoughnut(const double & radius_in, const double & radius_out) {
 
4840
Selector SelectorDoughnut(const double radius_in, const double radius_out) {
4297
4841
  return Selector(new SW_Doughnut(radius_in, radius_out));
4298
4842
}
4299
4843
class SW_Strip : public SW_WithReference {
4300
4844
public:
4301
 
  SW_Strip(const double &delta) : _delta(delta) {}
 
4845
  SW_Strip(const double delta) : _delta(delta) {}
4302
4846
  virtual SelectorWorker* copy(){ return new SW_Strip(*this);}
4303
4847
  virtual bool pass(const PseudoJet & jet) const {
4304
4848
    if (! _is_initialised)
4325
4869
protected:
4326
4870
  double _delta;
4327
4871
};
4328
 
Selector SelectorStrip(const double & half_width) {
 
4872
Selector SelectorStrip(const double half_width) {
4329
4873
  return Selector(new SW_Strip(half_width));
4330
4874
}
4331
4875
class SW_Rectangle : public SW_WithReference {
4332
4876
public:
4333
 
  SW_Rectangle(const double &delta_rap, const double &delta_phi)
 
4877
  SW_Rectangle(const double delta_rap, const double delta_phi)
4334
4878
    : _delta_rap(delta_rap),  _delta_phi(delta_phi) {}
4335
4879
  virtual SelectorWorker* copy(){ return new SW_Rectangle(*this);}
4336
4880
  virtual bool pass(const PseudoJet & jet) const {
4358
4902
protected:
4359
4903
  double _delta_rap, _delta_phi;
4360
4904
};
4361
 
Selector SelectorRectangle(const double & half_rap_width, const double & half_phi_width) {
 
4905
Selector SelectorRectangle(const double half_rap_width, const double half_phi_width) {
4362
4906
  return Selector(new SW_Rectangle(half_rap_width, half_phi_width));
4363
4907
}
4364
4908
class SW_PtFractionMin : public SW_WithReference {
4401
4945
  return *this;
4402
4946
}
4403
4947
FJCORE_END_NAMESPACE      // defined in fastjet/internal/base.hh
 
4948
#include <iomanip>
 
4949
using namespace std;
 
4950
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
4951
LazyTiling25::LazyTiling25(ClusterSequence & cs) :
 
4952
  _cs(cs), _jets(cs.jets())
 
4953
{
 
4954
#ifdef INSTRUMENT2
 
4955
  _ncall = 0; // gps tmp
 
4956
  _ncall_dtt = 0; // gps tmp
 
4957
#endif // INSTRUMENT2
 
4958
  _Rparam = cs.jet_def().R();
 
4959
  _R2 = _Rparam * _Rparam;
 
4960
  _invR2 = 1.0 / _R2;
 
4961
  _initialise_tiles();
 
4962
}
 
4963
void LazyTiling25::_initialise_tiles() {
 
4964
  double default_size = max(0.1,_Rparam)/2;
 
4965
  _tile_size_eta = default_size;
 
4966
  _n_tiles_phi   = max(5,int(floor(twopi/default_size)));
 
4967
  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
 
4968
#define _FASTJET_TILING25_USE_TILING_ANALYSIS_
 
4969
#ifdef  _FASTJET_TILING25_USE_TILING_ANALYSIS_
 
4970
  TilingExtent tiling_analysis(_cs);
 
4971
  _tiles_eta_min = tiling_analysis.minrap();
 
4972
  _tiles_eta_max = tiling_analysis.maxrap();
 
4973
#else // not _FASTJET_TILING25_USE_TILING_ANALYSIS_
 
4974
  _tiles_eta_min = 0.0;
 
4975
  _tiles_eta_max = 0.0;
 
4976
  const double maxrap = 7.0;
 
4977
  for(unsigned int i = 0; i < _jets.size(); i++) {
 
4978
    double eta = _jets[i].rap();
 
4979
    if (abs(eta) < maxrap) {
 
4980
      if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
 
4981
      if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
 
4982
    }
 
4983
  }
 
4984
#endif // _FASTJET_TILING25_USE_TILING_ANALYSIS_
 
4985
#define FASTJET_LAZY25_MIN3TILESY
 
4986
#ifdef FASTJET_LAZY25_MIN3TILESY
 
4987
   if (_tiles_eta_max - _tiles_eta_min < 3*_tile_size_eta) {
 
4988
     _tile_size_eta = (_tiles_eta_max - _tiles_eta_min)/3;
 
4989
     _tiles_ieta_min = 0;
 
4990
     _tiles_ieta_max = 2;
 
4991
     _tiles_eta_max -= _tile_size_eta;
 
4992
   } else {
 
4993
#endif //FASTJET_LAZY25_MIN3TILESY
 
4994
    _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
 
4995
    _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
 
4996
    _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
 
4997
    _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
 
4998
#ifdef FASTJET_LAZY25_MIN3TILESY
 
4999
   }
 
5000
#endif
 
5001
  _tile_half_size_eta = _tile_size_eta * 0.5;
 
5002
  _tile_half_size_phi = _tile_size_phi * 0.5;
 
5003
  vector<bool> use_periodic_delta_phi(_n_tiles_phi, false);
 
5004
  if (_n_tiles_phi <= 5) {
 
5005
    fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(), true);
 
5006
  } else {
 
5007
    use_periodic_delta_phi[0] = true;
 
5008
    use_periodic_delta_phi[1] = true;
 
5009
    use_periodic_delta_phi[_n_tiles_phi-2] = true;
 
5010
    use_periodic_delta_phi[_n_tiles_phi-1] = true;
 
5011
  }
 
5012
  _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
 
5013
  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
 
5014
    for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
 
5015
      Tile25 * tile = & _tiles[_tile_index(ieta,iphi)];
 
5016
      tile->head = NULL; // first element of tiles points to itself
 
5017
      tile->begin_tiles[0] =  tile;
 
5018
      Tile25 ** pptile = & (tile->begin_tiles[0]);
 
5019
      pptile++;
 
5020
      tile->surrounding_tiles = pptile;
 
5021
      if (ieta > _tiles_ieta_min) {
 
5022
        // with the itile subroutine, we can safely run tiles from
 
5023
        // idphi=-1 to idphi=+1, because it takes care of
 
5024
        // negative and positive boundaries
 
5025
        for (int idphi = -2; idphi <=+2; idphi++) {
 
5026
          *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
 
5027
          pptile++;
 
5028
        }       
 
5029
      }
 
5030
      if (ieta > _tiles_ieta_min + 1) {
 
5031
        // with the itile subroutine, we can safely run tiles from
 
5032
        // idphi=-1 to idphi=+1, because it takes care of
 
5033
        // negative and positive boundaries
 
5034
        for (int idphi = -2; idphi <= +2; idphi++) {
 
5035
          *pptile = & _tiles[_tile_index(ieta-2,iphi+idphi)];
 
5036
          pptile++;
 
5037
        }       
 
5038
      }
 
5039
      *pptile = & _tiles[_tile_index(ieta,iphi-1)];
 
5040
      pptile++;
 
5041
      *pptile = & _tiles[_tile_index(ieta,iphi-2)];
 
5042
      pptile++;
 
5043
      tile->RH_tiles = pptile;
 
5044
      *pptile = & _tiles[_tile_index(ieta,iphi+1)];
 
5045
      pptile++;
 
5046
      *pptile = & _tiles[_tile_index(ieta,iphi+2)];
 
5047
      pptile++;
 
5048
      if (ieta < _tiles_ieta_max) {
 
5049
        for (int idphi = -2; idphi <= +2; idphi++) {
 
5050
          *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
 
5051
          pptile++;
 
5052
        }       
 
5053
      }
 
5054
      if (ieta < _tiles_ieta_max - 1) {
 
5055
        for (int idphi = -2; idphi <= +2; idphi++) {
 
5056
          *pptile = & _tiles[_tile_index(ieta+2,iphi+idphi)];
 
5057
          pptile++;
 
5058
        }       
 
5059
      }
 
5060
      tile->end_tiles = pptile;
 
5061
      tile->tagged = false;
 
5062
      tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
 
5063
      tile->max_NN_dist = 0;
 
5064
      tile->eta_centre = (ieta-_tiles_ieta_min+0.5)*_tile_size_eta + _tiles_eta_min;
 
5065
      tile->phi_centre = (iphi+0.5)*_tile_size_phi;
 
5066
    }
 
5067
  }
 
5068
}
 
5069
int LazyTiling25::_tile_index(const double eta, const double phi) const {
 
5070
  int ieta, iphi;
 
5071
  if      (eta <= _tiles_eta_min) {ieta = 0;}
 
5072
  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
 
5073
  else {
 
5074
    ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
 
5075
    if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
 
5076
      ieta = _tiles_ieta_max-_tiles_ieta_min;} 
 
5077
  }
 
5078
  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
 
5079
  return (iphi + ieta * _n_tiles_phi);
 
5080
}
 
5081
inline void LazyTiling25::_tj_set_jetinfo( TiledJet * const jet,
 
5082
                                              const int _jets_index) {
 
5083
  _bj_set_jetinfo<>(jet, _jets_index);
 
5084
  jet->tile_index = _tile_index(jet->eta, jet->phi);
 
5085
  Tile25 * tile = &_tiles[jet->tile_index];
 
5086
  jet->previous   = NULL;
 
5087
  jet->next       = tile->head;
 
5088
  if (jet->next != NULL) {jet->next->previous = jet;}
 
5089
  tile->head      = jet;
 
5090
}
 
5091
void LazyTiling25::_bj_remove_from_tiles(TiledJet * const jet) {
 
5092
  Tile25 * tile = & _tiles[jet->tile_index];
 
5093
  if (jet->previous == NULL) {
 
5094
    tile->head = jet->next;
 
5095
  } else {
 
5096
    jet->previous->next = jet->next;
 
5097
  }
 
5098
  if (jet->next != NULL) {
 
5099
    jet->next->previous = jet->previous;
 
5100
  }
 
5101
}
 
5102
void LazyTiling25::_print_tiles(TiledJet * briefjets ) const {
 
5103
  for (vector<Tile25>::const_iterator tile = _tiles.begin(); 
 
5104
       tile < _tiles.end(); tile++) {
 
5105
    cout << "Tile " << tile - _tiles.begin()
 
5106
         << " at " << setw(10) << tile->eta_centre << "," << setw(10) << tile->phi_centre
 
5107
         << " = ";
 
5108
    vector<int> list;
 
5109
    for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
 
5110
      list.push_back(jetI-briefjets);
 
5111
    }
 
5112
    sort(list.begin(),list.end());
 
5113
    for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
 
5114
    cout <<"\n";
 
5115
  }
 
5116
}
 
5117
void LazyTiling25::_add_neighbours_to_tile_union(const int tile_index, 
 
5118
               vector<int> & tile_union, int & n_near_tiles) const {
 
5119
  for (Tile25 * const * near_tile = _tiles[tile_index].begin_tiles; 
 
5120
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
 
5121
    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5122
    n_near_tiles++;
 
5123
  }
 
5124
}
 
5125
inline void LazyTiling25::_add_untagged_neighbours_to_tile_union(
 
5126
               const int tile_index, 
 
5127
               vector<int> & tile_union, int & n_near_tiles)  {
 
5128
  for (Tile25 ** near_tile = _tiles[tile_index].begin_tiles; 
 
5129
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
 
5130
    if (! (*near_tile)->tagged) {
 
5131
      (*near_tile)->tagged = true;
 
5132
      tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5133
      n_near_tiles++;
 
5134
    }
 
5135
  }
 
5136
}
 
5137
inline void LazyTiling25::_add_untagged_neighbours_to_tile_union_using_max_info(
 
5138
               const TiledJet * jet, 
 
5139
               vector<int> & tile_union, int & n_near_tiles)  {
 
5140
  Tile25 & tile = _tiles[jet->tile_index];
 
5141
  for (Tile25 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
 
5142
    if ((*near_tile)->tagged) continue;
 
5143
    double dist = _distance_to_tile(jet, *near_tile) - tile_edge_security_margin;
 
5144
    if (dist > (*near_tile)->max_NN_dist) continue;
 
5145
    (*near_tile)->tagged = true;
 
5146
    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5147
    n_near_tiles++;
 
5148
  }
 
5149
}
 
5150
inline double LazyTiling25::_distance_to_tile(const TiledJet * bj, const Tile25 * tile) 
 
5151
#ifdef INSTRUMENT2
 
5152
   {
 
5153
  _ncall_dtt++; // GPS tmp
 
5154
#else
 
5155
  const {
 
5156
#endif // INSTRUMENT2
 
5157
  double deta;
 
5158
  if (_tiles[bj->tile_index].eta_centre == tile->eta_centre) deta = 0;
 
5159
  else   deta = std::abs(bj->eta - tile->eta_centre) - _tile_half_size_eta;
 
5160
  double dphi = std::abs(bj->phi - tile->phi_centre);
 
5161
  if (dphi > pi) dphi = twopi-dphi;
 
5162
  dphi -= _tile_half_size_phi;
 
5163
  if (dphi < 0) dphi = 0;
 
5164
  return dphi*dphi + deta*deta;
 
5165
}
 
5166
inline void LazyTiling25::_update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, vector<TiledJet *> & jets_for_minheap) {
 
5167
  double dist = _bj_dist(jetI,jetX);
 
5168
  if (dist < jetI->NN_dist) {
 
5169
    if (jetI != jetX) {
 
5170
      jetI->NN_dist = dist;
 
5171
      jetI->NN = jetX;
 
5172
      if (!jetI->minheap_update_needed()) {
 
5173
        jetI->label_minheap_update_needed();
 
5174
        jets_for_minheap.push_back(jetI);
 
5175
      }
 
5176
    }
 
5177
  }
 
5178
  if (dist < jetX->NN_dist) {
 
5179
    if (jetI != jetX) {
 
5180
      jetX->NN_dist = dist;
 
5181
      jetX->NN      = jetI;}
 
5182
  }
 
5183
}
 
5184
inline void LazyTiling25::_set_NN(TiledJet * jetI, 
 
5185
                              vector<TiledJet *> & jets_for_minheap) {
 
5186
  jetI->NN_dist = _R2;
 
5187
  jetI->NN      = NULL;
 
5188
  if (!jetI->minheap_update_needed()) {
 
5189
    jetI->label_minheap_update_needed();
 
5190
    jets_for_minheap.push_back(jetI);}
 
5191
  Tile25 * tile_ptr = &_tiles[jetI->tile_index];
 
5192
    for (Tile25 ** near_tile  = tile_ptr->begin_tiles; 
 
5193
         near_tile != tile_ptr->end_tiles; near_tile++) {
 
5194
      if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
 
5195
      for (TiledJet * jetJ  = (*near_tile)->head; 
 
5196
           jetJ != NULL; jetJ = jetJ->next) {
 
5197
        double dist = _bj_dist(jetI,jetJ);
 
5198
        if (dist < jetI->NN_dist && jetJ != jetI) {
 
5199
          jetI->NN_dist = dist; jetI->NN = jetJ;
 
5200
        }
 
5201
      }
 
5202
    }
 
5203
}
 
5204
void LazyTiling25::run() {
 
5205
  int n = _jets.size();
 
5206
  if (n == 0) return; 
 
5207
  TiledJet * briefjets = new TiledJet[n];
 
5208
  TiledJet * jetA = briefjets, * jetB;
 
5209
  TiledJet oldB = briefjets[0]; 
 
5210
  vector<int> tile_union(3*25);
 
5211
  for (int i = 0; i< n; i++) {
 
5212
    _tj_set_jetinfo(jetA, i);
 
5213
    jetA++; // move on to next entry of briefjets
 
5214
  }
 
5215
  TiledJet * head = briefjets; // a nicer way of naming start
 
5216
  vector<Tile25>::iterator tile;
 
5217
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5218
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5219
      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
 
5220
        double dist = _bj_dist_not_periodic(jetA,jetB);
 
5221
        if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5222
        if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5223
      }
 
5224
    }
 
5225
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5226
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
 
5227
    }
 
5228
  }
 
5229
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5230
    if (tile->use_periodic_delta_phi) {
 
5231
      for (Tile25 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
 
5232
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5233
          double dist_to_tile = _distance_to_tile(jetA, *RTile);
 
5234
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
 
5235
          bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
 
5236
          if (relevant_for_jetA || relevant_for_RTile) {
 
5237
            for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
 
5238
              double dist = _bj_dist(jetA,jetB);
 
5239
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5240
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5241
            }
 
5242
          } 
 
5243
        }
 
5244
      }
 
5245
    } else {
 
5246
      for (Tile25 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
 
5247
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5248
          double dist_to_tile = _distance_to_tile(jetA, *RTile);
 
5249
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
 
5250
          bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
 
5251
          if (relevant_for_jetA || relevant_for_RTile) {
 
5252
            for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
 
5253
              double dist = _bj_dist_not_periodic(jetA,jetB);
 
5254
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5255
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5256
            }
 
5257
          } 
 
5258
        }
 
5259
      }
 
5260
    }
 
5261
  }
 
5262
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5263
    tile->max_NN_dist = 0;
 
5264
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5265
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
 
5266
    }
 
5267
  }
 
5268
#ifdef INSTRUMENT2
 
5269
  cout << "intermediate ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
 
5270
#endif // INSTRUMENT2
 
5271
  vector<double> diJs(n);
 
5272
  for (int i = 0; i < n; i++) {
 
5273
    diJs[i] = _bj_diJ(&briefjets[i]);
 
5274
    briefjets[i].label_minheap_update_done();
 
5275
  }
 
5276
  MinHeap minheap(diJs);
 
5277
  vector<TiledJet *> jets_for_minheap;
 
5278
  jets_for_minheap.reserve(n); 
 
5279
  int history_location = n-1;
 
5280
  while (n > 0) {
 
5281
    double diJ_min = minheap.minval() *_invR2;
 
5282
    jetA = head + minheap.minloc();
 
5283
    history_location++;
 
5284
    jetB = jetA->NN;
 
5285
    if (jetB != NULL) {
 
5286
      if (jetA < jetB) {std::swap(jetA,jetB);}
 
5287
      int nn; // new jet index
 
5288
      _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
 
5289
      _bj_remove_from_tiles(jetA);
 
5290
      oldB = * jetB;  // take a copy because we will need it...
 
5291
      _bj_remove_from_tiles(jetB);
 
5292
      _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
 
5293
    } else {
 
5294
      _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
 
5295
      _bj_remove_from_tiles(jetA);
 
5296
    }
 
5297
    minheap.remove(jetA-head);
 
5298
    int n_near_tiles = 0;
 
5299
    if (jetB != NULL) {
 
5300
      Tile25 & jetB_tile = _tiles[jetB->tile_index];
 
5301
      for (Tile25 ** near_tile  = jetB_tile.begin_tiles; 
 
5302
                   near_tile != jetB_tile.end_tiles; near_tile++) {
 
5303
        double dist_to_tile = _distance_to_tile(jetB, *near_tile);
 
5304
        bool relevant_for_jetB  = dist_to_tile <= jetB->NN_dist;
 
5305
        bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
 
5306
        bool relevant = relevant_for_jetB || relevant_for_near_tile;
 
5307
        if (! relevant) continue;
 
5308
        tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5309
        (*near_tile)->tagged = true;
 
5310
        n_near_tiles++;
 
5311
        for (TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
 
5312
          if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
 
5313
          _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
 
5314
        }
 
5315
      }
 
5316
    }
 
5317
    int n_done_tiles = n_near_tiles;
 
5318
    _add_untagged_neighbours_to_tile_union_using_max_info(jetA, 
 
5319
                                           tile_union, n_near_tiles);
 
5320
    if (jetB != NULL) {
 
5321
        _add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
 
5322
                                                              tile_union,n_near_tiles);
 
5323
      jetB->label_minheap_update_needed();
 
5324
      jets_for_minheap.push_back(jetB);
 
5325
    }
 
5326
    for (int itile = 0; itile < n_done_tiles; itile++) {
 
5327
      _tiles[tile_union[itile]].tagged = false;
 
5328
    }
 
5329
    for (int itile = n_done_tiles; itile < n_near_tiles; itile++) {
 
5330
      Tile25 * tile_ptr = &_tiles[tile_union[itile]];
 
5331
      tile_ptr->tagged = false;
 
5332
      for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
 
5333
        if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
 
5334
          _set_NN(jetI, jets_for_minheap);
 
5335
        }
 
5336
      }
 
5337
    }
 
5338
    while (jets_for_minheap.size() > 0) {
 
5339
      TiledJet * jetI = jets_for_minheap.back(); 
 
5340
      jets_for_minheap.pop_back();
 
5341
      minheap.update(jetI-head, _bj_diJ(jetI));
 
5342
      jetI->label_minheap_update_done();
 
5343
      Tile25 & tile_I = _tiles[jetI->tile_index];
 
5344
      if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
 
5345
    }
 
5346
    n--;
 
5347
  }
 
5348
  delete[] briefjets;
 
5349
#ifdef INSTRUMENT2
 
5350
  cout << "ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
 
5351
#endif // INSTRUMENT2
 
5352
}
 
5353
FJCORE_END_NAMESPACE
 
5354
#include <iomanip>
 
5355
#include <limits>
 
5356
#include <cmath>
 
5357
using namespace std;
 
5358
#define _FASTJET_TILING2_USE_TILING_ANALYSIS_
 
5359
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
5360
LazyTiling9::LazyTiling9(ClusterSequence & cs) :
 
5361
  _cs(cs), _jets(cs.jets())
 
5362
{
 
5363
#ifdef INSTRUMENT2
 
5364
  _ncall = 0; // gps tmp
 
5365
  _ncall_dtt = 0; // gps tmp
 
5366
#endif // INSTRUMENT2
 
5367
  _Rparam = cs.jet_def().R();
 
5368
  _R2 = _Rparam * _Rparam;
 
5369
  _invR2 = 1.0 / _R2;
 
5370
  _initialise_tiles();
 
5371
}
 
5372
void LazyTiling9::_initialise_tiles() {
 
5373
  double default_size = max(0.1,_Rparam);
 
5374
  _tile_size_eta = default_size;
 
5375
  _n_tiles_phi   = max(3,int(floor(twopi/default_size)));
 
5376
  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
 
5377
#ifdef _FASTJET_TILING2_USE_TILING_ANALYSIS_
 
5378
  TilingExtent tiling_analysis(_cs);
 
5379
  _tiles_eta_min = tiling_analysis.minrap();
 
5380
  _tiles_eta_max = tiling_analysis.maxrap();
 
5381
#else
 
5382
  _tiles_eta_min = 0.0;
 
5383
  _tiles_eta_max = 0.0;
 
5384
  const double maxrap = 7.0;
 
5385
  for(unsigned int i = 0; i < _jets.size(); i++) {
 
5386
    double eta = _jets[i].rap();
 
5387
    if (abs(eta) < maxrap) {
 
5388
      if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
 
5389
      if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
 
5390
    }
 
5391
  }
 
5392
#endif
 
5393
#define FASTJET_LAZY9_MIN2TILESY
 
5394
#ifdef FASTJET_LAZY9_MIN2TILESY
 
5395
   if (_tiles_eta_max - _tiles_eta_min < 2*_tile_size_eta) {
 
5396
     _tile_size_eta = (_tiles_eta_max - _tiles_eta_min)/2;
 
5397
     _tiles_ieta_min = 0;
 
5398
     _tiles_ieta_max = 1;
 
5399
     _tiles_eta_max -= _tile_size_eta;
 
5400
   } else {
 
5401
#endif //FASTJET_LAZY9_MIN2TILESY
 
5402
  _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
 
5403
  _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
 
5404
  _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
 
5405
  _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
 
5406
#ifdef FASTJET_LAZY9_MIN2TILESY
 
5407
   }
 
5408
#endif
 
5409
  _tile_half_size_eta = _tile_size_eta * 0.5;
 
5410
  _tile_half_size_phi = _tile_size_phi * 0.5;
 
5411
  vector<bool> use_periodic_delta_phi(_n_tiles_phi, false);
 
5412
  if (_n_tiles_phi <= 3) {
 
5413
    fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(), true);
 
5414
  } else {
 
5415
    use_periodic_delta_phi[0] = true;
 
5416
    use_periodic_delta_phi[_n_tiles_phi-1] = true;
 
5417
  }
 
5418
  _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
 
5419
  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
 
5420
    for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
 
5421
      Tile2 * tile = & _tiles[_tile_index(ieta,iphi)];
 
5422
      tile->head = NULL; // first element of tiles points to itself
 
5423
      tile->begin_tiles[0] =  tile;
 
5424
      Tile2 ** pptile = & (tile->begin_tiles[0]);
 
5425
      pptile++;
 
5426
      tile->surrounding_tiles = pptile;
 
5427
      if (ieta > _tiles_ieta_min) {
 
5428
        // with the itile subroutine, we can safely run tiles from
 
5429
        // idphi=-1 to idphi=+1, because it takes care of
 
5430
        // negative and positive boundaries
 
5431
        for (int idphi = -1; idphi <=+1; idphi++) {
 
5432
          *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
 
5433
          pptile++;
 
5434
        }       
 
5435
      }
 
5436
      *pptile = & _tiles[_tile_index(ieta,iphi-1)];
 
5437
      pptile++;
 
5438
      tile->RH_tiles = pptile;
 
5439
      *pptile = & _tiles[_tile_index(ieta,iphi+1)];
 
5440
      pptile++;
 
5441
      if (ieta < _tiles_ieta_max) {
 
5442
        for (int idphi = -1; idphi <= +1; idphi++) {
 
5443
          *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
 
5444
          pptile++;
 
5445
        }       
 
5446
      }
 
5447
      tile->end_tiles = pptile;
 
5448
      tile->tagged = false;
 
5449
      tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
 
5450
      tile->max_NN_dist = 0;
 
5451
      tile->eta_centre = (ieta-_tiles_ieta_min+0.5)*_tile_size_eta + _tiles_eta_min;
 
5452
      tile->phi_centre = (iphi+0.5)*_tile_size_phi;
 
5453
    }
 
5454
  }
 
5455
}
 
5456
int LazyTiling9::_tile_index(const double eta, const double phi) const {
 
5457
  int ieta, iphi;
 
5458
  if      (eta <= _tiles_eta_min) {ieta = 0;}
 
5459
  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
 
5460
  else {
 
5461
    ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
 
5462
    if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
 
5463
      ieta = _tiles_ieta_max-_tiles_ieta_min;} 
 
5464
  }
 
5465
  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
 
5466
  return (iphi + ieta * _n_tiles_phi);
 
5467
}
 
5468
inline void LazyTiling9::_tj_set_jetinfo( TiledJet * const jet,
 
5469
                                              const int _jets_index) {
 
5470
  _bj_set_jetinfo<>(jet, _jets_index);
 
5471
  jet->tile_index = _tile_index(jet->eta, jet->phi);
 
5472
  Tile2 * tile = &_tiles[jet->tile_index];
 
5473
  jet->previous   = NULL;
 
5474
  jet->next       = tile->head;
 
5475
  if (jet->next != NULL) {jet->next->previous = jet;}
 
5476
  tile->head      = jet;
 
5477
}
 
5478
void LazyTiling9::_bj_remove_from_tiles(TiledJet * const jet) {
 
5479
  Tile2 * tile = & _tiles[jet->tile_index];
 
5480
  if (jet->previous == NULL) {
 
5481
    tile->head = jet->next;
 
5482
  } else {
 
5483
    jet->previous->next = jet->next;
 
5484
  }
 
5485
  if (jet->next != NULL) {
 
5486
    jet->next->previous = jet->previous;
 
5487
  }
 
5488
}
 
5489
void LazyTiling9::_print_tiles(TiledJet * briefjets ) const {
 
5490
  for (vector<Tile2>::const_iterator tile = _tiles.begin(); 
 
5491
       tile < _tiles.end(); tile++) {
 
5492
    cout << "Tile " << tile - _tiles.begin()<<" = ";
 
5493
    vector<int> list;
 
5494
    for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
 
5495
      list.push_back(jetI-briefjets);
 
5496
    }
 
5497
    sort(list.begin(),list.end());
 
5498
    for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
 
5499
    cout <<"\n";
 
5500
  }
 
5501
}
 
5502
void LazyTiling9::_add_neighbours_to_tile_union(const int tile_index, 
 
5503
               vector<int> & tile_union, int & n_near_tiles) const {
 
5504
  for (Tile2 * const * near_tile = _tiles[tile_index].begin_tiles; 
 
5505
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
 
5506
    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5507
    n_near_tiles++;
 
5508
  }
 
5509
}
 
5510
inline void LazyTiling9::_add_untagged_neighbours_to_tile_union(
 
5511
               const int tile_index, 
 
5512
               vector<int> & tile_union, int & n_near_tiles)  {
 
5513
  for (Tile2 ** near_tile = _tiles[tile_index].begin_tiles; 
 
5514
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
 
5515
    if (! (*near_tile)->tagged) {
 
5516
      (*near_tile)->tagged = true;
 
5517
      tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5518
      n_near_tiles++;
 
5519
    }
 
5520
  }
 
5521
}
 
5522
inline void LazyTiling9::_add_untagged_neighbours_to_tile_union_using_max_info(
 
5523
               const TiledJet * jet, 
 
5524
               vector<int> & tile_union, int & n_near_tiles)  {
 
5525
  Tile2 & tile = _tiles[jet->tile_index];
 
5526
  for (Tile2 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
 
5527
    if ((*near_tile)->tagged) continue;
 
5528
    double dist = _distance_to_tile(jet, *near_tile) - tile_edge_security_margin;
 
5529
    if (dist > (*near_tile)->max_NN_dist) continue;
 
5530
    (*near_tile)->tagged = true;
 
5531
    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5532
    n_near_tiles++;
 
5533
  }
 
5534
}
 
5535
inline double LazyTiling9::_distance_to_tile(const TiledJet * bj, const Tile2 * tile) 
 
5536
#ifdef INSTRUMENT2
 
5537
   {
 
5538
  _ncall_dtt++; // GPS tmp
 
5539
#else
 
5540
  const {
 
5541
#endif // INSTRUMENT2
 
5542
  double deta;
 
5543
  if (_tiles[bj->tile_index].eta_centre == tile->eta_centre) deta = 0;
 
5544
  else   deta = std::abs(bj->eta - tile->eta_centre) - _tile_half_size_eta;
 
5545
  double dphi = std::abs(bj->phi - tile->phi_centre);
 
5546
  if (dphi > pi) dphi = twopi-dphi;
 
5547
  dphi -= _tile_half_size_phi;
 
5548
  if (dphi < 0) dphi = 0;
 
5549
  return dphi*dphi + deta*deta;
 
5550
}
 
5551
inline void LazyTiling9::_update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, vector<TiledJet *> & jets_for_minheap) {
 
5552
  double dist = _bj_dist(jetI,jetX);
 
5553
  if (dist < jetI->NN_dist) {
 
5554
    if (jetI != jetX) {
 
5555
      jetI->NN_dist = dist;
 
5556
      jetI->NN = jetX;
 
5557
      if (!jetI->minheap_update_needed()) {
 
5558
        jetI->label_minheap_update_needed();
 
5559
        jets_for_minheap.push_back(jetI);
 
5560
      }
 
5561
    }
 
5562
  }
 
5563
  if (dist < jetX->NN_dist) {
 
5564
    if (jetI != jetX) {
 
5565
      jetX->NN_dist = dist;
 
5566
      jetX->NN      = jetI;}
 
5567
  }
 
5568
}
 
5569
inline void LazyTiling9::_set_NN(TiledJet * jetI, 
 
5570
                              vector<TiledJet *> & jets_for_minheap) {
 
5571
  jetI->NN_dist = _R2;
 
5572
  jetI->NN      = NULL;
 
5573
  if (!jetI->minheap_update_needed()) {
 
5574
    jetI->label_minheap_update_needed();
 
5575
    jets_for_minheap.push_back(jetI);}
 
5576
  Tile2 * tile_ptr = &_tiles[jetI->tile_index];
 
5577
    for (Tile2 ** near_tile  = tile_ptr->begin_tiles; 
 
5578
         near_tile != tile_ptr->end_tiles; near_tile++) {
 
5579
      if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
 
5580
      for (TiledJet * jetJ  = (*near_tile)->head; 
 
5581
           jetJ != NULL; jetJ = jetJ->next) {
 
5582
        double dist = _bj_dist(jetI,jetJ);
 
5583
        if (dist < jetI->NN_dist && jetJ != jetI) {
 
5584
          jetI->NN_dist = dist; jetI->NN = jetJ;
 
5585
        }
 
5586
      }
 
5587
    }
 
5588
}
 
5589
void LazyTiling9::run() {
 
5590
  int n = _jets.size();
 
5591
  if (n == 0) return; 
 
5592
  TiledJet * briefjets = new TiledJet[n];
 
5593
  TiledJet * jetA = briefjets, * jetB;
 
5594
  TiledJet oldB = briefjets[0]; 
 
5595
  vector<int> tile_union(3*n_tile_neighbours);
 
5596
  for (int i = 0; i< n; i++) {
 
5597
    _tj_set_jetinfo(jetA, i);
 
5598
    jetA++; // move on to next entry of briefjets
 
5599
  }
 
5600
  TiledJet * head = briefjets; // a nicer way of naming start
 
5601
  vector<Tile2>::iterator tile;
 
5602
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5603
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5604
      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
 
5605
        double dist = _bj_dist_not_periodic(jetA,jetB);
 
5606
        if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5607
        if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5608
      }
 
5609
    }
 
5610
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5611
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
 
5612
    }
 
5613
  }
 
5614
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5615
    if (tile->use_periodic_delta_phi) {
 
5616
      for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
 
5617
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5618
          double dist_to_tile = _distance_to_tile(jetA, *RTile);
 
5619
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
 
5620
          bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
 
5621
          if (relevant_for_jetA || relevant_for_RTile) {
 
5622
            for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
 
5623
              double dist = _bj_dist(jetA,jetB);
 
5624
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5625
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5626
            }
 
5627
          } 
 
5628
        }
 
5629
      }
 
5630
    } else {
 
5631
      for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
 
5632
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5633
          double dist_to_tile = _distance_to_tile(jetA, *RTile);
 
5634
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
 
5635
          bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
 
5636
          if (relevant_for_jetA || relevant_for_RTile) {
 
5637
            for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
 
5638
              double dist = _bj_dist_not_periodic(jetA,jetB);
 
5639
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5640
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5641
            }
 
5642
          } 
 
5643
        }
 
5644
      }
 
5645
    }
 
5646
  }
 
5647
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5648
    tile->max_NN_dist = 0;
 
5649
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5650
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
 
5651
    }
 
5652
  }
 
5653
#ifdef INSTRUMENT2
 
5654
  cout << "intermediate ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
 
5655
#endif // INSTRUMENT2
 
5656
  vector<double> diJs(n);
 
5657
  for (int i = 0; i < n; i++) {
 
5658
    diJs[i] = _bj_diJ(&briefjets[i]);
 
5659
    briefjets[i].label_minheap_update_done();
 
5660
  }
 
5661
  MinHeap minheap(diJs);
 
5662
  vector<TiledJet *> jets_for_minheap;
 
5663
  jets_for_minheap.reserve(n); 
 
5664
  int history_location = n-1;
 
5665
  while (n > 0) {
 
5666
    double diJ_min = minheap.minval() *_invR2;
 
5667
    jetA = head + minheap.minloc();
 
5668
    history_location++;
 
5669
    jetB = jetA->NN;
 
5670
    if (jetB != NULL) {
 
5671
      if (jetA < jetB) {std::swap(jetA,jetB);}
 
5672
      int nn; // new jet index
 
5673
      _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
 
5674
      _bj_remove_from_tiles(jetA);
 
5675
      oldB = * jetB;  // take a copy because we will need it...
 
5676
      _bj_remove_from_tiles(jetB);
 
5677
      _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
 
5678
    } else {
 
5679
      _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
 
5680
      _bj_remove_from_tiles(jetA);
 
5681
    }
 
5682
    minheap.remove(jetA-head);
 
5683
    int n_near_tiles = 0;
 
5684
    if (jetB != NULL) {
 
5685
      Tile2 & jetB_tile = _tiles[jetB->tile_index];
 
5686
      for (Tile2 ** near_tile  = jetB_tile.begin_tiles; 
 
5687
                   near_tile != jetB_tile.end_tiles; near_tile++) {
 
5688
        double dist_to_tile = _distance_to_tile(jetB, *near_tile);
 
5689
        bool relevant_for_jetB  = dist_to_tile <= jetB->NN_dist;
 
5690
        bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
 
5691
        bool relevant = relevant_for_jetB || relevant_for_near_tile;
 
5692
        if (! relevant) continue;
 
5693
        tile_union[n_near_tiles] = *near_tile - & _tiles[0];
 
5694
        (*near_tile)->tagged = true;
 
5695
        n_near_tiles++;
 
5696
        for (TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
 
5697
          if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
 
5698
          _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
 
5699
        }
 
5700
      }
 
5701
    }
 
5702
    int n_done_tiles = n_near_tiles;
 
5703
    _add_untagged_neighbours_to_tile_union_using_max_info(jetA, 
 
5704
                                           tile_union, n_near_tiles);
 
5705
    if (jetB != NULL) {
 
5706
        _add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
 
5707
                                                              tile_union,n_near_tiles);
 
5708
      jetB->label_minheap_update_needed();
 
5709
      jets_for_minheap.push_back(jetB);
 
5710
    }
 
5711
    for (int itile = 0; itile < n_done_tiles; itile++) {
 
5712
      _tiles[tile_union[itile]].tagged = false;
 
5713
    }
 
5714
    for (int itile = n_done_tiles; itile < n_near_tiles; itile++) {
 
5715
      Tile2 * tile_ptr = &_tiles[tile_union[itile]];
 
5716
      tile_ptr->tagged = false;
 
5717
      for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
 
5718
        if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
 
5719
          _set_NN(jetI, jets_for_minheap);
 
5720
        }
 
5721
      }
 
5722
    }
 
5723
    while (jets_for_minheap.size() > 0) {
 
5724
      TiledJet * jetI = jets_for_minheap.back(); 
 
5725
      jets_for_minheap.pop_back();
 
5726
      minheap.update(jetI-head, _bj_diJ(jetI));
 
5727
      jetI->label_minheap_update_done();
 
5728
      Tile2 & tile_I = _tiles[jetI->tile_index];
 
5729
      if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
 
5730
    }
 
5731
    n--;
 
5732
  }
 
5733
  delete[] briefjets;
 
5734
#ifdef INSTRUMENT2
 
5735
  cout << "ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
 
5736
#endif // INSTRUMENT2
 
5737
}
 
5738
FJCORE_END_NAMESPACE
 
5739
#include <iomanip>
 
5740
using namespace std;
 
5741
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
5742
LazyTiling9Alt::LazyTiling9Alt(ClusterSequence & cs) :
 
5743
  _cs(cs), _jets(cs.jets())
 
5744
{
 
5745
  _Rparam = cs.jet_def().R();
 
5746
  _R2 = _Rparam * _Rparam;
 
5747
  _invR2 = 1.0 / _R2;
 
5748
  _initialise_tiles();
 
5749
}
 
5750
void LazyTiling9Alt::_initialise_tiles() {
 
5751
  double default_size = max(0.1,_Rparam);
 
5752
  _tile_size_eta = default_size;
 
5753
  _n_tiles_phi   = max(3,int(floor(twopi/default_size)));
 
5754
  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
 
5755
  _tiles_eta_min = 0.0;
 
5756
  _tiles_eta_max = 0.0;
 
5757
  const double maxrap = 7.0;
 
5758
  for(unsigned int i = 0; i < _jets.size(); i++) {
 
5759
    double eta = _jets[i].rap();
 
5760
    if (abs(eta) < maxrap) {
 
5761
      if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
 
5762
      if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
 
5763
    }
 
5764
  }
 
5765
  _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
 
5766
  _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
 
5767
  _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
 
5768
  _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
 
5769
  _tile_half_size_eta = _tile_size_eta * 0.5;
 
5770
  _tile_half_size_phi = _tile_size_phi * 0.5;
 
5771
  vector<bool> use_periodic_delta_phi(_n_tiles_phi, false);
 
5772
  if (_n_tiles_phi <= 3) {
 
5773
    fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(), true);
 
5774
  } else {
 
5775
    use_periodic_delta_phi[0] = true;
 
5776
    use_periodic_delta_phi[_n_tiles_phi-1] = true;
 
5777
  }
 
5778
  _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
 
5779
  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
 
5780
    for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
 
5781
      Tile * tile = & _tiles[_tile_index(ieta,iphi)];
 
5782
      tile->head = NULL; // first element of tiles points to itself
 
5783
      tile->begin_tiles[0] =  Tile::TileFnPair(tile,&Tile::distance_to_centre);
 
5784
      Tile::TileFnPair * pptile = & (tile->begin_tiles[0]);
 
5785
      pptile++;
 
5786
      tile->surrounding_tiles = pptile;
 
5787
      if (ieta > _tiles_ieta_min) {
 
5788
        // with the itile subroutine, we can safely run tiles from
 
5789
        // idphi=-1 to idphi=+1, because it takes care of
 
5790
        // negative and positive boundaries
 
5791
        //for (int idphi = -1; idphi <=+1; idphi++) {
 
5792
        *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta-1,iphi-1)],
 
5793
                                   &Tile::distance_to_left_bottom);
 
5794
        pptile++;
 
5795
        *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta-1,iphi)],
 
5796
                                   &Tile::distance_to_left);
 
5797
        pptile++;
 
5798
        *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta-1,iphi+1)],
 
5799
                                   &Tile::distance_to_left_top);
 
5800
        pptile++;
 
5801
      }
 
5802
      *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta,iphi-1)], 
 
5803
                                 &Tile::distance_to_bottom);
 
5804
      pptile++;
 
5805
      tile->RH_tiles = pptile;
 
5806
      *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta,iphi+1)], 
 
5807
                                 &Tile::distance_to_top);
 
5808
      pptile++;
 
5809
      if (ieta < _tiles_ieta_max) {
 
5810
        //for (int idphi = -1; idphi <= +1; idphi++) {
 
5811
        //  *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
 
5812
        //  pptile++;
 
5813
        //}     
 
5814
        *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta+1,iphi-1)],
 
5815
                                   &Tile::distance_to_right_bottom);
 
5816
        pptile++;
 
5817
        *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta+1,iphi)],
 
5818
                                   &Tile::distance_to_right);
 
5819
        pptile++;
 
5820
        *pptile = Tile::TileFnPair(& _tiles[_tile_index(ieta+1,iphi+1)],
 
5821
                                   &Tile::distance_to_right_top);
 
5822
        pptile++;
 
5823
      }
 
5824
      tile->end_tiles = pptile;
 
5825
      tile->tagged = false;
 
5826
      tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
 
5827
      tile->max_NN_dist = 0;
 
5828
      tile->eta_min = ieta*_tile_size_eta;
 
5829
      tile->eta_max = (ieta+1)*_tile_size_eta;
 
5830
      tile->phi_min = iphi*_tile_size_phi;
 
5831
      tile->phi_max = (iphi+1)*_tile_size_phi;
 
5832
    }
 
5833
  }
 
5834
}
 
5835
int LazyTiling9Alt::_tile_index(const double eta, const double phi) const {
 
5836
  int ieta, iphi;
 
5837
  if      (eta <= _tiles_eta_min) {ieta = 0;}
 
5838
  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
 
5839
  else {
 
5840
    ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
 
5841
    if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
 
5842
      ieta = _tiles_ieta_max-_tiles_ieta_min;} 
 
5843
  }
 
5844
  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
 
5845
  return (iphi + ieta * _n_tiles_phi);
 
5846
}
 
5847
inline void LazyTiling9Alt::_tj_set_jetinfo( TiledJet * const jet,
 
5848
                                              const int _jets_index) {
 
5849
  _bj_set_jetinfo<>(jet, _jets_index);
 
5850
  jet->tile_index = _tile_index(jet->eta, jet->phi);
 
5851
  Tile * tile = &_tiles[jet->tile_index];
 
5852
  jet->previous   = NULL;
 
5853
  jet->next       = tile->head;
 
5854
  if (jet->next != NULL) {jet->next->previous = jet;}
 
5855
  tile->head      = jet;
 
5856
}
 
5857
void LazyTiling9Alt::_bj_remove_from_tiles(TiledJet * const jet) {
 
5858
  Tile * tile = & _tiles[jet->tile_index];
 
5859
  if (jet->previous == NULL) {
 
5860
    tile->head = jet->next;
 
5861
  } else {
 
5862
    jet->previous->next = jet->next;
 
5863
  }
 
5864
  if (jet->next != NULL) {
 
5865
    jet->next->previous = jet->previous;
 
5866
  }
 
5867
}
 
5868
void LazyTiling9Alt::_print_tiles(TiledJet * briefjets ) const {
 
5869
  for (vector<Tile>::const_iterator tile = _tiles.begin(); 
 
5870
       tile < _tiles.end(); tile++) {
 
5871
    cout << "Tile " << tile - _tiles.begin()<<" = ";
 
5872
    vector<int> list;
 
5873
    for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
 
5874
      list.push_back(jetI-briefjets);
 
5875
    }
 
5876
    sort(list.begin(),list.end());
 
5877
    for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
 
5878
    cout <<"\n";
 
5879
  }
 
5880
}
 
5881
void LazyTiling9Alt::_add_neighbours_to_tile_union(const int tile_index, 
 
5882
               vector<int> & tile_union, int & n_near_tiles) const {
 
5883
  for (Tile::TileFnPair const * near_tile = _tiles[tile_index].begin_tiles; 
 
5884
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
 
5885
    tile_union[n_near_tiles] = near_tile->first - & _tiles[0];
 
5886
    n_near_tiles++;
 
5887
  }
 
5888
}
 
5889
inline void LazyTiling9Alt::_add_untagged_neighbours_to_tile_union(
 
5890
               const int tile_index, 
 
5891
               vector<int> & tile_union, int & n_near_tiles)  {
 
5892
  for (Tile::TileFnPair * near_tile = _tiles[tile_index].begin_tiles; 
 
5893
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
 
5894
    if (! (near_tile->first)->tagged) {
 
5895
      (near_tile->first)->tagged = true;
 
5896
      tile_union[n_near_tiles] = near_tile->first - & _tiles[0];
 
5897
      n_near_tiles++;
 
5898
    }
 
5899
  }
 
5900
}
 
5901
inline void LazyTiling9Alt::_add_untagged_neighbours_to_tile_union_using_max_info(
 
5902
               const TiledJet * jet, 
 
5903
               vector<int> & tile_union, int & n_near_tiles)  {
 
5904
  Tile & tile = _tiles[jet->tile_index];
 
5905
  for (Tile::TileFnPair * near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
 
5906
    if ((near_tile->first)->tagged) continue;
 
5907
    double dist = (tile.*(near_tile->second))(jet) - tile_edge_security_margin;
 
5908
    if (dist > (near_tile->first)->max_NN_dist) continue;
 
5909
    (near_tile->first)->tagged = true;
 
5910
    tile_union[n_near_tiles] = near_tile->first - & _tiles[0];
 
5911
    n_near_tiles++;
 
5912
  }
 
5913
}
 
5914
ostream & operator<<(ostream & ostr, const TiledJet & jet) {
 
5915
  ostr << "j" << setw(3) << jet._jets_index << ":pt2,rap,phi=" ; ostr.flush();
 
5916
  ostr     << jet.kt2 << ","; ostr.flush();
 
5917
  ostr     << jet.eta << ","; ostr.flush();
 
5918
  ostr     << jet.phi; ostr.flush();
 
5919
  ostr     << ", tile=" << jet.tile_index; ostr.flush();
 
5920
  return ostr;
 
5921
}
 
5922
inline void LazyTiling9Alt::_update_jetX_jetI_NN(TiledJet * jetX, TiledJet * jetI, vector<TiledJet *> & jets_for_minheap) {
 
5923
  double dist = _bj_dist(jetI,jetX);
 
5924
  if (dist < jetI->NN_dist) {
 
5925
    if (jetI != jetX) {
 
5926
      jetI->NN_dist = dist;
 
5927
      jetI->NN = jetX;
 
5928
      if (!jetI->minheap_update_needed()) {
 
5929
        jetI->label_minheap_update_needed();
 
5930
        jets_for_minheap.push_back(jetI);
 
5931
      }
 
5932
    }
 
5933
  }
 
5934
  if (dist < jetX->NN_dist) {
 
5935
    if (jetI != jetX) {
 
5936
      jetX->NN_dist = dist;
 
5937
      jetX->NN      = jetI;}
 
5938
  }
 
5939
}
 
5940
inline void LazyTiling9Alt::_set_NN(TiledJet * jetI, 
 
5941
                            vector<TiledJet *> & jets_for_minheap) {
 
5942
  jetI->NN_dist = _R2;
 
5943
  jetI->NN      = NULL;
 
5944
  if (!jetI->minheap_update_needed()) {
 
5945
    jetI->label_minheap_update_needed();
 
5946
    jets_for_minheap.push_back(jetI);}
 
5947
  Tile * tile_ptr = &_tiles[jetI->tile_index];
 
5948
    for (Tile::TileFnPair * near_tile  = tile_ptr->begin_tiles; 
 
5949
         near_tile != tile_ptr->end_tiles; near_tile++) {
 
5950
      if (jetI->NN_dist < (tile_ptr->*(near_tile->second))(jetI)) continue;
 
5951
      for (TiledJet * jetJ  = (near_tile->first)->head; 
 
5952
           jetJ != NULL; jetJ = jetJ->next) {
 
5953
        double dist = _bj_dist(jetI,jetJ);
 
5954
        if (dist < jetI->NN_dist && jetJ != jetI) {
 
5955
          jetI->NN_dist = dist; jetI->NN = jetJ;
 
5956
        }
 
5957
      }
 
5958
    }
 
5959
}
 
5960
void LazyTiling9Alt::run() {
 
5961
  int n = _jets.size();
 
5962
  TiledJet * briefjets = new TiledJet[n];
 
5963
  TiledJet * jetA = briefjets, * jetB;
 
5964
  TiledJet oldB;
 
5965
  vector<int> tile_union(3*n_tile_neighbours);
 
5966
  for (int i = 0; i< n; i++) {
 
5967
    _tj_set_jetinfo(jetA, i);
 
5968
    jetA++; // move on to next entry of briefjets
 
5969
  }
 
5970
  TiledJet * head = briefjets; // a nicer way of naming start
 
5971
  vector<Tile>::iterator tile;
 
5972
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5973
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5974
      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
 
5975
        double dist = _bj_dist_not_periodic(jetA,jetB);
 
5976
        if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5977
        if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5978
      }
 
5979
    }
 
5980
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5981
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
 
5982
    }
 
5983
  }
 
5984
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
5985
    if (tile->use_periodic_delta_phi) {
 
5986
      for (Tile::TileFnPair * RTileFnPair = tile->RH_tiles; 
 
5987
           RTileFnPair != tile->end_tiles; RTileFnPair++) {
 
5988
        Tile *RTile = RTileFnPair->first;
 
5989
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
5990
          double dist_to_tile = ((*tile).*(RTileFnPair->second))(jetA);
 
5991
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
 
5992
          bool relevant_for_RTile = dist_to_tile <= RTile->max_NN_dist;
 
5993
          if (relevant_for_jetA || relevant_for_RTile) {
 
5994
            for (jetB = RTile->head; jetB != NULL; jetB = jetB->next) {
 
5995
              double dist = _bj_dist(jetA,jetB);
 
5996
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
5997
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
5998
            }
 
5999
          } 
 
6000
        }
 
6001
      }
 
6002
    } else {
 
6003
      for (Tile::TileFnPair* RTileFnPair = tile->RH_tiles;
 
6004
           RTileFnPair != tile->end_tiles; RTileFnPair++) {
 
6005
        Tile *RTile = RTileFnPair->first;
 
6006
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
6007
          double dist_to_tile = ((*tile).*(RTileFnPair->second))(jetA);
 
6008
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
 
6009
          bool relevant_for_RTile = dist_to_tile <= RTile->max_NN_dist;
 
6010
          if (relevant_for_jetA || relevant_for_RTile) {
 
6011
            for (jetB = RTile->head; jetB != NULL; jetB = jetB->next) {
 
6012
              double dist = _bj_dist_not_periodic(jetA,jetB);
 
6013
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
 
6014
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
 
6015
            }
 
6016
          } 
 
6017
        }
 
6018
      }
 
6019
    }
 
6020
  }
 
6021
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
 
6022
    tile->max_NN_dist = 0;
 
6023
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
 
6024
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
 
6025
    }
 
6026
  }
 
6027
  vector<double> diJs(n);
 
6028
  for (int i = 0; i < n; i++) {
 
6029
    diJs[i] = _bj_diJ(&briefjets[i]);
 
6030
    briefjets[i].label_minheap_update_done();
 
6031
  }
 
6032
  MinHeap minheap(diJs);
 
6033
  vector<TiledJet *> jets_for_minheap;
 
6034
  jets_for_minheap.reserve(n); 
 
6035
  int history_location = n-1;
 
6036
  while (n > 0) {
 
6037
    double diJ_min = minheap.minval() *_invR2;
 
6038
    jetA = head + minheap.minloc();
 
6039
    history_location++;
 
6040
    jetB = jetA->NN;
 
6041
    if (jetB != NULL) {
 
6042
      if (jetA < jetB) {std::swap(jetA,jetB);}
 
6043
      int nn; // new jet index
 
6044
      _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
 
6045
      _bj_remove_from_tiles(jetA);
 
6046
      oldB = * jetB;  // take a copy because we will need it...
 
6047
      _bj_remove_from_tiles(jetB);
 
6048
      _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
 
6049
    } else {
 
6050
      _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
 
6051
      _bj_remove_from_tiles(jetA);
 
6052
    }
 
6053
    minheap.remove(jetA-head);
 
6054
    int n_near_tiles = 0;
 
6055
    _add_untagged_neighbours_to_tile_union_using_max_info(jetA, 
 
6056
                                           tile_union, n_near_tiles);
 
6057
    if (jetB != NULL) {
 
6058
        _add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
 
6059
                                                              tile_union,n_near_tiles);
 
6060
      jetB->label_minheap_update_needed();
 
6061
      jets_for_minheap.push_back(jetB);
 
6062
    }
 
6063
    if (jetB != NULL) {
 
6064
      Tile & jetB_tile = _tiles[jetB->tile_index];
 
6065
      for (Tile::TileFnPair * near_tile_fn_pair  = jetB_tile.begin_tiles; 
 
6066
                   near_tile_fn_pair != jetB_tile.end_tiles; near_tile_fn_pair++) {
 
6067
        Tile * near_tile = near_tile_fn_pair->first;
 
6068
        double dist_to_tile = (jetB_tile.*(near_tile_fn_pair->second))(jetB);
 
6069
        bool relevant_for_jetB  = dist_to_tile <= jetB->NN_dist;
 
6070
        bool relevant_for_near_tile = dist_to_tile <= near_tile->max_NN_dist;
 
6071
        bool relevant = relevant_for_jetB || relevant_for_near_tile;
 
6072
        if (relevant) {
 
6073
          if (near_tile->tagged) {
 
6074
            for (TiledJet * jetI = near_tile->head; jetI != NULL; jetI = jetI->next) {
 
6075
              if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
 
6076
              _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
 
6077
            }
 
6078
          near_tile->tagged = false;
 
6079
          } else {
 
6080
            for (TiledJet * jetI = near_tile->head; jetI != NULL; jetI = jetI->next) {
 
6081
              _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
 
6082
            }
 
6083
          }
 
6084
        }
 
6085
        //     // -- Keep this old inline code for later speed tests
 
6086
      }
 
6087
    }
 
6088
    for (int itile = 0; itile < n_near_tiles; itile++) {
 
6089
      Tile * tile_ptr = &_tiles[tile_union[itile]];
 
6090
      if (!tile_ptr->tagged) continue; // because earlier loop may have undone the tag
 
6091
      tile_ptr->tagged = false;
 
6092
      for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
 
6093
        if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
 
6094
          _set_NN(jetI, jets_for_minheap);
 
6095
        }
 
6096
      }
 
6097
    }
 
6098
    while (jets_for_minheap.size() > 0) {
 
6099
      TiledJet * jetI = jets_for_minheap.back(); 
 
6100
      jets_for_minheap.pop_back();
 
6101
      minheap.update(jetI-head, _bj_diJ(jetI));
 
6102
      jetI->label_minheap_update_done();
 
6103
      Tile & tile_I = _tiles[jetI->tile_index];
 
6104
      if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
 
6105
    }
 
6106
    n--;
 
6107
  }
 
6108
  delete[] briefjets;
 
6109
}
 
6110
FJCORE_END_NAMESPACE
 
6111
#include <iomanip>
 
6112
#include <limits>
 
6113
#include <cmath>
 
6114
using namespace std;
 
6115
FJCORE_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
 
6116
TilingExtent::TilingExtent(ClusterSequence & cs) {
 
6117
  _determine_rapidity_extent(cs.jets());
 
6118
}
 
6119
void TilingExtent::_determine_rapidity_extent(const vector<PseudoJet> & particles) {
 
6120
  int nrap = 20; 
 
6121
  int nbins = 2*nrap;
 
6122
  vector<double> counts(nbins, 0);
 
6123
  _minrap =  numeric_limits<double>::max();
 
6124
  _maxrap = -numeric_limits<double>::max();
 
6125
  int ibin;
 
6126
  for (unsigned i = 0; i < particles.size(); i++) {
 
6127
    if (particles[i].E() == abs(particles[i].pz())) continue;
 
6128
    double rap = particles[i].rap();
 
6129
    if (rap < _minrap) _minrap = rap;
 
6130
    if (rap > _maxrap) _maxrap = rap;
 
6131
    ibin = int(rap+nrap); 
 
6132
    if (ibin < 0) ibin = 0;
 
6133
    if (ibin >= nbins) ibin = nbins - 1;
 
6134
    counts[ibin]++;
 
6135
  }
 
6136
  double max_in_bin = 0;
 
6137
  for (ibin = 0; ibin < nbins; ibin++) {
 
6138
    if (max_in_bin < counts[ibin]) max_in_bin = counts[ibin];
 
6139
  }
 
6140
  const double allowed_max_fraction = 0.25;
 
6141
  const double min_multiplicity = 4;
 
6142
  double allowed_max_cumul = floor(max(max_in_bin * allowed_max_fraction, min_multiplicity));
 
6143
  if (allowed_max_cumul > max_in_bin) allowed_max_cumul = max_in_bin;
 
6144
  double cumul_lo = 0;
 
6145
  _cumul2 = 0;
 
6146
  for (ibin = 0; ibin < nbins; ibin++) {
 
6147
    cumul_lo += counts[ibin];
 
6148
    if (cumul_lo >= allowed_max_cumul) {
 
6149
      double y = ibin-nrap;
 
6150
      if (y > _minrap) _minrap = y;
 
6151
      break;
 
6152
    }
 
6153
  }
 
6154
  assert(ibin != nbins); // internal consistency check that you found a bin
 
6155
  _cumul2 += cumul_lo*cumul_lo;
 
6156
  int ibin_lo = ibin;
 
6157
  double cumul_hi = 0;
 
6158
  for (ibin = nbins-1; ibin >= 0; ibin--) {
 
6159
    cumul_hi += counts[ibin];
 
6160
    if (cumul_hi >= allowed_max_cumul) {
 
6161
      double y = ibin-nrap+1; // +1 here is the rapidity bin width
 
6162
      if (y < _maxrap) _maxrap = y;
 
6163
      break;
 
6164
    }
 
6165
  }
 
6166
  assert(ibin >= 0); // internal consistency check that you found a bin
 
6167
  int ibin_hi = ibin;
 
6168
  assert(ibin_hi >= ibin_lo); 
 
6169
  if (ibin_hi == ibin_lo) {
 
6170
    _cumul2 = pow(double(cumul_lo + cumul_hi - counts[ibin_hi]), 2);
 
6171
  } else {
 
6172
    _cumul2 += cumul_hi*cumul_hi;
 
6173
    for (ibin = ibin_lo+1; ibin < ibin_hi; ibin++) {
 
6174
      _cumul2 += counts[ibin]*counts[ibin];
 
6175
    }
 
6176
  }
 
6177
}
 
6178
FJCORE_END_NAMESPACE