Fork me on GitHub

source: git/modules/PileUpMergerPythia8.cc@ 8c72f09

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8c72f09 was cab38f6, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

remove svn tags and fix formatting

  • Property mode set to 100644
File size: 4.4 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 P. 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/DelphesFormula.h"
32#include "classes/DelphesPileUpReader.h"
33
34#include "ExRootAnalysis/ExRootResult.h"
35#include "ExRootAnalysis/ExRootFilter.h"
36#include "ExRootAnalysis/ExRootClassifier.h"
37
38#include "Pythia.h"
39
40#include "TMath.h"
41#include "TString.h"
42#include "TFormula.h"
43#include "TRandom3.h"
44#include "TObjArray.h"
45#include "TDatabasePDG.h"
46#include "TLorentzVector.h"
47
48#include <algorithm>
49#include <stdexcept>
50#include <iostream>
51#include <sstream>
52
53using namespace std;
54
55//------------------------------------------------------------------------------
56
57PileUpMergerPythia8::PileUpMergerPythia8() :
58 fPythia(0), fItInputArray(0)
59{
60}
61
62//------------------------------------------------------------------------------
63
64PileUpMergerPythia8::~PileUpMergerPythia8()
65{
66}
67
68//------------------------------------------------------------------------------
69
70void PileUpMergerPythia8::Init()
71{
72 const char *fileName;
73
74 fMeanPileUp = GetDouble("MeanPileUp", 10);
75 fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
76
77 fPTMin = GetDouble("PTMin", 0.0);
78
79 fileName = GetString("ConfigFile", "MinBias.cmnd");
80 fPythia = new Pythia8::Pythia();
81 fPythia->readFile(fileName);
82
83 // import input array
84 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
85 fItInputArray = fInputArray->MakeIterator();
86
87 // create output arrays
88 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
89}
90
91//------------------------------------------------------------------------------
92
93void PileUpMergerPythia8::Finish()
94{
95 if(fPythia) delete fPythia;
96}
97
98//------------------------------------------------------------------------------
99
100void PileUpMergerPythia8::Process()
101{
102 TDatabasePDG *pdg = TDatabasePDG::Instance();
103 TParticlePDG *pdgParticle;
104 Int_t pid;
105 Float_t x, y, z, t;
106 Float_t px, py, pz, e;
107 Double_t dz, dphi;
108 Int_t poisson, event, i;
109 Candidate *candidate;
110 DelphesFactory *factory;
111
112 fItInputArray->Reset();
113 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
114 {
115 fOutputArray->Add(candidate);
116 }
117
118 factory = GetFactory();
119
120 poisson = gRandom->Poisson(fMeanPileUp);
121
122 for(event = 0; event < poisson; ++event)
123 {
124 while(!fPythia->next());
125
126 dz = gRandom->Gaus(0.0, fZVertexSpread);
127 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
128
129 for(i = 1; i < fPythia->event.size(); ++i)
130 {
131 Pythia8::Particle &particle = fPythia->event[i];
132
133 if(particle.status() != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
134
135 pid = particle.id();
136 px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.e();
137 x = particle.xProd(); y = particle.yProd(); z = particle.zProd(); t = particle.tProd();
138
139 candidate = factory->NewCandidate();
140
141 candidate->PID = pid;
142
143 candidate->Status = 1;
144
145 pdgParticle = pdg->GetParticle(pid);
146 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
147 candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
148
149 candidate->IsPU = 1;
150
151 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
152 candidate->Momentum.RotateZ(dphi);
153
154 candidate->Position.SetXYZT(x, y, z + dz, t);
155 candidate->Position.RotateZ(dphi);
156
157 fOutputArray->Add(candidate);
158 }
159 }
160}
161
162//------------------------------------------------------------------------------
163
Note: See TracBrowser for help on using the repository browser.