Fork me on GitHub

source: svn/trunk/modules/PileUpMergerPythia8.cc@ 1167

Last change on this file since 1167 was 1167, checked in by Pavel Demin, 11 years ago

add PileUpMergerPythia8

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 3.9 KB
Line 
1/** \class PileUpMergerPythia8
2 *
3 * Merges particles from pile-up sample into event
4 *
5 *
6 * $Date: 2013-07-03 23:04:05 +0000 (Wed, 03 Jul 2013) $
7 * $Revision: 1167 $
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 "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 fCMEnergy = GetDouble("CMEnergy", 14000.0);
65 fPTMin = GetDouble("PTMin", 0.5);
66
67 fRandomSeed = GetInt("RandomSeed", 1);
68 fBeamAID = GetInt("BeamAID", 2212);
69 fBeamBID = GetInt("BeamBID", 2212);
70
71 // Pythia8 pile-up initialization
72 fPythia = new Pythia8::Pythia();
73 fPythia->readString("SoftQCD:minBias = on");
74 fPythia->init(fBeamAID, fBeamBID, fCMEnergy);
75 fPythia->rndm.init(fRandomSeed);
76
77 // import input array
78 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
79 fItInputArray = fInputArray->MakeIterator();
80
81 // create output arrays
82 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
83}
84
85//------------------------------------------------------------------------------
86
87void PileUpMergerPythia8::Finish()
88{
89 if(fPythia) delete fPythia;
90}
91
92//------------------------------------------------------------------------------
93
94void PileUpMergerPythia8::Process()
95{
96 TDatabasePDG *pdg = TDatabasePDG::Instance();
97 TParticlePDG *pdgParticle;
98 Int_t pid;
99 Float_t x, y, z, t;
100 Float_t px, py, pz, e;
101 Double_t dz, dphi;
102 Int_t poisson, event, i;
103 Long64_t allEntries, entry;
104 Candidate *candidate;
105 DelphesFactory *factory;
106
107 fItInputArray->Reset();
108 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
109 {
110 fOutputArray->Add(candidate);
111 }
112
113 factory = GetFactory();
114
115 poisson = gRandom->Poisson(fMeanPileUp);
116
117 for(event = 0; event < poisson; ++event)
118 {
119 while(!fPythia->next());
120
121 dz = gRandom->Gaus(0.0, fZVertexSpread);
122 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
123
124 for(i = 0; i < fPythia->event.size(); ++i)
125 {
126 Pythia8::Particle &particle = pythia->event[i];
127
128 if(particle.status() != 1 || !particle.isVisible() || particle.pT() <= fPTMin) continue;
129
130 pid = particle.id();
131 px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.e();
132 x = particle.xProd(); y = particle.yProd(); z = particle.zProd(); t = particle.tProd();
133
134 candidate = factory->NewCandidate();
135
136 candidate->PID = pid;
137
138 candidate->Status = 1;
139
140 pdgParticle = pdg->GetParticle(pid);
141 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
142 candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
143
144 candidate->IsPU = 1;
145
146 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
147 candidate->Momentum.RotateZ(dphi);
148
149 candidate->Position.SetXYZT(x, y, z + dz, t);
150 candidate->Position.RotateZ(dphi);
151
152 fOutputArray->Add(candidate);
153 }
154 }
155}
156
157//------------------------------------------------------------------------------
158
Note: See TracBrowser for help on using the repository browser.