2 //--------------------------------------------------------------------------
5 // This software is part of the EvtGen package. If you use all or part
6 // of it, please give an appropriate acknowledgement.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 2011 University of Warwick, UK
11 // Module: EvtPythiaEngine
13 // Description: Interface to the Pytha 8 external generator
15 // Modification history:
17 // John Back April 2011 Module created
19 //------------------------------------------------------------------------
21 #include "EvtGenExternal/EvtPythiaEngine.hh"
23 #include "EvtGenBase/EvtPDL.hh"
24 #include "EvtGenBase/EvtDecayTable.hh"
25 #include "EvtGenBase/EvtSpinType.hh"
26 #include "EvtGenBase/EvtParticleFactory.hh"
27 #include "EvtGenBase/EvtReport.hh"
29 #include "EvtGenBase/EvtExtGeneratorCommandsTable.hh"
30 #include "EvtGenExternal/EvtPythia6CommandConverter.hh"
32 //#include "Pythia8/Event.h"
33 #include "pythia8175/include/Event.h"
40 EvtPythiaEngine::EvtPythiaEngine(std::string xmlDir, bool convertPhysCodes,
41 bool useEvtGenRandom) {
43 // Create two Pythia generators. One will be for generic
44 // Pythia decays in the decay.dec file. The other one will be to
45 // only decay aliased particles, which are in general "signal"
46 // decays different from those in the decay.dec file.
47 // Even though it is not possible to have two different particle
48 // versions in one Pythia generator, we can use two generators to
49 // get the required behaviour.
51 report(INFO,"EvtGen")<<"Creating generic Pythia generator"<<endl;
52 _genericPythiaGen = new Pythia8::Pythia(xmlDir);
53 _genericPartData = _genericPythiaGen->particleData;
55 report(INFO,"EvtGen")<<"Creating alias Pythia generator"<<endl;
56 _aliasPythiaGen = new Pythia8::Pythia(xmlDir);
57 _aliasPartData = _aliasPythiaGen->particleData;
59 _thePythiaGenerator = 0;
60 _daugPDGVector.clear(); _daugP4Vector.clear();
62 _convertPhysCodes = convertPhysCodes;
64 // Specify if we are going to use the random number generator (engine)
65 // from EvtGen for Pythia 8.
66 _useEvtGenRandom = useEvtGenRandom;
68 _evtgenRandom = new EvtPythiaRandom();
74 EvtPythiaEngine::~EvtPythiaEngine() {
76 delete _genericPythiaGen; _genericPythiaGen = 0;
77 delete _aliasPythiaGen; _aliasPythiaGen = 0;
79 delete _evtgenRandom; _evtgenRandom = 0;
81 _thePythiaGenerator = 0;
83 this->clearDaughterVectors();
84 this->clearPythiaModeMap();
88 void EvtPythiaEngine::clearDaughterVectors() {
89 _daugPDGVector.clear();
90 _daugP4Vector.clear();
93 void EvtPythiaEngine::clearPythiaModeMap() {
95 PythiaModeMap::iterator iter;
96 for (iter = _pythiaModeMap.begin(); iter != _pythiaModeMap.end(); ++iter) {
98 std::vector<int> modeVector = iter->second;
103 _pythiaModeMap.clear();
107 void EvtPythiaEngine::initialise() {
109 if (_initialised) {return;}
111 this->clearPythiaModeMap();
113 this->updateParticleLists();
115 // Hadron-level processes only (hadronized, string fragmentation and secondary decays).
116 // We do not want to generate the full pp or e+e- event structure etc..
117 _genericPythiaGen->readString("ProcessLevel:all = off");
118 _aliasPythiaGen->readString("ProcessLevel:all = off");
120 // Apply any other physics (or special particle) requirements/cuts etc..
121 this->updatePhysicsParameters();
123 // Set the random number generator
124 if (_useEvtGenRandom == true) {
126 _genericPythiaGen->setRndmEnginePtr(_evtgenRandom);
127 _aliasPythiaGen->setRndmEnginePtr(_evtgenRandom);
131 _genericPythiaGen->init();
132 _aliasPythiaGen->init();
138 bool EvtPythiaEngine::doDecay(EvtParticle* theParticle) {
140 // Store the mother particle within a Pythia8 Event object.
141 // Then do the hadron level decays.
142 // The EvtParticle must be a colour singlet (meson/baryon/lepton), i.e. not a gluon or quark
144 // We delete any daughters the particle may have, since we are asking Pythia
145 // to generate the decay anew. Also note that _any_ Pythia decay allowed for the particle
146 // will be generated and not the specific Pythia decay mode that EvtGen has already
147 // specified. This is necessary since we only want to initialise the Pythia decay table
148 // once; all Pythia branching fractions for a given mother particle are renormalised to sum to 1.0.
149 // In EvtGen decay.dec files, it may be the case that Pythia decays are only used
150 // for some of the particle decays (i.e. Pythia BF sum < 1.0). As we loop over many events,
151 // the total frequency for each Pythia decay mode will normalise correctly to what
152 // we wanted via the specifications made to the decay.dec file, even though event-by-event
153 // the EvtGen decay channel and the Pythia decay channel may be different.
155 if (_initialised == false) {this->initialise();}
157 if (theParticle == 0) {
158 report(INFO,"EvtGen")<<"Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay."<<endl;
162 // Delete EvtParticle daughters if they already exist
163 if (theParticle->getNDaug() != 0) {
164 bool keepChannel(false);
165 theParticle->deleteDaughters(keepChannel);
168 EvtId particleId = theParticle->getId();
169 int isAlias = particleId.isAlias();
171 // Choose the generator depending if we have an aliased (parent) particle or not
172 _thePythiaGenerator = _genericPythiaGen;
173 if (isAlias == 1) {_thePythiaGenerator = _aliasPythiaGen;}
175 //_thePythiaGenerator->settings.listChanged();
177 // Need to use the reference to the Pythia8::Event object,
178 // otherwise it will just return a new empty, default event object.
179 Pythia8::Event& theEvent = _thePythiaGenerator->event;
182 // Initialise the event to be the particle rest frame
183 int PDGCode = EvtPDL::getStdHep(particleId);
186 int colour(0), anticolour(0);
187 double px(0.0), py(0.0), pz(0.0);
188 double m0 = theParticle->mass();
191 theEvent.append(PDGCode, status, colour, anticolour, px, py, pz, E, m0);
193 // Generate the Pythia event
195 bool generatedEvent(false);
196 for (iTrial = 0; iTrial < 10; iTrial++) {
198 generatedEvent = _thePythiaGenerator->next();
199 if (generatedEvent) {break;}
205 if (generatedEvent) {
207 // Store the daughters for this particle from the Pythia decay tree.
208 // This is a recursive function that will continue looping through
209 // all available daughters until the first set of non-quark and non-gluon
210 // particles are encountered in the Pythia Event structure.
212 // First, clear up the internal vectors storing the daughter
213 // EvtId types and 4-momenta.
214 this->clearDaughterVectors();
216 // Now store the daughter info. Since this is a recursive function
217 // to loop through the full Pythia decay tree, we do not want to create
218 // EvtParticles here but in the next step.
219 this->storeDaughterInfo(theParticle, 1);
221 // Now create the EvtParticle daughters of the (parent) particle.
222 // We need to use the EvtParticle::makeDaughters function
223 // owing to the way EvtParticle stores parent-daughter information.
224 this->createDaughterEvtParticles(theParticle);
226 //theParticle->printTree();
227 //theEvent.list(true, true);
237 void EvtPythiaEngine::storeDaughterInfo(EvtParticle* theParticle, int startInt) {
239 Pythia8::Event& theEvent = _thePythiaGenerator->event;
241 std::vector<int> daugList = theEvent.daughterList(startInt);
243 std::vector<int>::iterator daugIter;
244 for (daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter) {
246 int daugInt = *daugIter;
248 // Ask if the daughter is a quark or gluon. If so, recursively search again.
249 Pythia8::Particle& daugParticle = theEvent[daugInt];
251 if (daugParticle.isQuark() || daugParticle.isGluon()) {
253 // Recursively search for correct daughter type
254 this->storeDaughterInfo(theParticle, daugInt);
258 // We have a daughter that is not a quark nor gluon particle.
259 // Make sure we are not double counting particles, since several quarks
260 // and gluons make one particle.
261 // Set the status flag for any "new" particle to say that we have stored it.
262 // Use status flag = 1000 (within the user allowed range for Pythia codes).
264 // Check that the status flag for the particle is not equal to 1000
265 int statusCode = daugParticle.status();
266 if (statusCode != 1000) {
268 int daugPDGInt = daugParticle.id();
270 double px = daugParticle.px();
271 double py = daugParticle.py();
272 double pz = daugParticle.pz();
273 double E = daugParticle.e();
274 EvtVector4R daughterP4(E, px, py, pz);
276 // Now store the EvtId and 4-momentum in the internal vectors
277 _daugPDGVector.push_back(daugPDGInt);
278 _daugP4Vector.push_back(daughterP4);
280 // Set the status flag for the Pythia particle to let us know
281 // that we have already considered it to avoid double counting.
282 daugParticle.status(1000);
284 } // Status code != 1000
292 void EvtPythiaEngine::createDaughterEvtParticles(EvtParticle* theParent) {
294 if (theParent == 0) {
295 report(INFO,"EvtGen")<<"Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null"<<endl;
299 // Get the list of Pythia decay modes defined for this particle id alias.
300 // It would be easier to just use the decay channel number that Pythia chose to use
301 // for the particle decay, but this is not accessible from the Pythia interface at present.
303 int nDaughters = _daugPDGVector.size();
304 std::vector<EvtId> daugAliasIdVect(0);
306 EvtId particleId = theParent->getId();
307 // Check to see if we have an anti-particle. If we do, charge conjugate the particle id to get the
308 // Pythia "alias" we can compare with the defined (particle) Pythia modes.
309 int PDGId = EvtPDL::getStdHep(particleId);
310 int aliasInt = particleId.getAlias();
311 int pythiaAliasInt(aliasInt);
314 // We have an anti-particle.
315 EvtId conjPartId = EvtPDL::chargeConj(particleId);
316 pythiaAliasInt = conjPartId.getAlias();
319 std::vector<int> pythiaModes = _pythiaModeMap[pythiaAliasInt];
321 // Loop over all available Pythia decay modes and find the channel that matches
322 // the daughter ids. Set each daughter id to also use the alias integer.
323 // This will then convert the Pythia generated channel to the EvtGen alias defined one.
325 std::vector<int>::iterator modeIter;
328 for (modeIter = pythiaModes.begin(); modeIter != pythiaModes.end(); ++modeIter) {
330 // Stop the loop if we have the right decay mode channel
331 if (gotMode) {break;}
333 int pythiaModeInt = *modeIter;
335 EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, pythiaModeInt);
337 if (decayModel != 0) {
339 int nModeDaug = decayModel->getNDaug();
341 // We need to make sure that the number of daughters match
342 if (nDaughters == nModeDaug) {
345 for (iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++) {
347 EvtId daugId = decayModel->getDaug(iModeDaug);
348 int daugPDGId = EvtPDL::getStdHep(daugId);
349 // Pythia has used the right PDG codes for this decay mode, even for conjugate modes
350 int pythiaPDGId = _daugPDGVector[iModeDaug];
352 if (daugPDGId == pythiaPDGId) {
353 daugAliasIdVect.push_back(daugId);
356 } // Loop over EvtGen mode daughters
358 int daugAliasSize = daugAliasIdVect.size();
359 if (daugAliasSize == nDaughters) {
360 // All daughter Id codes are accounted for. Set the flag to stop the loop.
363 // We do not have the correct daughter ordering. Clear the id vector
364 // and try another mode.
365 daugAliasIdVect.clear();
368 } // Same number of daughters
372 } // Loop over available Pythia modes
374 if (gotMode == false) {
376 // We did not find a match for the daughter aliases. Just use the normal PDG codes
377 // from the Pythia decay result
379 for (iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++) {
381 int daugPDGCode = _daugPDGVector[iPyDaug];
382 EvtId daugPyId = EvtPDL::evtIdFromStdHep(daugPDGCode);
383 daugAliasIdVect.push_back(daugPyId);
388 // Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised.
389 theParent->makeDaughters(nDaughters, daugAliasIdVect);
391 // Now set the 4-momenta of the daughters.
393 // Can use an iterator here, but we already had to use the vector size...
394 for (iDaug = 0; iDaug < nDaughters; iDaug++) {
396 EvtParticle* theDaughter = theParent->getDaug(iDaug);
398 // Set the correct 4-momentum for each daughter particle.
399 if (theDaughter != 0) {
400 EvtId theDaugId = daugAliasIdVect[iDaug];
401 const EvtVector4R theDaugP4 = _daugP4Vector[iDaug];
402 theDaughter->init(theDaugId, theDaugP4);
409 void EvtPythiaEngine::updateParticleLists() {
411 // Use the EvtGen decay table (decay/user.dec) to update particle entries
412 // for Pythia. Pythia 8 should use the latest PDG codes, so if the evt.pdl
413 // file is up to date, just let Pythia 8 find the particle properties
414 // knowing the PDG code integer. If we want to use evt.pdl for _all_
415 // particle properties, then we need to make sure that this is up to date,
416 // and modify the code in this class to read that data and use it...
417 // Using the PDG code only also avoids the need to convert EvtGen particle names
418 // to Pythia particle names.
420 // Loop over all entries in the EvtPDL particle data table.
421 // Aliases are added at the end with id numbers equal to the
422 // original particle, but with alias integer = lastPDLEntry+1 etc..
424 int nPDL = EvtPDL::entries();
426 // Reset the _addedPDGCodes map that keeps track
427 // of any new particles added to the Pythia input data stream
428 _addedPDGCodes.clear();
430 for (iPDL = 0; iPDL < nPDL; iPDL++) {
432 EvtId particleId = EvtPDL::getEntry(iPDL);
433 int aliasInt = particleId.getAlias();
435 // Check which particles have a Pythia decay defined.
436 // Get the list of all possible decays for the particle, using the alias integer.
437 // If the particle is not actually an alias, aliasInt = idInt.
439 // Should change isJetSet to isPythia eventually.
440 bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia(aliasInt);
442 if (hasPythiaDecays) {
444 int isAlias = particleId.isAlias();
446 int PDGCode = EvtPDL::getStdHep(particleId);
448 // Decide what generator to use depending on ether we have
449 // an alias particle or not
450 _thePythiaGenerator = _genericPythiaGen;
451 _theParticleData = _genericPartData;
453 _thePythiaGenerator = _aliasPythiaGen;
454 _theParticleData = _aliasPartData;
457 // Find the Pythia particle name given the standard PDG code integer
458 std::string dataName = _theParticleData.name(PDGCode);
459 bool alreadyStored(false);
460 if (_addedPDGCodes.find(abs(PDGCode)) != _addedPDGCodes.end()) {alreadyStored = true;}
462 if (dataName == " " && alreadyStored == false) {
464 // Particle and its antiparticle does not exist in the Pythia database.
465 // Create a new particle, then create the new decay modes.
466 this->createPythiaParticle(particleId, PDGCode);
470 // Particle exists in the Pythia database.
471 // Could update mass/lifetime values here. For now just use Pythia defaults.
475 // For the particle, create the Pythia decay modes.
476 // Update Pythia data tables.
477 this->updatePythiaDecayTable(particleId, aliasInt, PDGCode);
479 } // Loop over Pythia decays
481 } // Loop over EvtPDL entries
483 //report(INFO,"EvtGen")<<"Writing out changed generic Pythia decay list"<<endl;
484 //_genericPythiaGen->particleData.listChanged();
486 //report(INFO,"EvtGen")<<"Writing out changed alias Pythia decay list"<<endl;
487 //_aliasPythiaGen->particleData.listChanged();
491 void EvtPythiaEngine::updatePythiaDecayTable(EvtId& particleId, int aliasInt, int PDGCode) {
493 // Update the particle data table in Pythia.
494 // The tables store information about the allowed decay modes
495 // whre the PDGId for all particles must be positive; anti-particles are stored
496 // with the corresponding particle entry.
497 // Since we do not want to implement CP violation here, just use the same branching
498 // fractions for particle and anti-particle modes.
500 int nModes = EvtDecayTable::getInstance()->getNModes(aliasInt);
503 bool firstMode(true);
505 // Only process positive PDG codes.
506 if (PDGCode < 0) {return;}
508 // Keep track of which decay modes are Pythia decays for each aliasInt
509 std::vector<int> pythiaModes(0);
511 // Loop over the decay modes for this particle
512 for (iMode = 0; iMode < nModes; iMode++) {
514 EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, iMode);
516 if (decayModel != 0) {
518 int nDaug = decayModel->getNDaug();
520 // If the decay mode has no daughters, then that means that there will be
521 // no entries for any submode re-definitions for Pythia.
522 // This sometimes occurs for any mode using non-standard Pythia 6 codes.
523 // Do not refine the decay mode, i.e. accept the Pythia 8 default (if it exists).
526 // Check to see if we have a Pythia decay mode
527 std::string modelName = decayModel->getModelName();
529 if (modelName == "PYTHIA") {
531 // Keep track which decay mode is a Pythia one. We need this in order to
532 // reassign alias Id values for particles generated in the decay.
533 pythiaModes.push_back(iMode);
535 std::ostringstream oss;
536 oss.setf(std::ios::scientific);
537 // Write out the absolute value of the PDG code, since
538 // particles and anti-particles occupy the same part of the Pythia table.
542 // Create a new channel
543 oss <<":oneChannel = ";
547 oss <<":addChannel = ";
550 // Select all channels (particle and anti-particle).
551 // For CP violation, or different BFs for particle and anti-particle,
552 // use options 2 or 3 (not here).
554 oss << onMode << " ";
556 double BF = decayModel->getBranchingFraction();
559 // Need to convert the old Pythia physics mode integers with the new ones
560 // To do this, get the model argument and write a conversion method.
561 int modeInt = this->getModeInt(decayModel);
565 for (iDaug = 0; iDaug < nDaug; iDaug++) {
567 EvtId daugId = decayModel->getDaug(iDaug);
568 int daugPDG = EvtPDL::getStdHep(daugId);
569 oss << " " << daugPDG;
573 _thePythiaGenerator->readString(oss.str());
579 report(INFO,"EvtGen")<<"Warning in EvtPythiaEngine. Trying to redefine Pythia table for "
580 <<EvtPDL::name(particleId)<<" for a decay channel that has no daughters."
581 <<" Keeping Pythia default (if available)."<<endl;
587 report(INFO,"EvtGen")<<"Error in EvtPythiaEngine. decayModel is null for particle "
588 <<EvtPDL::name(particleId)<<" mode number "<<iMode<<endl;
594 _pythiaModeMap[aliasInt] = pythiaModes;
596 // Now, renormalise the decay branching fractions to sum to 1.0
597 std::ostringstream rescaleStr;
598 rescaleStr.setf(std::ios::scientific);
599 rescaleStr << PDGCode << ":rescaleBR = 1.0";
601 _thePythiaGenerator->readString(rescaleStr.str());
605 int EvtPythiaEngine::getModeInt(EvtDecayBase* decayModel) {
607 int tmpModeInt(0), modeInt(0);
609 if (decayModel != 0) {
611 int nVars = decayModel->getNArg();
612 // Just read the first integer, which specifies the Pythia decay model.
613 // Ignore any other values.
615 tmpModeInt = static_cast<int>(decayModel->getArg(0));
619 if (_convertPhysCodes) {
621 // Extra code to convert the old Pythia decay model integer MDME(ICC,2) to the new one.
622 // This should be removed eventually after updating decay.dec files to use
623 // the new convention.
625 if (tmpModeInt == 0) {
626 modeInt = 0; // phase-space
627 } else if (tmpModeInt == 1) {
628 modeInt = 1; // omega or phi -> 3pi
629 } else if (tmpModeInt == 2) {
630 modeInt = 11; // Dalitz decay
631 } else if (tmpModeInt == 3) {
632 modeInt = 2; // V -> PS PS
633 } else if (tmpModeInt == 4) {
634 modeInt = 92; // onium -> ggg or gg gamma
635 } else if (tmpModeInt == 11) {
636 modeInt = 42; // phase-space of hadrons from available quarks
637 } else if (tmpModeInt == 12) {
638 modeInt = 42; // phase-space for onia resonances
639 } else if (tmpModeInt == 13) {
640 modeInt = 43; // phase-space of at least 3 hadrons
641 } else if (tmpModeInt == 14) {
642 modeInt = 44; // phase-space of at least 4 hadrons
643 } else if (tmpModeInt == 15) {
644 modeInt = 45; // phase-space of at least 5 hadrons
645 } else if (tmpModeInt >= 22 && tmpModeInt <= 30) {
646 modeInt = tmpModeInt + 40; // phase space of hadrons with fixed multiplicity (modeInt - 60)
647 } else if (tmpModeInt == 31) {
648 modeInt = 42; // two or more quarks phase-space; one spectactor quark
649 } else if (tmpModeInt == 32) {
650 modeInt = 91; // qqbar or gg pair
651 } else if (tmpModeInt == 33) {
652 modeInt = 0; // triplet q X qbar, where X = gluon or colour singlet (superfluous, since g's are created anyway)
653 } else if (tmpModeInt == 41) {
654 modeInt = 21; // weak decay phase space, weighting nu_tau spectrum
655 } else if (tmpModeInt == 42) {
656 modeInt = 22; // weak decay V-A matrix element
657 } else if (tmpModeInt == 43) {
658 modeInt = 22; // weak decay V-A matrix element, quarks as jets (superfluous)
659 } else if (tmpModeInt == 44) {
660 modeInt = 22; // weak decay V-A matrix element, parton showers (superfluous)
661 } else if (tmpModeInt == 48) {
662 modeInt = 23; // weak decay V-A matrix element, at least 3 decay products
663 } else if (tmpModeInt == 50) {
664 modeInt = 0; // default behaviour
665 } else if (tmpModeInt == 51) {
666 modeInt = 0; // step threshold (channel switched off when mass daughters > mother mass
667 } else if (tmpModeInt == 52 || tmpModeInt == 53) {
668 modeInt = 0; // beta-factor threshold
669 } else if (tmpModeInt == 84) {
670 modeInt = 42; // unknown physics process - just use phase-space
671 } else if (tmpModeInt == 101) {
672 modeInt = 0; // continuation line
673 } else if (tmpModeInt == 102) {
674 modeInt = 0; // off mass shell particles.
676 report(INFO,"EvtGen")<<"Pythia mode integer "<<tmpModeInt
677 <<" is not recognised. Using phase-space model"<<endl;
678 modeInt = 0; // Use phase-space for anything else
683 // No need to convert the physics mode integer code
684 modeInt = tmpModeInt;
692 void EvtPythiaEngine::createPythiaParticle(EvtId& particleId, int PDGCode) {
694 // Use the EvtGen name, PDGId and other variables to define the new Pythia particle.
695 EvtId antiPartId = EvtPDL::chargeConj(particleId);
697 std::string aliasName = EvtPDL::name(particleId); // If not an alias, aliasName = normal name
698 std::string antiName = EvtPDL::name(antiPartId);
700 EvtSpinType::spintype spinType = EvtPDL::getSpinType(particleId);
701 int spin = EvtSpinType::getSpin2(spinType);
703 int charge = EvtPDL::chg3(particleId);
705 // Must set the correct colour type manually here, since the evt.pdl file
706 // does not store this information. This is required for quarks otherwise
707 // Pythia cannot generate the decay properly.
708 int PDGId = EvtPDL::getStdHep(particleId);
711 colour = 2; // gluons
712 } else if (PDGId <= 8 && PDGId > 0) {
713 colour = 1; // single quarks
716 double m0 = EvtPDL::getMeanMass(particleId);
717 double mWidth = EvtPDL::getWidth(particleId);
718 double mMin = EvtPDL::getMinMass(particleId);
719 double mMax = EvtPDL::getMaxMass(particleId);
721 double tau0 = EvtPDL::getctau(particleId);
723 std::ostringstream oss;
724 oss.setf(std::ios::scientific);
725 int absPDGCode = abs(PDGCode);
726 oss << absPDGCode << ":new = " << aliasName << " " << antiName << " "
727 << spin << " " << charge << " " << colour << " "
728 << m0 << " " << mWidth << " " << mMin << " " << mMax << " "
731 // Pass this information to Pythia
732 _thePythiaGenerator->readString(oss.str());
734 // Also store the absolute value of the PDG entry
735 // to keep track of which new particles have been added,
736 // which also automatically includes the anti-particle.
737 // We need to avoid creating new anti-particles when
738 // they already exist when the particle was added.
739 _addedPDGCodes[absPDGCode] = 1;
743 void EvtPythiaEngine::updatePhysicsParameters() {
745 // Update any more Pythia physics (or special particle) requirements/cuts etc..
746 // This should be used if any of the Pythia 6 parameters like JetSetPar MSTJ(i) = x
747 // are needed. Such commands will need to be implemented using the new interface
748 // pythiaGenerator->readString(cmd); Here cmd is a string telling Pythia 8
749 // what physics parameters to change. This will need to be done for the generic and
750 // alias generator pointers, as appropriate.
752 // Set the multiplicity level for hadronic weak decays
753 std::string multiWeakCut("ParticleDecays:multIncreaseWeak = 2.0");
754 _genericPythiaGen->readString(multiWeakCut);
755 _aliasPythiaGen->readString(multiWeakCut);
757 // Set the multiplicity level for all other decays
758 std::string multiCut("ParticleDecays:multIncrease = 4.5");
759 _genericPythiaGen->readString(multiCut);
760 _aliasPythiaGen->readString(multiCut);
762 //Now read in any custom configuration entered in the XML
763 GeneratorCommands commands = EvtExtGeneratorCommandsTable::getInstance()->getCommands("PYTHIA");
764 GeneratorCommands::iterator it = commands.begin();
766 for( ; it!=commands.end(); it++) {
768 Command command = *it;
769 std::vector<std::string> commandStrings;
771 if(command["VERSION"] == "PYTHIA6") {
772 report(INFO,"EvtGen")<<"Converting Pythia 6 command: "<<command["MODULE"]<<"("<<command["PARAM"]<<")="<<command["VALUE"]<<"..."<<endl;
773 commandStrings = convertPythia6Command(command);
774 } else if(command["VERSION"] == "PYTHIA8") {
775 commandStrings.push_back(command["MODULE"]+":"+command["PARAM"]+" = "+command["VALUE"]);
777 report(ERROR, "EvtGen") << "Pythia command received by EvtPythiaEngine has bad version:"<<endl;
778 report(ERROR, "EvtGen") << "Received "<<command["VERSION"]<<" but expected PYTHIA6 or PYTHIA8."<<endl;
779 report(ERROR, "EvtGen") << "The error is likely to be in EvtDecayTable.cpp"<<endl;
780 report(ERROR, "EvtGen") << "EvtGen will now abort."<<endl;
783 std::string generator = command["GENERATOR"];
784 if(generator == "GENERIC" || generator == "Generic" || generator == "generic" ||
785 generator == "BOTH" || generator == "Both" || generator == "both") {
786 std::vector<std::string>::iterator it2 = commandStrings.begin();
787 for( ; it2!=commandStrings.end(); it2++) {
788 report(INFO,"EvtGen")<<"Configuring generic Pythia generator: " << (*it2) << endl;
789 _genericPythiaGen->readString(*it2);
792 if(generator == "ALIAS" || generator == "Alias" || generator == "alias" ||
793 generator == "BOTH" || generator == "Both" || generator == "both") {
794 std::vector<std::string>::iterator it2 = commandStrings.begin();
795 for( ; it2!=commandStrings.end(); it2++) {
796 report(INFO,"EvtGen")<<"Configuring alias Pythia generator: " << (*it2) << endl;
797 _aliasPythiaGen->readString(*it2);