Line data Source code
1 : #ifdef EVTGEN_PYTHIA
2 : //--------------------------------------------------------------------------
3 : //
4 : // Environment:
5 : // This software is part of the EvtGen package. If you use all or part
6 : // of it, please give an appropriate acknowledgement.
7 : //
8 : // Copyright Information: See EvtGen/COPYRIGHT
9 : // Copyright (C) 2011 University of Warwick, UK
10 : //
11 : // Module: EvtPythiaEngine
12 : //
13 : // Description: Interface to the Pytha 8 external generator
14 : //
15 : // Modification history:
16 : //
17 : // John Back April 2011 Module created
18 : //
19 : //------------------------------------------------------------------------
20 :
21 : #include "EvtGenExternal/EvtPythiaEngine.hh"
22 :
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"
28 :
29 : #include "EvtGenBase/EvtExtGeneratorCommandsTable.hh"
30 : #include "EvtGenExternal/EvtPythia6CommandConverter.hh"
31 :
32 : #include "Pythia8/Event.h"
33 :
34 :
35 : #include <iostream>
36 : #include <sstream>
37 :
38 : using std::endl;
39 :
40 0 : EvtPythiaEngine::EvtPythiaEngine(std::string xmlDir, bool convertPhysCodes,
41 0 : bool useEvtGenRandom) {
42 :
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.
50 :
51 0 : report(Severity::Info,"EvtGen")<<"Creating generic Pythia generator"<<endl;
52 0 : _genericPythiaGen = new Pythia8::Pythia(xmlDir);
53 0 : _genericPartData = _genericPythiaGen->particleData;
54 :
55 0 : report(Severity::Info,"EvtGen")<<"Creating alias Pythia generator"<<endl;
56 0 : _aliasPythiaGen = new Pythia8::Pythia(xmlDir);
57 0 : _aliasPartData = _aliasPythiaGen->particleData;
58 :
59 0 : _thePythiaGenerator = 0;
60 0 : _daugPDGVector.clear(); _daugP4Vector.clear();
61 :
62 0 : _convertPhysCodes = convertPhysCodes;
63 :
64 : // Specify if we are going to use the random number generator (engine)
65 : // from EvtGen for Pythia 8.
66 0 : _useEvtGenRandom = useEvtGenRandom;
67 :
68 0 : _evtgenRandom = new EvtPythiaRandom();
69 :
70 0 : _initialised = false;
71 :
72 0 : }
73 :
74 0 : EvtPythiaEngine::~EvtPythiaEngine() {
75 :
76 0 : delete _genericPythiaGen; _genericPythiaGen = 0;
77 0 : delete _aliasPythiaGen; _aliasPythiaGen = 0;
78 :
79 0 : delete _evtgenRandom; _evtgenRandom = 0;
80 :
81 0 : _thePythiaGenerator = 0;
82 :
83 0 : this->clearDaughterVectors();
84 0 : this->clearPythiaModeMap();
85 :
86 0 : }
87 :
88 : void EvtPythiaEngine::clearDaughterVectors() {
89 0 : _daugPDGVector.clear();
90 0 : _daugP4Vector.clear();
91 0 : }
92 :
93 : void EvtPythiaEngine::clearPythiaModeMap() {
94 :
95 0 : PythiaModeMap::iterator iter;
96 0 : for (iter = _pythiaModeMap.begin(); iter != _pythiaModeMap.end(); ++iter) {
97 :
98 0 : std::vector<int> modeVector = iter->second;
99 0 : modeVector.clear();
100 :
101 0 : }
102 :
103 0 : _pythiaModeMap.clear();
104 :
105 0 : }
106 :
107 : void EvtPythiaEngine::initialise() {
108 :
109 0 : if (_initialised) {return;}
110 :
111 0 : this->clearPythiaModeMap();
112 :
113 0 : this->updateParticleLists();
114 :
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 0 : _genericPythiaGen->readString("ProcessLevel:all = off");
118 0 : _aliasPythiaGen->readString("ProcessLevel:all = off");
119 :
120 : // Apply any other physics (or special particle) requirements/cuts etc..
121 0 : this->updatePhysicsParameters();
122 :
123 : // Set the random number generator
124 0 : if (_useEvtGenRandom == true) {
125 :
126 0 : _genericPythiaGen->setRndmEnginePtr(_evtgenRandom);
127 0 : _aliasPythiaGen->setRndmEnginePtr(_evtgenRandom);
128 :
129 0 : }
130 :
131 0 : _genericPythiaGen->init();
132 0 : _aliasPythiaGen->init();
133 :
134 0 : _initialised = true;
135 :
136 0 : }
137 :
138 : bool EvtPythiaEngine::doDecay(EvtParticle* theParticle) {
139 :
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
143 :
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.
154 :
155 0 : if (_initialised == false) {this->initialise();}
156 :
157 0 : if (theParticle == 0) {
158 0 : report(Severity::Info,"EvtGen")<<"Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay."<<endl;
159 0 : return false;
160 : }
161 :
162 : // Delete EvtParticle daughters if they already exist
163 0 : if (theParticle->getNDaug() != 0) {
164 : bool keepChannel(false);
165 0 : theParticle->deleteDaughters(keepChannel);
166 0 : }
167 :
168 0 : EvtId particleId = theParticle->getId();
169 0 : int isAlias = particleId.isAlias();
170 :
171 : // Choose the generator depending if we have an aliased (parent) particle or not
172 0 : _thePythiaGenerator = _genericPythiaGen;
173 0 : if (isAlias == 1) {_thePythiaGenerator = _aliasPythiaGen;}
174 :
175 : //_thePythiaGenerator->settings.listChanged();
176 :
177 : // Need to use the reference to the Pythia8::Event object,
178 : // otherwise it will just return a new empty, default event object.
179 0 : Pythia8::Event& theEvent = _thePythiaGenerator->event;
180 0 : theEvent.reset();
181 :
182 : // Initialise the event to be the particle rest frame
183 0 : int PDGCode = EvtPDL::getStdHep(particleId);
184 :
185 : int status(1);
186 : int colour(0), anticolour(0);
187 : double px(0.0), py(0.0), pz(0.0);
188 0 : double m0 = theParticle->mass();
189 : double E = m0;
190 :
191 0 : theEvent.append(PDGCode, status, colour, anticolour, px, py, pz, E, m0);
192 :
193 : // Generate the Pythia event
194 : int iTrial(0);
195 : bool generatedEvent(false);
196 0 : for (iTrial = 0; iTrial < 10; iTrial++) {
197 :
198 0 : generatedEvent = _thePythiaGenerator->next();
199 0 : if (generatedEvent) {break;}
200 :
201 : }
202 :
203 : bool success(false);
204 :
205 0 : if (generatedEvent) {
206 :
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.
211 :
212 : // First, clear up the internal vectors storing the daughter
213 : // EvtId types and 4-momenta.
214 0 : this->clearDaughterVectors();
215 :
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 0 : this->storeDaughterInfo(theParticle, 1);
220 :
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 0 : this->createDaughterEvtParticles(theParticle);
225 :
226 : //theParticle->printTree();
227 : //theEvent.list(true, true);
228 :
229 : success = true;
230 :
231 0 : }
232 :
233 0 : return success;
234 :
235 0 : }
236 :
237 : void EvtPythiaEngine::storeDaughterInfo(EvtParticle* theParticle, int startInt) {
238 :
239 0 : Pythia8::Event& theEvent = _thePythiaGenerator->event;
240 :
241 0 : std::vector<int> daugList = theEvent.daughterList(startInt);
242 :
243 0 : std::vector<int>::iterator daugIter;
244 0 : for (daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter) {
245 :
246 0 : int daugInt = *daugIter;
247 :
248 : // Ask if the daughter is a quark or gluon. If so, recursively search again.
249 0 : Pythia8::Particle& daugParticle = theEvent[daugInt];
250 :
251 0 : if (daugParticle.isQuark() || daugParticle.isGluon()) {
252 :
253 : // Recursively search for correct daughter type
254 0 : this->storeDaughterInfo(theParticle, daugInt);
255 :
256 : } else {
257 :
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).
263 :
264 : // Check that the status flag for the particle is not equal to 1000
265 0 : int statusCode = daugParticle.status();
266 0 : if (statusCode != 1000) {
267 :
268 0 : int daugPDGInt = daugParticle.id();
269 :
270 0 : double px = daugParticle.px();
271 0 : double py = daugParticle.py();
272 0 : double pz = daugParticle.pz();
273 0 : double E = daugParticle.e();
274 0 : EvtVector4R daughterP4(E, px, py, pz);
275 :
276 : // Now store the EvtId and 4-momentum in the internal vectors
277 0 : _daugPDGVector.push_back(daugPDGInt);
278 0 : _daugP4Vector.push_back(daughterP4);
279 :
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 0 : daugParticle.status(1000);
283 :
284 0 : } // Status code != 1000
285 :
286 : }
287 :
288 : }
289 :
290 0 : }
291 :
292 : void EvtPythiaEngine::createDaughterEvtParticles(EvtParticle* theParent) {
293 :
294 0 : if (theParent == 0) {
295 0 : report(Severity::Info,"EvtGen")<<"Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null"<<endl;
296 0 : return;
297 : }
298 :
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.
302 :
303 0 : int nDaughters = _daugPDGVector.size();
304 0 : std::vector<EvtId> daugAliasIdVect(0);
305 :
306 0 : 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 0 : int PDGId = EvtPDL::getStdHep(particleId);
310 0 : int aliasInt = particleId.getAlias();
311 0 : int pythiaAliasInt(aliasInt);
312 :
313 0 : if (PDGId < 0) {
314 : // We have an anti-particle.
315 0 : EvtId conjPartId = EvtPDL::chargeConj(particleId);
316 0 : pythiaAliasInt = conjPartId.getAlias();
317 0 : }
318 :
319 0 : std::vector<int> pythiaModes = _pythiaModeMap[pythiaAliasInt];
320 :
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.
324 :
325 0 : std::vector<int>::iterator modeIter;
326 : bool gotMode(false);
327 :
328 0 : for (modeIter = pythiaModes.begin(); modeIter != pythiaModes.end(); ++modeIter) {
329 :
330 : // Stop the loop if we have the right decay mode channel
331 0 : if (gotMode) {break;}
332 :
333 0 : int pythiaModeInt = *modeIter;
334 :
335 0 : EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, pythiaModeInt);
336 :
337 0 : if (decayModel != 0) {
338 :
339 0 : int nModeDaug = decayModel->getNDaug();
340 :
341 : // We need to make sure that the number of daughters match
342 0 : if (nDaughters == nModeDaug) {
343 :
344 : int iModeDaug(0);
345 0 : for (iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++) {
346 :
347 0 : EvtId daugId = decayModel->getDaug(iModeDaug);
348 0 : int daugPDGId = EvtPDL::getStdHep(daugId);
349 : // Pythia has used the right PDG codes for this decay mode, even for conjugate modes
350 0 : int pythiaPDGId = _daugPDGVector[iModeDaug];
351 :
352 0 : if (daugPDGId == pythiaPDGId) {
353 0 : daugAliasIdVect.push_back(daugId);
354 : }
355 :
356 0 : } // Loop over EvtGen mode daughters
357 :
358 0 : int daugAliasSize = daugAliasIdVect.size();
359 0 : if (daugAliasSize == nDaughters) {
360 : // All daughter Id codes are accounted for. Set the flag to stop the loop.
361 : gotMode = true;
362 0 : } else {
363 : // We do not have the correct daughter ordering. Clear the id vector
364 : // and try another mode.
365 0 : daugAliasIdVect.clear();
366 : }
367 :
368 0 : } // Same number of daughters
369 :
370 0 : } // decayModel != 0
371 :
372 : } // Loop over available Pythia modes
373 :
374 0 : if (gotMode == false) {
375 :
376 : // We did not find a match for the daughter aliases. Just use the normal PDG codes
377 : // from the Pythia decay result
378 : int iPyDaug(0);
379 0 : for (iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++) {
380 :
381 0 : int daugPDGCode = _daugPDGVector[iPyDaug];
382 0 : EvtId daugPyId = EvtPDL::evtIdFromStdHep(daugPDGCode);
383 0 : daugAliasIdVect.push_back(daugPyId);
384 :
385 0 : }
386 0 : }
387 :
388 : // Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised.
389 0 : theParent->makeDaughters(nDaughters, daugAliasIdVect);
390 :
391 : // Now set the 4-momenta of the daughters.
392 : int iDaug(0);
393 : // Can use an iterator here, but we already had to use the vector size...
394 0 : for (iDaug = 0; iDaug < nDaughters; iDaug++) {
395 :
396 0 : EvtParticle* theDaughter = theParent->getDaug(iDaug);
397 :
398 : // Set the correct 4-momentum for each daughter particle.
399 0 : if (theDaughter != 0) {
400 0 : EvtId theDaugId = daugAliasIdVect[iDaug];
401 0 : const EvtVector4R theDaugP4 = _daugP4Vector[iDaug];
402 0 : theDaughter->init(theDaugId, theDaugP4);
403 0 : }
404 :
405 : }
406 :
407 0 : }
408 :
409 : void EvtPythiaEngine::updateParticleLists() {
410 :
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.
419 :
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..
423 : int iPDL;
424 0 : int nPDL = EvtPDL::entries();
425 :
426 : // Reset the _addedPDGCodes map that keeps track
427 : // of any new particles added to the Pythia input data stream
428 0 : _addedPDGCodes.clear();
429 :
430 0 : for (iPDL = 0; iPDL < nPDL; iPDL++) {
431 :
432 0 : EvtId particleId = EvtPDL::getEntry(iPDL);
433 0 : int aliasInt = particleId.getAlias();
434 :
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.
438 :
439 : // Should change isJetSet to isPythia eventually.
440 0 : bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia(aliasInt);
441 :
442 0 : if (hasPythiaDecays) {
443 :
444 0 : int isAlias = particleId.isAlias();
445 :
446 0 : int PDGCode = EvtPDL::getStdHep(particleId);
447 :
448 : // Decide what generator to use depending on ether we have
449 : // an alias particle or not
450 0 : _thePythiaGenerator = _genericPythiaGen;
451 0 : _theParticleData = _genericPartData;
452 0 : if (isAlias == 1) {
453 0 : _thePythiaGenerator = _aliasPythiaGen;
454 0 : _theParticleData = _aliasPartData;
455 0 : }
456 :
457 : // Find the Pythia particle name given the standard PDG code integer
458 0 : std::string dataName = _theParticleData.name(PDGCode);
459 : bool alreadyStored(false);
460 0 : if (_addedPDGCodes.find(abs(PDGCode)) != _addedPDGCodes.end()) {alreadyStored = true;}
461 :
462 0 : if (dataName == " " && alreadyStored == false) {
463 :
464 : // Particle and its antiparticle does not exist in the Pythia database.
465 : // Create a new particle, then create the new decay modes.
466 0 : this->createPythiaParticle(particleId, PDGCode);
467 :
468 : } else {
469 :
470 : // Particle exists in the Pythia database.
471 : // Could update mass/lifetime values here. For now just use Pythia defaults.
472 :
473 : }
474 :
475 : // For the particle, create the Pythia decay modes.
476 : // Update Pythia data tables.
477 0 : this->updatePythiaDecayTable(particleId, aliasInt, PDGCode);
478 :
479 0 : } // Loop over Pythia decays
480 :
481 0 : } // Loop over EvtPDL entries
482 :
483 : //report(Severity::Info,"EvtGen")<<"Writing out changed generic Pythia decay list"<<endl;
484 : //_genericPythiaGen->particleData.listChanged();
485 :
486 : //report(Severity::Info,"EvtGen")<<"Writing out changed alias Pythia decay list"<<endl;
487 : //_aliasPythiaGen->particleData.listChanged();
488 :
489 0 : }
490 :
491 : void EvtPythiaEngine::updatePythiaDecayTable(EvtId& particleId, int aliasInt, int PDGCode) {
492 :
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.
499 :
500 0 : int nModes = EvtDecayTable::getInstance()->getNModes(aliasInt);
501 0 : int iMode(0);
502 :
503 : bool firstMode(true);
504 :
505 : // Only process positive PDG codes.
506 0 : if (PDGCode < 0) {return;}
507 :
508 : // Keep track of which decay modes are Pythia decays for each aliasInt
509 0 : std::vector<int> pythiaModes(0);
510 :
511 : // Loop over the decay modes for this particle
512 0 : for (iMode = 0; iMode < nModes; iMode++) {
513 :
514 0 : EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, iMode);
515 :
516 0 : if (decayModel != 0) {
517 :
518 0 : int nDaug = decayModel->getNDaug();
519 :
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).
524 0 : if (nDaug > 0) {
525 :
526 : // Check to see if we have a Pythia decay mode
527 0 : std::string modelName = decayModel->getModelName();
528 :
529 0 : if (modelName == "PYTHIA") {
530 :
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 0 : pythiaModes.push_back(iMode);
534 :
535 0 : std::ostringstream oss;
536 0 : 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.
539 0 : oss << PDGCode;
540 :
541 0 : if (firstMode) {
542 : // Create a new channel
543 0 : oss <<":oneChannel = ";
544 : firstMode = false;
545 0 : } else {
546 : // Add the channel
547 0 : oss <<":addChannel = ";
548 : }
549 :
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).
553 : int onMode(1);
554 0 : oss << onMode << " ";
555 :
556 0 : double BF = decayModel->getBranchingFraction();
557 0 : oss << BF << " ";
558 :
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 0 : int modeInt = this->getModeInt(decayModel);
562 0 : oss << modeInt;
563 :
564 : int iDaug(0);
565 0 : for (iDaug = 0; iDaug < nDaug; iDaug++) {
566 :
567 0 : EvtId daugId = decayModel->getDaug(iDaug);
568 0 : int daugPDG = EvtPDL::getStdHep(daugId);
569 0 : oss << " " << daugPDG;
570 :
571 : } // Daughter list
572 :
573 0 : _thePythiaGenerator->readString(oss.str());
574 :
575 0 : } // is Pythia
576 :
577 0 : } else {
578 :
579 0 : report(Severity::Info,"EvtGen")<<"Warning in EvtPythiaEngine. Trying to redefine Pythia table for "
580 0 : <<EvtPDL::name(particleId)<<" for a decay channel that has no daughters."
581 0 : <<" Keeping Pythia default (if available)."<<endl;
582 :
583 : }
584 :
585 0 : } else {
586 :
587 0 : report(Severity::Info,"EvtGen")<<"Error in EvtPythiaEngine. decayModel is null for particle "
588 0 : <<EvtPDL::name(particleId)<<" mode number "<<iMode<<endl;
589 :
590 : }
591 :
592 : } // Loop over modes
593 :
594 0 : _pythiaModeMap[aliasInt] = pythiaModes;
595 :
596 : // Now, renormalise the decay branching fractions to sum to 1.0
597 0 : std::ostringstream rescaleStr;
598 0 : rescaleStr.setf(std::ios::scientific);
599 0 : rescaleStr << PDGCode << ":rescaleBR = 1.0";
600 :
601 0 : _thePythiaGenerator->readString(rescaleStr.str());
602 :
603 0 : }
604 :
605 : int EvtPythiaEngine::getModeInt(EvtDecayBase* decayModel) {
606 :
607 : int tmpModeInt(0), modeInt(0);
608 :
609 0 : if (decayModel != 0) {
610 :
611 0 : int nVars = decayModel->getNArg();
612 : // Just read the first integer, which specifies the Pythia decay model.
613 : // Ignore any other values.
614 0 : if (nVars > 0) {
615 0 : tmpModeInt = static_cast<int>(decayModel->getArg(0));
616 0 : }
617 0 : }
618 :
619 0 : if (_convertPhysCodes) {
620 :
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.
624 :
625 0 : if (tmpModeInt == 0) {
626 : modeInt = 0; // phase-space
627 0 : } else if (tmpModeInt == 1) {
628 : modeInt = 1; // omega or phi -> 3pi
629 0 : } else if (tmpModeInt == 2) {
630 : modeInt = 11; // Dalitz decay
631 0 : } else if (tmpModeInt == 3) {
632 : modeInt = 2; // V -> PS PS
633 0 : } else if (tmpModeInt == 4) {
634 : modeInt = 92; // onium -> ggg or gg gamma
635 0 : } else if (tmpModeInt == 11) {
636 : modeInt = 42; // phase-space of hadrons from available quarks
637 0 : } else if (tmpModeInt == 12) {
638 : modeInt = 42; // phase-space for onia resonances
639 0 : } else if (tmpModeInt == 13) {
640 : modeInt = 43; // phase-space of at least 3 hadrons
641 0 : } else if (tmpModeInt == 14) {
642 : modeInt = 44; // phase-space of at least 4 hadrons
643 0 : } else if (tmpModeInt == 15) {
644 : modeInt = 45; // phase-space of at least 5 hadrons
645 0 : } else if (tmpModeInt >= 22 && tmpModeInt <= 30) {
646 0 : modeInt = tmpModeInt + 40; // phase space of hadrons with fixed multiplicity (modeInt - 60)
647 0 : } else if (tmpModeInt == 31) {
648 : modeInt = 42; // two or more quarks phase-space; one spectactor quark
649 0 : } else if (tmpModeInt == 32) {
650 : modeInt = 91; // qqbar or gg pair
651 0 : } 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 0 : } else if (tmpModeInt == 41) {
654 : modeInt = 21; // weak decay phase space, weighting nu_tau spectrum
655 0 : } else if (tmpModeInt == 42) {
656 : modeInt = 22; // weak decay V-A matrix element
657 0 : } else if (tmpModeInt == 43) {
658 : modeInt = 22; // weak decay V-A matrix element, quarks as jets (superfluous)
659 0 : } else if (tmpModeInt == 44) {
660 : modeInt = 22; // weak decay V-A matrix element, parton showers (superfluous)
661 0 : } else if (tmpModeInt == 48) {
662 : modeInt = 23; // weak decay V-A matrix element, at least 3 decay products
663 0 : } else if (tmpModeInt == 50) {
664 : modeInt = 0; // default behaviour
665 0 : } else if (tmpModeInt == 51) {
666 : modeInt = 0; // step threshold (channel switched off when mass daughters > mother mass
667 0 : } else if (tmpModeInt == 52 || tmpModeInt == 53) {
668 : modeInt = 0; // beta-factor threshold
669 0 : } else if (tmpModeInt == 84) {
670 : modeInt = 42; // unknown physics process - just use phase-space
671 0 : } else if (tmpModeInt == 101) {
672 : modeInt = 0; // continuation line
673 0 : } else if (tmpModeInt == 102) {
674 : modeInt = 0; // off mass shell particles.
675 0 : } else {
676 0 : report(Severity::Info,"EvtGen")<<"Pythia mode integer "<<tmpModeInt
677 0 : <<" is not recognised. Using phase-space model"<<endl;
678 : modeInt = 0; // Use phase-space for anything else
679 : }
680 :
681 : } else {
682 :
683 : // No need to convert the physics mode integer code
684 : modeInt = tmpModeInt;
685 :
686 : }
687 :
688 0 : return modeInt;
689 :
690 : }
691 :
692 : void EvtPythiaEngine::createPythiaParticle(EvtId& particleId, int PDGCode) {
693 :
694 : // Use the EvtGen name, PDGId and other variables to define the new Pythia particle.
695 0 : EvtId antiPartId = EvtPDL::chargeConj(particleId);
696 :
697 0 : std::string aliasName = EvtPDL::name(particleId); // If not an alias, aliasName = normal name
698 0 : std::string antiName = EvtPDL::name(antiPartId);
699 :
700 0 : EvtSpinType::spintype spinType = EvtPDL::getSpinType(particleId);
701 0 : int spin = EvtSpinType::getSpin2(spinType);
702 :
703 0 : int charge = EvtPDL::chg3(particleId);
704 :
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 0 : int PDGId = EvtPDL::getStdHep(particleId);
709 : int colour(0);
710 0 : if (PDGId == 21) {
711 : colour = 2; // gluons
712 0 : } else if (PDGId <= 8 && PDGId > 0) {
713 : colour = 1; // single quarks
714 0 : }
715 :
716 0 : double m0 = EvtPDL::getMeanMass(particleId);
717 0 : double mWidth = EvtPDL::getWidth(particleId);
718 0 : double mMin = EvtPDL::getMinMass(particleId);
719 0 : double mMax = EvtPDL::getMaxMass(particleId);
720 :
721 0 : double tau0 = EvtPDL::getctau(particleId);
722 :
723 0 : std::ostringstream oss;
724 0 : oss.setf(std::ios::scientific);
725 0 : int absPDGCode = abs(PDGCode);
726 0 : oss << absPDGCode << ":new = " << aliasName << " " << antiName << " "
727 0 : << spin << " " << charge << " " << colour << " "
728 0 : << m0 << " " << mWidth << " " << mMin << " " << mMax << " "
729 0 : << tau0;
730 :
731 : // Pass this information to Pythia
732 0 : _thePythiaGenerator->readString(oss.str());
733 :
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 0 : _addedPDGCodes[absPDGCode] = 1;
740 :
741 0 : }
742 :
743 : void EvtPythiaEngine::updatePhysicsParameters() {
744 :
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.
751 :
752 : // Set the multiplicity level for hadronic weak decays
753 0 : std::string multiWeakCut("ParticleDecays:multIncreaseWeak = 2.0");
754 0 : _genericPythiaGen->readString(multiWeakCut);
755 0 : _aliasPythiaGen->readString(multiWeakCut);
756 :
757 : // Set the multiplicity level for all other decays
758 0 : std::string multiCut("ParticleDecays:multIncrease = 4.5");
759 0 : _genericPythiaGen->readString(multiCut);
760 0 : _aliasPythiaGen->readString(multiCut);
761 :
762 : //Now read in any custom configuration entered in the XML
763 0 : GeneratorCommands commands = EvtExtGeneratorCommandsTable::getInstance()->getCommands("PYTHIA");
764 0 : GeneratorCommands::iterator it = commands.begin();
765 :
766 0 : for( ; it!=commands.end(); it++) {
767 :
768 0 : Command command = *it;
769 0 : std::vector<std::string> commandStrings;
770 :
771 0 : if(command["VERSION"] == "PYTHIA6") {
772 0 : report(Severity::Info,"EvtGen")<<"Converting Pythia 6 command: "<<command["MODULE"]<<"("<<command["PARAM"]<<")="<<command["VALUE"]<<"..."<<endl;
773 0 : commandStrings = convertPythia6Command(command);
774 0 : } else if(command["VERSION"] == "PYTHIA8") {
775 0 : commandStrings.push_back(command["MODULE"]+":"+command["PARAM"]+" = "+command["VALUE"]);
776 : } else {
777 0 : report(Severity::Error, "EvtGen") << "Pythia command received by EvtPythiaEngine has bad version:"<<endl;
778 0 : report(Severity::Error, "EvtGen") << "Received "<<command["VERSION"]<<" but expected PYTHIA6 or PYTHIA8."<<endl;
779 0 : report(Severity::Error, "EvtGen") << "The error is likely to be in EvtDecayTable.cpp"<<endl;
780 0 : report(Severity::Error, "EvtGen") << "EvtGen will now abort."<<endl;
781 0 : ::abort();
782 : }
783 0 : std::string generator = command["GENERATOR"];
784 0 : if(generator == "GENERIC" || generator == "Generic" || generator == "generic" ||
785 0 : generator == "BOTH" || generator == "Both" || generator == "both") {
786 0 : std::vector<std::string>::iterator it2 = commandStrings.begin();
787 0 : for( ; it2!=commandStrings.end(); it2++) {
788 0 : report(Severity::Info,"EvtGen")<<"Configuring generic Pythia generator: " << (*it2) << endl;
789 0 : _genericPythiaGen->readString(*it2);
790 : }
791 0 : }
792 0 : if(generator == "ALIAS" || generator == "Alias" || generator == "alias" ||
793 0 : generator == "BOTH" || generator == "Both" || generator == "both") {
794 0 : std::vector<std::string>::iterator it2 = commandStrings.begin();
795 0 : for( ; it2!=commandStrings.end(); it2++) {
796 0 : report(Severity::Info,"EvtGen")<<"Configuring alias Pythia generator: " << (*it2) << endl;
797 0 : _aliasPythiaGen->readString(*it2);
798 : }
799 0 : }
800 0 : }
801 0 : }
802 :
803 : #endif
|
