Changes in modules/FastJetFinder.cc [341014c:35807af] in git

File:

• modules/FastJetFinder.cc (modified) (21 diffs)

modules/FastJetFinder.cc

@@ -17 +17 @@
  */
 
+
 /** \class FastJetFinder
  *
@@ -31 +32 @@
 #include "classes/DelphesFormula.h"
 
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
 #include "ExRootAnalysis/ExRootClassifier.h"
-#include "ExRootAnalysis/ExRootFilter.h"
-#include "ExRootAnalysis/ExRootResult.h"
-
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
 #include "TDatabasePDG.h"
-#include "TFormula.h"
 #include "TLorentzVector.h"
-#include "TMath.h"
-#include "TObjArray.h"
-#include "TRandom3.h"
-#include "TString.h"
 
 #include <algorithm>
+#include <stdexcept>
 #include <iostream>
 #include <sstream>
-#include <stdexcept>
 #include <vector>
 
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/JetDefinition.hh"
 #include "fastjet/ClusterSequence.hh"
+#include "fastjet/Selector.hh"
 #include "fastjet/ClusterSequenceArea.hh"
-#include "fastjet/JetDefinition.hh"
-#include "fastjet/PseudoJet.hh"
-#include "fastjet/Selector.hh"
 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
 
+#include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
+#include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
 #include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
-#include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
-#include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
-
-#include "fastjet/contribs/Nsubjettiness/ExtraRecombiners.hh"
+
+#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
 #include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
 #include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
-#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
+#include "fastjet/contribs/Nsubjettiness/ExtraRecombiners.hh"
 
 #include "fastjet/contribs/ValenciaPlugin/ValenciaPlugin.hh"
 
-#include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
 #include "fastjet/tools/Filter.hh"
 #include "fastjet/tools/Pruner.hh"
+#include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
 
 using namespace std;
 using namespace fastjet;
 using namespace fastjet::contrib;
+
 
 //------------------------------------------------------------------------------
@@ -81 +83 @@
   fDefinition(0), fAreaDefinition(0), fItInputArray(0)
 {
+
 }
 
@@ -87 +90 @@
 FastJetFinder::~FastJetFinder()
 {
+
 }
 
@@ -116 +120 @@
   fJetPTMin = GetDouble("JetPTMin", 10.0);
 
+ 
   //-- N(sub)jettiness parameters --
 
@@ -122 +127 @@
   fAxisMode = GetInt("AxisMode", 1);
   fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
-  fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8)
+  fN = GetInt("N", 2);                  // used only if Njettiness is used as jet clustering algo (case 8)
 
   //-- Exclusive clustering for e+e- collisions --
-
-  fNJets = GetInt("NJets", 2);
+  
+  fNJets = GetInt("NJets",2);
   fExclusiveClustering = GetBool("ExclusiveClustering", false);
 
   //-- Valencia Linear Collider algorithm
-
   fGamma = GetDouble("Gamma", 1.0);
   //fBeta parameter see above
-
+  
   fMeasureDef = new NormalizedMeasure(fBeta, fParameterR);
-
+   
   switch(fAxisMode)
   {
-  default:
-  case 1:
-    fAxesDef = new WTA_KT_Axes();
-    break;
-  case 2:
-    fAxesDef = new OnePass_WTA_KT_Axes();
-    break;
-  case 3:
-    fAxesDef = new KT_Axes();
-    break;
-  case 4:
-    fAxesDef = new OnePass_KT_Axes();
-  }
+    default:
+      case 1:
+        fAxesDef = new WTA_KT_Axes();
+        break;
+      case 2:
+        fAxesDef = new OnePass_WTA_KT_Axes();
+        break;
+      case 3:
+        fAxesDef = new KT_Axes();
+        break;
+      case 4:
+        fAxesDef = new OnePass_KT_Axes();
+   }
 
   //-- Trimming parameters --
-
+  
   fComputeTrimming = GetBool("ComputeTrimming", false);
   fRTrim = GetDouble("RTrim", 0.2);
   fPtFracTrim = GetDouble("PtFracTrim", 0.05);
+  
 
   //-- Pruning parameters --
-
+  
   fComputePruning = GetBool("ComputePruning", false);
   fZcutPrun = GetDouble("ZcutPrun", 0.1);
   fRcutPrun = GetDouble("RcutPrun", 0.5);
   fRPrun = GetDouble("RPrun", 0.8);
-
+ 
   //-- SoftDrop parameters --
-
-  fComputeSoftDrop = GetBool("ComputeSoftDrop", false);
-  fBetaSoftDrop = GetDouble("BetaSoftDrop", 0.0);
+  
+  fComputeSoftDrop     = GetBool("ComputeSoftDrop", false);
+  fBetaSoftDrop        = GetDouble("BetaSoftDrop", 0.0);
   fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
-  fR0SoftDrop = GetDouble("R0SoftDrop=", 0.8);
+  fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
+  
 
   // ---  Jet Area Parameters ---
-
   fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
   fComputeRho = GetBool("ComputeRho", false);
@@ -190 +195 @@
   switch(fAreaAlgorithm)
   {
-  default:
-  case 0:
-    fAreaDefinition = 0;
-    break;
-  case 1:
-    fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
-    break;
-  case 2:
-    fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
-    break;
-  case 3:
-    fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
-    break;
-  case 4:
-    fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
-    break;
-  case 5:
-    fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
-    break;
+    case 1:
+      fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
+      break;
+    case 2:
+      fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
+      break;
+    case 3:
+      fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
+      break;
+    case 4:
+      fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
+      break;
+    case 5:
+      fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
+      break;
+    default:
+    case 0:
+      fAreaDefinition = 0;
+      break;
   }
 
   switch(fJetAlgorithm)
   {
-  case 1:
-    plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
-    fDefinition = new JetDefinition(plugin);
-    break;
-  case 2:
-    plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
-    fDefinition = new JetDefinition(plugin);
-    break;
-  case 3:
-    plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
-    fDefinition = new JetDefinition(plugin);
-    break;
-  case 4:
-    fDefinition = new JetDefinition(kt_algorithm, fParameterR);
-    break;
-  case 5:
-    fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
-    break;
-  default:
-  case 6:
-    fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
-    break;
-  case 7:
-    recomb = new WinnerTakeAllRecombiner();
-    fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
-    break;
-  case 8:
-    fNjettinessPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
-    fDefinition = new JetDefinition(fNjettinessPlugin);
-    break;
+    case 1:
+      plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
+      fDefinition = new JetDefinition(plugin);
+      break;
+    case 2:
+      plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
+      fDefinition = new JetDefinition(plugin);
+      break;
+    case 3:
+      plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
+      fDefinition = new JetDefinition(plugin);
+      break;
+    case 4:
+      fDefinition = new JetDefinition(kt_algorithm, fParameterR);
+      break;
+    case 5:
+      fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
+      break;
+    default:
+    case 6:
+      fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
+      break;
+    case 7:
+      recomb = new WinnerTakeAllRecombiner();
+      fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
+      break;
+    case 8:
+      fNjettinessPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
+      fDefinition = new JetDefinition(fNjettinessPlugin);
+      break;
   case 9:
-    fValenciaPlugin = new ValenciaPlugin(fParameterR, fBeta, fGamma);
-    fDefinition = new JetDefinition(fValenciaPlugin);
-    break;
-  }
+      fValenciaPlugin = new ValenciaPlugin(fParameterR, fBeta, fGamma);
+      fDefinition = new JetDefinition(fValenciaPlugin);
+      break;
+
+  }
+
+   
 
   fPlugin = plugin;
@@ -262 +270 @@
 
     fEstimators.clear();
-    for(i = 0; i < size / 2; ++i)
-    {
-      etaMin = param[i * 2].GetDouble();
-      etaMax = param[i * 2 + 1].GetDouble();
+    for(i = 0; i < size/2; ++i)
+    {
+      etaMin = param[i*2].GetDouble();
+      etaMax = param[i*2 + 1].GetDouble();
       estimatorStruct.estimator = new JetMedianBackgroundEstimator(SelectorRapRange(etaMin, etaMax), *fDefinition, *fAreaDefinition);
       estimatorStruct.etaMin = etaMin;
@@ -282 +290 @@
   fOutputArray = ExportArray(GetString("OutputArray", "jets"));
   fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
-  fConstituentsOutputArray = ExportArray(GetString("ConstituentsOutputArray", "constituents"));
 }
 
@@ -289 +296 @@
 void FastJetFinder::Finish()
 {
-  vector<TEstimatorStruct>::iterator itEstimators;
+  vector< TEstimatorStruct >::iterator itEstimators;
 
   for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
@@ -299 +306 @@
   if(fDefinition) delete fDefinition;
   if(fAreaDefinition) delete fAreaDefinition;
-  if(fPlugin) delete static_cast<JetDefinition::Plugin *>(fPlugin);
-  if(fRecomb) delete static_cast<JetDefinition::Recombiner *>(fRecomb);
-  if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin *>(fNjettinessPlugin);
+  if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
+  if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
+  if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
   if(fAxesDef) delete fAxesDef;
   if(fMeasureDef) delete fMeasureDef;
-  if(fValenciaPlugin) delete static_cast<JetDefinition::Plugin *>(fValenciaPlugin);
+  if(fValenciaPlugin) delete static_cast<JetDefinition::Plugin*>(fValenciaPlugin);
+
 }
 
@@ -317 +325 @@
   Double_t time, timeWeight;
   Int_t number, ncharged, nneutrals;
-  Int_t charge;
+  Int_t charge; 
   Double_t rho = 0.0;
   PseudoJet jet, area;
   ClusterSequence *sequence;
-  vector<PseudoJet> inputList, outputList, subjets;
-  vector<PseudoJet>::iterator itInputList, itOutputList;
-  vector<TEstimatorStruct>::iterator itEstimators;
+  vector< PseudoJet > inputList, outputList, subjets;
+  vector< PseudoJet >::iterator itInputList, itOutputList;
+  vector< TEstimatorStruct >::iterator itEstimators;
   Double_t excl_ymerge23 = 0.0;
   Double_t excl_ymerge34 = 0.0;
   Double_t excl_ymerge45 = 0.0;
   Double_t excl_ymerge56 = 0.0;
-
+  
   DelphesFactory *factory = GetFactory();
 
@@ -336 +344 @@
   fItInputArray->Reset();
   number = 0;
-  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
+  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
     momentum = candidate->Momentum;
@@ -373 +381 @@
   outputList.clear();
 
+  
   if(fExclusiveClustering)
-  {
-    try
-    {
-      outputList = sorted_by_pt(sequence->exclusive_jets(fNJets));
-    }
-    catch(fastjet::Error)
-    {
-      outputList.clear();
-    }
-
-    excl_ymerge23 = sequence->exclusive_ymerge(2);
-    excl_ymerge34 = sequence->exclusive_ymerge(3);
-    excl_ymerge45 = sequence->exclusive_ymerge(4);
-    excl_ymerge56 = sequence->exclusive_ymerge(5);
-  }
+    {
+   try{
+     outputList = sorted_by_pt(sequence->exclusive_jets( fNJets ));
+   }
+   catch(fastjet::Error)
+     {
+       outputList.clear();
+     }
+
+      excl_ymerge23 = sequence->exclusive_ymerge( 2 );
+      excl_ymerge34 = sequence->exclusive_ymerge( 3 );
+      excl_ymerge45 = sequence->exclusive_ymerge( 4 );
+      excl_ymerge56 = sequence->exclusive_ymerge( 5 );
+    }
   else
-  {
-    outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
-  }
+    {
+      outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
+    }
 
   // loop over all jets and export them
   detaMax = 0.0;
   dphiMax = 0.0;
-
+  
   for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
   {
@@ -408 +416 @@
     if(fAreaDefinition) area = itOutputList->area_4vector();
 
+
+    
     candidate = factory->NewCandidate();
 
@@ -424 +434 @@
     {
       if(itInputList->user_index() < 0) continue;
-      constituent = static_cast<Candidate *>(fInputArray->At(itInputList->user_index()));
+      constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
 
       deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
@@ -431 +441 @@
       if(dphi > dphiMax) dphiMax = dphi;
 
-      if(constituent->Charge == 0)
-        nneutrals++;
-      else
-        ncharged++;
-
-      time += TMath::Sqrt(constituent->Momentum.E()) * (constituent->Position.T());
+      if(constituent->Charge == 0) nneutrals++;
+      else ncharged++;
+
+      time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
       timeWeight += TMath::Sqrt(constituent->Momentum.E());
 
       charge += constituent->Charge;
 
-      fConstituentsOutputArray->Add(constituent);
       candidate->AddCandidate(constituent);
     }
 
     candidate->Momentum = momentum;
-    candidate->Position.SetT(time / timeWeight);
+    candidate->Position.SetT(time/timeWeight);
     candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
 
     candidate->DeltaEta = detaMax;
     candidate->DeltaPhi = dphiMax;
-    candidate->Charge = charge;
+    candidate->Charge = charge; 
     candidate->NNeutrals = nneutrals;
     candidate->NCharged = ncharged;
+
 
     //for exclusive clustering, access y_n,n+1 as exclusive_ymerge (fNJets);
@@ -460 +468 @@
     candidate->ExclYmerge45 = excl_ymerge45;
     candidate->ExclYmerge56 = excl_ymerge56;
-
+    
     //------------------------------------
     // Trimming
@@ -468 +476 @@
     {
 
-      fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, fRTrim), fastjet::SelectorPtFractionMin(fPtFracTrim));
+      fastjet::Filter    trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim));
       fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
-
+      
       trimmed_jet = join(trimmed_jet.constituents());
-
+     
       candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
-
-      // four hardest subjets
+        
+      // four hardest subjets 
       subjets.clear();
       subjets = trimmed_jet.pieces();
       subjets = sorted_by_pt(subjets);
-
+      
       candidate->NSubJetsTrimmed = subjets.size();
 
-      for(size_t i = 0; i < subjets.size() and i < 4; i++)
+      for (size_t i = 0; i < subjets.size() and i < 4; i++)
       {
-        if(subjets.at(i).pt() < 0) continue;
-        candidate->TrimmedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
+            if(subjets.at(i).pt() < 0) continue ; 
+            candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
       }
     }
-
+    
+    
     //------------------------------------
     // Pruning
     //------------------------------------
-
+    
+    
     if(fComputePruning)
     {
 
-      fastjet::Pruner pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm, fRPrun), fZcutPrun, fRcutPrun);
+      fastjet::Pruner    pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm,fRPrun),fZcutPrun,fRcutPrun);
       fastjet::PseudoJet pruned_jet = pruner(*itOutputList);
 
       candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
-
-      // four hardest subjet
+         
+      // four hardest subjet 
       subjets.clear();
       subjets = pruned_jet.pieces();
       subjets = sorted_by_pt(subjets);
-
+      
       candidate->NSubJetsPruned = subjets.size();
 
-      for(size_t i = 0; i < subjets.size() and i < 4; i++)
+      for (size_t i = 0; i < subjets.size() and i < 4; i++)
       {
-        if(subjets.at(i).pt() < 0) continue;
-        candidate->PrunedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
+            if(subjets.at(i).pt() < 0) continue ; 
+            candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
       }
-    }
-
+
+    } 
+     
     //------------------------------------
     // SoftDrop
     //------------------------------------
-
+   
     if(fComputeSoftDrop)
     {
-
-      contrib::SoftDrop softDrop(fBetaSoftDrop, fSymmetryCutSoftDrop, fR0SoftDrop);
+    
+      contrib::SoftDrop softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop);
       fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
-
+      
       candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
-
-      // four hardest subjet
-
+        
+      // four hardest subjet 
+      
       subjets.clear();
-      subjets = softdrop_jet.pieces();
-      subjets = sorted_by_pt(subjets);
+      subjets    = softdrop_jet.pieces();
+      subjets    = sorted_by_pt(subjets);
       candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
 
       candidate->SoftDroppedJet = candidate->SoftDroppedP4[0];
 
-      for(size_t i = 0; i < subjets.size() and i < 4; i++)
+      for (size_t i = 0; i < subjets.size()  and i < 4; i++)
       {
-        if(subjets.at(i).pt() < 0) continue;
-        candidate->SoftDroppedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
-        if(i == 0) candidate->SoftDroppedSubJet1 = candidate->SoftDroppedP4[i + 1];
-        if(i == 1) candidate->SoftDroppedSubJet2 = candidate->SoftDroppedP4[i + 1];
+            if(subjets.at(i).pt() < 0) continue ; 
+            candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
+            if(i==0) candidate->SoftDroppedSubJet1 = candidate->SoftDroppedP4[i+1];
+            if(i==1) candidate->SoftDroppedSubJet2 = candidate->SoftDroppedP4[i+1];
       }
     }
-
+  
     // --- compute N-subjettiness with N = 1,2,3,4,5 ----
 
     if(fComputeNsubjettiness)
     {
-
+     
       Nsubjettiness nSub1(1, *fAxesDef, *fMeasureDef);
       Nsubjettiness nSub2(2, *fAxesDef, *fMeasureDef);
@@ -555 +566 @@
       Nsubjettiness nSub4(4, *fAxesDef, *fMeasureDef);
       Nsubjettiness nSub5(5, *fAxesDef, *fMeasureDef);
-
+     
       candidate->Tau[0] = nSub1(*itOutputList);
       candidate->Tau[1] = nSub2(*itOutputList);
@@ -561 +572 @@
       candidate->Tau[3] = nSub4(*itOutputList);
       candidate->Tau[4] = nSub5(*itOutputList);
+         
     }