Fork me on GitHub

source: git/readers/DelphesPythia8.cpp@ 19389b8

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

added matching in Pythia8

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