Belle II Software  release-05-01-25
EvtGenDecays.h
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * This file is a copy of Pythia8Plugins/EvtGen.h from Pythia 8 but *
4  * slightly modified for extended evtgen models *
5  * *
6  * Adapted by: Martin Ritter *
7  * *
8  **************************************************************************/
9 
10 // EvtGen.h is a part of the PYTHIA event generator.
11 // Copyright (C) 2016 Torbjorn Sjostrand.
12 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
13 // Please respect the MCnet Guidelines, see GUIDELINES for details.
14 // Author: Philip Ilten.
15 
16 // This file contains an EvtGen interface. HepMC and EvtGen must be enabled.
17 
18 #pragma once
19 
20 #include <generators/utilities/GeneratorConst.h>
21 #include "Pythia8/Pythia.h"
22 #include "EvtGen/EvtGen.hh"
23 #include "EvtGenBase/EvtRandomEngine.hh"
24 #include "EvtGenBase/EvtParticle.hh"
25 #include "EvtGenBase/EvtParticleFactory.hh"
26 #include "EvtGenBase/EvtPDL.hh"
27 #include "EvtGenBase/EvtDecayTable.hh"
28 #include "EvtGenBase/EvtParticleDecayList.hh"
29 #include "EvtGenBase/EvtDecayBase.hh"
30 #include "EvtGenExternal/EvtExternalGenList.hh"
31 
32 //==========================================================================
33 
35 class EvtGenRandom : public EvtRandomEngine {
36 
37 public:
38 
40  explicit EvtGenRandom(Pythia8::Rndm* rndmPtrIn) {rndmPtr = rndmPtrIn;}
41 
43  double random() {if (rndmPtr) return rndmPtr->flat(); else return -1.0;}
44 
46  Pythia8::Rndm* rndmPtr;
47 
48 };
49 
50 //==========================================================================
51 
94 class EvtGenDecays {
95 
96 public:
97 
99  EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, EvtGen* evtGen, bool limit = true);
100 
102  EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, std::string decayFile, std::string particleDataFile,
103  EvtExternalGenList* extPtrIn = 0, EvtAbsRadCorr* fsrPtrIn = 0,
104  int mixing = 1, bool xml = false, bool limit = true,
105  bool extUse = true, bool fsrUse = true);
106 
109  {
110  if (evtgen) delete evtgen;
111  if (extOwner && extPtr) delete extPtr;
112  if (fsrOwner && fsrPtr) delete fsrPtr;
113  }
114 
116  EvtGenDecays(const EvtGenDecays&) = delete;
117 
119  EvtGenDecays& operator=(const EvtGenDecays&) = delete;
120 
122  void getDecayLimits(bool limit);
123 
125  double decay();
126 
128  void exclude(int id) {excIds.insert(id);}
129 
131  void updatePythia();
132 
134  void updateEvtGen();
135 
137  void readDecayFile(std::string decayFile, bool xml = false)
138  {
139  evtgen->readUDecay(decayFile.c_str(), xml); updateData();
140  }
141 
142  bool extOwner;
143  bool fsrOwner;
144  EvtExternalGenList* extPtr;
145  EvtAbsRadCorr* fsrPtr;
146  std::list<EvtDecayBase*> models;
149  struct Signal {
150  int status;
151  EvtId egId;
152  std::vector<double> bfs;
153  std::vector<int> map;
154  EvtParticleDecayList modes;
155  };
157  std::map<int, Signal> signals;
158 
160  std::string signalSuffix;
161 
162 protected:
163 
165  void updateData(bool final = false);
166 
168  void updateEvent(Pythia8::Particle* pyPro, EvtParticle* egPro,
169  std::vector<int>* pySigs = 0, std::vector<EvtParticle*>* egSigs = 0,
170  std::vector<double>* bfs = 0, double* wgt = 0);
171 
173  bool checkVertex(Pythia8::Particle* pyPro);
174 
176  bool checkSignal(Pythia8::Particle* pyPro);
177 
179  bool checkOsc(EvtParticle* egPro);
180 
182  static const int NTRYDECAY = 1000;
183 
185  Pythia8::Pythia* pythiaPtr;
186 
189 
191  EvtGen* evtgen;
192 
193  std::set<int> incIds;
194  std::set<int> excIds;
197  bool updated;
198 
200  std::map<int, Signal>::iterator signal;
201 
206  double tau0Max;
207  double tauMax;
208  double rMax;
209  double xyMax;
210  double zMax;
211  bool limitTau0;
212  bool limitTau;
213  bool limitRadius;
215  bool limitDecay;
216 };
217 
218 //--------------------------------------------------------------------------
219 
220 // The constructor.
221 
222 // The EvtGenDecays object is associated with a single Pythia
223 // instance. This is to ensure a consistent random number generator
224 // across the two, as well as any updates to particle data, etc. Note
225 // that if multiple EvtGenDecays objects exist, that they will modify
226 // one anothers particle databases due to the design of EvtGen.
227 
228 // This constructor also sets all particles to be decayed by EvtGen as
229 // stable within Pythia. The parameters within Pythia used to check if
230 // a particle should be decayed, as described in the "Particle Decays"
231 // section of the Pythia manual, are set. Note that if the variable
232 // "limit" is set to "false", then no check will be made before
233 // decaying a particle with EvtGen.
234 
235 // The constructor is designed to have the exact same form as the
236 // EvtGen constructor except for these five differences.
237 // (1) The first variable is the pointer to the Pythia object.
238 // (2) The third last argument is a flag to limit decays based on the
239 // Pythia criteria (based on the particle decay vertex).
240 // (3) The second last argument is a flag if external models should be
241 // passed to EvtGen (default is true).
242 // (4) The last argument is a flag if an FSR model should be passed
243 // to EvtGen (default is true).
244 // (5) No random engine pointer is passed, as this is obtained from
245 // Pythia.
246 
247 // pythiaPtrIn: the pointer to the associated Pythia object.
248 // decayFile: the name of the decay file to pass to EvtGen.
249 // particleDataFile: the name of the particle data file to pass to EvtGen.
250 // extPtrIn: the optional EvtExternalGenList pointer, this must be
251 // be provided if fsrPtrIn is provided to avoid double
252 // initializations.
253 // fsrPtrIn: the EvtAbsRadCorr pointer to pass to EvtGen.
254 // mixing: the mixing type to pass to EvtGen.
255 // xml: flag to use XML files to pass to EvtGen.
256 // limit: flag to limit particle decays based on Pythia criteria.
257 // extUse: flag to use external models with EvtGen.
258 // fsrUse: flag to use radiative correction engine with EvtGen.
259 
260 
261 EvtGenDecays::EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, EvtGen* evtGen, bool limit):
262  extOwner(false), fsrOwner(false), extPtr(nullptr), fsrPtr(nullptr),
263  signalSuffix("_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(nullptr),
264  evtgen(evtGen), updated(false)
265 {
266  getDecayLimits(limit);
267 }
268 
269 EvtGenDecays::EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, std::string decayFile,
270  std::string particleDataFile, EvtExternalGenList* extPtrIn,
271  EvtAbsRadCorr* fsrPtrIn, int mixing, bool xml, bool limit,
272  bool extUse, bool fsrUse) : extPtr(extPtrIn), fsrPtr(fsrPtrIn),
273  signalSuffix("_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(&pythiaPtr->rndm),
274  updated(false)
275 {
276 
277  // Initialize EvtGen.
278  if (!extPtr && fsrPtr) {
279  if (pythiaPtr) pythiaPtr->info.errorMsg("Error in EvtGenDecays::"
280  "EvtGenDecays: extPtr is null but fsrPtr is provided");
281  return;
282  }
283  if (extPtr) extOwner = false;
284  else {extOwner = true; extPtr = new EvtExternalGenList();}
285  if (fsrPtr) fsrOwner = false;
286  else {fsrOwner = true; fsrPtr = extPtr->getPhotosModel();}
287  models = extPtr->getListOfModels();
288  evtgen = new EvtGen(decayFile.c_str(), particleDataFile.c_str(),
289  &rndm, fsrUse ? fsrPtr : 0, extUse ? &models : 0, mixing, xml);
290  getDecayLimits(limit);
291 }
292 
294 {
295  // Get the Pythia decay limits.
296  if (!pythiaPtr) return;
297  limitTau0 = pythiaPtr->settings.flag("ParticleDecays:limitTau0");
298  tau0Max = pythiaPtr->settings.parm("ParticleDecays:tau0Max");
299  limitTau = pythiaPtr->settings.flag("ParticleDecays:limitTau");
300  tauMax = pythiaPtr->settings.parm("ParticleDecays:tauMax");
301  limitRadius = pythiaPtr->settings.flag("ParticleDecays:limitRadius");
302  rMax = pythiaPtr->settings.parm("ParticleDecays:rMax");
303  limitCylinder = pythiaPtr->settings.flag("ParticleDecays:limitCylinder");
304  xyMax = pythiaPtr->settings.parm("ParticleDecays:xyMax");
305  zMax = pythiaPtr->settings.parm("ParticleDecays:zMax");
306  limitDecay = limit && (limitTau0 || limitTau ||
308 
309 }
310 
311 //--------------------------------------------------------------------------
312 
313 // Perform all decays and return the event weight.
314 
315 // All particles in the event record that can be decayed by EvtGen are
316 // decayed. If a particle is a signal particle, then this is stored in
317 // a vector of signal particles. A signal particle is only stored if
318 // its status is the same as the status provided in the signals map. A
319 // negative status in the signal map indicates that all statuses
320 // should be accepted. After all signal particles are identified, one
321 // is randomly chosen and decayed as signal. The remainder are decayed
322 // normally.
323 
324 // Forcing a signal decay changes the weight of an event from unity,
325 // and so the relative event weight is returned, given the forced
326 // signal decay. A weight of 0 indicates no signal in the event, while
327 // a weight of -1 indicates something is wrong, e.g. either the Pythia
328 // or EvtGen pointers are not available or the number of tries has
329 // been exceeded. For the event weight to be valid, one should not
330 // change the absolute branching fractions in the signal and inclusive
331 // definitions, but rather just remove the unwanted decay channels
332 // from the signal decay definition.
333 
335 {
336 
337  // Reset the signal and signal counters.
338  if (!pythiaPtr || !evtgen) return -1;
339  if (!updated) updateData(true);
340 
341  // Loop over all particles in the Pythia event.
342  Pythia8::Event& event = pythiaPtr->event;
343  std::vector<int> pySigs; std::vector<EvtParticle*> egSigs, egPrts;
344  std::vector<double> bfs; double wgt(1.);
345  for (int iPro = 0; iPro < event.size(); ++iPro) {
346 
347  // Check particle is final and can be decayed by EvtGen.
348  Pythia8::Particle* pyPro = &event[iPro];
349  if (!pyPro->isFinal()) continue;
350  if (incIds.find(pyPro->id()) == incIds.end()) continue;
351  if (pyPro->status() == 93 || pyPro->status() == 94) continue;
352 
353  // Decay the progenitor with EvtGen.
354  EvtParticle* egPro = EvtParticleFactory::particleFactory
355  (EvtPDL::evtIdFromStdHep(pyPro->id()),
356  EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz()));
357  egPrts.push_back(egPro);
358  egPro->setDiagonalSpinDensity();
359  evtgen->generateDecay(egPro);
360  pyPro->tau(egPro->getLifetime());
361  if (!checkVertex(pyPro)) continue;
362 
363  // Add oscillations to event record.
364  if (checkOsc(egPro))
365  updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
366 
367  // Undo decay if signal (duplicate to stop oscillations).
368  else if (checkSignal(pyPro)) {
369  pySigs.push_back(pyPro->index());
370  egSigs.push_back(egPro);
371  bfs.push_back(signal->second.bfs[0]);
372  wgt *= 1 - bfs.back();
373  egPro->deleteDaughters();
374  EvtParticle* egDau = EvtParticleFactory::particleFactory
375  (EvtPDL::evtIdFromStdHep(pyPro->id()),
376  EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz()));
377  egDau->addDaug(egPro);
378  egDau->setDiagonalSpinDensity();
379 
380  // If not signal, add to event record.
381  } else updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
382  }
383  if (pySigs.size() == 0) {
384  for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
385  egPrts[iPrt]->deleteTree();
386  return 0;
387  }
388 
389  // Determine the decays of the signal particles (signal or background).
390  std::vector<int> modes; int force(-1);
391  for (int iTry = 1; iTry <= NTRYDECAY; ++iTry) {
392  int n(0);
393 
394  modes.clear(); force = pythiaPtr->rndm.pick(bfs); n = 0;
395  for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
396  if (iSig == force) modes.push_back(0);
397  else modes.push_back(pythiaPtr->rndm.flat() > bfs[iSig]);
398  if (modes.back() == 0) ++n;
399  }
400  if (pythiaPtr->rndm.flat() <= 1.0 / n) break;
401  if (iTry == NTRYDECAY) {
402  wgt = 2.;
403  pythiaPtr->info.errorMsg("Warning in EvtGenDecays::decay: maximum "
404  "number of signal decay attempts exceeded");
405  }
406  }
407 
408  // Decay the signal particles and mark forced decay.
409  for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
410  EvtParticle* egSig = egSigs[iSig];
411  Pythia8::Particle* pySig = &event[pySigs[iSig]];
412  signal = signals.find(pySig->id());
413  if (egSig->getNDaug() > 0) egSig = egSig->getDaug(0);
414  if (modes[iSig] == 0) egSig->setId(signal->second.egId);
415  else {
416  signal->second.modes.getDecayModel(egSig);
417  egSig->setChannel(signal->second.map[egSig->getChannel()]);
418  }
419  if (iSig == force)
420  pySig->status(pySig->status() == 92 || pySig->status() == 94 ? 96 : 95);
421  evtgen->generateDecay(egSig);
422  updateEvent(pySig, egSigs[iSig]);
423  }
424 
425  // Delete all EvtGen particles and return weight.
426  for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
427  egPrts[iPrt]->deleteTree();
428  return 1. - wgt;
429 
430 }
431 
432 //--------------------------------------------------------------------------
433 
434 // Update the Pythia particle database from EvtGen.
435 
436 // Note that only the particle spin type, charge type, nominal mass,
437 // width, minimum mass, maximum mass, and nominal lifetime are
438 // set. The name string is not set.
439 
441 {
442  if (!pythiaPtr || !evtgen) return;
443  for (int entry = 0; entry < (int)EvtPDL::entries(); ++entry) {
444  EvtId egId = EvtPDL::getEntry(entry);
445  int pyId = EvtPDL::getStdHep(egId);
446  pythiaPtr->particleData.spinType(pyId, EvtPDL::getSpinType(egId));
447  pythiaPtr->particleData.chargeType(pyId, EvtPDL::chg3(egId));
448  pythiaPtr->particleData.m0(pyId, EvtPDL::getMass(egId));
449  pythiaPtr->particleData.mWidth(pyId, EvtPDL::getWidth(egId));
450  pythiaPtr->particleData.mMin(pyId, EvtPDL::getMinMass(egId));
451  pythiaPtr->particleData.mMax(pyId, EvtPDL::getMaxMass(egId));
452  pythiaPtr->particleData.tau0(pyId, EvtPDL::getctau(egId));
453  }
454 }
455 
456 //--------------------------------------------------------------------------
457 
458 // Update the EvtGen particle database from Pythia.
459 
460 // The particle update database between Pythia and EvtGen is not
461 // symmetric. Particularly, it is not possible to update the spin
462 // type, charge, or nominal lifetime in the EvtGen particle database.
463 
465 {
466  if (!pythiaPtr || !evtgen) return;
467  int pyId = pythiaPtr->particleData.nextId(0);
468  while (pyId != 0) {
469  EvtId egId = EvtPDL::evtIdFromStdHep(pyId);
470  EvtPDL::reSetMass(egId, pythiaPtr->particleData.m0(pyId));
471  EvtPDL::reSetWidth(egId, pythiaPtr->particleData.mWidth(pyId));
472  EvtPDL::reSetMassMin(egId, pythiaPtr->particleData.mMin(pyId));
473  EvtPDL::reSetMassMax(egId, pythiaPtr->particleData.mMax(pyId));
474  pyId = pythiaPtr->particleData.nextId(pyId);
475  }
476 }
477 
478 //--------------------------------------------------------------------------
479 
480 // Update the particles to decay with EvtGen, and the selected signals.
481 
482 // If final is false, then only signals are initialized in the signal
483 // map. Any particle or alias that ends with signalSuffix is taken as
484 // a signal particle. If final is true all particle entries in EvtGen
485 // are checked to see if they should be set stable in Pythia. If an
486 // EvtGen particle has no decay modes, then Pythia is still allowed to
487 // decay the particle. Additionally, the signal decay channels are
488 // turned off for the non-aliased signal particle.
489 
490 void EvtGenDecays::updateData(bool final)
491 {
492 
493  // Loop over the EvtGen entries.
494  if (!pythiaPtr) return;
495  EvtDecayTable* egTable = EvtDecayTable::getInstance();
496  if (!egTable) return;
497  for (int iEntry = 0; iEntry < (int)EvtPDL::entries(); ++iEntry) {
498  EvtId egId = EvtPDL::getEntry(iEntry);
499  int pyId = EvtPDL::getStdHep(egId);
500  if (egTable->getNModes(egId) == 0) continue;
501  if (excIds.find(pyId) != excIds.end()) continue;
502 
503  // Stop Pythia from decaying the particle and include in decay set.
504  if (final) {
505  incIds.insert(pyId);
506  pythiaPtr->particleData.mayDecay(pyId, false);
507  }
508 
509  // Check for signal.
510  std::string egName = EvtPDL::name(egId);
511  if (egName.size() <= signalSuffix.size() || egName.substr
512  (egName.size() - signalSuffix.size()) != signalSuffix) continue;
513  signal = signals.find(pyId);
514  if (signal == signals.end()) {
515  signals[pyId].status = -1;
516  signal = signals.find(pyId);
517  }
518  signal->second.egId = egId;
519  signal->second.bfs = std::vector<double>(2, 0);
520  if (!final) continue;
521 
522  // Get the signal and background decay modes.
523  std::vector<EvtParticleDecayList> egList = egTable->getDecayTable();
524  int sigIdx = egId.getAlias();
525  int bkgIdx = EvtPDL::evtIdFromStdHep(pyId).getAlias();
526  if (sigIdx > (int)egList.size() || bkgIdx > (int)egList.size()) continue;
527  EvtParticleDecayList sigModes = egList[sigIdx];
528  EvtParticleDecayList bkgModes = egList[bkgIdx];
529  EvtParticleDecayList allModes = egList[bkgIdx];
530  double sum(0);
531 
532  // Sum signal branching fractions.
533  for (int iMode = 0; iMode < sigModes.getNMode(); ++iMode) {
534  EvtDecayBase* mode = sigModes.getDecayModel(iMode);
535  if (!mode) continue;
536  signal->second.bfs[0] += mode->getBranchingFraction();
537  sum += mode->getBranchingFraction();
538  }
539 
540  // Sum remaining background branching fractions.
541  for (int iMode = 0; iMode < allModes.getNMode(); ++iMode) {
542  EvtDecayBase* mode = allModes.getDecayModel(iMode);
543  if (!mode) continue;
544  bool match(false);
545  for (int jMode = 0; jMode < sigModes.getNMode(); ++jMode)
546  if (mode->matchingDecay(*(sigModes.getDecayModel(jMode)))) {
547  match = true; break;
548  }
549  if (match) bkgModes.removeMode(mode);
550  else {
551  signal->second.map.push_back(iMode);
552  signal->second.bfs[1] += mode->getBranchingFraction();
553  sum += mode->getBranchingFraction();
554  }
555  }
556  signal->second.modes = bkgModes;
557  for (int iBf = 0; iBf < (int)signal->second.bfs.size(); ++iBf)
558  signal->second.bfs[iBf] /= sum;
559  }
560  if (final) updated = true;
561 
562 }
563 
564 //--------------------------------------------------------------------------
565 
566 // Update the Pythia event record with an EvtGen decay tree.
567 
568 // The production vertex of each particle (which can also be obtained
569 // in EvtGen via EvtParticle::get4Pos()) is set by the decay vertex of
570 // its mother, which in turn is calculated from the mother's
571 // lifetime. The status code 93 is used to indicate an external decay,
572 // while the status code 94 is used to indicate an oscillated external
573 // decay. If the progenitor has a single daughter with the same ID,
574 // this daughter is used as the progenitor. This is used to prevent
575 // double oscillations.
576 
577 // If the arguments after egPro are no NULL and a particle in the
578 // decay tree is a signal particle, the decay for this particle is
579 // removed and the particle is stored as a signal candidate in the
580 // pySigs and egSigs vectors, to be decayed later. However, if any of
581 // these arguments is NULL then the entire tree is written.
582 
583 void EvtGenDecays::updateEvent(Pythia8::Particle* pyPro, EvtParticle* egPro,
584  std::vector<int>* pySigs, std::vector<EvtParticle*>* egSigs,
585  std::vector<double>* bfs, double* wgt)
586 {
587 
588  // Set up the mother vector.
589  if (!pythiaPtr) return;
590  EvtParticle* egMom = egPro;
591  if (egPro->getNDaug() == 1 && egPro->getPDGId() ==
592  egPro->getDaug(0)->getPDGId()) egMom = egPro->getDaug(0);
593  Pythia8::Event& event = pythiaPtr->event;
594  std::vector< std::pair<EvtParticle*, int> >
595  moms(1, std::pair<EvtParticle*, int>(egMom, pyPro->index()));
596 
597  // Loop over the mothers.
598  while (moms.size() != 0) {
599 
600  // Check if particle can decay.
601  egMom = moms.back().first;
602  int iMom = moms.back().second;
603  Pythia8::Particle* pyMom = &event[iMom];
604  moms.pop_back();
605  if (!checkVertex(pyMom)) continue;
606  bool osc(checkOsc(egMom));
607 
608  // Set the children of the mother.
609  pyMom->daughters(event.size(), event.size() + egMom->getNDaug() - 1);
610  pyMom->statusNeg();
611  Pythia8::Vec4 vProd = pyMom->vDec();
612  for (int iDtr = 0 ; iDtr < (int)egMom->getNDaug(); ++iDtr) {
613  EvtParticle* egDtr = egMom->getDaug(iDtr);
614  int id = egDtr->getPDGId();
615  EvtVector4R p = egDtr->getP4Lab();
616  int idx = event.append(id, 93, iMom, 0, 0, 0, 0, 0, p.get(1),
617  p.get(2), p.get(3), p.get(0), egDtr->mass());
618  Pythia8::Particle* pyDtr = &event.back();
619  pyDtr->vProd(vProd);
620  pyDtr->tau(egDtr->getLifetime());
621  if (pySigs && egSigs && bfs && wgt && checkSignal(pyDtr)) {
622  pySigs->push_back(pyDtr->index());
623  egSigs->push_back(egDtr);
624  bfs->push_back(signal->second.bfs[0]);
625  (*wgt) *= 1 - bfs->back();
626  egDtr->deleteDaughters();
627  }
628  if (osc) pyDtr->status(94);
629  if (egDtr->getAttribute("FSR"))
630  pyDtr->status(Belle2::GeneratorConst::FSR_STATUS_CODE);
631  if (egDtr->getNDaug() > 0)
632  moms.push_back(std::pair<EvtParticle*, int>(egDtr, idx));
633  }
634  }
635 
636 }
637 
638 //--------------------------------------------------------------------------
639 
640 // Check if a particle can decay.
641 
642 // Modified slightly from ParticleDecays::checkVertex.
643 
644 bool EvtGenDecays::checkVertex(Pythia8::Particle* pyPro)
645 {
646  using ::Pythia8::pow2;
647 
648  if (!limitDecay) return true;
649  if (limitTau0 && pyPro->tau0() > tau0Max) return false;
650  if (limitTau && pyPro->tau() > tauMax) return false;
651  if (limitRadius && pow2(pyPro->xDec()) + pow2(pyPro->yDec())
652  + pow2(pyPro->zDec()) > pow2(rMax)) return false;
653  if (limitCylinder && (pow2(pyPro->xDec()) + pow2(pyPro->yDec())
654  > pow2(xyMax) || abs(pyPro->zDec()) > zMax)) return false;
655  return true;
656 }
657 
658 //--------------------------------------------------------------------------
659 
660 // Check if a particle is signal.
661 
662 bool EvtGenDecays::checkSignal(Pythia8::Particle* pyPro)
663 {
664  signal = signals.find(pyPro->id());
665  if (signal != signals.end() && (signal->second.status < 0 ||
666  signal->second.status == pyPro->status())) return true;
667  else return false;
668 }
669 
670 //--------------------------------------------------------------------------
671 
672 // Check if an EvtGen particle has oscillated.
673 
674 // The criteria defined here for oscillation is a single daughter but
675 // with a different ID from the mother.
676 
677 bool EvtGenDecays::checkOsc(EvtParticle* egPro)
678 {
679  if (egPro && egPro->getNDaug() == 1 &&
680  egPro->getPDGId() != egPro->getDaug(0)->getPDGId()) return true;
681  else return false;
682 }
683 
684 //==========================================================================
EvtGenDecays::rMax
double rMax
the pythia parameter ParticleDecays:rMax
Definition: EvtGenDecays.h:208
EvtGenDecays::readDecayFile
void readDecayFile(std::string decayFile, bool xml=false)
Read an EvtGen user decay file.
Definition: EvtGenDecays.h:137
EvtGenDecays::Signal::modes
EvtParticleDecayList modes
list of decay modes &
Definition: EvtGenDecays.h:154
EvtGenDecays::models
std::list< EvtDecayBase * > models
list of model pointers
Definition: EvtGenDecays.h:146
EvtGenDecays::checkSignal
bool checkSignal(Pythia8::Particle *pyPro)
Check if a particle is signal.
Definition: EvtGenDecays.h:662
EvtGenDecays::operator=
EvtGenDecays & operator=(const EvtGenDecays &)=delete
forbid assignment
EvtGenDecays::Signal::bfs
std::vector< double > bfs
branching fractions
Definition: EvtGenDecays.h:152
EvtGenDecays::zMax
double zMax
the pythia parameter ParticleDecays:zMax
Definition: EvtGenDecays.h:210
EvtGenDecays::updateEvent
void updateEvent(Pythia8::Particle *pyPro, EvtParticle *egPro, std::vector< int > *pySigs=0, std::vector< EvtParticle * > *egSigs=0, std::vector< double > *bfs=0, double *wgt=0)
Update the Pythia event record with an EvtGen decay tree.
Definition: EvtGenDecays.h:583
EvtGenDecays::signalSuffix
std::string signalSuffix
The suffix indicating an EvtGen particle or alias is signal.
Definition: EvtGenDecays.h:160
EvtGenDecays::extOwner
bool extOwner
bool has external model
Definition: EvtGenDecays.h:142
EvtGenRandom::random
double random()
Return a random number.
Definition: EvtGenDecays.h:50
EvtGenDecays
A class to perform decays via the external EvtGen decay program, see http://evtgen....
Definition: EvtGenDecays.h:94
EvtGenDecays::updatePythia
void updatePythia()
Update the Pythia particle database from EvtGen.
Definition: EvtGenDecays.h:440
EvtGenDecays::evtgen
EvtGen * evtgen
The EvtGen object.
Definition: EvtGenDecays.h:191
EvtGenDecays::Signal::map
std::vector< int > map
the map for statuses
Definition: EvtGenDecays.h:153
EvtGenRandom::EvtGenRandom
EvtGenRandom(Pythia8::Rndm *rndmPtrIn)
Constructor.
Definition: EvtGenDecays.h:47
EvtGenDecays::rndm
EvtGenRandom rndm
Random number wrapper for EvtGen.
Definition: EvtGenDecays.h:188
EvtGenDecays::NTRYDECAY
static const int NTRYDECAY
Number of times to try a decay sampling (constant).
Definition: EvtGenDecays.h:182
EvtGenDecays::EvtGenDecays
EvtGenDecays(Pythia8::Pythia *pythiaPtrIn, EvtGen *evtGen, bool limit=true)
Expert constructor which takes an already initialized EvtGen instance.
Definition: EvtGenDecays.h:261
EvtGenRandom
A class to wrap the Pythia random number generator for use by EvtGen.
Definition: EvtGenDecays.h:35
EvtGenDecays::tauMax
double tauMax
the pythia parameter ParticleDecays:tauMax
Definition: EvtGenDecays.h:207
EvtGenDecays::incIds
std::set< int > incIds
Set of particle IDs to include decays with EvtGen.
Definition: EvtGenDecays.h:193
EvtGenDecays::getDecayLimits
void getDecayLimits(bool limit)
Get the Decay limits from Pythia.
Definition: EvtGenDecays.h:293
EvtGenDecays::xyMax
double xyMax
the pythia parameter ParticleDecays:xyMax
Definition: EvtGenDecays.h:209
EvtGenDecays::excIds
std::set< int > excIds
Set of particle IDs to exclude decays with EvtGen.
Definition: EvtGenDecays.h:194
EvtGenDecays::signals
std::map< int, Signal > signals
signal particles and statuses
Definition: EvtGenDecays.h:157
EvtGenDecays::fsrPtr
EvtAbsRadCorr * fsrPtr
FSR model pointer.
Definition: EvtGenDecays.h:145
EvtGenDecays::tau0Max
double tau0Max
Parameters used to check if a particle should decay (as set via Pythia).
Definition: EvtGenDecays.h:206
EvtGenDecays::limitTau0
bool limitTau0
the pythia parameter ParticleDecays:limitTau0
Definition: EvtGenDecays.h:211
EvtGenDecays::updateData
void updateData(bool final=false)
Update the particles to decay with EvtGen, and the selected signals.
Definition: EvtGenDecays.h:490
EvtGenDecays::updated
bool updated
Flag whether the final particle update has been performed.
Definition: EvtGenDecays.h:197
EvtGenDecays::limitTau
bool limitTau
the pythia parameter ParticleDecays:limitTau
Definition: EvtGenDecays.h:212
EvtGenDecays::pythiaPtr
Pythia8::Pythia * pythiaPtr
The pointer to the associated Pythia object.
Definition: EvtGenDecays.h:185
EvtGenDecays::checkVertex
bool checkVertex(Pythia8::Particle *pyPro)
Check if a particle can decay.
Definition: EvtGenDecays.h:644
EvtGenDecays::checkOsc
bool checkOsc(EvtParticle *egPro)
Check if an EvtGen particle has oscillated.
Definition: EvtGenDecays.h:677
EvtGenDecays::signal
std::map< int, Signal >::iterator signal
The selected signal iterator.
Definition: EvtGenDecays.h:200
EvtGenDecays::updateEvtGen
void updateEvtGen()
Update the EvtGen particle database from Pythia.
Definition: EvtGenDecays.h:464
EvtGenDecays::limitRadius
bool limitRadius
the pythia parameter ParticleDecays:limitRadius
Definition: EvtGenDecays.h:213
EvtGenDecays::limitDecay
bool limitDecay
is the decay limited
Definition: EvtGenDecays.h:215
EvtGenRandom::rndmPtr
Pythia8::Rndm * rndmPtr
The random number pointer.
Definition: EvtGenDecays.h:53
EvtGenDecays::Signal::status
int status
status needed to consider a particle a signal
Definition: EvtGenDecays.h:150
EvtGenDecays::exclude
void exclude(int id)
Stop EvtGen decaying a particle.
Definition: EvtGenDecays.h:128
EvtGenDecays::Signal
Map of signal particle info.
Definition: EvtGenDecays.h:149
EvtGenDecays::Signal::egId
EvtId egId
the Lund ID
Definition: EvtGenDecays.h:151
EvtGenDecays::fsrOwner
bool fsrOwner
bool has FSR model
Definition: EvtGenDecays.h:143
EvtGenDecays::limitCylinder
bool limitCylinder
the pythia parameter ParticleDecays:limitCylinder
Definition: EvtGenDecays.h:214
EvtGenDecays::extPtr
EvtExternalGenList * extPtr
External model pointer.
Definition: EvtGenDecays.h:144
EvtGenDecays::~EvtGenDecays
~EvtGenDecays()
Destructor.
Definition: EvtGenDecays.h:108
EvtGenDecays::decay
double decay()
Perform all decays and return the event weight.
Definition: EvtGenDecays.h:334