Belle II Software development
FragmentationModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <generators/modules/fragmentation/FragmentationModule.h>
10
11#include <generators/evtgen/EvtGenInterface.h>
12#include <generators/utilities/GeneratorConst.h>
13
14#include <framework/gearbox/Unit.h>
15#include <framework/gearbox/Const.h>
16#include <framework/particledb/EvtGenDatabasePDG.h>
17#include <framework/utilities/FileSystem.h>
18
19#include <TRandom3.h>
20
21#include <mdst/dataobjects/MCParticleGraph.h>
22
23#include <framework/logging/Logger.h>
24#include <framework/utilities/IOIntercept.h>
25
26#include <string>
27
28using namespace std;
29using namespace Belle2;
30using namespace Pythia8;
31
32//-----------------------------------------------------------------
33// Register the Module
34//-----------------------------------------------------------------
35REG_MODULE(Fragmentation);
36
37//-----------------------------------------------------------------
38// Implementation
39//-----------------------------------------------------------------
41{
42 //Set module properties
43 setDescription("Fragmention of (u/d/s/c) quarks using PYTHIA8");
44
45 //Parameter definition
46 addParam("ParameterFile", m_parameterfile, "Input parameter file for PYTHIA",
47 FileSystem::findFile("generators/modules/fragmentation/data/pythia_belle2.dat"));
48 addParam("ListPYTHIAEvent", m_listEvent, "List event record of PYTHIA after hadronization", 0);
49 addParam("UseEvtGen", m_useEvtGen, "Use EvtGen for specific decays", 1);
50 addParam("DecFile", m_DecFile, "EvtGen decay file (DECAY.DEC)",
51 FileSystem::findFile("decfiles/dec/DECAY_BELLE2.DEC", true));
52 addParam("UserDecFile", m_UserDecFile, "User EvtGen decay file", std::string(""));
53 addParam("CoherentMixing", m_coherentMixing, "Decay the B0-B0bar coherently (should always be true)", true);
54 addParam("QuarkPairMotherParticle", m_quarkPairMotherParticle,
55 "PDG Code of the mother particle of the quark pair (default Z0 boson)", 23);
56 addParam("AdditionalPDGCodes", m_additionalPDGCodes, "Additional particles used by Pythia", {});
57
58 //initialize member variables
59 evtgen = 0;
60 nAdded = 0;
61 nQuarks = 0;
62 nVpho = 0;
63 nAll = 0;
64 nGood = 0;
65}
66
67
72
74{
75
76 // print internal pythia error statistics
78 statLogCapture.start();
79 m_Pythia->stat();
80 statLogCapture.finish();
81
82 if (nAll != nGood) {
83 double ratio = 0.; //ratio of good over all events
84 if (nAll) ratio = 100.0 * nGood / nAll;
85
86 B2WARNING("Not all events could be fragmented: " << nAll - nGood << " events failed.");
87 B2WARNING("Total number of events: " << nAll << ", of these fragmented: " << nGood << ", success-ratio (should be >97%): " << ratio
88 << "%");
89 B2WARNING("Please contact the generator librarian if the success ratio is below 97%.");
90 B2WARNING("Please treat the success-ratio as correction of the effective cross section due to unphysical events.");
91 }
92}
93
94//-----------------------------------------------------------------
95// Initialize
96//-----------------------------------------------------------------
98{
99 m_mcparticles.isRequired(m_particleList);
100
101 B2DEBUG(150, "Initialize PYTHIA8");
102
103 // Generator and the shorthand m_PythiaEvent = pythia->event are declared in .h file
104 // A simple way to collect all the changes is to store the parameter values in a separate file,
105 // with one line per change. This should be done between the creation of the Pythia object
106 // and the init call for it.
107
109 initLogCapture.start();
110
111 m_Pythia = new Pythia;
112 m_PythiaEvent = &m_Pythia->event;
113 (*m_PythiaEvent) = 0;
114
115 // Load EvtGen particle data.
118
119 // Switch off ProcessLevel
120 m_Pythia->readString("ProcessLevel:all = off");
121
122 // Read the PYTHIA input file, overrides parameters
123 if (!m_Pythia->readFile(m_parameterfile))
124 B2FATAL("Cannot read Pythia parameter file.");
125
126 // Set framework generator
127 FragmentationRndm* fragRndm = new FragmentationRndm();
128 m_Pythia->setRndmEnginePtr(fragRndm);
129
130 // Initialize PYTHIA
131 m_Pythia->init();
132
133 // Set EvtGen (after m_Pythia->init())
134 evtgen = 0;
135
136 if (m_useEvtGen) {
137 B2INFO("Using PYTHIA EvtGen Interface");
138 const std::string defaultDecFile = FileSystem::findFile("decfiles/dec/DECAY_BELLE2.DEC", true);
139 if (m_DecFile.empty()) {
140 B2ERROR("No global decay file defined, please make sure the parameter 'DecFile' is set correctly");
141 return;
142 }
143 if (defaultDecFile.empty()) {
144 B2WARNING("Cannot find default decay file");
145 } else if (defaultDecFile != m_DecFile) {
146 B2INFO("Using non-standard DECAY file \"" << m_DecFile << "\"");
147 }
148 evtgen = new EvtGenDecays(m_Pythia, evtGen);
149 // Set signal particle suffix.
150 // Signal particle names must end in "sig" in user decay file.
151 evtgen->signalSuffix = "sig";
152 // Read the user decay file.
153 evtgen->readDecayFile(m_UserDecFile);
154 // Workaround for Pythia bug. It is the only way to call
155 // EvtGenDecays::updateData(true) to disable particle decays
156 // for all particles from EvtGen decay table. Thus, EvtGen generation
157 // has to be called before generation of the first Pythia event.
158 // Since Pythia event is currently empty, it actually only updates
159 // the particle properties (exactly what is necessary).
160 evtgen->decay();
161 }
162
163 // List variable(s) that differ from their defaults
164 m_Pythia->settings.listChanged();
165
166 initLogCapture.finish();
167}
168
169//-----------------------------------------------------------------
170// Event
171//-----------------------------------------------------------------
173{
174 // Reset the indices of the graph
175 mcParticleGraph.clear();
177
178 // Reset PYTHIA event record to allow for new event
179 IOIntercept::OutputToLogMessages resetLogCapture("EvtGen", LogConfig::c_Debug, LogConfig::c_Error, 100, 100);
180 resetLogCapture.start();
181 m_PythiaEvent->reset();
182 resetLogCapture.finish();
183
184 // Reset counter for added quarks and vphos
185 nAdded = 0;
186 nQuarks = 0;
187 nVpho = 0;
188
189 // Store which MCParticle index belongs to which Pythia index
190 std::map<int, int> indexPYTHIA;
191
192 // Store which Pythia index belongs to which MCParticle index
193 std::map<int, int> indexMCGraph;
194
195 // Store position of the quark (mother of hadronized final state)
196 int quarkPosition = 0;
197
198 // Loop over all particles to find the quark pair
199 int nPart = m_mcparticles.getEntries();
200 for (int iPart = 0; iPart < nPart; iPart++) {
201 MCParticle* currParticle = m_mcparticles[iPart];
202
203 //returns quark id if it finds a quark, zero otherwise
204 //increments a counter for the number of found quarks
205 int pythiaIndex = addParticleToPYTHIA(*currParticle);
206
207 if (pythiaIndex != 0) {
208 indexPYTHIA[nAdded] = iPart;
209 if (pythiaIndex > 0) quarkPosition = iPart;
210 }
211 }
212
213 // Check needed if virtual exchange boson and two quarks are present
214 if (nQuarks != 2) {
215 B2FATAL("Invalid number of quarks: " << nQuarks << " (should be 2)!");
216 }
217
218 if (nVpho != 1) {
219 B2WARNING("No virtual exchange particle given, no PYTHIA FSR in Decays");
220 } else {
221 // Adding QCD and QED FSR
222 m_Pythia->forceTimeShower(2, 3, 20.00);
223 }
224
225 // Check needed if event is energetically possible
226 // ...
227
228 // Do the fragmentation using PYTHIA
229 setReturnValue(1); //return value is 1...
230 nAll = nAll + 1;
231
233 eventLogCapture.start();
234 int success = m_Pythia->next();
235 eventLogCapture.finish();
236
237 if (!success) {
239 listLogCapture.start();
240 m_PythiaEvent->list();
241 listLogCapture.finish();
242
243 setReturnValue(-1); //return value becomes -1 if trials were not successful
244 } else {
245 nGood = nGood + 1;
246 }
247
248 // use evtgen to perform the decay
249 if (m_useEvtGen) {
251 decayLogCapture.start();
252 evtgen->decay();
253 decayLogCapture.finish();
254 }
255
256 // Loop over the PYTHIA list and assign the mother-daughter relation
257 // Might not work if the mother appear below the daughter in the event record
258 for (int iPythiaPart = 0; iPythiaPart < m_Pythia->event.size(); ++iPythiaPart) {
259 auto oldindex = indexPYTHIA.find(iPythiaPart);
260
261 //skip "system" particle generated by PYTHIA
262 if (m_Pythia->event[iPythiaPart].id() == 90) continue;
263
264 if (oldindex == end(indexPYTHIA)) {
265 // --> new particle
266
267 // Add to particle grapg
268 int position = mcParticleGraph.size();
269 mcParticleGraph.addParticle();
270 indexMCGraph[iPythiaPart] = position;
271
273
274 // from PYTHIA manual: "<1: an empty entry, with no meaningful information and
275 // therefore to be skipped unconditionally (should not occur in PYTHIA)"
276 if (m_Pythia->event[iPythiaPart].statusHepMC() < 1) continue;
277
278 // Set PDG code
279 p->setPDG(m_Pythia->event[iPythiaPart].id());
280
281 // Set four vector
282 ROOT::Math::PxPyPzEVector p4(m_Pythia->event[iPythiaPart].px(), m_Pythia->event[iPythiaPart].py(),
283 m_Pythia->event[iPythiaPart].pz(),
284 m_Pythia->event[iPythiaPart].e());
285 p->set4Vector(p4);
286 p->setMass(m_Pythia->event[iPythiaPart].m());
287
288 // Set vertex
289 p->setProductionVertex(m_Pythia->event[iPythiaPart].xProd() * Unit::mm, m_Pythia->event[iPythiaPart].yProd() * Unit::mm,
290 m_Pythia->event[iPythiaPart].zProd() * Unit::mm);
291 p->setProductionTime(m_Pythia->event[iPythiaPart].tProd() * Unit::mm / Const::speedOfLight);
292 p->setValidVertex(true);
293
294 // Set all(!) particles from the generator to primary
296
297 // Set FSR flag from PYTHIA TimeShower:QEDshowerByQ
298 if (m_Pythia->event[iPythiaPart].status() == 51 && m_Pythia->event[iPythiaPart].id() == 22) {
300 }
301
302 // Set PHOTOS flag from PYTHIA-EvtGen
303 if (m_Pythia->event[iPythiaPart].status() == GeneratorConst::FSR_STATUS_CODE && m_Pythia->event[iPythiaPart].id() == 22) {
305 }
306
307 // Set stable at generator level
308 if (m_Pythia->event[iPythiaPart].statusHepMC() == 1) {
310 }
311
312 //set mother
313 const int motherid = m_Pythia->event[iPythiaPart].mother1();
314
315 //check if mother exists in indexMCGraph
316 auto motherindex = indexMCGraph.find(motherid);
317
318 if (motherindex != end(indexMCGraph)) {
319 int motheridingraph = indexMCGraph[motherid];
320 MCParticleGraph::GraphParticle* q = &mcParticleGraph[motheridingraph];
321 p->comesFrom(*q);
322 } else {
323 // Particle has no mother from PYTHIA, add quark as mother
325 p->comesFrom(*q);
326 }
327
328 } else {
329 // particle is already in the graph
330 // modify here if needed
331 // ...
332 }
333 }
334
335 // Print original PYTHIA list
336 if (m_listEvent) m_PythiaEvent->list();
337
338 // Create new ParticleGraph
339 mcParticleGraph.generateList(m_particleList,
341
342}
343
344//-----------------------------------------------------------------
345// addParticleToPYTHIA
346//-----------------------------------------------------------------
348{
349 //get PDG code
350 const int id = mcParticle.getPDG();
351
352 //check that this particle is a quark or a virtual photon(-like)
353 bool isVPho = false;
354 bool isQuark = false;
355 if (abs(id) >= 1 && abs(id) <= 5) isQuark = true;
356 if (id == m_quarkPairMotherParticle) isVPho = true;
357
358 if (!(isVPho || isQuark)) return 0;
359
360 //check that there is no daughter for the quarks
361 if (isQuark && mcParticle.getDaughters().size()) {
362 B2WARNING("Quark already has a daughter!");
363 return 0;
364 }
365
366 //get some basic kinematics
367 const double mass = mcParticle.getMass();
368 const ROOT::Math::XYZVector& p = mcParticle.getMomentum();
369 const double energy = sqrt(mass * mass + p.Mag2());
370
371 //add this (anti)quark to the m_PythiaEvent
372 if (id == m_quarkPairMotherParticle) {
373 m_PythiaEvent->append(m_quarkPairMotherParticle, -22, 0, 0, 2, 3, 0, 0, p.X(), p.Y(), p.Z(), energy, mass);
374 nVpho++;
375 } else if (id > 0) {
376 m_PythiaEvent->append(id, 23, 1, 0, 0, 0, 101, 0, p.X(), p.Y(), p.Z(), energy, mass, 20.0);
377 nQuarks++;
378 } else if (id < 0) {
379 m_PythiaEvent->append(id, 23, 1, 0, 0, 0, 0, 101, p.X(), p.Y(), p.Z(), energy, mass, 20.0);
380 nQuarks++;
381 }
382
383 nAdded++;
384
385 return id;
386}
387
389{
390 Pythia8::ParticleData* particleData = &pythia->particleData;
391 int nParticles = EvtPDL::entries();
392 for (int i = 0; i < nParticles; ++i) {
393 EvtId evtgenParticle = EvtPDL::getEntry(i);
394 /*
395 * Pythia uses absolute value of the PDG code, thus, only positive codes
396 * are necessary. Only leptons (11 <= pdg <= 20) and hadrons (pdg > 100)
397 * are updated.
398 */
399 int pdg = EvtPDL::getStdHep(evtgenParticle);
400 if (pdg <= 10 || (pdg > 20 && pdg <= 100))
401 continue;
402 EvtId evtgenAntiParticle = EvtPDL::chargeConj(evtgenParticle);
403 if (particleData->isParticle(pdg)) {
404 particleData->setAll(pdg,
405 EvtPDL::name(evtgenParticle),
406 EvtPDL::name(evtgenAntiParticle),
407 EvtPDL::getSpinType(evtgenParticle),
408 EvtPDL::chg3(evtgenParticle),
409 // colType == 0 for uncolored particles.
410 0,
411 EvtPDL::getMeanMass(evtgenParticle),
412 EvtPDL::getWidth(evtgenParticle),
413 EvtPDL::getMinMass(evtgenParticle),
414 EvtPDL::getMaxMass(evtgenParticle),
415 EvtPDL::getctau(evtgenParticle));
416 } else if (std::find(m_additionalPDGCodes.begin(), m_additionalPDGCodes.end(), pdg) != m_additionalPDGCodes.end()) {
417 particleData->addParticle(
418 pdg,
419 EvtPDL::name(evtgenParticle),
420 EvtPDL::name(evtgenAntiParticle),
421 EvtPDL::getSpinType(evtgenParticle),
422 EvtPDL::chg3(evtgenParticle),
423 // colType == 0 for uncolored particles.
424 0,
425 EvtPDL::getMeanMass(evtgenParticle),
426 EvtPDL::getWidth(evtgenParticle),
427 EvtPDL::getMinMass(evtgenParticle),
428 EvtPDL::getMaxMass(evtgenParticle),
429 EvtPDL::getctau(evtgenParticle));
430 }
431 }
432}
433
434//-----------------------------------------------------------------
435// random generator for PYTHIA
436//-----------------------------------------------------------------
437FragmentationRndm::FragmentationRndm() : Pythia8::RndmEngine()
438{
439
440}
441
443{
444 double value = gRandom->Rndm();
445 return value;
446}
static const double speedOfLight
[cm/ns]
Definition Const.h:695
static EvtGen * createEvtGen(const std::string &decayFileName, bool coherentMixing)
Create and initialize an EvtGen instance:
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
EvtGenDecays * evtgen
EvtGen decay engine inside PYTHIA8.
int nGood
number of events with successful fragmentation.
void loadEvtGenParticleData(Pythia8::Pythia *pythia)
Load EvtGen particle data.
bool m_coherentMixing
decay the B0-B0bar coherently.
int nVpho
number of virtual exchange particles.
virtual void initialize() override
Initialize the module.
virtual void event() override
Event method (process events)
Pythia8::Pythia * m_Pythia
Pythia generator.
StoreArray< MCParticle > m_mcparticles
store array for the MCParticles
virtual ~FragmentationModule()
Destructor.
virtual void terminate() override
terminate the module
Pythia8::Event * m_PythiaEvent
Pythia event.
int addParticleToPYTHIA(const MCParticle &mcParticle)
Add particle to Pythia event.
std::string m_DecFile
EvtGen decay file.
int m_useEvtGen
use EvtGen for some decays.
int nAll
number of events created.
std::string m_UserDecFile
User EvtGen decay file.
std::string m_particleList
The name of the MCParticle collection.
MCParticleGraph mcParticleGraph
An instance of the MCParticle graph.
std::string m_parameterfile
Module parameters.
std::vector< int > m_additionalPDGCodes
Additional particles used in Pythia.
int nAdded
number of added particles.
int m_quarkPairMotherParticle
PDG Code of the mother particle of the quark pair.
int m_listEvent
list event generated by PYTHIA.
Minimal class for external random generator to be used in PYTHIA.
double flat()
flat random generator.
bool start()
Start intercepting the output.
Capture stdout and stderr and convert into log messages.
bool finish()
Finish the capture and emit the message if output has appeard on stdout or stderr.
@ c_Error
Error: for things that went wrong and have to be fixed.
Definition LogConfig.h:30
@ c_Debug
Debug: for code development.
Definition LogConfig.h:26
@ c_Warning
Warning: for potential problems that the user should pay attention to.
Definition LogConfig.h:29
Class to represent Particle data in graph.
@ c_checkCyclic
Check for cyclic dependencies.
@ c_clearParticles
Clear the particle list before adding the graph.
@ c_setDecayInfo
Set decay time and vertex.
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
@ c_IsFSRPhoton
bit 7: Particle is from final state radiation
Definition MCParticle.h:61
@ c_IsPHOTOSPhoton
bit 8: Particle is an radiative photon from PHOTOS
Definition MCParticle.h:63
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition MCParticle.h:49
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition MCParticle.cc:52
float getMass() const
Return the particle mass in GeV.
Definition MCParticle.h:124
int getPDG() const
Return PDG code of particle.
Definition MCParticle.h:101
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition MCParticle.h:187
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Module()
Constructor.
Definition Module.cc:30
void setReturnValue(int value)
Sets the return value for this module as integer.
Definition Module.cc:220
static const double mm
[millimeters]
Definition Unit.h:70
A class to perform decays via the external EvtGen decay program, see http://evtgen....
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.