1 | #include <stdexcept>
|
---|
2 | #include <iostream>
|
---|
3 | #include <sstream>
|
---|
4 |
|
---|
5 | #include <signal.h>
|
---|
6 |
|
---|
7 | #include "Pythia8/Pythia.h"
|
---|
8 |
|
---|
9 | #include "TROOT.h"
|
---|
10 | #include "TApplication.h"
|
---|
11 |
|
---|
12 | #include "TFile.h"
|
---|
13 | #include "TObjArray.h"
|
---|
14 | #include "TStopwatch.h"
|
---|
15 | #include "TDatabasePDG.h"
|
---|
16 | #include "TParticlePDG.h"
|
---|
17 | #include "TLorentzVector.h"
|
---|
18 |
|
---|
19 | #include "modules/Delphes.h"
|
---|
20 | #include "classes/DelphesClasses.h"
|
---|
21 | #include "classes/DelphesFactory.h"
|
---|
22 |
|
---|
23 | #include "ExRootAnalysis/ExRootTreeWriter.h"
|
---|
24 | #include "ExRootAnalysis/ExRootTreeBranch.h"
|
---|
25 | #include "ExRootAnalysis/ExRootProgressBar.h"
|
---|
26 |
|
---|
27 | using namespace std;
|
---|
28 |
|
---|
29 | //---------------------------------------------------------------------------
|
---|
30 |
|
---|
31 | void ConvertInput(Long64_t eventCounter, Pythia8::Pythia *pythia,
|
---|
32 | ExRootTreeBranch *branch, DelphesFactory *factory,
|
---|
33 | TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray,
|
---|
34 | TStopwatch *readStopWatch, TStopwatch *procStopWatch)
|
---|
35 | {
|
---|
36 | int i;
|
---|
37 |
|
---|
38 | HepMCEvent *element;
|
---|
39 | Candidate *candidate;
|
---|
40 | TDatabasePDG *pdg;
|
---|
41 | TParticlePDG *pdgParticle;
|
---|
42 | Int_t pdgCode;
|
---|
43 |
|
---|
44 | Int_t pid, status;
|
---|
45 | Double_t px, py, pz, e, mass;
|
---|
46 | Double_t x, y, z, t;
|
---|
47 |
|
---|
48 | // event information
|
---|
49 | element = static_cast<HepMCEvent *>(branch->NewEntry());
|
---|
50 |
|
---|
51 | element->Number = eventCounter;
|
---|
52 |
|
---|
53 | element->ProcessID = pythia->info.code();
|
---|
54 | element->MPI = 1;
|
---|
55 | element->Weight = pythia->info.weight();
|
---|
56 | element->Scale = pythia->info.QRen();
|
---|
57 | element->AlphaQED = pythia->info.alphaEM();
|
---|
58 | element->AlphaQCD = pythia->info.alphaS();
|
---|
59 |
|
---|
60 | element->ID1 = pythia->info.id1();
|
---|
61 | element->ID2 = pythia->info.id2();
|
---|
62 | element->X1 = pythia->info.x1();
|
---|
63 | element->X2 = pythia->info.x2();
|
---|
64 | element->ScalePDF = pythia->info.QFac();
|
---|
65 | element->PDF1 = pythia->info.pdf1();
|
---|
66 | element->PDF2 = pythia->info.pdf2();
|
---|
67 |
|
---|
68 | element->ReadTime = readStopWatch->RealTime();
|
---|
69 | element->ProcTime = procStopWatch->RealTime();
|
---|
70 |
|
---|
71 | pdg = TDatabasePDG::Instance();
|
---|
72 |
|
---|
73 | for(i = 0; i < pythia->event.size(); ++i)
|
---|
74 | {
|
---|
75 | Pythia8::Particle &particle = pythia->event[i];
|
---|
76 |
|
---|
77 | pid = particle.id();
|
---|
78 | status = pythia->event.statusHepMC(i);
|
---|
79 | px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.e(); mass = particle.m();
|
---|
80 | x = particle.xProd(); y = particle.yProd(); z = particle.zProd(); t = particle.tProd();
|
---|
81 |
|
---|
82 | candidate = factory->NewCandidate();
|
---|
83 |
|
---|
84 | candidate->PID = pid;
|
---|
85 | pdgCode = TMath::Abs(candidate->PID);
|
---|
86 |
|
---|
87 | candidate->Status = status;
|
---|
88 |
|
---|
89 | candidate->M1 = particle.mother1() - 1;
|
---|
90 | candidate->M2 = particle.mother2() - 1;
|
---|
91 |
|
---|
92 | candidate->D1 = particle.daughter1() - 1;
|
---|
93 | candidate->D2 = particle.daughter2() - 1;
|
---|
94 |
|
---|
95 | pdgParticle = pdg->GetParticle(pid);
|
---|
96 | candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
|
---|
97 | candidate->Mass = mass;
|
---|
98 |
|
---|
99 | candidate->Momentum.SetPxPyPzE(px, py, pz, e);
|
---|
100 |
|
---|
101 | candidate->Position.SetXYZT(x, y, z, t);
|
---|
102 |
|
---|
103 | allParticleOutputArray->Add(candidate);
|
---|
104 |
|
---|
105 | if(!pdgParticle) continue;
|
---|
106 |
|
---|
107 | if(status == 1)
|
---|
108 | {
|
---|
109 | stableParticleOutputArray->Add(candidate);
|
---|
110 | }
|
---|
111 | else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
|
---|
112 | {
|
---|
113 | partonOutputArray->Add(candidate);
|
---|
114 | }
|
---|
115 | }
|
---|
116 | }
|
---|
117 |
|
---|
118 | //---------------------------------------------------------------------------
|
---|
119 |
|
---|
120 | static bool interrupted = false;
|
---|
121 |
|
---|
122 | void SignalHandler(int sig)
|
---|
123 | {
|
---|
124 | interrupted = true;
|
---|
125 | }
|
---|
126 |
|
---|
127 | //---------------------------------------------------------------------------
|
---|
128 |
|
---|
129 | int main(int argc, char *argv[])
|
---|
130 | {
|
---|
131 | char appName[] = "DelphesPythia8";
|
---|
132 | stringstream message;
|
---|
133 | TFile *outputFile = 0;
|
---|
134 | TStopwatch readStopWatch, procStopWatch;
|
---|
135 | ExRootTreeWriter *treeWriter = 0;
|
---|
136 | ExRootTreeBranch *branchEvent = 0;
|
---|
137 | ExRootConfReader *confReader = 0;
|
---|
138 | Delphes *modularDelphes = 0;
|
---|
139 | DelphesFactory *factory = 0;
|
---|
140 | TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
|
---|
141 | Long64_t eventCounter, errorCounter;
|
---|
142 | Long64_t numberOfEvents, timesAllowErrors;
|
---|
143 |
|
---|
144 | Pythia8::Pythia *pythia = 0;
|
---|
145 |
|
---|
146 | if(argc != 4)
|
---|
147 | {
|
---|
148 | cout << " Usage: " << appName << " config_file" << " pythia_card" << " output_file" << endl;
|
---|
149 | cout << " config_file - configuration file in Tcl format," << endl;
|
---|
150 | cout << " pythia_card - Pythia8 configuration file," << endl;
|
---|
151 | cout << " output_file - output file in ROOT format." << endl;
|
---|
152 | return 1;
|
---|
153 | }
|
---|
154 |
|
---|
155 | signal(SIGINT, SignalHandler);
|
---|
156 |
|
---|
157 | gROOT->SetBatch();
|
---|
158 |
|
---|
159 | int appargc = 1;
|
---|
160 | char *appargv[] = {appName};
|
---|
161 | TApplication app(appName, &appargc, appargv);
|
---|
162 |
|
---|
163 | try
|
---|
164 | {
|
---|
165 | outputFile = TFile::Open(argv[3], "CREATE");
|
---|
166 |
|
---|
167 | if(outputFile == NULL)
|
---|
168 | {
|
---|
169 | message << "can't create output file " << argv[3];
|
---|
170 | throw runtime_error(message.str());
|
---|
171 | }
|
---|
172 |
|
---|
173 | treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
|
---|
174 |
|
---|
175 | branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
|
---|
176 |
|
---|
177 | confReader = new ExRootConfReader;
|
---|
178 | confReader->ReadFile(argv[1]);
|
---|
179 |
|
---|
180 | modularDelphes = new Delphes("Delphes");
|
---|
181 | modularDelphes->SetConfReader(confReader);
|
---|
182 | modularDelphes->SetTreeWriter(treeWriter);
|
---|
183 |
|
---|
184 | factory = modularDelphes->GetFactory();
|
---|
185 | allParticleOutputArray = modularDelphes->ExportArray("allParticles");
|
---|
186 | stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
|
---|
187 | partonOutputArray = modularDelphes->ExportArray("partons");
|
---|
188 |
|
---|
189 | modularDelphes->InitTask();
|
---|
190 |
|
---|
191 | // Initialize pythia
|
---|
192 | pythia = new Pythia8::Pythia;
|
---|
193 |
|
---|
194 | if(pythia == NULL)
|
---|
195 | {
|
---|
196 | throw runtime_error("can't create Pythia instance");
|
---|
197 | }
|
---|
198 |
|
---|
199 | // Read in commands from configuration file
|
---|
200 | pythia->readFile(argv[2]);
|
---|
201 |
|
---|
202 | // Extract settings to be used in the main program
|
---|
203 | numberOfEvents = pythia->mode("Main:numberOfEvents");
|
---|
204 | timesAllowErrors = pythia->mode("Main:timesAllowErrors");
|
---|
205 |
|
---|
206 | pythia->init();
|
---|
207 |
|
---|
208 | // ExRootProgressBar progressBar(numberOfEvents - 1);
|
---|
209 | ExRootProgressBar progressBar(-1);
|
---|
210 |
|
---|
211 | // Loop over all events
|
---|
212 | errorCounter = 0;
|
---|
213 | treeWriter->Clear();
|
---|
214 | modularDelphes->Clear();
|
---|
215 | readStopWatch.Start();
|
---|
216 | for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
|
---|
217 | {
|
---|
218 | if(!pythia->next())
|
---|
219 | {
|
---|
220 | // If failure because reached end of file then exit event loop
|
---|
221 | if (pythia->info.atEndOfFile())
|
---|
222 | {
|
---|
223 | cerr << "Aborted since reached end of Les Houches Event File" << endl;
|
---|
224 | break;
|
---|
225 | }
|
---|
226 |
|
---|
227 | // First few failures write off as "acceptable" errors, then quit
|
---|
228 | if (++errorCounter < timesAllowErrors) continue;
|
---|
229 | cerr << "Event generation aborted prematurely, owing to error!" << endl;
|
---|
230 | break;
|
---|
231 | }
|
---|
232 |
|
---|
233 | readStopWatch.Stop();
|
---|
234 |
|
---|
235 | procStopWatch.Start();
|
---|
236 | ConvertInput(eventCounter, pythia, branchEvent, factory,
|
---|
237 | allParticleOutputArray, stableParticleOutputArray, partonOutputArray,
|
---|
238 | &readStopWatch, &procStopWatch);
|
---|
239 | modularDelphes->ProcessTask();
|
---|
240 | procStopWatch.Stop();
|
---|
241 |
|
---|
242 | treeWriter->Fill();
|
---|
243 |
|
---|
244 | treeWriter->Clear();
|
---|
245 | modularDelphes->Clear();
|
---|
246 |
|
---|
247 | readStopWatch.Start();
|
---|
248 | progressBar.Update(eventCounter, eventCounter);
|
---|
249 | }
|
---|
250 |
|
---|
251 | progressBar.Update(eventCounter, eventCounter, kTRUE);
|
---|
252 | progressBar.Finish();
|
---|
253 |
|
---|
254 | pythia->statistics();
|
---|
255 |
|
---|
256 | modularDelphes->FinishTask();
|
---|
257 | treeWriter->Write();
|
---|
258 |
|
---|
259 | cout << "** Exiting..." << endl;
|
---|
260 |
|
---|
261 | delete pythia;
|
---|
262 | delete modularDelphes;
|
---|
263 | delete confReader;
|
---|
264 | delete treeWriter;
|
---|
265 | delete outputFile;
|
---|
266 |
|
---|
267 | return 0;
|
---|
268 | }
|
---|
269 | catch(runtime_error &e)
|
---|
270 | {
|
---|
271 | if(treeWriter) delete treeWriter;
|
---|
272 | if(outputFile) delete outputFile;
|
---|
273 | cerr << "** ERROR: " << e.what() << endl;
|
---|
274 | return 1;
|
---|
275 | }
|
---|
276 | }
|
---|