Fork me on GitHub

Changeset 9687203 in git for modules/FastJetFinder.cc


Ignore:
Timestamp:
Apr 16, 2014, 4:32:53 PM (11 years ago)
Author:
mselvaggi <mselvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
82575a3
Parents:
a0431dc
Message:

added nsubjettiness

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/FastJetFinder.cc

    ra0431dc r9687203  
    4747#include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
    4848
     49#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
     50#include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
     51#include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
     52#include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
     53
    4954using namespace std;
    5055using namespace fastjet;
     56using namespace fastjet::contrib;
     57
    5158
    5259//------------------------------------------------------------------------------
    5360
    5461FastJetFinder::FastJetFinder() :
    55   fPlugin(0), fDefinition(0), fAreaDefinition(0), fItInputArray(0)
     62  fPlugin(0), fDefinition(0), fAreaDefinition(0), fItInputArray(0), fRecomb(0), fNjettinessPlugin(0)
    5663{
    5764
     
    6976void FastJetFinder::Init()
    7077{
     78 
    7179  JetDefinition::Plugin *plugin = NULL;
    72 
     80  JetDefinition::Recombiner *recomb = NULL;
     81  NjettinessPlugin *njet_plugin = NULL;
     82 
    7383  // read eta ranges
    7484
     
    99109  fJetPTMin = GetDouble("JetPTMin", 10.0);
    100110
     111  //-- N(sub)jettiness parameters --
     112 
     113  fComputeNsubjettiness = GetBool("ComputeNsubjettiness", false);
     114  fBeta = GetDouble("Beta", 1.0);
     115  fAxisMode = GetInt("AxisMode", 1);
     116  fRcutOff = GetDouble("RcutOff", 0.8); //used only if Njettiness is used as jet clustering algo (case 8)
     117  fN = GetInt("N", 2);                  //used only if Njettiness is used as jet clustering algo (case 8)
     118 
    101119  // ---  Jet Area Parameters ---
    102120  fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
    103121  fComputeRho = GetBool("ComputeRho", false);
     122 
    104123  // - ghost based areas -
    105124  fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
     
    109128  fPtScatter = GetDouble("PtScatter", 0.1);
    110129  fMeanGhostPt = GetDouble("MeanGhostPt", 1.0E-100);
     130 
    111131  // - voronoi based areas -
    112132  fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
     
    159179      fDefinition = new fastjet::JetDefinition(fastjet::antikt_algorithm, fParameterR);
    160180      break;
    161   }
     181    case 7:
     182      recomb = new fastjet::contrib::WinnerTakeAllRecombiner();
     183      fDefinition = new fastjet::JetDefinition(fastjet::antikt_algorithm, fParameterR, recomb, Best);
     184      break;
     185    case 8:
     186      njet_plugin = new fastjet::contrib::NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
     187      fDefinition = new fastjet::JetDefinition(njet_plugin);
     188      break;
     189  }
     190 
    162191
    163192  fPlugin = plugin;
    164 
     193  fRecomb = recomb;
     194  fNjettinessPlugin = njet_plugin;
     195 
    165196  ClusterSequence::print_banner();
    166197
     
    184215  if(fAreaDefinition) delete fAreaDefinition;
    185216  if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
     217  if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
     218  if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
    186219}
    187220
     
    248281  outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
    249282
     283
    250284  // loop over all jets and export them
    251285  detaMax = 0.0;
     
    265299    inputList.clear();
    266300    inputList = sequence->constituents(*itOutputList);
     301     
    267302    for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
    268303    {
     
    289324    candidate->DeltaPhi = dphiMax;
    290325
     326    // --- compute N-subjettiness with N = 1,2,3,4,5 ----
     327   
     328    if(fComputeNsubjettiness)
     329    {
     330      Njettiness::AxesMode axisMode;
     331     
     332      if (fAxisMode == 1) axisMode = Njettiness::wta_kt_axes;
     333      if (fAxisMode == 2) axisMode = Njettiness::onepass_wta_kt_axes;
     334      if (fAxisMode == 3) axisMode = Njettiness::kt_axes;
     335      if (fAxisMode == 4) axisMode = Njettiness::onepass_kt_axes;
     336       
     337      Njettiness::MeasureMode measureMode = Njettiness::unnormalized_measure;
     338     
     339      Nsubjettiness nSub1(1, axisMode, measureMode, fBeta);
     340      Nsubjettiness nSub2(2, axisMode, measureMode, fBeta);
     341      Nsubjettiness nSub3(3, axisMode, measureMode, fBeta);
     342      Nsubjettiness nSub4(4, axisMode, measureMode, fBeta);
     343      Nsubjettiness nSub5(5, axisMode, measureMode, fBeta);
     344   
     345      candidate -> Tau1 = nSub1(*itOutputList);
     346      candidate -> Tau2 = nSub2(*itOutputList);
     347      candidate -> Tau3 = nSub3(*itOutputList);
     348      candidate -> Tau4 = nSub4(*itOutputList);
     349      candidate -> Tau5 = nSub5(*itOutputList);
     350    }
     351
     352
    291353    fOutputArray->Add(candidate);
    292354  }
Note: See TracChangeset for help on using the changeset viewer.