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

File:

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

modules/FastJetFinder.cc

@@ -17 +17 @@
  */
 
-
 /** \class FastJetFinder
  *
@@ -32 +31 @@
 #include "classes/DelphesFormula.h"
 
+#include "ExRootAnalysis/ExRootClassifier.h"
+#include "ExRootAnalysis/ExRootFilter.h"
 #include "ExRootAnalysis/ExRootResult.h"
-#include "ExRootAnalysis/ExRootFilter.h"
-#include "ExRootAnalysis/ExRootClassifier.h"
-
+
+#include "TDatabasePDG.h"
+#include "TFormula.h"
+#include "TLorentzVector.h"
 #include "TMath.h"
+#include "TObjArray.h"
+#include "TRandom3.h"
 #include "TString.h"
-#include "TFormula.h"
-#include "TRandom3.h"
-#include "TObjArray.h"
-#include "TDatabasePDG.h"
-#include "TLorentzVector.h"
 
 #include <algorithm>
-#include <stdexcept>
 #include <iostream>
 #include <sstream>
+#include <stdexcept>
 #include <vector>
 
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/JetDefinition.hh"
 #include "fastjet/PseudoJet.hh"
-#include "fastjet/JetDefinition.hh"
-#include "fastjet/ClusterSequence.hh"
 #include "fastjet/Selector.hh"
-#include "fastjet/ClusterSequenceArea.hh"
 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
 
+#include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
+#include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
 #include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
-#include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
-#include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
-
-#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
+
+#include "fastjet/contribs/Nsubjettiness/ExtraRecombiners.hh"
 #include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
 #include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
-#include "fastjet/contribs/Nsubjettiness/ExtraRecombiners.hh"
+#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.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;
-
 
 //------------------------------------------------------------------------------
@@ -83 +81 @@
   fDefinition(0), fAreaDefinition(0), fItInputArray(0)
 {
-
 }
 
@@ -90 +87 @@
 FastJetFinder::~FastJetFinder()
 {
-
 }
 
@@ -120 +116 @@
   fJetPTMin = GetDouble("JetPTMin", 10.0);
 
- 
   //-- N(sub)jettiness parameters --
 
@@ -127 +122 @@
   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);
@@ -195 +190 @@
   switch(fAreaAlgorithm)
   {
-    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;
+  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;
   }
 
   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;
@@ -270 +262 @@
 
     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;
@@ -290 +282 @@
   fOutputArray = ExportArray(GetString("OutputArray", "jets"));
   fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
+  fConstituentsOutputArray = ExportArray(GetString("ConstituentsOutputArray", "constituents"));
 }
 
@@ -296 +289 @@
 void FastJetFinder::Finish()
 {
-  vector< TEstimatorStruct >::iterator itEstimators;
+  vector<TEstimatorStruct>::iterator itEstimators;
 
   for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
@@ -306 +299 @@
   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);
 }
 
@@ -325 +317 @@
   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();
 
@@ -344 +336 @@
   fItInputArray->Reset();
   number = 0;
-  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
+  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
   {
     momentum = candidate->Momentum;
@@ -381 +373 @@
   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)
   {
@@ -416 +408 @@
     if(fAreaDefinition) area = itOutputList->area_4vector();
 
-
-    
     candidate = factory->NewCandidate();
 
@@ -434 +424 @@
     {
       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());
@@ -441 +431 @@
       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);
@@ -468 +460 @@
     candidate->ExclYmerge45 = excl_ymerge45;
     candidate->ExclYmerge56 = excl_ymerge56;
-    
+
     //------------------------------------
     // Trimming
@@ -476 +468 @@
     {
 
-      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);
@@ -566 +555 @@
       Nsubjettiness nSub4(4, *fAxesDef, *fMeasureDef);
       Nsubjettiness nSub5(5, *fAxesDef, *fMeasureDef);
-     
+
       candidate->Tau[0] = nSub1(*itOutputList);
       candidate->Tau[1] = nSub2(*itOutputList);
@@ -572 +561 @@
       candidate->Tau[3] = nSub4(*itOutputList);
       candidate->Tau[4] = nSub5(*itOutputList);
-         
     }