Changes in modules/PileUpMergerPythia8.cc [2d494a6:cab38f6] in git

File:

• modules/PileUpMergerPythia8.cc (modified) (8 diffs)

modules/PileUpMergerPythia8.cc

@@ -21 +21 @@
  *  Merges particles from pile-up sample into event
  *
- *  \author M. Selvaggi - UCL, Louvain-la-Neuve
+ *  \author P. Selvaggi - UCL, Louvain-la-Neuve
  *
  */
@@ -29 +29 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesTF2.h"
+#include "classes/DelphesFormula.h"
 #include "classes/DelphesPileUpReader.h"
 
@@ -56 +56 @@
 
 PileUpMergerPythia8::PileUpMergerPythia8() :
-  fFunction(0), fPythia(0), fItInputArray(0)
+  fPythia(0), fItInputArray(0)
 {
-  fFunction = new DelphesTF2;
 }
 
@@ -65 +64 @@
 PileUpMergerPythia8::~PileUpMergerPythia8()
 {
-  delete fFunction;
 }
 
@@ -74 +72 @@
   const char *fileName;
 
-  fPileUpDistribution = GetInt("PileUpDistribution", 0);
-
   fMeanPileUp  = GetDouble("MeanPileUp", 10);
-
-  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);
+  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
 
   fPTMin = GetDouble("PTMin", 0.0);
-
-  fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
-  fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
 
   fileName = GetString("ConfigFile", "MinBias.cmnd");
@@ -100 +86 @@
 
   // create output arrays
-  fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
-  fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
+  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
 }
 
@@ -117 +102 @@
   TDatabasePDG *pdg = TDatabasePDG::Instance();
   TParticlePDG *pdgParticle;
-  Int_t pid, status;
-  Float_t x, y, z, t, vx, vy;
+  Int_t pid;
+  Float_t x, y, z, t;
   Float_t px, py, pz, e;
-  Double_t dz, dphi, dt;
-  Int_t numberOfEvents, event, numberOfParticles, i;
-  Candidate *candidate, *vertex;
+  Double_t dz, dphi;
+  Int_t poisson, event, i;
+  Candidate *candidate;
   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())))
   {
-    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;
+    fOutputArray->Add(candidate);
   }
 
   factory = GetFactory();
 
-  vertex = factory->NewCandidate();
-  vertex->Position.SetXYZT(vx, vy, dz, dt);
-  fVertexOutputArray->Add(vertex);
+  poisson = gRandom->Poisson(fMeanPileUp);
 
-  // --- 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)
+  for(event = 0; event < poisson; ++event)
   {
     while(!fPythia->next());
 
-   // --- 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
-
+    dz = gRandom->Gaus(0.0, fZVertexSpread);
     dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
 
-    vx = 0.0;
-    vy = 0.0;
-    numberOfParticles = fPythia->event.size();
-    for(i = 1; i < numberOfParticles; ++i)
+    for(i = 1; i < fPythia->event.size(); ++i)
     {
       Pythia8::Particle &particle = fPythia->event[i];
 
-      status = particle.statusHepMC();
-
-      if(status != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
+      if(particle.status() != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
 
       pid = particle.id();
@@ -219 +152 @@
       candidate->Momentum.RotateZ(dphi);
 
-      x -= fInputBeamSpotX;
-      y -= fInputBeamSpotY;
-      candidate->Position.SetXYZT(x, y, z + dz, t + dt);
+      candidate->Position.SetXYZT(x, y, z + dz, t);
       candidate->Position.RotateZ(dphi);
-      candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
 
-      vx += candidate->Position.X();
-      vy += candidate->Position.Y();
-
-      fParticleOutputArray->Add(candidate);
+      fOutputArray->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);
   }
 }
 
 //------------------------------------------------------------------------------
+