Changeset d77b51d in git for modules

Timestamp:

Sep 29, 2015, 2:08:10 PM (9 years ago)

Author:

Michele Selvaggi <michele.selvaggi@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

a98c7ef

Parents:

d870fc5 (diff), 06ec139 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge remote-tracking branch 'upstream/master'

Location:

modules

Files:

• AngularSmearing.cc (modified) (2 diffs)
• BTagging.cc (modified) (11 diffs)
• BTagging.h (modified) (2 diffs)
• CMakeLists.txt (modified) (2 diffs)
• Calorimeter.cc (modified) (10 diffs)
• Calorimeter.h (modified) (2 diffs)
• Efficiency.cc (modified) (2 diffs)
• EnergyScale.cc (modified) (1 diff)
• EnergySmearing.cc (modified) (3 diffs)
• FastJetFinder.cc (modified) (14 diffs)
• FastJetFinder.h (modified) (4 diffs)
• FastJetLinkDef.h (modified) (2 diffs)
• IdentificationMap.cc (modified) (6 diffs)
• IdentificationMap.h (modified) (1 diff)
• ImpactParameterSmearing.cc (modified) (4 diffs)
• Isolation.cc (modified) (4 diffs)
• JetFakeParticle.cc (added)
• JetFakeParticle.h (added)
• JetFlavorAssociation.cc (added)
• JetFlavorAssociation.h (added)
• JetPileUpSubtractor.cc (modified) (1 diff)
• ModulesLinkDef.h (modified) (6 diffs)
• MomentumSmearing.cc (modified) (2 diffs)
• PdgCodeFilter.cc (modified) (3 diffs)
• PdgCodeFilter.h (modified) (1 diff)
• PhotonConversions.cc (added)
• PhotonConversions.h (added)
• PileUpJetID.cc (modified) (8 diffs)
• PileUpJetID.h (modified) (5 diffs)
• PileUpMerger.cc (modified) (7 diffs)
• PileUpMerger.h (modified) (1 diff)
• PileUpMergerPythia8.cc (modified) (8 diffs)
• PileUpMergerPythia8.h (modified) (3 diffs)
• RunPUPPI.cc (added)
• RunPUPPI.h (added)
• SimpleCalorimeter.cc (modified) (4 diffs)
• SimpleCalorimeter.h (modified) (2 diffs)
• TaggingParticlesSkimmer.cc (added)
• TaggingParticlesSkimmer.h (added)
• TauTagging.cc (modified) (7 diffs)
• TauTagging.h (modified) (2 diffs)
• TrackPileUpSubtractor.cc (modified) (2 diffs)
• TrackPileUpSubtractor.h (modified) (1 diff)
• TreeWriter.cc (modified) (10 diffs)
• UniqueObjectFinder.cc (modified) (5 diffs)
• UniqueObjectFinder.h (modified) (2 diffs)

modules/AngularSmearing.cc

@@ -100 +100 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi;
+  Double_t pt, eta, phi, e;
 
   fItInputArray->Reset();
@@ -110 +110 @@
     phi = candidatePosition.Phi();
     pt = candidateMomentum.Pt();
+    e = candidateMomentum.E();
 
     // apply smearing formula for eta,phi
 
-    eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta));
-    phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta));
+    eta = gRandom->Gaus(eta, fFormulaEta->Eval(pt, eta, phi, e));
+    phi = gRandom->Gaus(phi, fFormulaPhi->Eval(pt, eta, phi, e));
     
     if(pt <= 0.0) continue;

modules/BTagging.cc

@@ -22 +22 @@
  *  Determines origin of jet,
  *  applies b-tagging efficiency (miss identification rate) formulas
- *  and sets b-tagging flags 
+ *  and sets b-tagging flags
  *
  *  \author P. Demin - UCL, Louvain-la-Neuve
@@ -34 +34 @@
 #include "classes/DelphesFormula.h"
 
-#include "ExRootAnalysis/ExRootResult.h"
-#include "ExRootAnalysis/ExRootFilter.h"
-#include "ExRootAnalysis/ExRootClassifier.h"
-
 #include "TMath.h"
 #include "TString.h"
@@ -46 +42 @@
 #include "TLorentzVector.h"
 
-#include <algorithm> 
+#include <algorithm>
 #include <stdexcept>
 #include <iostream>
@@ -55 +51 @@
 //------------------------------------------------------------------------------
 
-class BTaggingPartonClassifier : public ExRootClassifier
+BTagging::BTagging() :
+  fItJetInputArray(0)
 {
-public:
-
-  BTaggingPartonClassifier() {}
-
-  Int_t GetCategory(TObject *object);
-  
-  Double_t fEtaMax, fPTMin;
-};
-
-//------------------------------------------------------------------------------
-
-Int_t BTaggingPartonClassifier::GetCategory(TObject *object)
-{
-  Candidate *parton = static_cast<Candidate*>(object);
-  const TLorentzVector &momentum = parton->Momentum;
-  Int_t pdgCode;
-
-  if(momentum.Pt() <= fPTMin || TMath::Abs(momentum.Eta()) > fEtaMax) return -1;
-  
-  pdgCode = TMath::Abs(parton->PID);
-  if(pdgCode != 21 && pdgCode > 5) return -1;
-
-  return 0;
-}
-
-//------------------------------------------------------------------------------
-
-BTagging::BTagging() :
-  fClassifier(0), fFilter(0),
-  fItPartonInputArray(0), fItJetInputArray(0)
-{
-  fClassifier = new BTaggingPartonClassifier;
 }
 
@@ -95 +60 @@
 BTagging::~BTagging()
 {
-  if(fClassifier) delete fClassifier;
 }
 
@@ -109 +73 @@
   fBitNumber = GetInt("BitNumber", 0);
 
-  fDeltaR = GetDouble("DeltaR", 0.5);
-
-  fClassifier->fPTMin = GetDouble("PartonPTMin", 1.0);
-  fClassifier->fEtaMax = GetDouble("PartonEtaMax", 2.5);
-
   // read efficiency formulas
   param = GetParam("EfficiencyFormula");
   size = param.GetSize();
-  
+
   fEfficiencyMap.clear();
   for(i = 0; i < size/2; ++i)
@@ -139 +98 @@
   // import input array(s)
 
-  fPartonInputArray = ImportArray(GetString("PartonInputArray", "Delphes/partons"));
-  fItPartonInputArray = fPartonInputArray->MakeIterator();
-
-  fFilter = new ExRootFilter(fPartonInputArray);
-  
   fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
   fItJetInputArray = fJetInputArray->MakeIterator();
@@ -155 +109 @@
   DelphesFormula *formula;
 
-  if(fFilter) delete fFilter;
   if(fItJetInputArray) delete fItJetInputArray;
-  if(fItPartonInputArray) delete fItPartonInputArray;
 
   for(itEfficiencyMap = fEfficiencyMap.begin(); itEfficiencyMap != fEfficiencyMap.end(); ++itEfficiencyMap)
@@ -170 +122 @@
 void BTagging::Process()
 {
-  Candidate *jet, *parton;
-  Double_t pt, eta, phi;
-  TObjArray *partonArray;
+  Candidate *jet;
+  Double_t pt, eta, phi, e;
   map< Int_t, DelphesFormula * >::iterator itEfficiencyMap;
   DelphesFormula *formula;
-  Int_t pdgCode, pdgCodeMax;
 
-  // select quark and gluons
-  fFilter->Reset();
-  partonArray = fFilter->GetSubArray(fClassifier, 0);
-  
-  if(partonArray == 0) return;
-
-  TIter itPartonArray(partonArray);
-  
   // loop over all input jets
   fItJetInputArray->Reset();
@@ -190 +132 @@
   {
     const TLorentzVector &jetMomentum = jet->Momentum;
-    pdgCodeMax = -1;
     eta = jetMomentum.Eta();
     phi = jetMomentum.Phi();
     pt = jetMomentum.Pt();
+    e = jetMomentum.E();
 
-    // loop over all input partons
-    itPartonArray.Reset();
-    while((parton = static_cast<Candidate*>(itPartonArray.Next())))
-    {
-      pdgCode = TMath::Abs(parton->PID);
-      if(pdgCode == 21) pdgCode = 0;
-      if(jetMomentum.DeltaR(parton->Momentum) <= fDeltaR)
-      {
-        if(pdgCodeMax < pdgCode) pdgCodeMax = pdgCode;
-      }
-    }
-    if(pdgCodeMax == 0) pdgCodeMax = 21;
-    if(pdgCodeMax == -1) pdgCodeMax = 0;
-
-    // find an efficency formula
-    itEfficiencyMap = fEfficiencyMap.find(pdgCodeMax);
+    // find an efficiency formula
+    itEfficiencyMap = fEfficiencyMap.find(jet->Flavor);
     if(itEfficiencyMap == fEfficiencyMap.end())
     {
@@ -217 +145 @@
     formula = itEfficiencyMap->second;
 
-    // apply an efficency formula
-    jet->BTag |= (gRandom->Uniform() <= formula->Eval(pt, eta)) << fBitNumber;
+    // apply an efficiency formula
+    jet->BTag |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber;
+
+    // find an efficiency formula for algo flavor definition
+    itEfficiencyMap = fEfficiencyMap.find(jet->FlavorAlgo);
+    if(itEfficiencyMap == fEfficiencyMap.end())
+    {
+      itEfficiencyMap = fEfficiencyMap.find(0);
+    }
+    formula = itEfficiencyMap->second;
+
+    // apply an efficiency formula
+    jet->BTagAlgo |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber;
+
+    // find an efficiency formula for phys flavor definition
+    itEfficiencyMap = fEfficiencyMap.find(jet->FlavorPhys);
+    if(itEfficiencyMap == fEfficiencyMap.end())
+    {
+      itEfficiencyMap = fEfficiencyMap.find(0);
+    }
+    formula = itEfficiencyMap->second;
+
+    // apply an efficiency formula
+    jet->BTagPhys |= (gRandom->Uniform() <= formula->Eval(pt, eta, phi, e)) << fBitNumber;
   }
 }

modules/BTagging.h

@@ -37 +37 @@
 class DelphesFormula;
 
-class ExRootFilter;
-class BTaggingPartonClassifier;
-
 class BTagging: public DelphesModule
 {
@@ -55 +52 @@
   Int_t fBitNumber;
 
-  Double_t fDeltaR;
-
 #if !defined(__CINT__) && !defined(__CLING__)
   std::map< Int_t, DelphesFormula * > fEfficiencyMap; //!
 #endif
 
-  BTaggingPartonClassifier *fClassifier; //!
-  
-  ExRootFilter *fFilter;
-
-  TIterator *fItPartonInputArray; //!
-  
   TIterator *fItJetInputArray; //!
-
-  const TObjArray *fPartonInputArray; //!
   
   const TObjArray *fJetInputArray; //!

modules/CMakeLists.txt

@@ -7 +7 @@
 file(GLOB sources *.cc)
 file(GLOB headers *.h)
+list(REMOVE_ITEM headers ${CMAKE_CURRENT_SOURCE_DIR}/FastJetLinkDef.h)
 list(REMOVE_ITEM headers ${CMAKE_CURRENT_SOURCE_DIR}/ModulesLinkDef.h)
 list(REMOVE_ITEM headers ${CMAKE_CURRENT_SOURCE_DIR}/Pythia8LinkDef.h)
 
+DELPHES_GENERATE_DICTIONARY(FastJetDict ${headers} LINKDEF FastJetLinkDef.h)
 DELPHES_GENERATE_DICTIONARY(ModulesDict ${headers} LINKDEF ModulesLinkDef.h)
 
@@ -15 +17 @@
 list(REMOVE_ITEM sources ${CMAKE_CURRENT_SOURCE_DIR}/PileUpMergerPythia8.cc)
 
-add_library(modules OBJECT ${sources} ModulesDict.cxx)
+add_library(modules OBJECT ${sources} FastJetDict.cxx ModulesDict.cxx)

modules/Calorimeter.cc

@@ -142 +142 @@
   }
 
-/*
-  TFractionMap::iterator itFractionMap;
-  for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
-  {
-    cout << itFractionMap->first << "   " << itFractionMap->second.first  << "   " << itFractionMap->second.second << endl;
-  }
-*/
+  // read min E value for timing measurement in ECAL
+  fTimingEnergyMin = GetDouble("TimingEnergyMin",4.);
+  // For timing
+  // So far this flag needs to be false
+  // Curved extrapolation not supported
+  fElectronsFromTrack = false;
 
   // read min E value for towers to be saved
@@ -158 +157 @@
 
   // switch on or off the dithering of the center of calorimeter towers
-  fDitherTowerCenter = GetBool("DitherTowerCenter", true);
+  fSmearTowerCenter = GetBool("SmearTowerCenter", true);
 
   // read resolution formulas
@@ -356 +355 @@
       fTrackHCalEnergy = 0.0;
 
-      fTowerECalTime = 0.0;
-      fTowerHCalTime = 0.0;
-
-      fTrackECalTime = 0.0;
-      fTrackHCalTime = 0.0;
-
-      fTowerECalTimeWeight = 0.0;
-      fTowerHCalTimeWeight = 0.0;
-
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -380 +370 @@
       position = track->Position;
 
-
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -387 +376 @@
       fTrackHCalEnergy += hcalEnergy;
 
-      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
-      fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy);
+      if(ecalEnergy > fTimingEnergyMin && fTower)
+      {
+        if(fElectronsFromTrack)
+        {
+          fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, track->Position.T()));
+        }
+      }
 
       fTowerTrackArray->Add(track);
@@ -412 +403 @@
     fTowerHCalEnergy += hcalEnergy;
 
-    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
-    fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy);
-
+    if(ecalEnergy > fTimingEnergyMin && fTower)
+    {
+      if (abs(particle->PID) != 11 || !fElectronsFromTrack)
+      {
+        fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
+      }
+    }
 
     fTower->AddCandidate(particle);
@@ -434 +426 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
-  Double_t ecalTime, hcalTime, time;
+  Float_t weight, sumWeightedTime, sumWeight;
 
   if(!fTower) return;
@@ -444 +436 @@
   hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
 
-  ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
-  hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
-
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
@@ -454 +443 @@
 
   energy = ecalEnergy + hcalEnergy;
-  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
-
-  if(fDitherTowerCenter)
+
+  if(fSmearTowerCenter)
   {
     eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
@@ -469 +457 @@
   pt = energy / TMath::CosH(eta);
 
-  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
+  // Time calculation for tower
+  fTower->NTimeHits = 0;
+  sumWeightedTime = 0.0;
+  sumWeight = 0.0;
+
+  for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
+  {
+    weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
+    sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
+    sumWeight += weight;
+    fTower->NTimeHits++;
+  }
+
+  if(sumWeight > 0.0)
+  {
+    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
+  }
+  else
+  {
+    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
+  }
+
+
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;

modules/Calorimeter.h

@@ -60 +60 @@
   Double_t fTrackECalEnergy, fTrackHCalEnergy;
 
-  Double_t fTowerECalTime, fTowerHCalTime;
-  Double_t fTrackECalTime, fTrackHCalTime;
-
-  Double_t fTowerECalTimeWeight, fTowerHCalTimeWeight;
-  Double_t fTrackECalTimeWeight, fTrackHCalTimeWeight;
+  Double_t fTimingEnergyMin;
+  Bool_t fElectronsFromTrack;
 
   Int_t fTowerTrackHits, fTowerPhotonHits;
@@ -74 +71 @@
   Double_t fHCalEnergySignificanceMin;
 
-  Bool_t fDitherTowerCenter;
+  Bool_t fSmearTowerCenter;
 
   TFractionMap fFractionMap; //!

modules/Efficiency.cc

@@ -96 +96 @@
 { 
   Candidate *candidate;
-  Double_t pt, eta, phi;
+  Double_t pt, eta, phi, e;
 
   fItInputArray->Reset();
@@ -106 +106 @@
     phi = candidatePosition.Phi();
     pt = candidateMomentum.Pt();
+    e = candidateMomentum.E();
 
     // apply an efficency formula
-    if(gRandom->Uniform() > fFormula->Eval(pt, eta)) continue;
+    if(gRandom->Uniform() > fFormula->Eval(pt, eta, phi, e)) continue;
     
     fOutputArray->Add(candidate);

modules/EnergyScale.cc

@@ -104 +104 @@
     momentum = candidate->Momentum;
 
-    scale = fFormula->Eval(momentum.Pt(), momentum.Eta());
+    scale = fFormula->Eval(momentum.Pt(), momentum.Eta(), momentum.Phi(), momentum.E());
 
     if(scale > 0.0) momentum *= scale;

modules/EnergySmearing.cc

@@ -96 +96 @@
 {
   Candidate *candidate, *mother;
-  Double_t energy, eta, phi;
+  Double_t pt, energy, eta, phi;
 
   fItInputArray->Reset();
@@ -103 +103 @@
     const TLorentzVector &candidatePosition = candidate->Position;
     const TLorentzVector &candidateMomentum = candidate->Momentum;
+    
+    pt = candidatePosition.Pt();
     eta = candidatePosition.Eta();
     phi = candidatePosition.Phi();
@@ -108 +110 @@
  
     // apply smearing formula
-    energy = gRandom->Gaus(energy, fFormula->Eval(0.0, eta, 0.0, energy));
+    energy = gRandom->Gaus(energy, fFormula->Eval(pt, eta, phi, energy));
      
     if(energy <= 0.0) continue;

modules/FastJetFinder.cc

@@ -66 +66 @@
 #include "fastjet/contribs/Nsubjettiness/WinnerTakeAllRecombiner.hh"
 
+#include "fastjet/tools/Filter.hh"
+#include "fastjet/tools/Pruner.hh"
+#include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
+
 using namespace std;
 using namespace fastjet;
@@ -92 +96 @@
   JetDefinition::Plugin *plugin = 0;
   JetDefinition::Recombiner *recomb = 0;
-  NjettinessPlugin *njetPlugin = 0;
-
-  // read eta ranges
-
-  ExRootConfParam param = GetParam("RhoEtaRange");
+  ExRootConfParam param;
   Long_t i, size;
-
-  fEtaRangeMap.clear();
-  size = param.GetSize();
-  for(i = 0; i < size/2; ++i)
-  {
-    fEtaRangeMap[param[i*2].GetDouble()] = param[i*2 + 1].GetDouble();
-  }
+  Double_t etaMin, etaMax;
+  TEstimatorStruct estimatorStruct;
 
   // define algorithm
@@ -130 +125 @@
   fN = GetInt("N", 2);                  // used only if Njettiness is used as jet clustering algo (case 8)
 
+  //-- 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);
@@ -197 +214 @@
       break;
     case 8:
-      njetPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
-      fDefinition = new JetDefinition(njetPlugin);
+      fNjettinessPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
+      fDefinition = new JetDefinition(fNjettinessPlugin);
       break;
   }
@@ -204 +221 @@
   fPlugin = plugin;
   fRecomb = recomb;
-  fNjettinessPlugin = njetPlugin;
 
   ClusterSequence::print_banner();
+
+  if(fComputeRho && fAreaDefinition)
+  {
+    // read eta ranges
+
+    param = GetParam("RhoEtaRange");
+    size = param.GetSize();
+
+    fEstimators.clear();
+    for(i = 0; i < size/2; ++i)
+    {
+      etaMin = param[i*2].GetDouble();
+      etaMax = param[i*2 + 1].GetDouble();
+      estimatorStruct.estimator = new JetMedianBackgroundEstimator(SelectorEtaRange(etaMin, etaMax), *fDefinition, *fAreaDefinition);
+      estimatorStruct.etaMin = etaMin;
+      estimatorStruct.etaMax = etaMax;
+      fEstimators.push_back(estimatorStruct);
+    }
+  }
 
   // import input array
@@ -223 +258 @@
 void FastJetFinder::Finish()
 {
+  vector< TEstimatorStruct >::iterator itEstimators;
+
+  for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
+  {
+    if(itEstimators->estimator) delete itEstimators->estimator;
+  }
+
   if(fItInputArray) delete fItInputArray;
   if(fDefinition) delete fDefinition;
@@ -241 +283 @@
   Double_t time, timeWeight;
   Int_t number;
+  Int_t charge; 
   Double_t rho = 0.0;
   PseudoJet jet, area;
-  vector<PseudoJet> inputList, outputList;
   ClusterSequence *sequence;
-  map< Double_t, Double_t >::iterator itEtaRangeMap;
+  vector< PseudoJet > inputList, outputList, subjets;
+  vector< PseudoJet >::iterator itInputList, itOutputList;
+  vector< TEstimatorStruct >::iterator itEstimators;
 
   DelphesFactory *factory = GetFactory();
@@ -276 +320 @@
   if(fComputeRho && fAreaDefinition)
   {
-    for(itEtaRangeMap = fEtaRangeMap.begin(); itEtaRangeMap != fEtaRangeMap.end(); ++itEtaRangeMap)
-    {
-      Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
-      JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
-      estimator.set_particles(inputList);
-      rho = estimator.rho();
+    for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
+    {
+      itEstimators->estimator->set_particles(inputList);
+      rho = itEstimators->estimator->rho();
 
       candidate = factory->NewCandidate();
       candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
-      candidate->Edges[0] = itEtaRangeMap->first;
-      candidate->Edges[1] = itEtaRangeMap->second;
+      candidate->Edges[0] = itEstimators->etaMin;
+      candidate->Edges[1] = itEstimators->etaMax;
       fRhoOutputArray->Add(candidate);
     }
@@ -298 +340 @@
   detaMax = 0.0;
   dphiMax = 0.0;
-  vector<PseudoJet>::iterator itInputList, itOutputList;
   for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
   {
@@ -314 +355 @@
     timeWeight = 0.0;
 
+    charge = 0;
+
     inputList.clear();
     inputList = sequence->constituents(*itOutputList);
@@ -319 +362 @@
     for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
     {
+      if(itInputList->user_index() < 0) continue;
       constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
 
@@ -329 +373 @@
       timeWeight += TMath::Sqrt(constituent->Momentum.E());
 
+      charge += constituent->Charge;
+
       candidate->AddCandidate(constituent);
     }
@@ -338 +384 @@
     candidate->DeltaEta = detaMax;
     candidate->DeltaPhi = dphiMax;
-
+    candidate->Charge = charge; 
+    //------------------------------------
+    // Trimming
+    //------------------------------------
+
+     if(fComputeTrimming)
+    {
+
+      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 
+      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 ; 
+        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::PseudoJet pruned_jet = pruner(*itOutputList);
+
+      candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
+         
+      // 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 ; 
+        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 
+      
+      subjets.clear();
+      subjets    = softdrop_jet.pieces();
+      subjets    = sorted_by_pt(subjets);
+      candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
+
+      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());
+      }
+    }
+  
     // --- compute N-subjettiness with N = 1,2,3,4,5 ----
 
@@ -377 +503 @@
     }
 
-
     fOutputArray->Add(candidate);
   }

modules/FastJetFinder.h

@@ -30 +30 @@
 #include "classes/DelphesModule.h"
 
-#include <map>
+#include <vector>
 
 class TObjArray;
@@ -38 +38 @@
   class JetDefinition;
   class AreaDefinition;
-  class Selector;
+  class JetMedianBackgroundEstimator;
   namespace contrib {
     class NjettinessPlugin;
@@ -83 +83 @@
   Int_t fN ;
 
+  //-- Trimming parameters --
+  
+  Bool_t fComputeTrimming;
+  Double_t fRTrim;
+  Double_t fPtFracTrim;
+  
+  //-- Pruning parameters --
+
+  Bool_t fComputePruning;
+  Double_t fZcutPrun;
+  Double_t fRcutPrun;
+  Double_t fRPrun;
+
+  //-- SoftDrop parameters --
+
+  Bool_t fComputeSoftDrop;
+  Double_t fBetaSoftDrop;
+  Double_t fSymmetryCutSoftDrop;
+  Double_t fR0SoftDrop;
+
   // --- FastJet Area method --------
 
@@ -100 +120 @@
   Double_t fEffectiveRfact;
 
-  std::map< Double_t, Double_t > fEtaRangeMap; //!
+#if !defined(__CINT__) && !defined(__CLING__)
+  struct TEstimatorStruct
+  {
+    fastjet::JetMedianBackgroundEstimator *estimator;
+    Double_t etaMin, etaMax;
+  };
+
+  std::vector< TEstimatorStruct > fEstimators; //!
+#endif
 
   TIterator *fItInputArray; //!

modules/FastJetLinkDef.h

@@ -28 +28 @@
 #include "modules/FastJetFinder.h"
 #include "modules/FastJetGridMedianEstimator.h"
+#include "modules/RunPUPPI.h"
 
 #ifdef __CINT__
@@ -37 +38 @@
 #pragma link C++ class FastJetFinder+;
 #pragma link C++ class FastJetGridMedianEstimator+;
+#pragma link C++ class RunPUPPI+;
 
 #endif

modules/IdentificationMap.cc

@@ -69 +69 @@
 void IdentificationMap::Init()
 {
-  // read efficiency formula
-
-
   TMisIDMap::iterator itEfficiencyMap;
   ExRootConfParam param;
@@ -87 +84 @@
     formula->Compile(param[i*3 + 2].GetString());
     pdg = param[i*3].GetInt();
-    fEfficiencyMap.insert(make_pair(pdg,make_pair(param[i*3 + 1].GetInt(),formula)));
-
-   // cout<<param[i*3].GetInt()<<","<<param[i*3+1].GetInt()<<","<<param[i*3 + 2].GetString()<<endl;
-
+    fEfficiencyMap.insert(make_pair(pdg, make_pair(param[i*3 + 1].GetInt(), formula)));
   }
 
@@ -100 +94 @@
     formula->Compile("1.0");
 
-    fEfficiencyMap.insert(make_pair(0,make_pair(0,formula)));
+    fEfficiencyMap.insert(make_pair(0, make_pair(0, formula)));
   }
 
@@ -126 +120 @@
     if(formula) delete formula;
   }
-
 }
 
@@ -134 +127 @@
 {
   Candidate *candidate;
-  Double_t pt, eta, phi;
+  Double_t pt, eta, phi, e;
   TMisIDMap::iterator itEfficiencyMap;
   pair <TMisIDMap::iterator, TMisIDMap::iterator> range;
   DelphesFormula *formula;
-  Int_t pdgIn, pdgOut, charge;
+  Int_t pdgCodeIn, pdgCodeOut, charge;
 
-  Double_t P, Pi;
-
- // cout<<"------------ New Event ------------"<<endl;
+  Double_t p, r, total;
 
   fItInputArray->Reset();
@@ -152 +143 @@
     phi = candidatePosition.Phi();
     pt = candidateMomentum.Pt();
-    pdgIn = candidate->PID;
+    e = candidateMomentum.E();
+   
+    pdgCodeIn = candidate->PID;
     charge = candidate->Charge;
 
-   // cout<<"------------ New Candidate ------------"<<endl;
-   // cout<<candidate->PID<<"   "<<pt<<","<<eta<<","<<phi<<endl;
+    // first check that PID of this particle is specified in the map
+    // otherwise, look for PID = 0
 
-    P = 1.0;
+    itEfficiencyMap = fEfficiencyMap.find(pdgCodeIn);
 
-    //first check that PID of this particle is specified in cfg, if not set look for PID=0
-
-    itEfficiencyMap = fEfficiencyMap.find(pdgIn);
-
-    range = fEfficiencyMap.equal_range(pdgIn);
-    if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgIn);
+    range = fEfficiencyMap.equal_range(pdgCodeIn);
+    if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgCodeIn);
     if(range.first == range.second) range = fEfficiencyMap.equal_range(0);
 
-    //loop over submap for this pid
-    for (TMisIDMap::iterator it=range.first; it!=range.second; ++it)
+    r = gRandom->Uniform();
+    total = 0.0;
+
+    // loop over sub-map for this PID
+    for(TMisIDMap::iterator it = range.first; it != range.second; ++it)
     {
+      formula = (it->second).second;
+      pdgCodeOut = (it->second).first;
 
-      formula = (it->second).second;
-      pdgOut = (it->second).first;
+      p = formula->Eval(pt, eta, phi, e);
 
-      Pi = formula->Eval(pt, eta);
-
-      // check that sum of probabilities does not exceed 1.
-      P = (P - Pi)/P;
-
-      if( P < 0.0 ) continue;
-      else
+      if(total <= r && r < total + p)
       {
-
-       //randomly assign a PID to particle according to map
-       Double_t rndm = gRandom->Uniform();
-
-       if(rndm > P)
-       {
-         //change PID of particle
-         if(pdgOut != 0) candidate->PID = charge*pdgOut;
-         fOutputArray->Add(candidate);
-         break;
-       }
+        // change PID of particle
+        if(pdgCodeOut != 0) candidate->PID = charge*pdgCodeOut;
+        fOutputArray->Add(candidate);
+        break;
       }
 
+      total += p;
     }
-
-   }
-
+  }
 }
 

modules/IdentificationMap.h

@@ -49 +49 @@
 private:
 
-  typedef std::multimap< Int_t, std::pair<Int_t , DelphesFormula * > > TMisIDMap; //!
+  typedef std::multimap< Int_t, std::pair< Int_t, DelphesFormula * > > TMisIDMap; //!
 
-  TMisIDMap fEfficiencyMap;
+  TMisIDMap fEfficiencyMap; //!
 
   TIterator *fItInputArray; //!

modules/ImpactParameterSmearing.cc

@@ -97 +97 @@
   Candidate *candidate, *particle, *mother;
   Double_t xd, yd, zd, dxy, sx, sy, sz, ddxy;
-  Double_t pt, eta, px, py;
+  Double_t pt, eta, px, py, phi, e;
 
   fItInputArray->Reset();
@@ -110 +110 @@
     eta = candidateMomentum.Eta();
     pt = candidateMomentum.Pt();
+    phi = candidateMomentum.Phi();
+    e = candidateMomentum.E();
+    
     px = candidateMomentum.Px();
     py = candidateMomentum.Py();
@@ -119 +122 @@
 
     // calculate smeared values
-    sx = gRandom->Gaus(0.0, fFormula->Eval(pt, eta));
-    sy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta));
-    sz = gRandom->Gaus(0.0, fFormula->Eval(pt, eta));
+    sx = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
+    sy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
+    sz = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
 
     xd += sx;
@@ -130 +133 @@
     dxy = (xd*py - yd*px)/pt;
 
-    ddxy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta));
+    ddxy = gRandom->Gaus(0.0, fFormula->Eval(pt, eta, phi, e));
 
     // fill smeared values in candidate

modules/Isolation.cc

@@ -25 +25 @@
  *  the transverse momenta fraction within (PTRatioMin, PTRatioMax].
  *
- *  \author P. Demin - UCL, Louvain-la-Neuve
+ *  \author P. Demin, M. Selvaggi, R. Gerosa - UCL, Louvain-la-Neuve
  *
  */
@@ -153 +153 @@
   Candidate *candidate, *isolation, *object;
   TObjArray *isolationArray;
-  Double_t sum, ratio;
+  Double_t sumCharged, sumNeutral, sumAllParticles, sumChargedPU, sumDBeta, ratioDBeta, sumRhoCorr, ratioRhoCorr;
   Int_t counter;
   Double_t eta = 0.0;
   Double_t rho = 0.0;
-
-  if(fRhoInputArray && fRhoInputArray->GetEntriesFast() > 0)
-  {
-    candidate = static_cast<Candidate*>(fRhoInputArray->At(0));
-    rho = candidate->Momentum.Pt();
-  }
 
   // select isolation objects
@@ -178 +172 @@
     const TLorentzVector &candidateMomentum = candidate->Momentum;
     eta = TMath::Abs(candidateMomentum.Eta());
-
-    // loop over all input tracks
-    sum = 0.0;
-    counter = 0;
-    itIsolationArray.Reset();
-    while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
-    {
-      const TLorentzVector &isolationMomentum = isolation->Momentum;
-
-      if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
-         candidate->GetUniqueID() != isolation->GetUniqueID())
-      {
-        sum += isolationMomentum.Pt();
-        ++counter;
-      }
-    }
 
     // find rho
@@ -209 +187 @@
     }
 
-    // correct sum for pile-up contamination
-    sum = sum - rho*fDeltaRMax*fDeltaRMax*TMath::Pi();
-
-    ratio = sum/candidateMomentum.Pt();
-    if((fUsePTSum && sum > fPTSumMax) || ratio > fPTRatioMax) continue;
-
+    // loop over all input tracks
+    
+    sumNeutral = 0.0;
+    sumCharged = 0.0;
+    sumChargedPU = 0.0;
+    sumAllParticles = 0.0;
+   
+    counter = 0;
+    itIsolationArray.Reset();
+    
+    while((isolation = static_cast<Candidate*>(itIsolationArray.Next())))
+    {
+      const TLorentzVector &isolationMomentum = isolation->Momentum;
+
+      if(candidateMomentum.DeltaR(isolationMomentum) <= fDeltaRMax &&
+         candidate->GetUniqueID() != isolation->GetUniqueID())
+      {
+        sumAllParticles += isolationMomentum.Pt();
+        if(isolation->Charge !=0) 
+        { 
+          sumCharged += isolationMomentum.Pt();
+          if(isolation->IsRecoPU != 0) sumChargedPU += isolationMomentum.Pt();
+        } 
+        else
+        {
+          sumNeutral += isolationMomentum.Pt();
+        }
+        ++counter;
+      }
+    }
+
+    // find rho
+    rho = 0.0;
+    if(fRhoInputArray)
+    {
+      fItRhoInputArray->Reset();
+      while((object = static_cast<Candidate*>(fItRhoInputArray->Next())))
+      {
+        if(eta >= object->Edges[0] && eta < object->Edges[1])
+        {
+          rho = object->Momentum.Pt();
+        }
+      }
+    }
+
+     // correct sum for pile-up contamination
+    sumDBeta = sumCharged + TMath::Max(sumNeutral-0.5*sumChargedPU,0.0);
+    sumRhoCorr = sumCharged + TMath::Max(sumNeutral-TMath::Max(rho,0.0)*fDeltaRMax*fDeltaRMax*TMath::Pi(),0.0);
+    ratioDBeta = sumDBeta/candidateMomentum.Pt();
+    ratioRhoCorr = sumRhoCorr/candidateMomentum.Pt();
+    
+    candidate->IsolationVar = ratioDBeta;
+    candidate->IsolationVarRhoCorr = ratioRhoCorr;
+    candidate->SumPtCharged = sumCharged;
+    candidate->SumPtNeutral = sumNeutral;
+    candidate->SumPtChargedPU = sumChargedPU;
+    candidate->SumPt = sumAllParticles;
+
+    if((fUsePTSum && sumDBeta > fPTSumMax) || (!fUsePTSum && ratioDBeta > fPTRatioMax)) continue;
     fOutputArray->Add(candidate);
   }

modules/JetPileUpSubtractor.cc

@@ -109 +109 @@
     momentum = candidate->Momentum;
     area = candidate->Area;
-    eta = TMath::Abs(momentum.Eta());
+    eta = momentum.Eta();
 
     // find rho

modules/ModulesLinkDef.h

@@ -29 +29 @@
 
 #include "modules/AngularSmearing.h"
+#include "modules/PhotonConversions.h"
 #include "modules/ParticlePropagator.h"
 #include "modules/Efficiency.h"
@@ -50 +51 @@
 #include "modules/JetPileUpSubtractor.h"
 #include "modules/TrackPileUpSubtractor.h"
+#include "modules/TaggingParticlesSkimmer.h"
 #include "modules/PileUpJetID.h"
 #include "modules/ConstituentFilter.h"
@@ -57 +59 @@
 #include "modules/Weighter.h"
 #include "modules/Hector.h"
+#include "modules/JetFlavorAssociation.h"
+#include "modules/JetFakeParticle.h"
 #include "modules/ExampleModule.h"
 
@@ -68 +72 @@
 
 #pragma link C++ class AngularSmearing+;
+#pragma link C++ class PhotonConversions+;
 #pragma link C++ class ParticlePropagator+;
 #pragma link C++ class Efficiency+;
@@ -89 +94 @@
 #pragma link C++ class JetPileUpSubtractor+;
 #pragma link C++ class TrackPileUpSubtractor+;
+#pragma link C++ class TaggingParticlesSkimmer+;
 #pragma link C++ class PileUpJetID+;
 #pragma link C++ class ConstituentFilter+;
@@ -96 +102 @@
 #pragma link C++ class Weighter+;
 #pragma link C++ class Hector+;
+#pragma link C++ class JetFlavorAssociation+;
+#pragma link C++ class JetFakeParticle+;
 #pragma link C++ class ExampleModule+;
 

modules/MomentumSmearing.cc

@@ -96 +96 @@
 {
   Candidate *candidate, *mother;
-  Double_t pt, eta, phi;
+  Double_t pt, eta, phi, e;
 
   fItInputArray->Reset();
@@ -106 +106 @@
     phi = candidatePosition.Phi();
     pt = candidateMomentum.Pt();
+    e = candidateMomentum.E();
 
     // apply smearing formula
-    pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta) * pt);
+    pt = gRandom->Gaus(pt, fFormula->Eval(pt, eta, phi, e) * pt);
     
     if(pt <= 0.0) continue;

modules/PdgCodeFilter.cc

@@ -74 +74 @@
   fPTMin = GetDouble("PTMin", 0.0);
 
+  fInvertPdg = GetBool("InvertPdg", false);
+
+  fRequireStatus = GetBool("RequireStatus", false);
+  fStatus = GetInt("Status", 1);
+
   // import input array
   fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
@@ -109 +114 @@
   Double_t pt;
 
+  const Bool_t requireStatus = fRequireStatus;
+  const Bool_t invertPdg = fInvertPdg;
+  const int reqStatus = fStatus;
+
   fItInputArray->Reset();
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
@@ -116 +125 @@
     pt = candidateMomentum.Pt();
 
+    if(pt < fPTMin) continue;
+    if(requireStatus && (candidate->Status != reqStatus)) continue;
+
     pass = kTRUE;
-
-    if(pt < fPTMin) pass = kFALSE;
     if(find(fPdgCodes.begin(), fPdgCodes.end(), pdgCode) != fPdgCodes.end()) pass = kFALSE;
 
+    if (invertPdg) pass = !pass;
     if(pass) fOutputArray->Add(candidate);
   }

modules/PdgCodeFilter.h

@@ -50 +50 @@
 
   Double_t fPTMin; //!
+  Bool_t fInvertPdg; //!
+  Bool_t fRequireStatus; //!
+  Int_t fStatus; //!
 
   std::vector<Int_t> fPdgCodes;

modules/PileUpJetID.cc

@@ -1 @@
-/*
- *  Delphes: a framework for fast simulation of a generic collider experiment
- *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
- *
- *  This program is free software: you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation, either version 3 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
 
 /** \class PileUpJetID
  *
- *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
+ *  CMS PileUp Jet ID Variables
  *
- *  \author S. Zenz, December 2013
+ *  \author S. Zenz
  *
  */
@@ -41 +23 @@
 #include "TRandom3.h"
 #include "TObjArray.h"
-#include "TDatabasePDG.h"
+//#include "TDatabasePDG.h"
 #include "TLorentzVector.h"
 
@@ -54 +36 @@
 
 PileUpJetID::PileUpJetID() :
-  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
+  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
 {
 
@@ -70 +52 @@
 void PileUpJetID::Init()
 {
+
   fJetPTMin = GetDouble("JetPTMin", 20.0);
   fParameterR = GetDouble("ParameterR", 0.5);
   fUseConstituents = GetInt("UseConstituents", 0);
 
+  fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
+  fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
+  fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
+  fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
+  fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
+  fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
+  fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
+
   fAverageEachTower = false; // for timing
 
@@ -81 +72 @@
   fItJetInputArray = fJetInputArray->MakeIterator();
 
-  fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
+  fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
   fItTrackInputArray = fTrackInputArray->MakeIterator();
 
-  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
+  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
   fItNeutralInputArray = fNeutralInputArray->MakeIterator();
 
-  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
-  fItVertexInputArray = fVertexInputArray->MakeIterator();
-
-  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
-
   // create output array(s)
 
   fOutputArray = ExportArray(GetString("OutputArray", "jets"));
+
+  fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
+
 }
 
@@ -101 +90 @@
 void PileUpJetID::Finish()
 {
+  //  cout << "In finish" << endl;
+
   if(fItJetInputArray) delete fItJetInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
   if(fItNeutralInputArray) delete fItNeutralInputArray;
-  if(fItVertexInputArray) delete fItVertexInputArray;
+
 }
 
@@ -113 +104 @@
   Candidate *candidate, *constituent;
   TLorentzVector momentum, area;
-  Int_t i, nc, nn;
-  Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq;
-  Double_t dr, pt, pt_ann[5];
-  Double_t zvtx = 0.0;
-
-  Candidate *track;
-
-  // find z position of primary vertex
-
-  fItVertexInputArray->Reset();
-  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
-  {
-    if(!candidate->IsPU)
-    {
-      zvtx = candidate->Position.Z();
-      break;
-    }
-  }
+
+  Candidate *trk;
 
   // loop over all input candidates
@@ -139 +114 @@
     area = candidate->Area;
 
-    sumpt = 0.0;
-    sumptch = 0.0;
-    sumptchpv = 0.0;
-    sumptchpu = 0.0;
-    sumdrsqptsq = 0.0;
-    sumptsq = 0.0;
-    nc = 0;
-    nn = 0;
-
-    for(i = 0; i < 5; ++i)
-    {
-      pt_ann[i] = 0.0;
-    }
-
-    if(fUseConstituents)
-    {
+    float sumT0 = 0.;
+    float sumT1 = 0.;
+    float sumT10 = 0.;
+    float sumT20 = 0.;
+    float sumT30 = 0.;
+    float sumT40 = 0.;
+    float sumWeightsForT = 0.;
+    candidate->NTimeHits = 0;
+
+    float sumpt = 0.;
+    float sumptch = 0.;
+    float sumptchpv = 0.;
+    float sumptchpu = 0.;
+    float sumdrsqptsq = 0.;
+    float sumptsq = 0.;
+    int nc = 0;
+    int nn = 0;
+    float pt_ann[5];
+
+    for (int i = 0 ; i < 5 ; i++) {
+      pt_ann[i] = 0.;
+    }
+
+    if (fUseConstituents) {
       TIter itConstituents(candidate->GetCandidates());
-      while((constituent = static_cast<Candidate*>(itConstituents.Next())))
-      {
-        pt = constituent->Momentum.Pt();
-        dr = candidate->Momentum.DeltaR(constituent->Momentum);
-        sumpt += pt;
-        sumdrsqptsq += dr*dr*pt*pt;
-        sumptsq += pt*pt;
-        if(constituent->Charge == 0)
-        {
-          // neutrals
-          ++nn;
-        }
-        else
-        {
-          // charged
-          if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution)
-          {
-            sumptchpu += pt;
-          }
-          else
-          {
-            sumptchpv += pt;
-          }
-          sumptch += pt;
-          ++nc;
-        }
-        for(i = 0; i < 5; ++i)
-        {
-          if(dr > 0.1*i && dr < 0.1*(i + 1))
-          {
-            pt_ann[i] += pt;
-          }
-        }
-      }
-    }
-    else
-    {
+      while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
+        float pt = constituent->Momentum.Pt();
+        float dr = candidate->Momentum.DeltaR(constituent->Momentum);
+        //      cout << " There exists a constituent with dr=" << dr << endl;
+        sumpt += pt;
+        sumdrsqptsq += dr*dr*pt*pt;
+        sumptsq += pt*pt;
+        if (constituent->Charge == 0) {
+          nn++;
+        } else {
+          if (constituent->IsRecoPU) {
+            sumptchpu += pt;
+          } else {
+            sumptchpv += pt;
+          }
+          sumptch += pt;
+          nc++;
+        }
+        for (int i = 0 ; i < 5 ; i++) {
+          if (dr > 0.1*i && dr < 0.1*(i+1)) {
+            pt_ann[i] += pt;
+          }
+        }
+        float tow_sumT = 0;
+        float tow_sumW = 0;
+        for (int i = 0 ; i < constituent->ECalEnergyTimePairs.size() ; i++) {
+          float w = TMath::Sqrt(constituent->ECalEnergyTimePairs[i].first);
+          if (fAverageEachTower) {
+            tow_sumT += w*constituent->ECalEnergyTimePairs[i].second;
+            tow_sumW += w;
+          } else {
+            sumT0 += w*constituent->ECalEnergyTimePairs[i].second;
+            sumT1 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.001);
+            sumT10 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.010);
+            sumT20 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.020);
+            sumT30 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.030);
+            sumT40 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.040);
+            sumWeightsForT += w;
+            candidate->NTimeHits++;
+          }
+        }
+        if (fAverageEachTower && tow_sumW > 0.) {
+          sumT0 += tow_sumT;
+          sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001);
+          sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010);
+          sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020);
+          sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030);
+          sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040);
+          sumWeightsForT += tow_sumW;
+          candidate->NTimeHits++;
+        }
+      }
+    } else {
       // Not using constituents, using dr
       fItTrackInputArray->Reset();
-      while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
-      {
-        if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR)
-        {
-          pt = track->Momentum.Pt();
-          sumpt += pt;
-          sumptch += pt;
-          if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution)
-          {
-            sumptchpu += pt;
-          }
-          else
-          {
-            sumptchpv += pt;
-          }
-          dr = candidate->Momentum.DeltaR(track->Momentum);
-          sumdrsqptsq += dr*dr*pt*pt;
-          sumptsq += pt*pt;
-          nc++;
-          for(i = 0; i < 5; ++i)
-          {
-            if(dr > 0.1*i && dr < 0.1*(i + 1))
-            {
-              pt_ann[i] += pt;
-            }
-          }
-        }
-      }
-
+      while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
+        if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
+          float pt = trk->Momentum.Pt();
+          sumpt += pt;
+          sumptch += pt;
+          if (trk->IsRecoPU) {
+            sumptchpu += pt;
+          } else {
+            sumptchpv += pt;
+          }
+          float dr = candidate->Momentum.DeltaR(trk->Momentum);
+          sumdrsqptsq += dr*dr*pt*pt;
+          sumptsq += pt*pt;
+          nc++;
+          for (int i = 0 ; i < 5 ; i++) {
+            if (dr > 0.1*i && dr < 0.1*(i+1)) {
+              pt_ann[i] += pt;
+            }
+          }
+        }
+      }
       fItNeutralInputArray->Reset();
-      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next())))
-      {
-        if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
-        {
-          pt = constituent->Momentum.Pt();
-          sumpt += pt;
-          dr = candidate->Momentum.DeltaR(constituent->Momentum);
-          sumdrsqptsq += dr*dr*pt*pt;
-          sumptsq += pt*pt;
-          nn++;
-          for(i = 0; i < 5; ++i)
-          {
-            if(dr > 0.1*i && dr < 0.1*(i + 1))
-            {
-              pt_ann[i] += pt;
-            }
-          }
-        }
-      }
-    }
-
-    if(sumptch > 0.0)
-    {
-      candidate->Beta = sumptchpu/sumptch;
-      candidate->BetaStar = sumptchpv/sumptch;
-    }
-    else
-    {
-      candidate->Beta = -999.0;
-      candidate->BetaStar = -999.0;
-    }
-    if(sumptsq > 0.0)
-    {
+      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
+        if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
+          float pt = constituent->Momentum.Pt();
+          sumpt += pt;
+          float dr = candidate->Momentum.DeltaR(constituent->Momentum);
+          sumdrsqptsq += dr*dr*pt*pt;
+          sumptsq += pt*pt;
+          nn++;
+          for (int i = 0 ; i < 5 ; i++) {
+            if (dr > 0.1*i && dr < 0.1*(i+1)) {
+            pt_ann[i] += pt;
+            }
+          }
+        }
+      }
+    }
+
+    if (sumptch > 0.) {
+      candidate->Beta = sumptchpv/sumptch;
+      candidate->BetaStar = sumptchpu/sumptch;
+    } else {
+      candidate->Beta = -999.;
+      candidate->BetaStar = -999.;
+    }
+    if (sumptsq > 0.) {
       candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
-    }
-    else
-    {
-      candidate->MeanSqDeltaR = -999.0;
+    } else {
+      candidate->MeanSqDeltaR = -999.;
     }
     candidate->NCharged = nc;
     candidate->NNeutrals = nn;
-    if(sumpt > 0.0)
-    {
+    if (sumpt > 0.) {
       candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
-      for(i = 0; i < 5; ++i)
-      {
+      for (int i = 0 ; i < 5 ; i++) {
         candidate->FracPt[i] = pt_ann[i]/sumpt;
       }
-    }
-    else
-    {
-      candidate->PTD = -999.0;
-      for(i = 0; i < 5; ++i)
-      {
-        candidate->FracPt[i] = -999.0;
+    } else {
+      candidate->PTD = -999.;
+      for (int i = 0 ; i < 5 ; i++) {
+        candidate->FracPt[i] = -999.;
       }
     }
 
     fOutputArray->Add(candidate);
+
+    // New stuff
+    /*
+    fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
+    fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
+    fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
+    fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
+    fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
+    */
+
+    bool passId = false;
+    if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) {
+      if (fabs(candidate->Momentum.Eta())<1.5) {
+        passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
+      } else if (fabs(candidate->Momentum.Eta())<4.0) {
+        passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
+      } else {
+        passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
+      }
+    }
+
+    //    cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt() 
+    //   << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
+
+    if (passId) {
+      if (fUseConstituents) {
+        TIter itConstituents(candidate->GetCandidates());
+        while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
+          if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) {
+            fNeutralsInPassingJets->Add(constituent);
+            //      cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
+          }
+        }
+      } else { // use DeltaR
+        fItNeutralInputArray->Reset();
+        while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
+          if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) {
+            fNeutralsInPassingJets->Add(constituent);
+            //            cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
+          }
+        }
+      }
+    }
+
+
   }
 }
 
 //------------------------------------------------------------------------------
-

modules/PileUpJetID.h

@@ -1 @@
-/*
- *  Delphes: a framework for fast simulation of a generic collider experiment
- *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
- *
- *  This program is free software: you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation, either version 3 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
 #ifndef PileUpJetID_h
 #define PileUpJetID_h
@@ -22 +4 @@
 /** \class PileUpJetID
  *
- *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
+ *  CMS PileUp Jet ID Variables
  *
- *  \author S. Zenz, December 2013
+ *  \author S. Zenz
  *
  */
@@ -34 +16 @@
 
 class TObjArray;
+class DelphesFormula;
 
 class PileUpJetID: public DelphesModule
@@ -51 +34 @@
   Double_t fParameterR;
 
+  Double_t fMeanSqDeltaRMaxBarrel; // |eta| < 1.5
+  Double_t fBetaMinBarrel; // |eta| < 2.5
+  Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < |eta| < 4.0
+  Double_t fBetaMinEndcap; // 1.5 < |eta| < 4.0
+  Double_t fMeanSqDeltaRMaxForward; // |eta| > 4.0
+
+  Double_t fNeutralPTMin;
+  Double_t fJetPTMinForNeutrals;
+
+  /*
+JAY
+---
+
+|Eta|<1.5
+
+meanSqDeltaR betaStar SigEff BgdEff
+0.13 0.92 96% 8%
+0.13 0.95 97% 16%
+0.13 0.97 98% 27%
+
+|Eta|>1.5
+
+meanSqDeltaR betaStar SigEff BgdEff
+0.14 0.91 95% 15%
+0.14 0.94 97% 19%
+0.14 0.97 98% 29%
+
+BRYAN
+-----
+
+Barrel (MeanSqDR, Beta, sig eff, bg eff): 
+0.10, 0.08, 90%, 8%
+0.11, 0.12, 90%, 6%
+0.13, 0.16, 89%, 5%
+
+Endcap (MeanSqDR, Beta, sig eff, bg eff):
+0.07, 0.06, 89%, 4%
+0.08, 0.08, 92%, 6%
+0.09, 0.08, 95%, 10%
+0.10, 0.08, 97%, 13%
+
+SETH GUESSES FOR |eta| > 4.0
+----------------------------
+
+MeanSqDeltaR 
+0.07 
+0.10 
+0.14 
+0.2
+  */
+
   // If set to true, may have weird results for PFCHS
   // If set to false, uses everything within dR < fParameterR even if in other jets &c.
   // Results should be very similar for PF
-  Int_t fUseConstituents;
+  Int_t fUseConstituents; 
 
   Bool_t fAverageEachTower;
@@ -62 +96 @@
   const TObjArray *fJetInputArray; //!
 
-  const TObjArray *fTrackInputArray; //!
-  const TObjArray *fNeutralInputArray; //!
+  const TObjArray *fTrackInputArray; // SCZ
+  const TObjArray *fNeutralInputArray; 
 
-  TIterator *fItTrackInputArray; //!
-  TIterator *fItNeutralInputArray; //!
+  TIterator *fItTrackInputArray; // SCZ
+  TIterator *fItNeutralInputArray; // SCZ
 
   TObjArray *fOutputArray; //!
+  TObjArray *fNeutralsInPassingJets; // SCZ
 
-  TIterator *fItVertexInputArray; //!
-  const TObjArray *fVertexInputArray; //!
 
-  Double_t fZVertexResolution;
-
-  ClassDef(PileUpJetID, 1)
+  ClassDef(PileUpJetID, 2)
 };
 
 #endif
-

modules/PileUpMerger.cc

@@ -80 +80 @@
   fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
 
+  fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
+  fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
+  fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
+  fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
+
   // read vertex smearing formula
 
@@ -111 +116 @@
   TParticlePDG *pdgParticle;
   Int_t pid;
-  Float_t x, y, z, t;
+  Float_t x, y, z, t, vx, vy;
   Float_t px, py, pz, e;
   Double_t dz, dphi, dt;
-  Int_t numberOfEvents, event;
+  Int_t numberOfEvents, event, numberOfParticles;
   Long64_t allEntries, entry;
-  Candidate *candidate, *vertexcandidate;
+  Candidate *candidate, *vertex;
   DelphesFactory *factory;
 
@@ -123 +128 @@
   fItInputArray->Reset();
 
-  // --- Deal with Primary vertex first  ------
+  // --- Deal with primary vertex first  ------
 
   fFunction->GetRandom2(dz, dt);
@@ -129 +134 @@
   dt *= c_light*1.0E3; // necessary in order to make t in mm/c
   dz *= 1.0E3; // necessary in order to make z in mm
-
+  vx = 0.0;
+  vy = 0.0;
+  numberOfParticles = fInputArray->GetEntriesFast();
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
+    vx += candidate->Position.X();
+    vy += candidate->Position.Y();
     z = candidate->Position.Z();
     t = candidate->Position.T();
@@ -139 +148 @@
   }
 
+  if(numberOfParticles > 0)
+  {
+    vx /= numberOfParticles;
+    vy /= numberOfParticles;
+  }
+
   factory = GetFactory();
 
-  vertexcandidate = factory->NewCandidate();
-  vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
-  fVertexOutputArray->Add(vertexcandidate);
+  vertex = factory->NewCandidate();
+  vertex->Position.SetXYZT(vx, vy, dz, dt);
+  fVertexOutputArray->Add(vertex);
 
   // --- Then with pile-up vertices  ------
@@ -181 +196 @@
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    vertexcandidate = factory->NewCandidate();
-    vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
-    vertexcandidate->IsPU = 1;
-
-    fVertexOutputArray->Add(vertexcandidate);
-
+    vx = 0.0;
+    vy = 0.0;
+    numberOfParticles = 0;
     while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
     {
@@ -204 +216 @@
       candidate->Momentum.RotateZ(dphi);
 
+      x -= fInputBeamSpotX;
+      y -= fInputBeamSpotY;
       candidate->Position.SetXYZT(x, y, z + dz, t + dt);
       candidate->Position.RotateZ(dphi);
+      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
+
+      vx += candidate->Position.X();
+      vy += candidate->Position.Y();
+      ++numberOfParticles;
 
       fParticleOutputArray->Add(candidate);
     }
-  }
-}
-
-//------------------------------------------------------------------------------
-
+
+    if(numberOfParticles > 0)
+    {
+      vx /= numberOfParticles;
+      vy /= numberOfParticles;
+    }
+
+    vertex = factory->NewCandidate();
+    vertex->Position.SetXYZT(vx, vy, dz, dt);
+    vertex->IsPU = 1;
+
+    fVertexOutputArray->Add(vertex);
+  }
+}
+
+//------------------------------------------------------------------------------

modules/PileUpMerger.h

@@ -53 +53 @@
   Double_t fTVertexSpread;
 
+  Double_t fInputBeamSpotX;
+  Double_t fInputBeamSpotY;
+  Double_t fOutputBeamSpotX;
+  Double_t fOutputBeamSpotY;
+
   DelphesTF2 *fFunction; //!
 

modules/PileUpMergerPythia8.cc

@@ -21 +21 @@
  *  Merges particles from pile-up sample into event
  *
- *  \author P. Selvaggi - UCL, Louvain-la-Neuve
+ *  \author M. Selvaggi - UCL, Louvain-la-Neuve
  *
  */
@@ -29 +29 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesFormula.h"
+#include "classes/DelphesTF2.h"
 #include "classes/DelphesPileUpReader.h"
 
@@ -56 +56 @@
 
 PileUpMergerPythia8::PileUpMergerPythia8() :
-  fPythia(0), fItInputArray(0)
-{
+  fFunction(0), fPythia(0), fItInputArray(0)
+{
+  fFunction = new DelphesTF2;
 }
 
@@ -64 +65 @@
 PileUpMergerPythia8::~PileUpMergerPythia8()
 {
+  delete fFunction;
 }
 
@@ -72 +74 @@
   const char *fileName;
 
+  fPileUpDistribution = GetInt("PileUpDistribution", 0);
+
   fMeanPileUp  = GetDouble("MeanPileUp", 10);
-  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
+
+  fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
+  fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
+
+  fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
+  fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
+  fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
+  fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
 
   fPTMin = GetDouble("PTMin", 0.0);
+
+  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
+  fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
 
   fileName = GetString("ConfigFile", "MinBias.cmnd");
@@ -86 +100 @@
 
   // create output arrays
-  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
+  fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
+  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
 }
 
@@ -102 +117 @@
   TDatabasePDG *pdg = TDatabasePDG::Instance();
   TParticlePDG *pdgParticle;
-  Int_t pid;
-  Float_t x, y, z, t;
+  Int_t pid, status;
+  Float_t x, y, z, t, vx, vy;
   Float_t px, py, pz, e;
-  Double_t dz, dphi;
-  Int_t poisson, event, i;
-  Candidate *candidate;
+  Double_t dz, dphi, dt;
+  Int_t numberOfEvents, event, numberOfParticles, i;
+  Candidate *candidate, *vertex;
   DelphesFactory *factory;
 
+  const Double_t c_light = 2.99792458E8;
+
   fItInputArray->Reset();
+
+  // --- Deal with primary vertex first  ------
+
+  fFunction->GetRandom2(dz, dt);
+
+  dt *= c_light*1.0E3; // necessary in order to make t in mm/c
+  dz *= 1.0E3; // necessary in order to make z in mm
+  vx = 0.0;
+  vy = 0.0;
+  numberOfParticles = fInputArray->GetEntriesFast();
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
   {
-    fOutputArray->Add(candidate);
+    vx += candidate->Position.X();
+    vy += candidate->Position.Y();
+    z = candidate->Position.Z();
+    t = candidate->Position.T();
+    candidate->Position.SetZ(z + dz);
+    candidate->Position.SetT(t + dt);
+    fParticleOutputArray->Add(candidate);
+  }
+
+  if(numberOfParticles > 0)
+  {
+    vx /= numberOfParticles;
+    vy /= numberOfParticles;
   }
 
   factory = GetFactory();
 
-  poisson = gRandom->Poisson(fMeanPileUp);
-
-  for(event = 0; event < poisson; ++event)
+  vertex = factory->NewCandidate();
+  vertex->Position.SetXYZT(vx, vy, dz, dt);
+  fVertexOutputArray->Add(vertex);
+
+  // --- Then with pile-up vertices  ------
+
+  switch(fPileUpDistribution)
+  {
+    case 0:
+      numberOfEvents = gRandom->Poisson(fMeanPileUp);
+      break;
+    case 1:
+      numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
+      break;
+    default:
+      numberOfEvents = gRandom->Poisson(fMeanPileUp);
+      break;
+  }
+
+  for(event = 0; event < numberOfEvents; ++event)
   {
     while(!fPythia->next());
 
-    dz = gRandom->Gaus(0.0, fZVertexSpread);
+   // --- Pile-up vertex smearing
+
+    fFunction->GetRandom2(dz, dt);
+
+    dt *= c_light*1.0E3; // necessary in order to make t in mm/c
+    dz *= 1.0E3; // necessary in order to make z in mm
+
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    for(i = 1; i < fPythia->event.size(); ++i)
+    vx = 0.0;
+    vy = 0.0;
+    numberOfParticles = fPythia->event.size();
+    for(i = 1; i < numberOfParticles; ++i)
     {
       Pythia8::Particle &particle = fPythia->event[i];
 
-      if(particle.status() != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
+      status = particle.statusHepMC();
+
+      if(status != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
 
       pid = particle.id();
@@ -152 +219 @@
       candidate->Momentum.RotateZ(dphi);
 
-      candidate->Position.SetXYZT(x, y, z + dz, t);
+      x -= fInputBeamSpotX;
+      y -= fInputBeamSpotY;
+      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
       candidate->Position.RotateZ(dphi);
-
-      fOutputArray->Add(candidate);
+      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
+
+      vx += candidate->Position.X();
+      vy += candidate->Position.Y();
+
+      fParticleOutputArray->Add(candidate);
     }
-  }
-}
-
-//------------------------------------------------------------------------------
-
+
+    if(numberOfParticles > 0)
+    {
+      vx /= numberOfParticles;
+      vy /= numberOfParticles;
+    }
+
+    vertex = factory->NewCandidate();
+    vertex->Position.SetXYZT(vx, vy, dz, dt);
+    vertex->IsPU = 1;
+
+    fVertexOutputArray->Add(vertex);
+  }
+}
+
+//------------------------------------------------------------------------------

modules/PileUpMergerPythia8.h

@@ -31 +31 @@
 
 class TObjArray;
+class DelphesTF2;
 
 namespace Pythia8
@@ -50 +51 @@
 private:
 
+  Int_t fPileUpDistribution;
   Double_t fMeanPileUp;
+
   Double_t fZVertexSpread;
+  Double_t fTVertexSpread;
+
+  Double_t fInputBeamSpotX;
+  Double_t fInputBeamSpotY;
+  Double_t fOutputBeamSpotX;
+  Double_t fOutputBeamSpotY;
+
   Double_t fPTMin;
+
+  DelphesTF2 *fFunction; //!
 
   Pythia8::Pythia *fPythia; //!
@@ -60 +72 @@
   const TObjArray *fInputArray; //!
 
-  TObjArray *fOutputArray; //!
+  TObjArray *fParticleOutputArray; //!
+  TObjArray *fVertexOutputArray; //!
 
   ClassDef(PileUpMergerPythia8, 1)

modules/SimpleCalorimeter.cc

@@ -150 +150 @@
   fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
 
+  // flag that says if current calo is Ecal of Hcal (will then fill correct values of Eem and Ehad)
+  fIsEcal = GetBool("IsEcal", false);
+
   // switch on or off the dithering of the center of calorimeter towers
-  fDitherTowerCenter = GetBool("DitherTowerCenter", true);
+  fSmearTowerCenter = GetBool("SmearTowerCenter", true);
 
   // read resolution formulas
@@ -409 +412 @@
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
 
-  if(fDitherTowerCenter)
+  if(fSmearTowerCenter)
   {
     eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
@@ -424 +427 @@
   fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+
+  fTower->Eem = (!fIsEcal) ? 0 : energy;
+  fTower->Ehad = (fIsEcal) ? 0 : energy; 
 
   fTower->Edges[0] = fTowerEdges[0];
@@ -447 +453 @@
     pt = energy / TMath::CosH(eta);
 
+    tower->Eem = (!fIsEcal) ? 0 : energy;
+    tower->Ehad = (fIsEcal) ? 0 : energy; 
+
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     fEFlowTowerOutputArray->Add(tower);

modules/SimpleCalorimeter.h

@@ +1 @@
+
 /*
  *  Delphes: a framework for fast simulation of a generic collider experiment
@@ -72 +73 @@
   Double_t fEnergySignificanceMin;
 
-  Bool_t fDitherTowerCenter;
+  Bool_t fSmearTowerCenter;
+
+  Bool_t fIsEcal; //! 
 
   TFractionMap fFractionMap; //!

modules/TauTagging.cc

@@ -33 +33 @@
 #include "classes/DelphesFactory.h"
 #include "classes/DelphesFormula.h"
-
-#include "ExRootAnalysis/ExRootResult.h"
-#include "ExRootAnalysis/ExRootFilter.h"
-#include "ExRootAnalysis/ExRootClassifier.h"
 
 #include "TMath.h"
@@ -53 +49 @@
 using namespace std;
 
-//------------------------------------------------------------------------------
-
-class TauTaggingPartonClassifier : public ExRootClassifier
-{
-public:
-
-  TauTaggingPartonClassifier(const TObjArray *array);
-
-  Int_t GetCategory(TObject *object);
-
-  Double_t fEtaMax, fPTMin;
-
-  const TObjArray *fParticleInputArray;
-};
 
 //------------------------------------------------------------------------------
@@ -79 +61 @@
 {
   Candidate *tau = static_cast<Candidate *>(object);
-  Candidate *daughter = 0;
+  Candidate *daughter1 = 0;
+  Candidate *daughter2 = 0;
+  
   const TLorentzVector &momentum = tau->Momentum;
-  Int_t pdgCode, i;
+  Int_t pdgCode, i, j;
 
   pdgCode = TMath::Abs(tau->PID);
@@ -90 +74 @@
   if(tau->D1 < 0) return -1;
 
+  if(tau->D2 < tau->D1) return -1;
+
   if(tau->D1 >= fParticleInputArray->GetEntriesFast() ||
      tau->D2 >= fParticleInputArray->GetEntriesFast())
@@ -98 +84 @@
   for(i = tau->D1; i <= tau->D2; ++i)
   {
-    daughter = static_cast<Candidate *>(fParticleInputArray->At(i));
-    pdgCode = TMath::Abs(daughter->PID);
-    if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15 || pdgCode == 24) return -1;
+    daughter1 = static_cast<Candidate *>(fParticleInputArray->At(i));
+    pdgCode = TMath::Abs(daughter1->PID);
+    if(pdgCode == 11 || pdgCode == 13 || pdgCode == 15) return -1;
+    else if(pdgCode == 24)
+    {
+     if(daughter1->D1 < 0) return -1;
+     for(j = daughter1->D1; j <= daughter1->D2; ++j)
+     {
+       daughter2 = static_cast<Candidate*>(fParticleInputArray->At(j));
+       pdgCode = TMath::Abs(daughter2->PID);
+       if(pdgCode == 11 || pdgCode == 13) return -1;
+     }
+        
+    }
   }
 
@@ -206 +203 @@
   tauArray = fFilter->GetSubArray(fClassifier, 0);
 
-  if(tauArray == 0) return;
-
-  TIter itTauArray(tauArray);
-
   // loop over all input jets
   fItJetInputArray->Reset();
@@ -222 +215 @@
 
     // loop over all input taus
-    itTauArray.Reset();
-    while((tau = static_cast<Candidate *>(itTauArray.Next())))
-    {
-      if(tau->D1 < 0) continue;
-
-      if(tau->D1 >= fParticleInputArray->GetEntriesFast() ||
-         tau->D2 >= fParticleInputArray->GetEntriesFast())
+    if(tauArray){
+      TIter itTauArray(tauArray);
+      while((tau = static_cast<Candidate *>(itTauArray.Next())))
       {
-        throw runtime_error("tau's daughter index is greater than the ParticleInputArray size");
-      }
-
-      tauMomentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
-
-      for(i = tau->D1; i <= tau->D2; ++i)
-      {
-        daughter = static_cast<Candidate *>(fParticleInputArray->At(i));
-        if(TMath::Abs(daughter->PID) == 16) continue;
-        tauMomentum += daughter->Momentum;
-      }
-
-      if(jetMomentum.DeltaR(tauMomentum) <= fDeltaR)
-      {
-        pdgCode = 15;
-        charge = tau->Charge;
+        if(tau->D1 < 0) continue;
+
+        if(tau->D1 >= fParticleInputArray->GetEntriesFast() ||
+           tau->D2 >= fParticleInputArray->GetEntriesFast())
+        {
+          throw runtime_error("tau's daughter index is greater than the ParticleInputArray size");
+        }
+
+        tauMomentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
+
+        for(i = tau->D1; i <= tau->D2; ++i)
+        {
+          daughter = static_cast<Candidate *>(fParticleInputArray->At(i));
+          if(TMath::Abs(daughter->PID) == 16) continue;
+          tauMomentum += daughter->Momentum;
+        }
+
+        if(jetMomentum.DeltaR(tauMomentum) <= fDeltaR)
+        {
+          pdgCode = 15;
+          charge = tau->Charge;
+        }
       }
     }

modules/TauTagging.h

@@ -31 +31 @@
 
 #include "classes/DelphesModule.h"
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
 
 #include <map>
@@ -76 +79 @@
 };
 
+
+//------------------------------------------------------------------------------
+
+class TauTaggingPartonClassifier : public ExRootClassifier
+{
+public:
+
+  TauTaggingPartonClassifier(const TObjArray *array);
+
+  Int_t GetCategory(TObject *object);
+
+  Double_t fEtaMax, fPTMin;
+
+  const TObjArray *fParticleInputArray;
+};
+
+
 #endif

modules/TrackPileUpSubtractor.cc

@@ -74 +74 @@
   fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
 
+  fPTMin = GetDouble("PTMin", 0.);
   // import arrays with output from other modules
-
+   
   ExRootConfParam param = GetParam("InputArray");
   Long_t i, size;
@@ -147 +148 @@
       // assume perfect pile-up subtraction for tracks outside fZVertexResolution
       
-      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue;
-
-      array->Add(candidate);
+      if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) candidate->IsRecoPU = 1;
+      else 
+      {
+         candidate->IsRecoPU = 0;
+         if( candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
+      }
     }
   }

modules/TrackPileUpSubtractor.h

@@ -50 +50 @@
   Double_t fZVertexResolution;
 
+  Double_t fPTMin; 
+
   std::map< TIterator *, TObjArray * > fInputMap; //!
 

modules/TreeWriter.cc

@@ -291 +291 @@
     entry->ZOuter = position.Z();
     entry->TOuter = position.T()*1.0E-3/c_light;
-    
+
     entry->Dxy = candidate->Dxy;
     entry->SDxy = candidate->SDxy ;
@@ -297 +297 @@
     entry->Yd = candidate->Yd;
     entry->Zd = candidate->Zd;
-   
+
     const TLorentzVector &momentum = candidate->Momentum;
 
@@ -362 +362 @@
 
     entry->T = position.T()*1.0E-3/c_light;
+    entry->NTimeHits = candidate->NTimeHits;
 
     FillParticles(candidate, &entry->Particles);
@@ -403 +404 @@
     entry->T = position.T()*1.0E-3/c_light;
 
+    // Isolation variables
+
+    entry->IsolationVar = candidate->IsolationVar;
+    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
+    entry->SumPtCharged = candidate->SumPtCharged ;
+    entry->SumPtNeutral = candidate->SumPtNeutral ;
+    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
+    entry->SumPt = candidate->SumPt ;
+
     entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9;
 
@@ -442 +452 @@
     entry->T = position.T()*1.0E-3/c_light;
 
+    // Isolation variables
+
+    entry->IsolationVar = candidate->IsolationVar;
+    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
+    entry->SumPtCharged = candidate->SumPtCharged ;
+    entry->SumPtNeutral = candidate->SumPtNeutral ;
+    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
+    entry->SumPt = candidate->SumPt ;
+
+
     entry->Charge = candidate->Charge;
 
@@ -469 +489 @@
     const TLorentzVector &momentum = candidate->Momentum;
     const TLorentzVector &position = candidate->Position;
-
 
     pt = momentum.Pt();
@@ -488 +507 @@
     entry->T = position.T()*1.0E-3/c_light;
 
+    // Isolation variables
+
+    entry->IsolationVar = candidate->IsolationVar;
+    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr ;
+    entry->SumPtCharged = candidate->SumPtCharged ;
+    entry->SumPtNeutral = candidate->SumPtNeutral ;
+    entry->SumPtChargedPU = candidate->SumPtChargedPU ;
+    entry->SumPt = candidate->SumPt ;
+
     entry->Charge = candidate->Charge;
 
@@ -504 +532 @@
   Double_t ecalEnergy, hcalEnergy;
   const Double_t c_light = 2.99792458E8;
+  Int_t i;
 
   array->Sort();
@@ -532 +561 @@
     entry->Mass = momentum.M();
 
+    entry->Area = candidate->Area;
+
     entry->DeltaEta = candidate->DeltaEta;
     entry->DeltaPhi = candidate->DeltaPhi;
 
+    entry->Flavor = candidate->Flavor;
+    entry->FlavorAlgo = candidate->FlavorAlgo;
+    entry->FlavorPhys = candidate->FlavorPhys;
+
     entry->BTag = candidate->BTag;
+
+    entry->BTagAlgo = candidate->BTagAlgo;
+    entry->BTagPhys = candidate->BTagPhys;
+
     entry->TauTag = candidate->TauTag;
 
@@ -561 +600 @@
     entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
     entry->PTD = candidate->PTD;
-    entry->FracPt[0] = candidate->FracPt[0];
-    entry->FracPt[1] = candidate->FracPt[1];
-    entry->FracPt[2] = candidate->FracPt[2];
-    entry->FracPt[3] = candidate->FracPt[3];
-    entry->FracPt[4] = candidate->FracPt[4];
-   
-    //--- N-subjettiness variables ----
-    
-    entry->Tau1 = candidate->Tau[0];
-    entry->Tau2 = candidate->Tau[1];
-    entry->Tau3 = candidate->Tau[2];
-    entry->Tau4 = candidate->Tau[3];
-    entry->Tau5 = candidate->Tau[4];
-   
+
+    //--- Sub-structure variables ----
+
+    entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed;
+    entry->NSubJetsPruned = candidate->NSubJetsPruned;
+    entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped;
+
+    for(i = 0; i < 5; i++)
+    {
+      entry->FracPt[i] = candidate -> FracPt[i];
+      entry->Tau[i] = candidate -> Tau[i];
+      entry->TrimmedP4[i] = candidate -> TrimmedP4[i];
+      entry->PrunedP4[i] = candidate -> PrunedP4[i];
+      entry->SoftDroppedP4[i] = candidate -> SoftDroppedP4[i];
+    }
+
     FillParticles(candidate, &entry->Particles);
   }

modules/UniqueObjectFinder.cc

@@ -74 +74 @@
   TIterator *iterator;
 
+  fInputMap.clear();
+
   size = param.GetSize();
   for(i = 0; i < size/2; ++i)
@@ -80 +82 @@
     iterator = array->MakeIterator();
 
-    fInputMap[iterator] = ExportArray(param[i*2 + 1].GetString());
+    fInputMap.push_back(make_pair(iterator, ExportArray(param[i*2 + 1].GetString())));
   }
 }
@@ -88 +90 @@
 void UniqueObjectFinder::Finish()
 {
-  map< TIterator *, TObjArray * >::iterator itInputMap;
+  vector< pair< TIterator *, TObjArray * > >::iterator itInputMap;
   TIterator *iterator;
 
@@ -104 +106 @@
 {
   Candidate *candidate;
-  map< TIterator *, TObjArray * >::iterator itInputMap;
+  vector< pair< TIterator *, TObjArray * > >::iterator itInputMap;
   TIterator *iterator;
   TObjArray *array;
@@ -128 +130 @@
 //------------------------------------------------------------------------------
 
-Bool_t UniqueObjectFinder::Unique(Candidate *candidate, map< TIterator *, TObjArray * >::iterator itInputMap)
+Bool_t UniqueObjectFinder::Unique(Candidate *candidate, vector< pair< TIterator *, TObjArray * > >::iterator itInputMap)
 {
   Candidate *previousCandidate;
-  map< TIterator *, TObjArray * >::iterator previousItInputMap;
+  vector< pair< TIterator *, TObjArray * > >::iterator previousItInputMap;
   TObjArray *array;
 

modules/UniqueObjectFinder.h

@@ -30 +30 @@
 #include "classes/DelphesModule.h"
 
-#include <map>
+#include <vector>
+#include <utility>
 
 class TIterator;
@@ -49 +50 @@
 private:
 
-  Bool_t Unique(Candidate *candidate, std::map< TIterator *, TObjArray * >::iterator itInputMap);
+  Bool_t Unique(Candidate *candidate, std::vector< std::pair< TIterator *, TObjArray * > >::iterator itInputMap);
 
-  std::map< TIterator *, TObjArray * > fInputMap; //!
+  std::vector< std::pair< TIterator *, TObjArray * > > fInputMap; //!
 
   ClassDef(UniqueObjectFinder, 1)