Fork me on GitHub

source: git/modules/PileUpMerger.cc@ 8b04b31

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8b04b31 was 7c0fcd5, checked in by Pavel Demin <demin@…>, 10 years ago

delete duplicate license file and prepend GPLv3 header to all source code files

  • Property mode set to 100644
File size: 5.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 PileUpMerger
20 *
21 * Merges particles from pile-up sample into event
22 *
23 *
24 * $Date: 2013-02-12 15:13:59 +0100 (Tue, 12 Feb 2013) $
25 * $Revision: 907 $
26 *
27 *
28 * \author M. Selvaggi - UCL, Louvain-la-Neuve
29 *
30 */
31
32#include "modules/PileUpMerger.h"
33
34#include "classes/DelphesClasses.h"
35#include "classes/DelphesFactory.h"
36#include "classes/DelphesTF2.h"
37#include "classes/DelphesPileUpReader.h"
38
39#include "ExRootAnalysis/ExRootResult.h"
40#include "ExRootAnalysis/ExRootFilter.h"
41#include "ExRootAnalysis/ExRootClassifier.h"
42
43#include "TMath.h"
44#include "TString.h"
45#include "TFormula.h"
46#include "TRandom3.h"
47#include "TObjArray.h"
48#include "TDatabasePDG.h"
49#include "TLorentzVector.h"
50
51#include <algorithm>
52#include <stdexcept>
53#include <iostream>
54#include <sstream>
55
56using namespace std;
57
58//------------------------------------------------------------------------------
59
60PileUpMerger::PileUpMerger() :
61 fFunction(0), fReader(0), fItInputArray(0)
62{
63 fFunction = new DelphesTF2;
64}
65
66
67//------------------------------------------------------------------------------
68
69PileUpMerger::~PileUpMerger()
70{
71 delete fFunction;
72}
73
74//------------------------------------------------------------------------------
75
76void PileUpMerger::Init()
77{
78 const char *fileName;
79
80 fPileUpDistribution = GetInt("PileUpDistribution", 0);
81
82 fMeanPileUp = GetDouble("MeanPileUp", 10);
83
84 fZVertexSpread = GetDouble("ZVertexSpread", 0.15);
85 fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09);
86
87 // read vertex smearing formula
88
89 fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));
90 fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);
91
92 fileName = GetString("PileUpFile", "MinBias.pileup");
93 fReader = new DelphesPileUpReader(fileName);
94
95 // import input array
96 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
97 fItInputArray = fInputArray->MakeIterator();
98
99 // create output arrays
100 fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles"));
101 fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
102}
103
104//------------------------------------------------------------------------------
105
106void PileUpMerger::Finish()
107{
108 if(fReader) delete fReader;
109}
110
111//------------------------------------------------------------------------------
112
113void PileUpMerger::Process()
114{
115 TDatabasePDG *pdg = TDatabasePDG::Instance();
116 TParticlePDG *pdgParticle;
117 Int_t pid;
118 Float_t x, y, z, t;
119 Float_t px, py, pz, e;
120 Double_t dz, dphi, dt;
121 Int_t numberOfEvents, event;
122 Long64_t allEntries, entry;
123 Candidate *candidate, *vertexcandidate;
124 DelphesFactory *factory;
125
126 const Double_t c_light = 2.99792458E8;
127
128 fItInputArray->Reset();
129
130 // --- Deal with Primary vertex first ------
131
132 fFunction->GetRandom2(dz, dt);
133
134 dt *= c_light*1.0E3; // necessary in order to make t in mm/c
135 dz *= 1.0E3; // necessary in order to make z in mm
136
137 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
138 {
139 z = candidate->Position.Z();
140 t = candidate->Position.T();
141 candidate->Position.SetZ(z + dz);
142 candidate->Position.SetT(t + dt);
143 fParticleOutputArray->Add(candidate);
144 }
145
146 factory = GetFactory();
147
148 vertexcandidate = factory->NewCandidate();
149 vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
150 fVertexOutputArray->Add(vertexcandidate);
151
152 // --- Then with pile-up vertices ------
153
154 switch(fPileUpDistribution)
155 {
156 case 0:
157 numberOfEvents = gRandom->Poisson(fMeanPileUp);
158 break;
159 case 1:
160 numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1);
161 break;
162 default:
163 numberOfEvents = gRandom->Poisson(fMeanPileUp);
164 break;
165 }
166
167 allEntries = fReader->GetEntries();
168
169 for(event = 0; event < numberOfEvents; ++event)
170 {
171 do
172 {
173 entry = TMath::Nint(gRandom->Rndm()*allEntries);
174 }
175 while(entry >= allEntries);
176
177 fReader->ReadEntry(entry);
178
179 // --- Pile-up vertex smearing
180
181 fFunction->GetRandom2(dz, dt);
182
183 dt *= c_light*1.0E3; // necessary in order to make t in mm/c
184 dz *= 1.0E3; // necessary in order to make z in mm
185
186 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
187
188 vertexcandidate = factory->NewCandidate();
189 vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt);
190 vertexcandidate->IsPU = 1;
191
192 fVertexOutputArray->Add(vertexcandidate);
193
194 while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e))
195 {
196 candidate = factory->NewCandidate();
197
198 candidate->PID = pid;
199
200 candidate->Status = 1;
201
202 pdgParticle = pdg->GetParticle(pid);
203 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
204 candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
205
206 candidate->IsPU = 1;
207
208 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
209 candidate->Momentum.RotateZ(dphi);
210
211 candidate->Position.SetXYZT(x, y, z + dz, t + dt);
212 candidate->Position.RotateZ(dphi);
213
214 fParticleOutputArray->Add(candidate);
215 }
216 }
217}
218
219//------------------------------------------------------------------------------
220
Note: See TracBrowser for help on using the repository browser.