Fork me on GitHub

Ignore:
Timestamp:
Sep 3, 2014, 3:18:54 PM (10 years ago)
Author:
Pavel Demin <demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
be2222c
Parents:
5b5a56b
Message:

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

File:
1 edited

Legend:

Unmodified
Added
Removed

  • external/fastjet/contribs/Nsubjettiness/Njettiness.cc

    r5b5a56b r35cdc46  
    55//  Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
    66//
     7//  $Id: Njettiness.cc 677 2014-06-12 18:56:46Z jthaler $
    78//----------------------------------------------------------------------
    89// This file is part of FastJet contrib.
    … …  
    3536///////
    3637
    37 // Helper function to correlate one pass minimization with appropriate measure
    38 void Njettiness::setOnePassAxesFinder(MeasureMode measure_mode, AxesFinder* startingFinder, double beta, double Rcutoff) {
    39    if (measure_mode == normalized_measure || measure_mode == unnormalized_measure || measure_mode == normalized_cutoff_measure || measure_mode == unnormalized_cutoff_measure) {
    40       _axesFinder = new AxesFinderFromOnePassMinimization(startingFinder, beta, Rcutoff);
    41    }
    42    else if (measure_mode == geometric_measure || measure_mode == geometric_cutoff_measure) {
    43       _axesFinder = new AxesFinderFromGeometricMinimization(startingFinder, beta, Rcutoff);
    44    }
    45    else {
    46       std::cerr << "Minimization only set up for normalized_measure, unnormalized_measure, normalized_cutoff_measure, unnormalized_cutoff_measure, geometric_measure, geometric_cutoff_measure" << std::endl;
    47       exit(1); }
    48 }
    49 
    50 // Parsing needed for constructor to set AxesFinder and MeasureFunction
    51 // All of the parameter handling is here, and checking that number of parameters is correct.
    52 void Njettiness::setMeasureFunctionandAxesFinder(AxesMode axes_mode, MeasureMode measure_mode, double para1, double para2, double para3, double para4) {
     38Njettiness::Njettiness(const AxesDefinition & axes_def, const MeasureDefinition & measure_def)
     39: _axes_def(axes_def.create()), _measure_def(measure_def.create()) {
     40   setMeasureFunctionAndAxesFinder();  // call helper function to do the hard work
     41}
     42
     43Njettiness::Njettiness(AxesMode axes_mode, const MeasureDefinition & measure_def)
     44: _axes_def(createAxesDef(axes_mode)), _measure_def(measure_def.create()) {
     45   setMeasureFunctionAndAxesFinder();  // call helper function to do the hard work
     46}
     47   
     48// Convert from MeasureMode enum to MeasureDefinition
     49// This returns a pointer that will be claimed by a SharedPtr
     50MeasureDefinition* Njettiness::createMeasureDef(MeasureMode measure_mode, int num_para, double para1, double para2, double para3) const {
    5351
    5452   // definition of maximum Rcutoff for non-cutoff measures, changed later by other measures
    5553   double Rcutoff = std::numeric_limits<double>::max();  //large number
    56    // Most (but all measures have some kind of beta value)
    57    double beta = NAN;
     54   // Most (but not all) measures have some kind of beta value
     55   double beta = std::numeric_limits<double>::quiet_NaN();
    5856   // The normalized measures have an R0 value.
    59    double R0 = NAN;
    60 
     57   double R0 = std::numeric_limits<double>::quiet_NaN();
     58   
    6159   // Find the MeasureFunction and set the parameters.
    6260   switch (measure_mode) {
    … …  
    6462         beta = para1;
    6563         R0 = para2;
    66          if(correctParameterCount(2, para1, para2, para3, para4))
    67             _measureFunction = new DefaultNormalizedMeasure(beta, R0, Rcutoff); //normalized_measure requires 2 parameters, beta and R0
    68          else {
    69             std::cerr << "normalized_measure needs 2 parameters (beta and R0)" << std::endl;
    70             exit(1); }
     64         if(num_para == 2) {
     65            return new NormalizedMeasure(beta,R0);
     66         } else {
     67            throw Error("normalized_measure needs 2 parameters (beta and R0)");
     68         }
    7169         break;
    7270      case unnormalized_measure:
    7371         beta = para1;
    74          if(correctParameterCount(1, para1, para2, para3, para4))
    75             _measureFunction = new DefaultUnnormalizedMeasure(beta, Rcutoff); //unnormalized_measure requires 1 parameter, beta
    76          else {
    77             std::cerr << "unnormalized_measure needs 1 parameter (beta)" << std::endl;
    78             exit(1); }
     72         if(num_para == 1) {
     73            return new UnnormalizedMeasure(beta);
     74         } else {
     75            throw Error("unnormalized_measure needs 1 parameter (beta)");
     76         }
    7977         break;
    8078      case geometric_measure:
    8179         beta = para1;
    82          if(correctParameterCount(1, para1, para2, para3, para4))
    83             _measureFunction = new GeometricMeasure(beta,Rcutoff); //geometric_measure requires 1 parameter, beta
    84          else {
    85             std::cerr << "geometric_measure needs 1 parameter (beta)" << std::endl;
    86             exit(1); }
     80         if (num_para == 1) {
     81            return new GeometricMeasure(beta);
     82         } else {
     83            throw Error("geometric_measure needs 1 parameter (beta)");
     84         }
    8785         break;
    8886      case normalized_cutoff_measure:
    … …  
    9088         R0 = para2;
    9189         Rcutoff = para3; //Rcutoff parameter is 3rd parameter in normalized_cutoff_measure
    92          if(correctParameterCount(3, para1, para2, para3, para4))
    93             _measureFunction = new DefaultNormalizedMeasure(beta, R0, Rcutoff); //normalized_cutoff_measure requires 3 parameters, beta, R0, and Rcutoff
    94          else {
    95             std::cerr << "normalized_cutoff_measure has 3 parameters (beta, R0, Rcutoff)" << std::endl;
    96             exit(1); }
     90         if (num_para == 3) {
     91            return new NormalizedCutoffMeasure(beta,R0,Rcutoff);
     92         } else {
     93            throw Error("normalized_cutoff_measure has 3 parameters (beta, R0, Rcutoff)");
     94         }
    9795         break;
    9896      case unnormalized_cutoff_measure:
    9997         beta = para1;
    10098         Rcutoff = para2; //Rcutoff parameter is 2nd parameter in normalized_cutoff_measure
    101          if (correctParameterCount(2, para1, para2, para3, para4))
    102             _measureFunction = new DefaultUnnormalizedMeasure(beta, Rcutoff); //unnormalized_cutoff_measure requires 2 parameters, beta and Rcutoff
    103          else {
    104             std::cerr << "unnormalized_cutoff_measure has 2 parameters (beta, Rcutoff)" << std::endl;
    105             exit(1); }
     99         if (num_para == 2) {
     100            return new UnnormalizedCutoffMeasure(beta,Rcutoff);
     101         } else {
     102            throw Error("unnormalized_cutoff_measure has 2 parameters (beta, Rcutoff)");
     103         }
    106104         break;
    107105      case geometric_cutoff_measure:
    108106         beta = para1;
    109107         Rcutoff = para2; //Rcutoff parameter is 2nd parameter in geometric_cutoff_measure
    110          if(correctParameterCount(2, para1, para2, para3, para4))
    111             _measureFunction = new GeometricMeasure(beta,Rcutoff); //geometric_cutoff_measure requires 2 parameters, beta and Rcutoff
    112          else {
    113             std::cerr << "geometric_cutoff_measure has 2 parameters (beta,Rcutoff)" << std::endl;
    114             exit(1); }
     108         if(num_para == 2) {
     109           return new GeometricCutoffMeasure(beta,Rcutoff);
     110         } else {
     111            throw Error("geometric_cutoff_measure has 2 parameters (beta, Rcutoff)");
     112         }
    115113         break;
    116114      default:
    117115         assert(false);
    118116         break;
    119    }   
    120 
    121    // Choose which AxesFinder from user input.
    122    // Uses setOnePassAxesFinder helpful function to use beta and Rcutoff values about (if needed)
     117   }
     118   return NULL;
     119}
     120
     121// Convert from AxesMode enum to AxesDefinition
     122// This returns a pointer that will be claimed by a SharedPtr
     123AxesDefinition* Njettiness::createAxesDef(Njettiness::AxesMode axes_mode) const {
     124   
    123125   switch (axes_mode) {
    124126      case wta_kt_axes:
    125          _axesFinder = new AxesFinderFromWTA_KT();
    126          break;
     127         return new WTA_KT_Axes();
    127128      case wta_ca_axes:
    128          _axesFinder = new AxesFinderFromWTA_CA();
    129          break;
     129         return new WTA_CA_Axes();
    130130      case kt_axes:
    131          _axesFinder = new AxesFinderFromKT();
    132          break;
     131         return new KT_Axes();
    133132      case ca_axes:
    134          _axesFinder = new AxesFinderFromCA();
    135          break;
     133         return new CA_Axes();
    136134      case antikt_0p2_axes:
    137          _axesFinder = new AxesFinderFromAntiKT(0.2);     
    138          break;
     135         return new AntiKT_Axes(0.2);
    139136      case onepass_wta_kt_axes:
    140          setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_KT(), beta, Rcutoff);
    141          break;
     137         return new OnePass_WTA_KT_Axes();
    142138      case onepass_wta_ca_axes:
    143          setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_CA(), beta, Rcutoff);
    144          break;
     139         return new OnePass_WTA_CA_Axes();
    145140      case onepass_kt_axes:
    146          setOnePassAxesFinder(measure_mode, new AxesFinderFromKT(), beta, Rcutoff);
    147          break;
     141         return new OnePass_KT_Axes();
    148142      case onepass_ca_axes:
    149          setOnePassAxesFinder(measure_mode, new AxesFinderFromCA(), beta, Rcutoff);
    150          break;
     143         return new OnePass_CA_Axes();
    151144      case onepass_antikt_0p2_axes:
    152          setOnePassAxesFinder(measure_mode, new AxesFinderFromAntiKT(0.2), beta, Rcutoff);
    153          break;
     145         return new OnePass_AntiKT_Axes(0.2);
    154146      case onepass_manual_axes:
    155          setOnePassAxesFinder(measure_mode, new AxesFinderFromUserInput(), beta, Rcutoff);
    156          break;
    157       case min_axes: //full minimization is not defined for geometric_measure.
    158          if (measure_mode == normalized_measure || measure_mode == unnormalized_measure || measure_mode == normalized_cutoff_measure || measure_mode == unnormalized_cutoff_measure)
    159             //Defaults to 100 iteration to find minimum
    160             _axesFinder = new AxesFinderFromKmeansMinimization(new AxesFinderFromKT(), beta, Rcutoff, 100);
    161          else {
    162             std::cerr << "Multi-pass minimization only set up for normalized_measure, unnormalized_measure, normalized_cutoff_measure, unnormalized_cutoff_measure." << std::endl;
    163             exit(1);
    164          }
    165          break;
     147         return new OnePass_Manual_Axes();
     148      case min_axes:
     149         return new MultiPass_Axes(100);
    166150      case manual_axes:
    167          _axesFinder = new AxesFinderFromUserInput();
    168          break;
    169 // These options have been commented out because they have not been fully tested
    170 //      case wta2_kt_axes: // option for alpha = 2 added
    171 //         _axesFinder = new AxesFinderFromWTA2_KT();
    172 //         break;
    173 //      case wta2_ca_axes: // option for alpha = 2 added
    174 //         _axesFinder = new AxesFinderFromWTA2_CA();
    175 //         break;
    176 //      case onepass_wta2_kt_axes: // option for alpha = 2 added
    177 //         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA2_KT(), beta, Rcutoff);
    178 //         break;
    179 //      case onepass_wta2_ca_axes: // option for alpha = 2 added
    180 //         setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA2_CA(), beta, Rcutoff);
    181 //         break;
     151         return new Manual_Axes();
    182152      default:
    183153         assert(false);
    184          break;
    185       }   
    186 
     154         return NULL;
     155   }
     156}
     157
     158   
     159// Parsing needed for constructor to set AxesFinder and MeasureFunction
     160// All of the parameter handling is here, and checking that number of parameters is correct.
     161void Njettiness::setMeasureFunctionAndAxesFinder() {
     162   // Get the correct MeasureFunction and AxesFinders
     163   _measureFunction.reset(_measure_def->createMeasureFunction());
     164   _startingAxesFinder.reset(_axes_def->createStartingAxesFinder(*_measure_def));
     165   _finishingAxesFinder.reset(_axes_def->createFinishingAxesFinder(*_measure_def));
    187166}
    188167
    189168// setAxes for Manual mode
    190 void Njettiness::setAxes(std::vector<fastjet::PseudoJet> myAxes) {
    191    if (_current_axes_mode == manual_axes || _current_axes_mode == onepass_manual_axes) {
     169void Njettiness::setAxes(const std::vector<fastjet::PseudoJet> & myAxes) {
     170   if (_axes_def->supportsManualAxes()) {
    192171      _currentAxes = myAxes;
    193    }
    194    else {
    195       std::cerr << "You can only use setAxes if using manual_axes or onepass_manual_axes measure mode" << std::endl;
    196       exit(1);
     172   } else {
     173      throw Error("You can only use setAxes for manual AxesDefinitions");
    197174   }
    198175}
    … …  
    200177// Calculates and returns all TauComponents that user would want.
    201178// This information is stored in _current_tau_components for later access as well.
    202 TauComponents Njettiness::getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) {
     179TauComponents Njettiness::getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) const {
    203180   if (inputJets.size() <= n_jets) {  //if not enough particles, return zero
    204181      _currentAxes = inputJets;
    … …  
    206183      _current_tau_components = TauComponents();
    207184      _seedAxes = _currentAxes;
     185      _currentJets = _currentAxes;
     186      _currentBeam = PseudoJet(0.0,0.0,0.0,0.0);
    208187   } else {
    209       _currentAxes = _axesFinder->getAxes(n_jets,inputJets,_currentAxes); // sets current Axes
    210       _seedAxes = _axesFinder->seedAxes(); // sets seed Axes (if one pass minimization was used)
    211       _current_tau_components = _measureFunction->result(inputJets, _currentAxes);  // sets current Tau Values
     188
     189      _seedAxes = _startingAxesFinder->getAxes(n_jets,inputJets,_currentAxes); //sets starting point for minimization
     190      if (_finishingAxesFinder) {
     191         _currentAxes = _finishingAxesFinder->getAxes(n_jets,inputJets,_seedAxes);
     192      } else {
     193         _currentAxes = _seedAxes;
     194      }
     195     
     196      // Find partition and store information
     197      // (jet information in _currentJets, beam in _currentBeam)
     198      _currentJets = _measureFunction->get_partition(inputJets,_currentAxes,&_currentBeam);
     199     
     200      // Find tau value and store information
     201      _current_tau_components = _measureFunction->result_from_partition(_currentJets, _currentAxes,&_currentBeam);  // sets current Tau Values
    212202   }
    213203   return _current_tau_components;
    … …  
    219209// Each vector element is a list of ints corresponding to the indices in
    220210// particles of the particles belonging to that jet.
    221 // TODO:  Consider moving to MeasureFunction
    222 std::vector<std::list<int> > Njettiness::getPartition(const std::vector<fastjet::PseudoJet> & particles) {
    223    std::vector<std::list<int> > partitions(_currentAxes.size());
    224 
    225    for (unsigned i = 0; i < particles.size(); i++) {
    226      
    227       int j_min = -1;
    228       // find minimum distance
    229       double minR = std::numeric_limits<double>::max();  //large number
    230       for (unsigned j = 0; j < _currentAxes.size(); j++) {
    231          double tempR = _measureFunction->jet_distance_squared(particles[i],_currentAxes[j]); // delta R distance
    232          if (tempR < minR) {
    233             minR = tempR;
    234             j_min = j;
    235          }
    236       }
    237       if (_measureFunction->do_cluster(particles[i],_currentAxes[j_min])) partitions[j_min].push_back(i);
    238    }
    239    return partitions;
    240 }
    241 
    242 // Having found axes, assign each particle in particles to an axis, and return a set of jets.
    243 // Each jet is the sum of particles closest to an axis (Njet = Naxes).
    244 // TODO:  Consider moving to MeasureFunction
    245 std::vector<fastjet::PseudoJet> Njettiness::getJets(const std::vector<fastjet::PseudoJet> & particles) {
    246    
    247    std::vector<fastjet::PseudoJet> jets(_currentAxes.size());
    248 
    249    std::vector<std::list<int> > partition = getPartition(particles);
    250    for (unsigned j = 0; j < partition.size(); ++j) {
    251       std::list<int>::const_iterator it, itE;
    252       for (it = partition[j].begin(), itE = partition[j].end(); it != itE; ++it) {
    253          jets[j] += particles[*it];
    254       }
    255    }
    256    return jets;
    257 }
    258 
     211std::vector<std::list<int> > Njettiness::getPartitionList(const std::vector<fastjet::PseudoJet> & particles) const {
     212   // core code is in MeasureFunction
     213   return _measureFunction->get_partition_list(particles,_currentAxes);
     214}
     215
     216   
    259217} // namespace contrib
    260218
Note: See TracChangeset for help on using the changeset viewer.