Fork me on GitHub

source: git/modules/PileUpMerger.cc@ f6b9fec

ImprovedOutputFile Timing dual_readout llp
Last change on this file since f6b9fec was 0ed696f, checked in by pavel <pavel@…>, 11 years ago

fix formatting in PileUpMerger

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