Fork me on GitHub

source: git/readers/DelphesPythia8.cpp@ 93bf9c6

Timing
Last change on this file since 93bf9c6 was 2c81caa, checked in by Michele Selvaggi <michele.selvaggi@…>, 5 years ago

adapted pythia reader to rhadrons

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