Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/PileUpMergerPythia8.cc

    r2d494a6 rdacd9c5  
    2929#include "classes/DelphesClasses.h"
    3030#include "classes/DelphesFactory.h"
    31 #include "classes/DelphesTF2.h"
     31#include "classes/DelphesFormula.h"
    3232#include "classes/DelphesPileUpReader.h"
    3333
     
    5656
    5757PileUpMergerPythia8::PileUpMergerPythia8() :
    58   fFunction(0), fPythia(0), fItInputArray(0)
     58  fPythia(0), fItInputArray(0)
    5959{
    60   fFunction = new DelphesTF2;
    6160}
    6261
     
    6564PileUpMergerPythia8::~PileUpMergerPythia8()
    6665{
    67   delete fFunction;
    6866}
    6967
     
    7472  const char *fileName;
    7573
    76   fPileUpDistribution = GetInt("PileUpDistribution", 0);
    77 
    7874  fMeanPileUp  = GetDouble("MeanPileUp", 10);
    79 
    80   fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
    81   fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
    82 
    83   fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0);
    84   fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0);
    85   fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0);
    86   fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0);
     75  fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
    8776
    8877  fPTMin = GetDouble("PTMin", 0.0);
    89 
    90   fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
    91   fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
    9278
    9379  fileName = GetString("ConfigFile", "MinBias.cmnd");
     
    10086
    10187  // create output arrays
    102   fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
    103   fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
     88  fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
    10489}
    10590
     
    118103  TParticlePDG *pdgParticle;
    119104  Int_t pid, status;
    120   Float_t x, y, z, t, vx, vy;
     105  Float_t x, y, z, t;
    121106  Float_t px, py, pz, e;
    122   Double_t dz, dphi, dt;
    123   Int_t numberOfEvents, event, numberOfParticles, i;
    124   Candidate *candidate, *vertex;
     107  Double_t dz, dphi;
     108  Int_t poisson, event, i;
     109  Candidate *candidate;
    125110  DelphesFactory *factory;
    126111
    127   const Double_t c_light = 2.99792458E8;
    128 
    129112  fItInputArray->Reset();
    130 
    131   // --- Deal with primary vertex first  ------
    132 
    133   fFunction->GetRandom2(dz, dt);
    134 
    135   dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    136   dz *= 1.0E3; // necessary in order to make z in mm
    137   vx = 0.0;
    138   vy = 0.0;
    139   numberOfParticles = fInputArray->GetEntriesFast();
    140113  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    141114  {
    142     vx += candidate->Position.X();
    143     vy += candidate->Position.Y();
    144     z = candidate->Position.Z();
    145     t = candidate->Position.T();
    146     candidate->Position.SetZ(z + dz);
    147     candidate->Position.SetT(t + dt);
    148     fParticleOutputArray->Add(candidate);
    149   }
    150 
    151   if(numberOfParticles > 0)
    152   {
    153     vx /= numberOfParticles;
    154     vy /= numberOfParticles;
     115    fOutputArray->Add(candidate);
    155116  }
    156117
    157118  factory = GetFactory();
    158119
    159   vertex = factory->NewCandidate();
    160   vertex->Position.SetXYZT(vx, vy, dz, dt);
    161   fVertexOutputArray->Add(vertex);
     120  poisson = gRandom->Poisson(fMeanPileUp);
    162121
    163   // --- Then with pile-up vertices  ------
    164 
    165   switch(fPileUpDistribution)
    166   {
    167     case 0:
    168       numberOfEvents = gRandom->Poisson(fMeanPileUp);
    169       break;
    170     case 1:
    171       numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
    172       break;
    173     default:
    174       numberOfEvents = gRandom->Poisson(fMeanPileUp);
    175       break;
    176   }
    177 
    178   for(event = 0; event < numberOfEvents; ++event)
     122  for(event = 0; event < poisson; ++event)
    179123  {
    180124    while(!fPythia->next());
    181125
    182    // --- Pile-up vertex smearing
    183 
    184     fFunction->GetRandom2(dz, dt);
    185 
    186     dt *= c_light*1.0E3; // necessary in order to make t in mm/c
    187     dz *= 1.0E3; // necessary in order to make z in mm
    188 
     126    dz = gRandom->Gaus(0.0, fZVertexSpread);
    189127    dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    190128
    191     vx = 0.0;
    192     vy = 0.0;
    193     numberOfParticles = fPythia->event.size();
    194     for(i = 1; i < numberOfParticles; ++i)
     129    for(i = 1; i < fPythia->event.size(); ++i)
    195130    {
    196131      Pythia8::Particle &particle = fPythia->event[i];
     
    208143      candidate->PID = pid;
    209144
    210       candidate->Status = 1;
     145      candidate->Status = status;
    211146
    212147      pdgParticle = pdg->GetParticle(pid);
     
    219154      candidate->Momentum.RotateZ(dphi);
    220155
    221       x -= fInputBeamSpotX;
    222       y -= fInputBeamSpotY;
    223       candidate->Position.SetXYZT(x, y, z + dz, t + dt);
     156      candidate->Position.SetXYZT(x, y, z + dz, t);
    224157      candidate->Position.RotateZ(dphi);
    225       candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
    226158
    227       vx += candidate->Position.X();
    228       vy += candidate->Position.Y();
    229 
    230       fParticleOutputArray->Add(candidate);
     159      fOutputArray->Add(candidate);
    231160    }
    232 
    233     if(numberOfParticles > 0)
    234     {
    235       vx /= numberOfParticles;
    236       vy /= numberOfParticles;
    237     }
    238 
    239     vertex = factory->NewCandidate();
    240     vertex->Position.SetXYZT(vx, vy, dz, dt);
    241     vertex->IsPU = 1;
    242 
    243     fVertexOutputArray->Add(vertex);
    244161  }
    245162}
    246163
    247164//------------------------------------------------------------------------------
     165
Note: See TracChangeset for help on using the changeset viewer.