Fork me on GitHub

source: svn/trunk/modules/PileUpMergerPythia8.cc

Last change on this file was 1395, checked in by Pavel Demin, 10 years ago

synchronize with GitHub

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 3.6 KB
Line 
1/** \class PileUpMergerPythia8
2 *
3 * Merges particles from pile-up sample into event
4 *
5 *
6 * $Date: 2014-07-08 11:59:26 +0000 (Tue, 08 Jul 2014) $
7 * $Revision: 1395 $
8 *
9 *
10 * \author P. Selvaggi - UCL, Louvain-la-Neuve
11 *
12 */
13
14#include "modules/PileUpMergerPythia8.h"
15
16#include "classes/DelphesClasses.h"
17#include "classes/DelphesFactory.h"
18#include "classes/DelphesFormula.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 "Pythia8/Pythia.h"
26
27#include "TMath.h"
28#include "TString.h"
29#include "TFormula.h"
30#include "TRandom3.h"
31#include "TObjArray.h"
32#include "TDatabasePDG.h"
33#include "TLorentzVector.h"
34
35#include <algorithm>
36#include <stdexcept>
37#include <iostream>
38#include <sstream>
39
40using namespace std;
41
42//------------------------------------------------------------------------------
43
44PileUpMergerPythia8::PileUpMergerPythia8() :
45 fPythia(0), fItInputArray(0)
46{
47}
48
49//------------------------------------------------------------------------------
50
51PileUpMergerPythia8::~PileUpMergerPythia8()
52{
53}
54
55//------------------------------------------------------------------------------
56
57void PileUpMergerPythia8::Init()
58{
59 const char *fileName;
60
61 fMeanPileUp = GetDouble("MeanPileUp", 10);
62 fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3;
63
64 fPTMin = GetDouble("PTMin", 0.0);
65
66 fileName = GetString("ConfigFile", "MinBias.cmnd");
67 fPythia = new Pythia8::Pythia();
68 fPythia->readFile(fileName);
69
70 // import input array
71 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
72 fItInputArray = fInputArray->MakeIterator();
73
74 // create output arrays
75 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
76}
77
78//------------------------------------------------------------------------------
79
80void PileUpMergerPythia8::Finish()
81{
82 if(fPythia) delete fPythia;
83}
84
85//------------------------------------------------------------------------------
86
87void PileUpMergerPythia8::Process()
88{
89 TDatabasePDG *pdg = TDatabasePDG::Instance();
90 TParticlePDG *pdgParticle;
91 Int_t pid;
92 Float_t x, y, z, t;
93 Float_t px, py, pz, e;
94 Double_t dz, dphi;
95 Int_t poisson, event, i;
96 Candidate *candidate;
97 DelphesFactory *factory;
98
99 fItInputArray->Reset();
100 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
101 {
102 fOutputArray->Add(candidate);
103 }
104
105 factory = GetFactory();
106
107 poisson = gRandom->Poisson(fMeanPileUp);
108
109 for(event = 0; event < poisson; ++event)
110 {
111 while(!fPythia->next());
112
113 dz = gRandom->Gaus(0.0, fZVertexSpread);
114 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
115
116 for(i = 0; i < fPythia->event.size(); ++i)
117 {
118 Pythia8::Particle &particle = fPythia->event[i];
119
120 if(particle.status() != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
121
122 pid = particle.id();
123 px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.e();
124 x = particle.xProd(); y = particle.yProd(); z = particle.zProd(); t = particle.tProd();
125
126 candidate = factory->NewCandidate();
127
128 candidate->PID = pid;
129
130 candidate->Status = 1;
131
132 pdgParticle = pdg->GetParticle(pid);
133 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
134 candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
135
136 candidate->IsPU = 1;
137
138 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
139 candidate->Momentum.RotateZ(dphi);
140
141 candidate->Position.SetXYZT(x, y, z + dz, t);
142 candidate->Position.RotateZ(dphi);
143
144 fOutputArray->Add(candidate);
145 }
146 }
147}
148
149//------------------------------------------------------------------------------
150
Note: See TracBrowser for help on using the repository browser.