Fork me on GitHub

Ignore:
Timestamp:
Dec 9, 2014, 1:27:13 PM (10 years ago)
Author:
Michele <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
37deb3b, 9e991f8
Parents:
f6b6ee7 (diff), e7e90df (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'TestFastJet310b1'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • external/fastjet/contribs/Nsubjettiness/AxesFinder.hh

    rf6b6ee7 r49234af  
    55//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
    66//
     7//  $Id: AxesFinder.hh 678 2014-06-12 20:43:03Z jthaler $
    78//----------------------------------------------------------------------
    89// This file is part of FastJet contrib.
    … …  
    5354class AxesFinder {
    5455   
    55 protected:
    56    AxesFinder* _startingFinder; // storing a possible starting finder if needed
    57    std::vector<fastjet::PseudoJet> _seedAxes;
    58    
    59    AxesFinder(AxesFinder* startingFinder = NULL) : _startingFinder(startingFinder) {}
    60    
    61 public:
    62    virtual ~AxesFinder(){
    63       if (_startingFinder) delete _startingFinder;  //TODO: Convert to smart pointers to avoid this.
    64    }
    65    
    66    // Allow setting of seedAxes from a starting finder
    67    std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector<fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
    68       if (_startingFinder) {
    69          _seedAxes = _startingFinder->getAxes(n_jets,inputs,currentAxes);
    70          return getBetterAxes(n_jets,inputs,_seedAxes);
    71       } else {
    72          _seedAxes = getBetterAxes(n_jets,inputs,currentAxes);
    73          return _seedAxes;
    74       }
    75    }
    76    
    77    // say what the current seed axes are
    78    std::vector<fastjet::PseudoJet> seedAxes() const {
    79       return _seedAxes;
    80    }
    81    
    82    // This function should be overloaded, and updates the seedAxes
    83    virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector<fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& seedAxes) = 0;
    84    
    85 };
    86 
     56public:
     57   
     58   // This function should be overloaded, and updates the seedAxes to return new axes
     59   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
     60                                                   const std::vector<fastjet::PseudoJet>& inputs,
     61                                                   const std::vector<fastjet::PseudoJet>& seedAxes) const = 0;
     62   // convenient shorthand for squaring
     63   static inline double sq(double x) {return x*x;}
     64
     65   //virtual destructor
     66   virtual ~AxesFinder(){}
     67   
     68};
     69
     70   
    8771//------------------------------------------------------------------------
    8872/// \class AxesFinderFromExclusiveJetDefinition
    … …  
    9074// with different jet algorithms.
    9175class AxesFinderFromExclusiveJetDefinition : public AxesFinder {
    92 
    93    private:
    94       fastjet::JetDefinition _def;
    95    
    96    public:
    97       AxesFinderFromExclusiveJetDefinition(fastjet::JetDefinition def) : _def(def) {}
    98      
    99       virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
    100          fastjet::ClusterSequence jet_clust_seq(inputs, _def);
    101          return jet_clust_seq.exclusive_jets(n_jets);
    102       }
     76   
     77public:
     78   AxesFinderFromExclusiveJetDefinition(fastjet::JetDefinition def)
     79   : _def(def) {}
     80   
     81   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
     82                                                   const std::vector <fastjet::PseudoJet> & inputs,
     83                                                   const std::vector<fastjet::PseudoJet>& /*seedAxes*/) const {
     84      fastjet::ClusterSequence jet_clust_seq(inputs, _def);
     85      return jet_clust_seq.exclusive_jets(n_jets);
     86   }
     87   
     88private:
     89   fastjet::JetDefinition _def;
     90
    10391};
    10492
    … …  
    10896// winner take all recombination scheme.
    10997class AxesFinderFromWTA_KT : public AxesFinderFromExclusiveJetDefinition {
    110    private:
    111       const WinnerTakeAllRecombiner *recomb;
    112    public:
    113       AxesFinderFromWTA_KT() : AxesFinderFromExclusiveJetDefinition(
    114          fastjet::JetDefinition(fastjet::kt_algorithm,
    115          fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
    116          recomb = new WinnerTakeAllRecombiner(),
    117          fastjet::Best)) {}
    118       ~AxesFinderFromWTA_KT() {delete recomb;}
    119    };
     98
     99public:
     100   AxesFinderFromWTA_KT()
     101   : AxesFinderFromExclusiveJetDefinition(
     102      fastjet::JetDefinition(fastjet::kt_algorithm,
     103      fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
     104      &_recomb,
     105      fastjet::Best)) {}
     106   
     107private:
     108   const WinnerTakeAllRecombiner _recomb;
     109
     110};
    120111   
    121112//------------------------------------------------------------------------
    … …  
    124115// winner take all recombination scheme.
    125116class AxesFinderFromWTA_CA : public AxesFinderFromExclusiveJetDefinition {
    126    private:
    127       const WinnerTakeAllRecombiner *recomb;
    128    public:
    129       AxesFinderFromWTA_CA() : AxesFinderFromExclusiveJetDefinition(
    130          fastjet::JetDefinition(fastjet::cambridge_algorithm,
    131          fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
    132          recomb = new WinnerTakeAllRecombiner(),
    133          fastjet::Best)) {}
    134       ~AxesFinderFromWTA_CA() {delete recomb;}
    135 };
    136 
    137 //  The following classes are for testing, and are commented out for initial release
    138 //
    139 ////------------------------------------------------------------------------
    140 ///// \class AxesFinderFromWTA2_KT
    141 //// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a
    142 //// winner take all recombination scheme with alpha = 2.
    143 //class AxesFinderFromWTA2_KT : public AxesFinderFromExclusiveJetDefinition {
    144 //   private:
    145 //      const WinnerTakeAllRecombiner *recomb;
    146 //   public:
    147 //      AxesFinderFromWTA2_KT() : AxesFinderFromExclusiveJetDefinition(
    148 //         fastjet::JetDefinition(fastjet::kt_algorithm,
    149 //         fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
    150 //         recomb = new WinnerTakeAllRecombiner(2), // uses alpha = 2 here
    151 //         fastjet::Best)) {}
    152 //      ~AxesFinderFromWTA2_KT() {delete recomb;}
    153 //   };
    154 //   
    155 ////------------------------------------------------------------------------
    156 ///// \class AxesFinderFromWTA2_CA
    157 //// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a
    158 //// winner take all recombination scheme with alpha = 2.
    159 //class AxesFinderFromWTA2_CA : public AxesFinderFromExclusiveJetDefinition {
    160 //   private:
    161 //      const WinnerTakeAllRecombiner *recomb;
    162 //   public:
    163 //      AxesFinderFromWTA2_CA() : AxesFinderFromExclusiveJetDefinition(
    164 //         fastjet::JetDefinition(fastjet::cambridge_algorithm,
    165 //         fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
    166 //         recomb = new WinnerTakeAllRecombiner(2), //uses alpha = 2 here
    167 //         fastjet::Best)) {}
    168 //      ~AxesFinderFromWTA2_CA() {delete recomb;}
    169 //};
     117public:
     118   AxesFinderFromWTA_CA()
     119   : AxesFinderFromExclusiveJetDefinition(
     120      fastjet::JetDefinition(fastjet::cambridge_algorithm,
     121      fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
     122      &_recomb,
     123      fastjet::Best)) {}
     124   
     125private:
     126   const WinnerTakeAllRecombiner _recomb;
     127};
     128
    170129
    171130//------------------------------------------------------------------------
    … …  
    174133// E_scheme recombination.
    175134class AxesFinderFromKT : public AxesFinderFromExclusiveJetDefinition {
    176    public:
    177       AxesFinderFromKT() : AxesFinderFromExclusiveJetDefinition(
    178          fastjet::JetDefinition(fastjet::kt_algorithm,
    179          fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
    180          fastjet::E_scheme,
    181          fastjet::Best)) {}
     135public:
     136   AxesFinderFromKT()
     137   : AxesFinderFromExclusiveJetDefinition(
     138      fastjet::JetDefinition(fastjet::kt_algorithm,
     139      fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
     140      fastjet::E_scheme,
     141      fastjet::Best)) {}
    182142};
    183143
    … …  
    187147// E_scheme recombination.
    188148class AxesFinderFromCA : public AxesFinderFromExclusiveJetDefinition {
    189    public:
    190       AxesFinderFromCA() : AxesFinderFromExclusiveJetDefinition(
    191          fastjet::JetDefinition(fastjet::cambridge_algorithm,
    192                                 fastjet::JetDefinition::max_allowable_R,  //maximum jet radius constant
    193          fastjet::E_scheme,
    194          fastjet::Best)) {}
     149public:
     150   AxesFinderFromCA()
     151   : AxesFinderFromExclusiveJetDefinition(
     152      fastjet::JetDefinition(fastjet::cambridge_algorithm,
     153                             fastjet::JetDefinition::max_allowable_R,  //maximum jet radius constant
     154      fastjet::E_scheme,
     155      fastjet::Best)) {}
    195156};
    196157
    … …  
    201162// This can be implemented with different jet algorithms.
    202163class AxesFinderFromHardestJetDefinition : public AxesFinder {
    203 
    204    private:
    205       fastjet::JetDefinition _def;
    206    
    207    public:
    208       AxesFinderFromHardestJetDefinition(fastjet::JetDefinition def) : _def(def) {}
    209      
    210       virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
    211          fastjet::ClusterSequence jet_clust_seq(inputs, _def);
    212          std::vector<fastjet::PseudoJet> myJets = sorted_by_pt(jet_clust_seq.inclusive_jets());
    213          myJets.resize(n_jets);  // only keep n hardest
    214          return myJets;
    215       }     
     164public:
     165   AxesFinderFromHardestJetDefinition(fastjet::JetDefinition def)
     166   : _def(def) {}
     167   
     168   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
     169                                                   const std::vector <fastjet::PseudoJet> & inputs,
     170                                                   const std::vector<fastjet::PseudoJet>& /*seedAxes*/) const {
     171      fastjet::ClusterSequence jet_clust_seq(inputs, _def);
     172      std::vector<fastjet::PseudoJet> myJets = sorted_by_pt(jet_clust_seq.inclusive_jets());
     173      myJets.resize(n_jets);  // only keep n hardest
     174      return myJets;
     175   }
     176   
     177private:
     178   fastjet::JetDefinition _def;
    216179};
    217180
    … …  
    221184// to an anti kT algorithm and E_scheme.
    222185class AxesFinderFromAntiKT : public AxesFinderFromHardestJetDefinition {
    223    public:
    224       AxesFinderFromAntiKT(double R0) : AxesFinderFromHardestJetDefinition(fastjet::JetDefinition(fastjet::antikt_algorithm,R0,fastjet::E_scheme,fastjet::Best)) {}
     186public:
     187   AxesFinderFromAntiKT(double R0)
     188   : AxesFinderFromHardestJetDefinition(
     189      fastjet::JetDefinition(fastjet::antikt_algorithm,
     190                             R0,fastjet::E_scheme,fastjet::Best)) {}
    225191};
    226192
    … …  
    231197class AxesFinderFromUserInput : public AxesFinder {
    232198
    233    public:
    234       AxesFinderFromUserInput() {}
    235      
    236       virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
    237          assert(currentAxes.size() == (unsigned int) n_jets);
    238          return currentAxes;
    239       }
     199public:
     200   AxesFinderFromUserInput() {}
     201   
     202   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & /*inputs*/, const std::vector<fastjet::PseudoJet>& currentAxes) const {
     203      assert(currentAxes.size() == (unsigned int) n_jets);
     204      (void)(n_jets);  // adding this line to fix unused-parameter warning
     205      return currentAxes;
     206   }
    240207};
    241208
    … …  
    249216class AxesFinderFromOnePassMinimization : public AxesFinder {
    250217
    251    private:
    252       double _precision;  // Desired precision in axes alignment
    253       int _halt;  // maximum number of steps per iteration
    254      
    255       double _beta;
    256       double _Rcutoff;
    257      
    258       DefaultUnnormalizedMeasure _measureFunction;
    259    
    260    public:
    261 
    262       // From a startingFinder, try to minimize the unnormalized_measure
    263       AxesFinderFromOnePassMinimization(AxesFinder* startingFinder, double beta, double Rcutoff)
    264          : AxesFinder(startingFinder),
    265            _precision(0.0001), //hard coded for now
    266            _halt(1000), //hard coded for now
    267            _beta(beta),
    268            _Rcutoff(Rcutoff),
    269            _measureFunction(beta, Rcutoff)
    270            {}
    271    
    272       virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);
    273 
    274       template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
    275                                   const std::vector <fastjet::PseudoJet> & inputJets);
    276    
    277       std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
    278                                       const std::vector <fastjet::PseudoJet> & inputJets);
     218public:
     219
     220   // From a startingFinder, try to minimize the unnormalized_measure
     221   AxesFinderFromOnePassMinimization(double beta, double Rcutoff)
     222      : _precision(0.0001), //hard coded for now
     223        _halt(1000), //hard coded for now
     224        _beta(beta),
     225        _Rcutoff(Rcutoff),
     226        _measureFunction(beta, Rcutoff)
     227        {}
     228
     229   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets,
     230                                                   const std::vector <fastjet::PseudoJet> & inputJets,
     231                                                   const std::vector<fastjet::PseudoJet>& currentAxes) const;
     232   
     233private:
     234   double _precision;  // Desired precision in axes alignment
     235   int _halt;  // maximum number of steps per iteration
     236   
     237   double _beta;
     238   double _Rcutoff;
     239   
     240   DefaultUnnormalizedMeasureFunction _measureFunction;
     241   
     242   template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
     243                                                              const std::vector <fastjet::PseudoJet> & inputJets) const;
     244   
     245   std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
     246                                         const std::vector <fastjet::PseudoJet> & inputJets) const;
    279247
    280248};
    … …  
    288256class AxesFinderFromKmeansMinimization : public AxesFinder{
    289257
    290    private:
    291       int _n_iterations;   // Number of iterations to run  (0 for no minimization, 1 for one-pass, >>1 for global minimum)
    292       double _noise_range; // noise range for random initialization
    293    
    294       DefaultUnnormalizedMeasure _measureFunction; //function to test whether minimum is reached
    295    
    296       AxesFinderFromOnePassMinimization _onePassFinder;  //one pass finder for minimization
    297 
    298       PseudoJet jiggle(const PseudoJet& axis);
    299    
    300    public:
    301       AxesFinderFromKmeansMinimization(AxesFinder *startingFinder, double beta, double Rcutoff, int n_iterations) :
    302          AxesFinder(startingFinder),
    303          _n_iterations(n_iterations),
    304          _noise_range(1.0), // hard coded for the time being
    305          _measureFunction(beta, Rcutoff),
    306          _onePassFinder(NULL, beta, Rcutoff)
    307          {}
    308 
    309       virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);
    310 
     258public:
     259   AxesFinderFromKmeansMinimization(double beta, double Rcutoff, int n_iterations)
     260   :  _n_iterations(n_iterations),
     261      _noise_range(1.0), // hard coded for the time being
     262      _measureFunction(beta, Rcutoff),
     263      _onePassFinder(beta, Rcutoff)
     264      {}
     265
     266   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes) const;
     267   
     268private:
     269   int _n_iterations;   // Number of iterations to run  (0 for no minimization, 1 for one-pass, >>1 for global minimum)
     270   double _noise_range; // noise range for random initialization
     271   
     272   DefaultUnnormalizedMeasureFunction _measureFunction; //function to test whether minimum is reached
     273   
     274   AxesFinderFromOnePassMinimization _onePassFinder;  //one pass finder that is repeatedly called
     275   
     276   PseudoJet jiggle(const PseudoJet& axis) const;
    311277};
    312278
    … …  
    317283class AxesFinderFromGeometricMinimization : public AxesFinder {
    318284
    319    private:
    320       MeasureFunction* _function;
    321       double _Rcutoff;
    322       double _nAttempts;
    323       double _accuracy;
    324 
    325    
    326    public:
    327       AxesFinderFromGeometricMinimization(AxesFinder* startingFinder, double beta, double Rcutoff) : AxesFinder(startingFinder), _Rcutoff(Rcutoff) {
    328          if (beta != 2.0) {
    329             std::cerr << "Geometric minimization is currently only defined for beta = 2.0." << std::endl;
    330             exit(1);
    331          }
    332          
    333          _nAttempts = 100;
    334          _accuracy = 0.000000001;
    335          _function = new GeometricMeasure(beta,_Rcutoff);
     285public:
     286   AxesFinderFromGeometricMinimization(double beta, double Rcutoff)
     287   :  _nAttempts(100),
     288      _accuracy(0.000000001),
     289      _function(beta,Rcutoff)
     290   {
     291      if (beta != 2.0) {
     292         throw Error("Geometric minimization is currently only defined for beta = 2.0.");
    336293      }
    337 
    338       ~AxesFinderFromGeometricMinimization() {
    339          delete _function;
    340       }
    341    
    342       virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes);
     294   }
     295
     296   virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes) const;
     297
     298private:
     299   double _nAttempts;
     300   double _accuracy;
     301   GeometricMeasureFunction _function;
     302
     303
    343304};
    344305
    … …  
    348309// in order to better facilitate calculations.
    349310class LightLikeAxis {
    350 private:
    351    double _rap, _phi, _weight, _mom;
    352    
    353    double DistanceSq(double rap2, double phi2) const {
    354       double rap1 = _rap;
    355       double phi1 = _phi;
    356      
    357       double distRap = rap1-rap2;
    358       double distPhi = std::fabs(phi1-phi2);
    359       if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
    360       return sq(distRap) + sq(distPhi);
    361    }
    362    
    363    double Distance(double rap2, double phi2) const {
    364       return std::sqrt(DistanceSq(rap2,phi2));
    365    }
    366    
    367    
     311
    368312public:
    369313   LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
    … …  
    399343   }
    400344
     345private:
     346   double _rap, _phi, _weight, _mom;
     347   
     348   double DistanceSq(double rap2, double phi2) const {
     349      double rap1 = _rap;
     350      double phi1 = _phi;
     351     
     352      double distRap = rap1-rap2;
     353      double distPhi = std::fabs(phi1-phi2);
     354      if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
     355      return distRap*distRap + distPhi*distPhi;
     356   }
     357   
     358   double Distance(double rap2, double phi2) const {
     359      return std::sqrt(DistanceSq(rap2,phi2));
     360   }
     361   
     362   
    401363};
    402364
Note: See TracChangeset for help on using the changeset viewer.