Fork me on GitHub

source: git/readers/DelphesPythia8.cpp@ e576ef60

ImprovedOutputFile Timing dual_readout llp
Last change on this file since e576ef60 was 6934875, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

changed into flat generation in E, included neutron, charged pions, c-jets in particle gun

  • Property mode set to 100644
File size: 12.5 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// Single-particle gun. The particle must be a colour singlet.
150// Input: flavour, energy, direction (theta, phi).
151// If theta < 0 then random choice over solid angle.
152// Optional final argument to put particle at rest => E = m.
153// from pythia8 example 21
154
155void fillParticle(int id, double p_max, double eta_max,
156 Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
157{
158 // Reset event record to allow for new event.
159 event.reset();
160
161 // Generate uniform pt and eta.
162 double pt, eta, phi, pp, ee, mm;
163
164 //pmin = 0.1 GeV for single particles
165 pp = pow(10, - 1.0 + (log10(p_max) + 1.0) * rndm.flat());
166 eta = (2.0 * rndm.flat() - 1.0) * eta_max;
167 phi = 2.0 * M_PI * rndm.flat();
168 mm = pdt.mSel(id);
169 ee = Pythia8::sqrtpos(pp*pp + mm*mm);
170 pt = pp / cosh(eta);
171
172 // Store the particle in the event record.
173 event.append(id, 1, 0, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
174}
175
176void fillPartons(int id, double p_max, double eta_max,
177 Pythia8::Event &event, Pythia8::ParticleData &pdt, Pythia8::Rndm &rndm)
178{
179
180 // Reset event record to allow for new event.
181 event.reset();
182
183 // Generate uniform pt and eta.
184 double pt, eta, phi, pp, ee, mm;
185
186 //pmin = 1 GeV for jets
187 pp = pow(10, log10(p_max) * rndm.flat());
188 eta = (2.0 * rndm.flat() - 1.0) * eta_max;
189 phi = 2.0 * M_PI * rndm.flat();
190 mm = pdt.mSel(id);
191 ee = Pythia8::sqrtpos(pp*pp + mm*mm);
192 pt = pp / cosh(eta);
193
194 if( (id == 4 || id == 5) && pt < 10.0) return;
195
196 if(id == 21)
197 {
198 event.append(21, 23, 101, 102, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee);
199 event.append(21, 23, 102, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee);
200 }
201 else
202 {
203 event.append(id, 23, 101, 0, pt * cos(phi), pt * sin(phi), pt * sinh(eta), ee, mm);
204 event.append(-id, 23, 0, 101, -pt * cos(phi), -pt * sin(phi), -pt * sinh(eta), ee, mm);
205 }
206}
207
208
209//---------------------------------------------------------------------------
210
211int main(int argc, char *argv[])
212{
213 char appName[] = "DelphesPythia8";
214 stringstream message;
215 FILE *inputFile = 0;
216 TFile *outputFile = 0;
217 TStopwatch readStopWatch, procStopWatch;
218 ExRootTreeWriter *treeWriter = 0;
219 ExRootTreeBranch *branchEvent = 0;
220 ExRootTreeBranch *branchEventLHEF = 0, *branchWeightLHEF = 0;
221 ExRootConfReader *confReader = 0;
222 Delphes *modularDelphes = 0;
223 DelphesFactory *factory = 0;
224 TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
225 TObjArray *stableParticleOutputArrayLHEF = 0, *allParticleOutputArrayLHEF = 0, *partonOutputArrayLHEF = 0;
226 DelphesLHEFReader *reader = 0;
227 Long64_t eventCounter, errorCounter;
228 Long64_t numberOfEvents, timesAllowErrors;
229
230 Pythia8::Pythia *pythia = 0;
231
232 // for matching
233 Pythia8::CombineMatchingInput *combined = 0;
234 Pythia8::UserHooks* matching = 0;
235
236 if(argc != 4)
237 {
238 cout << " Usage: " << appName << " config_file" << " pythia_card" << " output_file" << endl;
239 cout << " config_file - configuration file in Tcl format," << endl;
240 cout << " pythia_card - Pythia8 configuration file," << endl;
241 cout << " output_file - output file in ROOT format." << endl;
242 return 1;
243 }
244
245 signal(SIGINT, SignalHandler);
246
247 gROOT->SetBatch();
248
249 int appargc = 1;
250 char *appargv[] = {appName};
251 TApplication app(appName, &appargc, appargv);
252
253 try
254 {
255 outputFile = TFile::Open(argv[3], "CREATE");
256
257 if(outputFile == NULL)
258 {
259 message << "can't create output file " << argv[3];
260 throw runtime_error(message.str());
261 }
262
263 treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
264
265 branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
266
267 confReader = new ExRootConfReader;
268 confReader->ReadFile(argv[1]);
269
270 modularDelphes = new Delphes("Delphes");
271 modularDelphes->SetConfReader(confReader);
272 modularDelphes->SetTreeWriter(treeWriter);
273
274 factory = modularDelphes->GetFactory();
275 allParticleOutputArray = modularDelphes->ExportArray("allParticles");
276 stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
277 partonOutputArray = modularDelphes->ExportArray("partons");
278
279 // Initialize Pythia
280 pythia = new Pythia8::Pythia;
281
282 // jet matching
283 matching = combined->getHook(*pythia);
284 if (!matching)
285 {
286 throw runtime_error("can't do matching");
287 }
288 pythia->setUserHooksPtr(matching);
289
290
291 if(pythia == NULL)
292 {
293 throw runtime_error("can't create Pythia instance");
294 }
295
296 // Read in commands from configuration file
297 if(!pythia->readFile(argv[2]))
298 {
299 message << "can't read Pythia8 configuration file " << argv[2] << endl;
300 throw runtime_error(message.str());
301 }
302
303 // Extract settings to be used in the main program
304 numberOfEvents = pythia->mode("Main:numberOfEvents");
305 timesAllowErrors = pythia->mode("Main:timesAllowErrors");
306
307 // Check if particle gun
308 if (!pythia->flag("Main:spareFlag1"))
309 {
310 inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
311 if(inputFile)
312 {
313 reader = new DelphesLHEFReader;
314 reader->SetInputFile(inputFile);
315
316 branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
317 branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
318
319 allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
320 stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
321 partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
322 }
323 }
324
325 modularDelphes->InitTask();
326
327 pythia->init();
328
329 // ExRootProgressBar progressBar(numberOfEvents - 1);
330 ExRootProgressBar progressBar(-1);
331
332 // Loop over all events
333 errorCounter = 0;
334 treeWriter->Clear();
335 modularDelphes->Clear();
336 readStopWatch.Start();
337 for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
338 {
339 while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
340 stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
341
342 if (pythia->flag("Main:spareFlag1"))
343 {
344 if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 15 || pythia->mode("Main:spareMode1") == 22 || pythia->mode("Main:spareMode1") == 211 || pythia->mode("Main:spareMode1") == 2112)
345 {
346 fillParticle(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
347 }
348 else
349 {
350 fillPartons(pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->parm("Main:spareParm2"), pythia->event, pythia->particleData, pythia->rndm);
351 }
352 }
353
354 if(!pythia->next())
355 {
356 // If failure because reached end of file then exit event loop
357 if(pythia->info.atEndOfFile())
358 {
359 cerr << "Aborted since reached end of Les Houches Event File" << endl;
360 break;
361 }
362
363 // First few failures write off as "acceptable" errors, then quit
364 if(++errorCounter > timesAllowErrors)
365 {
366 cerr << "Event generation aborted prematurely, owing to error!" << endl;
367 break;
368 }
369
370 modularDelphes->Clear();
371 reader->Clear();
372 continue;
373 }
374
375 readStopWatch.Stop();
376
377 procStopWatch.Start();
378 ConvertInput(eventCounter, pythia, branchEvent, factory,
379 allParticleOutputArray, stableParticleOutputArray, partonOutputArray,
380 &readStopWatch, &procStopWatch);
381 modularDelphes->ProcessTask();
382 procStopWatch.Stop();
383
384 if(reader)
385 {
386 reader->AnalyzeEvent(branchEventLHEF, eventCounter, &readStopWatch, &procStopWatch);
387 reader->AnalyzeWeight(branchWeightLHEF);
388 }
389
390 treeWriter->Fill();
391
392 treeWriter->Clear();
393 modularDelphes->Clear();
394 if(reader) reader->Clear();
395
396 readStopWatch.Start();
397 progressBar.Update(eventCounter, eventCounter);
398 }
399
400 progressBar.Update(eventCounter, eventCounter, kTRUE);
401 progressBar.Finish();
402
403 pythia->stat();
404
405 modularDelphes->FinishTask();
406 treeWriter->Write();
407
408 cout << "** Exiting..." << endl;
409
410 delete reader;
411 delete pythia;
412 delete modularDelphes;
413 delete confReader;
414 delete treeWriter;
415 delete outputFile;
416
417 return 0;
418 }
419 catch(runtime_error &e)
420 {
421 if(treeWriter) delete treeWriter;
422 if(outputFile) delete outputFile;
423 cerr << "** ERROR: " << e.what() << endl;
424 return 1;
425 }
426}
Note: See TracBrowser for help on using the repository browser.