Fork me on GitHub

source: git/readers/DelphesPythia8.cpp@ bdf9e2e

Last change on this file since bdf9e2e was 78d1846, checked in by Pavel Demin <pavel-demin@…>, 8 years ago

use fillParticle for all particles that aren't partons

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