Index: examples/CMakeLists.txt
===================================================================
--- examples/CMakeLists.txt	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ examples/CMakeLists.txt	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -1,5 +1,5 @@
 include_directories(
   ${CMAKE_SOURCE_DIR}
-  ${DelphesExternals_INCLUDE_DIR} 
+  ${DelphesExternals_INCLUDE_DIR}
 )
 
@@ -17,2 +17,6 @@
 # take all other relevant files and put them into examples/
 install(FILES ${macros} DESTINATION examples)
+install(DIRECTORY ExternalFastJet DESTINATION examples)
+if(PYTHIA8_FOUND)
+  install(DIRECTORY Pythia8 DESTINATION examples)
+endif()
Index: modules/CMakeLists.txt
===================================================================
--- modules/CMakeLists.txt	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/CMakeLists.txt	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -28,6 +28,11 @@
 if (NOT ${ROOT_VERSION} VERSION_LESS "6.0.0")
   install(FILES
-      ${PROJECT_BINARY_DIR}/modules/libModulesDict_rdict.pcm
-      ${PROJECT_BINARY_DIR}/modules/libFastJetDict_rdict.pcm
+    ${PROJECT_BINARY_DIR}/modules/libModulesDict_rdict.pcm
+    ${PROJECT_BINARY_DIR}/modules/libFastJetDict_rdict.pcm
+    DESTINATION lib)
+  if(PYTHIA8_FOUND)
+    install(FILES
+      ${PROJECT_BINARY_DIR}/modules/libPythia8Dict_rdict.pcm
       DESTINATION lib)
+  endif()
 endif()
Index: modules/FastJetFinder.cc
===================================================================
--- modules/FastJetFinder.cc	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/FastJetFinder.cc	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -120,5 +120,4 @@
   fJetPTMin = GetDouble("JetPTMin", 10.0);
 
- 
   //-- N(sub)jettiness parameters --
 
@@ -130,53 +129,54 @@
 
   //-- Exclusive clustering for e+e- collisions --
-  
+
   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();
-   }
+    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);
   fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
   fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
-  
 
   // ---  Jet Area Parameters ---
+
   fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
   fComputeRho = GetBool("ComputeRho", false);
@@ -195,23 +195,23 @@
   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;
+    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;
   }
 
@@ -248,5 +248,5 @@
       fDefinition = new JetDefinition(fNjettinessPlugin);
       break;
-  case 9:
+    case 9:
       fValenciaPlugin = new ValenciaPlugin(fParameterR, fBeta, fGamma);
       fDefinition = new JetDefinition(fValenciaPlugin);
@@ -255,5 +255,5 @@
   }
 
-   
+
 
   fPlugin = plugin;
@@ -290,4 +290,5 @@
   fOutputArray = ExportArray(GetString("OutputArray", "jets"));
   fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
+  fConstituentsOutputArray = ExportArray(GetString("ConstituentsOutputArray", "constituents"));
 }
 
@@ -325,5 +326,5 @@
   Double_t time, timeWeight;
   Int_t number, ncharged, nneutrals;
-  Int_t charge; 
+  Int_t charge;
   Double_t rho = 0.0;
   PseudoJet jet, area;
@@ -336,5 +337,5 @@
   Double_t excl_ymerge45 = 0.0;
   Double_t excl_ymerge56 = 0.0;
-  
+
   DelphesFactory *factory = GetFactory();
 
@@ -381,5 +382,4 @@
   outputList.clear();
 
-  
   if(fExclusiveClustering)
     {
@@ -405,5 +405,5 @@
   detaMax = 0.0;
   dphiMax = 0.0;
-  
+
   for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
   {
@@ -416,6 +416,4 @@
     if(fAreaDefinition) area = itOutputList->area_4vector();
 
-
-    
     candidate = factory->NewCandidate();
 
@@ -449,4 +447,5 @@
       charge += constituent->Charge;
 
+      fConstituentsOutputArray->Add(constituent);
       candidate->AddCandidate(constituent);
     }
@@ -458,5 +457,5 @@
     candidate->DeltaEta = detaMax;
     candidate->DeltaPhi = dphiMax;
-    candidate->Charge = charge; 
+    candidate->Charge = charge;
     candidate->NNeutrals = nneutrals;
     candidate->NCharged = ncharged;
@@ -468,5 +467,5 @@
     candidate->ExclYmerge45 = excl_ymerge45;
     candidate->ExclYmerge56 = excl_ymerge56;
-    
+
     //------------------------------------
     // Trimming
@@ -478,29 +477,29 @@
       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++)
       {
-	    if(subjets.at(i).pt() < 0) continue ; 
+	    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)
     {
@@ -510,34 +509,34 @@
 
       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++)
       {
-	    if(subjets.at(i).pt() < 0) continue ; 
+	    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);
       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();
@@ -549,5 +548,5 @@
       for (size_t i = 0; i < subjets.size()  and i < 4; i++)
       {
-	    if(subjets.at(i).pt() < 0) continue ; 
+	    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];
@@ -555,10 +554,10 @@
       }
     }
-  
+
     // --- compute N-subjettiness with N = 1,2,3,4,5 ----
 
     if(fComputeNsubjettiness)
     {
-     
+
       Nsubjettiness nSub1(1, *fAxesDef, *fMeasureDef);
       Nsubjettiness nSub2(2, *fAxesDef, *fMeasureDef);
@@ -566,5 +565,5 @@
       Nsubjettiness nSub4(4, *fAxesDef, *fMeasureDef);
       Nsubjettiness nSub5(5, *fAxesDef, *fMeasureDef);
-     
+
       candidate->Tau[0] = nSub1(*itOutputList);
       candidate->Tau[1] = nSub2(*itOutputList);
@@ -572,5 +571,5 @@
       candidate->Tau[3] = nSub4(*itOutputList);
       candidate->Tau[4] = nSub5(*itOutputList);
-         
+
     }
 
Index: modules/FastJetFinder.h
===================================================================
--- modules/FastJetFinder.h	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/FastJetFinder.h	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -42,6 +42,6 @@
     class NjettinessPlugin;
     class ValenciaPlugin;
-    class AxesDefinition; 
-    class MeasureDefinition;    
+    class AxesDefinition;
+    class MeasureDefinition;
   }
 }
@@ -63,6 +63,6 @@
   void *fRecomb; //!
 
-  fastjet::contrib::AxesDefinition *fAxesDef; 
-  fastjet::contrib::MeasureDefinition *fMeasureDef;    
+  fastjet::contrib::AxesDefinition *fAxesDef;
+  fastjet::contrib::MeasureDefinition *fMeasureDef;
 
   fastjet::contrib::NjettinessPlugin *fNjettinessPlugin; //!
@@ -90,5 +90,5 @@
   //-- Valencia Linear Collider algorithm
   Double_t fGamma;
-  
+
   //-- N (sub)jettiness parameters --
 
@@ -100,9 +100,9 @@
 
   //-- Trimming parameters --
-  
+
   Bool_t fComputeTrimming;
   Double_t fRTrim;
   Double_t fPtFracTrim;
-  
+
   //-- Pruning parameters --
 
@@ -152,4 +152,5 @@
   TObjArray *fOutputArray; //!
   TObjArray *fRhoOutputArray; //!
+  TObjArray *fConstituentsOutputArray; //!
 
   ClassDef(FastJetFinder, 1)
Index: modules/PdgCodeFilter.cc
===================================================================
--- modules/PdgCodeFilter.cc	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/PdgCodeFilter.cc	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -76,4 +76,7 @@
   fInvert = GetBool("Invert", false);
 
+  // no pileup 
+  fRequireNotPileup = GetBool("RequireNotPileup", false);
+
   fRequireStatus = GetBool("RequireStatus", false);
   fStatus = GetInt("Status", 1);
@@ -127,4 +130,5 @@
     if(fRequireStatus && (candidate->Status != fStatus)) continue;
     if(fRequireCharge && (candidate->Charge != fCharge)) continue;
+    if(fRequireNotPileup && (candidate->IsPU   >0     )) continue;
 
     pass = kTRUE;
Index: modules/PdgCodeFilter.h
===================================================================
--- modules/PdgCodeFilter.h	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/PdgCodeFilter.h	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -55,5 +55,5 @@
   Bool_t fRequireCharge; //!
   Int_t fCharge; //!
-
+  Bool_t fRequireNotPileup; //!
 
   std::vector<Int_t> fPdgCodes;
Index: modules/StatusPidFilter.cc
===================================================================
--- modules/StatusPidFilter.cc	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/StatusPidFilter.cc	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -153,6 +153,8 @@
 {
   // PT threshold
-
   fPTMin = GetDouble("PTMin", 0.5);
+
+  // keep or remove pileup particles
+  fRequireNotPileup = GetBool("RequireNotPileup", false);
 
   // import input array
@@ -227,4 +229,7 @@
     if(!pass || (candidate->Momentum.Pt() < fPTMin && !(is_b_hadron || is_b_quark || is_tau_daughter || is_W_daughter)) ) continue;
 
+    // not pileup particles
+    if(fRequireNotPileup && (candidate->IsPU >0)) continue;
+
     fOutputArray->Add(candidate);
   }
Index: modules/StatusPidFilter.h
===================================================================
--- modules/StatusPidFilter.h	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/StatusPidFilter.h	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -51,4 +51,6 @@
   Double_t fPTMin; //!
 
+  Bool_t fRequireNotPileup; //!
+
   TIterator *fItInputArray; //!
 
Index: modules/UniqueObjectFinder.cc
===================================================================
--- modules/UniqueObjectFinder.cc	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/UniqueObjectFinder.cc	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -67,4 +67,7 @@
 void UniqueObjectFinder::Init()
 {
+  // use GetUniqueID algorithm to find unique objects (faster than the default Overlaps method)
+  fUseUniqueID = GetBool("UseUniqueID", false);
+
   // import arrays with output from other modules
 
@@ -146,7 +149,17 @@
     while((previousCandidate = static_cast<Candidate*>(iterator.Next())))
     {
-      if(candidate->Overlaps(previousCandidate))
+      if(fUseUniqueID)
       {
-        return kFALSE;
+        if(candidate->GetUniqueID() == previousCandidate->GetUniqueID())
+        {
+          return kFALSE;
+        }
+      }
+      else
+      {
+        if(candidate->Overlaps(previousCandidate))
+        {
+          return kFALSE;
+        }
       }
     }
Index: modules/UniqueObjectFinder.h
===================================================================
--- modules/UniqueObjectFinder.h	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ modules/UniqueObjectFinder.h	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -50,4 +50,6 @@
 private:
 
+  Bool_t fUseUniqueID;
+
   Bool_t Unique(Candidate *candidate, std::vector< std::pair< TIterator *, TObjArray * > >::iterator itInputMap);
 
Index: readers/DelphesProMC.cpp
===================================================================
--- readers/DelphesProMC.cpp	(revision dd64cff451c2f6c84f259f0c660b8cd87e641503)
+++ readers/DelphesProMC.cpp	(revision 8a11cdcd860e1b3a14e0892c626695b66aa6d896)
@@ -125,4 +125,7 @@
     candidate->Status = status;
 
+    candidate->IsPU=0;
+    if (mutableParticles->barcode(i)>0) candidate->IsPU=1; // pileup particle
+
     candidate->M1 = mutableParticles->mother1(i);
     candidate->M2 = mutableParticles->mother2(i);
@@ -231,4 +234,8 @@
       cout << "** Reading " << argv[i] << endl;
 
+      // use 64 bit 
+      //inputFile = new ProMCBook(argv[i], "r", true);
+
+      //use 32 bit zip (faster but limitted to 64k events)
       inputFile = new ProMCBook(argv[i], "r");
 
