Fork me on GitHub

source: git/modules/PileUpMergerPythia8.cc@ 77249aa

Timing
Last change on this file since 77249aa was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

  • Property mode set to 100644
File size: 6.8 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/** \class PileUpMergerPythia8
20 *
21 * Merges particles from pile-up sample into event
22 *
23 * \author M. Selvaggi - UCL, Louvain-la-Neuve
24 *
25 */
26
27#include "modules/PileUpMergerPythia8.h"
28
29#include "classes/DelphesClasses.h"
30#include "classes/DelphesFactory.h"
31#include "classes/DelphesPileUpReader.h"
32#include "classes/DelphesTF2.h"
33
34#include "ExRootAnalysis/ExRootClassifier.h"
35#include "ExRootAnalysis/ExRootFilter.h"
36#include "ExRootAnalysis/ExRootResult.h"
37
38#include "Pythia.h"
39
40#include "TDatabasePDG.h"
41#include "TFormula.h"
42#include "TLorentzVector.h"
43#include "TMath.h"
44#include "TObjArray.h"
45#include "TRandom3.h"
46#include "TString.h"
47
48#include <algorithm>
49#include <iostream>
50#include <sstream>
51#include <stdexcept>
52
53using namespace std;
54
55//------------------------------------------------------------------------------
56
57PileUpMergerPythia8::PileUpMergerPythia8() :
58 fFunction(0), fPythia(0), fItInputArray(0)
59{
60 fFunction = new DelphesTF2;
61}
62
63//------------------------------------------------------------------------------
64
65PileUpMergerPythia8::~PileUpMergerPythia8()
66{
67 delete fFunction;
68}
69
70//------------------------------------------------------------------------------
71
72void PileUpMergerPythia8::Init()
73{
74 const char *fileName;
75
76 fPileUpDistribution = GetInt("PileUpDistribution", 0);
77
78 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);
87
88 fPTMin = GetDouble("PTMin", 0.0);
89
90 fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
91 fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
92
93 fileName = GetString("ConfigFile", "MinBias.cmnd");
94 fPythia = new Pythia8::Pythia();
95 fPythia->readFile(fileName);
96
97 // import input array
98 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
99 fItInputArray = fInputArray->MakeIterator();
100
101 // create output arrays
102 fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
103 fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
104}
105
106//------------------------------------------------------------------------------
107
108void PileUpMergerPythia8::Finish()
109{
110 if(fPythia) delete fPythia;
111}
112
113//------------------------------------------------------------------------------
114
115void PileUpMergerPythia8::Process()
116{
117 TDatabasePDG *pdg = TDatabasePDG::Instance();
118 TParticlePDG *pdgParticle;
119 Int_t pid, status;
120 Float_t x, y, z, t, vx, vy;
121 Float_t px, py, pz, e;
122 Double_t dz, dphi, dt;
123 Int_t numberOfEvents, event, numberOfParticles, i;
124 Candidate *candidate, *vertex;
125 DelphesFactory *factory;
126
127 const Double_t c_light = 2.99792458E8;
128
129 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();
140 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
141 {
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;
155 }
156
157 factory = GetFactory();
158
159 vertex = factory->NewCandidate();
160 vertex->Position.SetXYZT(vx, vy, dz, dt);
161 fVertexOutputArray->Add(vertex);
162
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)
179 {
180 while(!fPythia->next())
181 ;
182
183 // --- Pile-up vertex smearing
184
185 fFunction->GetRandom2(dz, dt);
186
187 dt *= c_light * 1.0E3; // necessary in order to make t in mm/c
188 dz *= 1.0E3; // necessary in order to make z in mm
189
190 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
191
192 vx = 0.0;
193 vy = 0.0;
194 numberOfParticles = fPythia->event.size();
195 for(i = 1; i < numberOfParticles; ++i)
196 {
197 Pythia8::Particle &particle = fPythia->event[i];
198
199 status = particle.statusHepMC();
200
201 if(status != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
202
203 pid = particle.id();
204 px = particle.px();
205 py = particle.py();
206 pz = particle.pz();
207 e = particle.e();
208 x = particle.xProd();
209 y = particle.yProd();
210 z = particle.zProd();
211 t = particle.tProd();
212
213 candidate = factory->NewCandidate();
214
215 candidate->PID = pid;
216
217 candidate->Status = 1;
218
219 pdgParticle = pdg->GetParticle(pid);
220 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
221 candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
222
223 candidate->IsPU = 1;
224
225 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
226 candidate->Momentum.RotateZ(dphi);
227
228 x -= fInputBeamSpotX;
229 y -= fInputBeamSpotY;
230 candidate->Position.SetXYZT(x, y, z + dz, t + dt);
231 candidate->Position.RotateZ(dphi);
232 candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);
233
234 vx += candidate->Position.X();
235 vy += candidate->Position.Y();
236
237 fParticleOutputArray->Add(candidate);
238 }
239
240 if(numberOfParticles > 0)
241 {
242 vx /= numberOfParticles;
243 vy /= numberOfParticles;
244 }
245
246 vertex = factory->NewCandidate();
247 vertex->Position.SetXYZT(vx, vy, dz, dt);
248 vertex->IsPU = 1;
249
250 fVertexOutputArray->Add(vertex);
251 }
252}
253
254//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.