Fork me on GitHub

source: git/modules/PileUpMergerPythia8.cc@ 94dcf80

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 94dcf80 was c466459, checked in by pavel <pavel@…>, 11 years ago

fix genMakefile and PileUpMergerPythia8

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