Fork me on GitHub

source: git/readers/DelphesPythia8.cpp@ bf2c695

Last change on this file since bf2c695 was 5eb063e, checked in by Roberto Preghenella <preghenella@…>, 4 years ago

Temporary fix UserHooks API changed in Pythia 8.3 series

  • Property mode set to 100644
File size: 12.7 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#include <iostream>
20#include <sstream>
21#include <stdexcept>
22#include <string>
23
24#include <signal.h>
25
26#include "Pythia.h"
27#include "Pythia8Plugins/CombineMatchingInput.h"
28
29#include "TApplication.h"
30#include "TROOT.h"
31
32#include "TDatabasePDG.h"
33#include "TFile.h"
34#include "TLorentzVector.h"
35#include "TObjArray.h"
36#include "TParticlePDG.h"
37#include "TStopwatch.h"
38
39#include "classes/DelphesClasses.h"
40#include "classes/DelphesFactory.h"
41#include "classes/DelphesLHEFReader.h"
42#include "modules/Delphes.h"
43
44#include "ExRootAnalysis/ExRootProgressBar.h"
45#include "ExRootAnalysis/ExRootTreeBranch.h"
46#include "ExRootAnalysis/ExRootTreeWriter.h"
47
48using namespace std;
49
50//---------------------------------------------------------------------------
51
52void ConvertInput(Long64_t eventCounter, Pythia8::Pythia *pythia,
53 ExRootTreeBranch *branch, DelphesFactory *factory,
54 TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray,
55 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
56{
57 int i;
58
59 HepMCEvent *element;
60 Candidate *candidate;
61 TDatabasePDG *pdg;
62 TParticlePDG *pdgParticle;
63 Int_t pdgCode;
64
65 Int_t pid, status;
66 Double_t px, py, pz, e, mass;
67 Double_t x, y, z, t;
68
69 // event information
70 element = static_cast<HepMCEvent *>(branch->NewEntry());
71
72 element->Number = eventCounter;
73
74 element->ProcessID = pythia->info.code();
75 element->MPI = 1;
76 element->Weight = pythia->info.weight();
77 element->Scale = pythia->info.QRen();
78 element->AlphaQED = pythia->info.alphaEM();
79 element->AlphaQCD = pythia->info.alphaS();
80
81 element->ID1 = pythia->info.id1();
82 element->ID2 = pythia->info.id2();
83 element->X1 = pythia->info.x1();
84 element->X2 = pythia->info.x2();
85 element->ScalePDF = pythia->info.QFac();
86 element->PDF1 = pythia->info.pdf1();
87 element->PDF2 = pythia->info.pdf2();
88
89 element->ReadTime = readStopWatch->RealTime();
90 element->ProcTime = procStopWatch->RealTime();
91
92 pdg = TDatabasePDG::Instance();
93
94 for(i = 1; i < pythia->event.size(); ++i)
95 {
96 Pythia8::Particle &particle = pythia->event[i];
97
98 pid = particle.id();
99 status = particle.statusHepMC();
100 px = particle.px();
101 py = particle.py();
102 pz = particle.pz();
103 e = particle.e();
104 mass = particle.m();
105 x = particle.xProd();
106 y = particle.yProd();
107 z = particle.zProd();
108 t = particle.tProd();
109
110 candidate = factory->NewCandidate();
111
112 candidate->PID = pid;
113 pdgCode = TMath::Abs(candidate->PID);
114
115 candidate->Status = status;
116
117 candidate->M1 = particle.mother1() - 1;
118 candidate->M2 = particle.mother2() - 1;
119
120 candidate->D1 = particle.daughter1() - 1;
121 candidate->D2 = particle.daughter2() - 1;
122
123 pdgParticle = pdg->GetParticle(pid);
124 candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
125 candidate->Mass = mass;
126
127 candidate->Momentum.SetPxPyPzE(px, py, pz, e);
128
129 candidate->Position.SetXYZT(x, y, z, t);
130
131 allParticleOutputArray->Add(candidate);
132
133 if(!pdgParticle) continue;
134
135 if(status == 1)
136 {
137 stableParticleOutputArray->Add(candidate);
138 }
139 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
140 {
141 partonOutputArray->Add(candidate);
142 }
143 }
144}
145
146//---------------------------------------------------------------------------
147
148static bool interrupted = false;
149
150void SignalHandler(int sig)
151{
152 interrupted = true;
153}
154
155//---------------------------------------------------------------------------
156
157/*
158Single-particle gun. The particle must be a colour singlet.
159Input: flavour, energy, direction (theta, phi).
160If theta < 0 then random choice over solid angle.
161Optional final argument to put particle at rest => E = m.
162from pythia8 example 21
163*/
164
165void fillParticle(int id, double pMax, double etaMax,
166 Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
167{
168 // Reset event record to allow for new event.
169 event.reset();
170
171 // Generate uniform pt and eta.
172 double pt, eta, phi, pp, ee, mm;
173
174 // pMin = 0.1 GeV for single particles
175 pp = pow(10, -1.0 + (log10(pMax) + 1.0) * rndm.flat());
176 eta = (2.0 * rndm.flat() - 1.0) * etaMax;
177 phi = 2.0 * M_PI * rndm.flat();
178 mm = pdt.mSel(id);
179 ee = Pythia8::sqrtpos(pp * pp + mm * mm);
180 pt = pp / cosh(eta);
181
182 // Store the particle in the event record.
183 event.append(id, 1, 0, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
184}
185
186//---------------------------------------------------------------------------
187
188void fillPartons(int id, double pMax, double etaMax,
189 Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
190{
191 // Reset event record to allow for new event.
192 event.reset();
193
194 // Generate uniform pt and eta.
195 double pt, eta, phi, pp, ee, mm;
196
197 // pMin = 1 GeV for jets
198 pp = pow(10, log10(pMax) * rndm.flat());
199 eta = (2.0 * rndm.flat() - 1.0) * etaMax;
200 phi = 2.0 * M_PI * rndm.flat();
201 mm = pdt.mSel(id);
202 ee = Pythia8::sqrtpos(pp * pp + mm * mm);
203 pt = pp / cosh(eta);
204
205 if((id == 4 || id == 5) && pt < 10.0) return;
206
207 if(id == 21)
208 {
209 event.append(21, 23, 101, 102, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee);
210 event.append(21, 23, 102, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee);
211 }
212 else
213 {
214 event.append(id, 23, 101, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
215 event.append(-id, 23, 0, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee, mm);
216 }
217}
218
219//---------------------------------------------------------------------------
220
221int main(int argc, char *argv[])
222{
223 char appName[] = "DelphesPythia8";
224 stringstream message;
225 FILE *inputFile = 0;
226 TFile *outputFile = 0;
227 TStopwatch readStopWatch, procStopWatch;
228 ExRootTreeWriter *treeWriter = 0;
229 ExRootTreeBranch *branchEvent = 0;
230 ExRootTreeBranch *branchEventLHEF = 0, *branchWeightLHEF = 0;
231 ExRootConfReader *confReader = 0;
232 Delphes *modularDelphes = 0;
233 DelphesFactory *factory = 0;
234 TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
235 TObjArray *stableParticleOutputArrayLHEF = 0, *allParticleOutputArrayLHEF = 0, *partonOutputArrayLHEF = 0;
236 DelphesLHEFReader *reader = 0;
237 Long64_t eventCounter, errorCounter;
238 Long64_t numberOfEvents, timesAllowErrors;
239 Bool_t spareFlag1;
240 Int_t spareMode1;
241 Double_t spareParm1, spareParm2;
242
243 Pythia8::Pythia *pythia = 0;
244
245 // for matching
246 Pythia8::CombineMatchingInput *combined = 0;
247 Pythia8::UserHooks *matching = 0;
248
249 if(argc != 4)
250 {
251 cout << " Usage: " << appName << " config_file"
252 << " pythia_card"
253 << " output_file" << endl;
254 cout << " config_file - configuration file in Tcl format," << endl;
255 cout << " pythia_card - Pythia8 configuration file," << endl;
256 cout << " output_file - output file in ROOT format." << endl;
257 return 1;
258 }
259
260 signal(SIGINT, SignalHandler);
261
262 gROOT->SetBatch();
263
264 int appargc = 1;
265 char *appargv[] = {appName};
266 TApplication app(appName, &appargc, appargv);
267
268 try
269 {
270 outputFile = TFile::Open(argv[3], "CREATE");
271
272 if(outputFile == NULL)
273 {
274 message << "can't create output file " << argv[3];
275 throw runtime_error(message.str());
276 }
277
278 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
279
280 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
281
282 confReader = new ExRootConfReader;
283 confReader->ReadFile(argv[1]);
284
285 modularDelphes = new Delphes("Delphes");
286 modularDelphes->SetConfReader(confReader);
287 modularDelphes->SetTreeWriter(treeWriter);
288
289 factory = modularDelphes->GetFactory();
290 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
291 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
292 partonOutputArray = modularDelphes->ExportArray("partons");
293
294 // Initialize Pythia
295 pythia = new Pythia8::Pythia;
296
297 // jet matching
298#if PYTHIA_VERSION_INTEGER < 8300
299 matching = combined->getHook(*pythia);
300 if(!matching)
301 {
302 throw runtime_error("can't do matching");
303 }
304 pythia->setUserHooksPtr(matching);
305#endif
306
307 if(pythia == NULL)
308 {
309 throw runtime_error("can't create Pythia instance");
310 }
311
312 // Read in commands from configuration file
313 if(!pythia->readFile(argv[2]))
314 {
315 message << "can't read Pythia8 configuration file " << argv[2] << endl;
316 throw runtime_error(message.str());
317 }
318
319 // Extract settings to be used in the main program
320 numberOfEvents = pythia->mode("Main:numberOfEvents");
321 timesAllowErrors = pythia->mode("Main:timesAllowErrors");
322
323 spareFlag1 = pythia->flag("Main:spareFlag1");
324 spareMode1 = pythia->mode("Main:spareMode1");
325 spareParm1 = pythia->parm("Main:spareParm1");
326 spareParm2 = pythia->parm("Main:spareParm2");
327
328 // Check if particle gun
329 if(!spareFlag1)
330 {
331 inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
332 if(inputFile)
333 {
334 reader = new DelphesLHEFReader;
335 reader->SetInputFile(inputFile);
336
337 branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
338 branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
339
340 allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
341 stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
342 partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
343 }
344 }
345
346 modularDelphes->InitTask();
347
348 pythia->init();
349
350 // ExRootProgressBar progressBar(numberOfEvents - 1);
351 ExRootProgressBar progressBar(-1);
352
353 // Loop over all events
354 errorCounter = 0;
355 treeWriter->Clear();
356 modularDelphes->Clear();
357 readStopWatch.Start();
358 for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
359 {
360 while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF, stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady())
361 ;
362
363 if(spareFlag1)
364 {
365 if((spareMode1 >= 1 && spareMode1 <= 5) || spareMode1 == 21)
366 {
367 fillPartons(spareMode1, spareParm1, spareParm2, pythia->event, pythia->particleData, pythia->rndm);
368 }
369 else
370 {
371 fillParticle(spareMode1, spareParm1, spareParm2, pythia->event, pythia->particleData, pythia->rndm);
372 }
373 }
374
375 if(!pythia->next())
376 {
377 // If failure because reached end of file then exit event loop
378 if(pythia->info.atEndOfFile())
379 {
380 cerr << "Aborted since reached end of Les Houches Event File" << endl;
381 break;
382 }
383
384 // First few failures write off as "acceptable" errors, then quit
385 if(++errorCounter > timesAllowErrors)
386 {
387 cerr << "Event generation aborted prematurely, owing to error!" << endl;
388 break;
389 }
390
391 modularDelphes->Clear();
392 reader->Clear();
393 continue;
394 }
395
396 readStopWatch.Stop();
397
398 procStopWatch.Start();
399 ConvertInput(eventCounter, pythia, branchEvent, factory,
400 allParticleOutputArray, stableParticleOutputArray, partonOutputArray,
401 &readStopWatch, &procStopWatch);
402 modularDelphes->ProcessTask();
403 procStopWatch.Stop();
404
405 if(reader)
406 {
407 reader->AnalyzeEvent(branchEventLHEF, eventCounter, &readStopWatch, &procStopWatch);
408 reader->AnalyzeWeight(branchWeightLHEF);
409 }
410
411 treeWriter->Fill();
412
413 treeWriter->Clear();
414 modularDelphes->Clear();
415 if(reader) reader->Clear();
416
417 readStopWatch.Start();
418 progressBar.Update(eventCounter, eventCounter);
419 }
420
421 progressBar.Update(eventCounter, eventCounter, kTRUE);
422 progressBar.Finish();
423
424 pythia->stat();
425
426 modularDelphes->FinishTask();
427 treeWriter->Write();
428
429 cout << "** Exiting..." << endl;
430
431 delete reader;
432 delete pythia;
433 delete modularDelphes;
434 delete confReader;
435 delete treeWriter;
436 delete outputFile;
437
438 return 0;
439 }
440 catch(runtime_error &e)
441 {
442 if(treeWriter) delete treeWriter;
443 if(outputFile) delete outputFile;
444 cerr << "** ERROR: " << e.what() << endl;
445 return 1;
446 }
447}
Note: See TracBrowser for help on using the repository browser.