Index: modules/PileUpMerger.cc
===================================================================
--- modules/PileUpMerger.cc	(revision b62c2da882823815cc5aab13f1901f5277f125a0)
+++ modules/PileUpMerger.cc	(revision 2d494a6c69b4a113bb7bbef0cb3ba494fd2aa400)
@@ -80,4 +80,9 @@
   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,10 +116,10 @@
   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,5 +128,5 @@
   fItInputArray->Reset();
 
-  // --- Deal with Primary vertex first  ------
+  // --- Deal with primary vertex first  ------
 
   fFunction->GetRandom2(dz, dt);
@@ -129,7 +134,11 @@
   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,9 +148,15 @@
   }
 
+  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,10 +196,7 @@
     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,12 +216,30 @@
       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);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/PileUpMerger.h
===================================================================
--- modules/PileUpMerger.h	(revision b62c2da882823815cc5aab13f1901f5277f125a0)
+++ modules/PileUpMerger.h	(revision 2d494a6c69b4a113bb7bbef0cb3ba494fd2aa400)
@@ -53,4 +53,9 @@
   Double_t fTVertexSpread;
 
+  Double_t fInputBeamSpotX;
+  Double_t fInputBeamSpotY;
+  Double_t fOutputBeamSpotX;
+  Double_t fOutputBeamSpotY;
+
   DelphesTF2 *fFunction; //!
 
Index: modules/PileUpMergerPythia8.cc
===================================================================
--- modules/PileUpMergerPythia8.cc	(revision b62c2da882823815cc5aab13f1901f5277f125a0)
+++ modules/PileUpMergerPythia8.cc	(revision 2d494a6c69b4a113bb7bbef0cb3ba494fd2aa400)
@@ -29,5 +29,5 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesFormula.h"
+#include "classes/DelphesTF2.h"
 #include "classes/DelphesPileUpReader.h"
 
@@ -56,6 +56,7 @@
 
 PileUpMergerPythia8::PileUpMergerPythia8() :
-  fPythia(0), fItInputArray(0)
-{
+  fFunction(0), fPythia(0), fItInputArray(0)
+{
+  fFunction = new DelphesTF2;
 }
 
@@ -64,4 +65,5 @@
 PileUpMergerPythia8::~PileUpMergerPythia8()
 {
+  delete fFunction;
 }
 
@@ -72,8 +74,20 @@
   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,5 +100,6 @@
 
   // create output arrays
-  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
+  fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
+  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
 }
 
@@ -103,29 +118,79 @@
   TParticlePDG *pdgParticle;
   Int_t pid, status;
-  Float_t x, y, z, t;
+  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];
@@ -143,5 +208,5 @@
       candidate->PID = pid;
 
-      candidate->Status = status;
+      candidate->Status = 1;
 
       pdgParticle = pdg->GetParticle(pid);
@@ -154,12 +219,29 @@
       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);
+  }
+}
+
+//------------------------------------------------------------------------------
Index: modules/PileUpMergerPythia8.h
===================================================================
--- modules/PileUpMergerPythia8.h	(revision b62c2da882823815cc5aab13f1901f5277f125a0)
+++ modules/PileUpMergerPythia8.h	(revision 2d494a6c69b4a113bb7bbef0cb3ba494fd2aa400)
@@ -31,4 +31,5 @@
 
 class TObjArray;
+class DelphesTF2;
 
 namespace Pythia8
@@ -50,7 +51,18 @@
 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,5 +72,6 @@
   const TObjArray *fInputArray; //!
 
-  TObjArray *fOutputArray; //!
+  TObjArray *fParticleOutputArray; //!
+  TObjArray *fVertexOutputArray; //!
 
   ClassDef(PileUpMergerPythia8, 1)
