]>
Commit | Line | Data |
---|---|---|
0ca57c2f | 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 | #include "pythia8175/include/Event.h" | |
34 | ||
35 | #include <iostream> | |
36 | #include <sstream> | |
37 | ||
38 | using std::endl; | |
39 | ||
40 | EvtPythiaEngine::EvtPythiaEngine(std::string xmlDir, bool convertPhysCodes, | |
41 | 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 | report(INFO,"EvtGen")<<"Creating generic Pythia generator"<<endl; | |
52 | _genericPythiaGen = new Pythia8::Pythia(xmlDir); | |
53 | _genericPartData = _genericPythiaGen->particleData; | |
54 | ||
55 | report(INFO,"EvtGen")<<"Creating alias Pythia generator"<<endl; | |
56 | _aliasPythiaGen = new Pythia8::Pythia(xmlDir); | |
57 | _aliasPartData = _aliasPythiaGen->particleData; | |
58 | ||
59 | _thePythiaGenerator = 0; | |
60 | _daugPDGVector.clear(); _daugP4Vector.clear(); | |
61 | ||
62 | _convertPhysCodes = convertPhysCodes; | |
63 | ||
64 | // Specify if we are going to use the random number generator (engine) | |
65 | // from EvtGen for Pythia 8. | |
66 | _useEvtGenRandom = useEvtGenRandom; | |
67 | ||
68 | _evtgenRandom = new EvtPythiaRandom(); | |
69 | ||
70 | _initialised = false; | |
71 | ||
72 | } | |
73 | ||
74 | EvtPythiaEngine::~EvtPythiaEngine() { | |
75 | ||
76 | delete _genericPythiaGen; _genericPythiaGen = 0; | |
77 | delete _aliasPythiaGen; _aliasPythiaGen = 0; | |
78 | ||
79 | delete _evtgenRandom; _evtgenRandom = 0; | |
80 | ||
81 | _thePythiaGenerator = 0; | |
82 | ||
83 | this->clearDaughterVectors(); | |
84 | this->clearPythiaModeMap(); | |
85 | ||
86 | } | |
87 | ||
88 | void EvtPythiaEngine::clearDaughterVectors() { | |
89 | _daugPDGVector.clear(); | |
90 | _daugP4Vector.clear(); | |
91 | } | |
92 | ||
93 | void EvtPythiaEngine::clearPythiaModeMap() { | |
94 | ||
95 | PythiaModeMap::iterator iter; | |
96 | for (iter = _pythiaModeMap.begin(); iter != _pythiaModeMap.end(); ++iter) { | |
97 | ||
98 | std::vector<int> modeVector = iter->second; | |
99 | modeVector.clear(); | |
100 | ||
101 | } | |
102 | ||
103 | _pythiaModeMap.clear(); | |
104 | ||
105 | } | |
106 | ||
107 | void EvtPythiaEngine::initialise() { | |
108 | ||
109 | if (_initialised) {return;} | |
110 | ||
111 | this->clearPythiaModeMap(); | |
112 | ||
113 | 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 | _genericPythiaGen->readString("ProcessLevel:all = off"); | |
118 | _aliasPythiaGen->readString("ProcessLevel:all = off"); | |
119 | ||
120 | // Apply any other physics (or special particle) requirements/cuts etc.. | |
121 | this->updatePhysicsParameters(); | |
122 | ||
123 | // Set the random number generator | |
124 | if (_useEvtGenRandom == true) { | |
125 | ||
126 | _genericPythiaGen->setRndmEnginePtr(_evtgenRandom); | |
127 | _aliasPythiaGen->setRndmEnginePtr(_evtgenRandom); | |
128 | ||
129 | } | |
130 | ||
131 | _genericPythiaGen->init(); | |
132 | _aliasPythiaGen->init(); | |
133 | ||
134 | _initialised = true; | |
135 | ||
136 | } | |
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 | if (_initialised == false) {this->initialise();} | |
156 | ||
157 | if (theParticle == 0) { | |
158 | report(INFO,"EvtGen")<<"Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay."<<endl; | |
159 | return false; | |
160 | } | |
161 | ||
162 | // Delete EvtParticle daughters if they already exist | |
163 | if (theParticle->getNDaug() != 0) { | |
164 | bool keepChannel(false); | |
165 | theParticle->deleteDaughters(keepChannel); | |
166 | } | |
167 | ||
168 | EvtId particleId = theParticle->getId(); | |
169 | int isAlias = particleId.isAlias(); | |
170 | ||
171 | // Choose the generator depending if we have an aliased (parent) particle or not | |
172 | _thePythiaGenerator = _genericPythiaGen; | |
173 | 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 | Pythia8::Event& theEvent = _thePythiaGenerator->event; | |
180 | theEvent.reset(); | |
181 | ||
182 | // Initialise the event to be the particle rest frame | |
183 | 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 | double m0 = theParticle->mass(); | |
189 | double E = m0; | |
190 | ||
191 | 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 | for (iTrial = 0; iTrial < 10; iTrial++) { | |
197 | ||
198 | generatedEvent = _thePythiaGenerator->next(); | |
199 | if (generatedEvent) {break;} | |
200 | ||
201 | } | |
202 | ||
203 | bool success(false); | |
204 | ||
205 | 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 | 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 | 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 | this->createDaughterEvtParticles(theParticle); | |
225 | ||
226 | //theParticle->printTree(); | |
227 | //theEvent.list(true, true); | |
228 | ||
229 | success = true; | |
230 | ||
231 | } | |
232 | ||
233 | return success; | |
234 | ||
235 | } | |
236 | ||
237 | void EvtPythiaEngine::storeDaughterInfo(EvtParticle* theParticle, int startInt) { | |
238 | ||
239 | Pythia8::Event& theEvent = _thePythiaGenerator->event; | |
240 | ||
241 | std::vector<int> daugList = theEvent.daughterList(startInt); | |
242 | ||
243 | std::vector<int>::iterator daugIter; | |
244 | for (daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter) { | |
245 | ||
246 | int daugInt = *daugIter; | |
247 | ||
248 | // Ask if the daughter is a quark or gluon. If so, recursively search again. | |
249 | Pythia8::Particle& daugParticle = theEvent[daugInt]; | |
250 | ||
251 | if (daugParticle.isQuark() || daugParticle.isGluon()) { | |
252 | ||
253 | // Recursively search for correct daughter type | |
254 | 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 | int statusCode = daugParticle.status(); | |
266 | if (statusCode != 1000) { | |
267 | ||
268 | int daugPDGInt = daugParticle.id(); | |
269 | ||
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); | |
275 | ||
276 | // Now store the EvtId and 4-momentum in the internal vectors | |
277 | _daugPDGVector.push_back(daugPDGInt); | |
278 | _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 | daugParticle.status(1000); | |
283 | ||
284 | } // Status code != 1000 | |
285 | ||
286 | } | |
287 | ||
288 | } | |
289 | ||
290 | } | |
291 | ||
292 | void EvtPythiaEngine::createDaughterEvtParticles(EvtParticle* theParent) { | |
293 | ||
294 | if (theParent == 0) { | |
295 | report(INFO,"EvtGen")<<"Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null"<<endl; | |
296 | 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 | int nDaughters = _daugPDGVector.size(); | |
304 | std::vector<EvtId> daugAliasIdVect(0); | |
305 | ||
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); | |
312 | ||
313 | if (PDGId < 0) { | |
314 | // We have an anti-particle. | |
315 | EvtId conjPartId = EvtPDL::chargeConj(particleId); | |
316 | pythiaAliasInt = conjPartId.getAlias(); | |
317 | } | |
318 | ||
319 | 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 | std::vector<int>::iterator modeIter; | |
326 | bool gotMode(false); | |
327 | ||
328 | for (modeIter = pythiaModes.begin(); modeIter != pythiaModes.end(); ++modeIter) { | |
329 | ||
330 | // Stop the loop if we have the right decay mode channel | |
331 | if (gotMode) {break;} | |
332 | ||
333 | int pythiaModeInt = *modeIter; | |
334 | ||
335 | EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, pythiaModeInt); | |
336 | ||
337 | if (decayModel != 0) { | |
338 | ||
339 | int nModeDaug = decayModel->getNDaug(); | |
340 | ||
341 | // We need to make sure that the number of daughters match | |
342 | if (nDaughters == nModeDaug) { | |
343 | ||
344 | int iModeDaug(0); | |
345 | for (iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++) { | |
346 | ||
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]; | |
351 | ||
352 | if (daugPDGId == pythiaPDGId) { | |
353 | daugAliasIdVect.push_back(daugId); | |
354 | } | |
355 | ||
356 | } // Loop over EvtGen mode daughters | |
357 | ||
358 | int daugAliasSize = daugAliasIdVect.size(); | |
359 | if (daugAliasSize == nDaughters) { | |
360 | // All daughter Id codes are accounted for. Set the flag to stop the loop. | |
361 | gotMode = true; | |
362 | } else { | |
363 | // We do not have the correct daughter ordering. Clear the id vector | |
364 | // and try another mode. | |
365 | daugAliasIdVect.clear(); | |
366 | } | |
367 | ||
368 | } // Same number of daughters | |
369 | ||
370 | } // decayModel != 0 | |
371 | ||
372 | } // Loop over available Pythia modes | |
373 | ||
374 | 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 | for (iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++) { | |
380 | ||
381 | int daugPDGCode = _daugPDGVector[iPyDaug]; | |
382 | EvtId daugPyId = EvtPDL::evtIdFromStdHep(daugPDGCode); | |
383 | daugAliasIdVect.push_back(daugPyId); | |
384 | ||
385 | } | |
386 | } | |
387 | ||
388 | // Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised. | |
389 | 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 | for (iDaug = 0; iDaug < nDaughters; iDaug++) { | |
395 | ||
396 | EvtParticle* theDaughter = theParent->getDaug(iDaug); | |
397 | ||
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); | |
403 | } | |
404 | ||
405 | } | |
406 | ||
407 | } | |
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 | 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 | _addedPDGCodes.clear(); | |
429 | ||
430 | for (iPDL = 0; iPDL < nPDL; iPDL++) { | |
431 | ||
432 | EvtId particleId = EvtPDL::getEntry(iPDL); | |
433 | 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 | bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia(aliasInt); | |
441 | ||
442 | if (hasPythiaDecays) { | |
443 | ||
444 | int isAlias = particleId.isAlias(); | |
445 | ||
446 | int PDGCode = EvtPDL::getStdHep(particleId); | |
447 | ||
448 | // Decide what generator to use depending on ether we have | |
449 | // an alias particle or not | |
450 | _thePythiaGenerator = _genericPythiaGen; | |
451 | _theParticleData = _genericPartData; | |
452 | if (isAlias == 1) { | |
453 | _thePythiaGenerator = _aliasPythiaGen; | |
454 | _theParticleData = _aliasPartData; | |
455 | } | |
456 | ||
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;} | |
461 | ||
462 | 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 | 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 | this->updatePythiaDecayTable(particleId, aliasInt, PDGCode); | |
478 | ||
479 | } // Loop over Pythia decays | |
480 | ||
481 | } // Loop over EvtPDL entries | |
482 | ||
483 | //report(INFO,"EvtGen")<<"Writing out changed generic Pythia decay list"<<endl; | |
484 | //_genericPythiaGen->particleData.listChanged(); | |
485 | ||
486 | //report(INFO,"EvtGen")<<"Writing out changed alias Pythia decay list"<<endl; | |
487 | //_aliasPythiaGen->particleData.listChanged(); | |
488 | ||
489 | } | |
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 | int nModes = EvtDecayTable::getInstance()->getNModes(aliasInt); | |
501 | int iMode(0); | |
502 | ||
503 | bool firstMode(true); | |
504 | ||
505 | // Only process positive PDG codes. | |
506 | if (PDGCode < 0) {return;} | |
507 | ||
508 | // Keep track of which decay modes are Pythia decays for each aliasInt | |
509 | std::vector<int> pythiaModes(0); | |
510 | ||
511 | // Loop over the decay modes for this particle | |
512 | for (iMode = 0; iMode < nModes; iMode++) { | |
513 | ||
514 | EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, iMode); | |
515 | ||
516 | if (decayModel != 0) { | |
517 | ||
518 | 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 | if (nDaug > 0) { | |
525 | ||
526 | // Check to see if we have a Pythia decay mode | |
527 | std::string modelName = decayModel->getModelName(); | |
528 | ||
529 | 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 | pythiaModes.push_back(iMode); | |
534 | ||
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. | |
539 | oss << PDGCode; | |
540 | ||
541 | if (firstMode) { | |
542 | // Create a new channel | |
543 | oss <<":oneChannel = "; | |
544 | firstMode = false; | |
545 | } else { | |
546 | // Add the channel | |
547 | 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 | oss << onMode << " "; | |
555 | ||
556 | double BF = decayModel->getBranchingFraction(); | |
557 | 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 | int modeInt = this->getModeInt(decayModel); | |
562 | oss << modeInt; | |
563 | ||
564 | int iDaug(0); | |
565 | for (iDaug = 0; iDaug < nDaug; iDaug++) { | |
566 | ||
567 | EvtId daugId = decayModel->getDaug(iDaug); | |
568 | int daugPDG = EvtPDL::getStdHep(daugId); | |
569 | oss << " " << daugPDG; | |
570 | ||
571 | } // Daughter list | |
572 | ||
573 | _thePythiaGenerator->readString(oss.str()); | |
574 | ||
575 | } // is Pythia | |
576 | ||
577 | } else { | |
578 | ||
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; | |
582 | ||
583 | } | |
584 | ||
585 | } else { | |
586 | ||
587 | report(INFO,"EvtGen")<<"Error in EvtPythiaEngine. decayModel is null for particle " | |
588 | <<EvtPDL::name(particleId)<<" mode number "<<iMode<<endl; | |
589 | ||
590 | } | |
591 | ||
592 | } // Loop over modes | |
593 | ||
594 | _pythiaModeMap[aliasInt] = pythiaModes; | |
595 | ||
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"; | |
600 | ||
601 | _thePythiaGenerator->readString(rescaleStr.str()); | |
602 | ||
603 | } | |
604 | ||
605 | int EvtPythiaEngine::getModeInt(EvtDecayBase* decayModel) { | |
606 | ||
607 | int tmpModeInt(0), modeInt(0); | |
608 | ||
609 | if (decayModel != 0) { | |
610 | ||
611 | int nVars = decayModel->getNArg(); | |
612 | // Just read the first integer, which specifies the Pythia decay model. | |
613 | // Ignore any other values. | |
614 | if (nVars > 0) { | |
615 | tmpModeInt = static_cast<int>(decayModel->getArg(0)); | |
616 | } | |
617 | } | |
618 | ||
619 | 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 | 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. | |
675 | } else { | |
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 | |
679 | } | |
680 | ||
681 | } else { | |
682 | ||
683 | // No need to convert the physics mode integer code | |
684 | modeInt = tmpModeInt; | |
685 | ||
686 | } | |
687 | ||
688 | 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 | EvtId antiPartId = EvtPDL::chargeConj(particleId); | |
696 | ||
697 | std::string aliasName = EvtPDL::name(particleId); // If not an alias, aliasName = normal name | |
698 | std::string antiName = EvtPDL::name(antiPartId); | |
699 | ||
700 | EvtSpinType::spintype spinType = EvtPDL::getSpinType(particleId); | |
701 | int spin = EvtSpinType::getSpin2(spinType); | |
702 | ||
703 | 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 | int PDGId = EvtPDL::getStdHep(particleId); | |
709 | int colour(0); | |
710 | if (PDGId == 21) { | |
711 | colour = 2; // gluons | |
712 | } else if (PDGId <= 8 && PDGId > 0) { | |
713 | colour = 1; // single quarks | |
714 | } | |
715 | ||
716 | double m0 = EvtPDL::getMeanMass(particleId); | |
717 | double mWidth = EvtPDL::getWidth(particleId); | |
718 | double mMin = EvtPDL::getMinMass(particleId); | |
719 | double mMax = EvtPDL::getMaxMass(particleId); | |
720 | ||
721 | double tau0 = EvtPDL::getctau(particleId); | |
722 | ||
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 << " " | |
729 | << tau0; | |
730 | ||
731 | // Pass this information to Pythia | |
732 | _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 | _addedPDGCodes[absPDGCode] = 1; | |
740 | ||
741 | } | |
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 | std::string multiWeakCut("ParticleDecays:multIncreaseWeak = 2.0"); | |
754 | _genericPythiaGen->readString(multiWeakCut); | |
755 | _aliasPythiaGen->readString(multiWeakCut); | |
756 | ||
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); | |
761 | ||
762 | //Now read in any custom configuration entered in the XML | |
763 | GeneratorCommands commands = EvtExtGeneratorCommandsTable::getInstance()->getCommands("PYTHIA"); | |
764 | GeneratorCommands::iterator it = commands.begin(); | |
765 | ||
766 | for( ; it!=commands.end(); it++) { | |
767 | ||
768 | Command command = *it; | |
769 | std::vector<std::string> commandStrings; | |
770 | ||
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"]); | |
776 | } else { | |
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; | |
781 | ::abort(); | |
782 | } | |
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); | |
790 | } | |
791 | } | |
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); | |
798 | } | |
799 | } | |
800 | } | |
801 | } | |
802 | ||
803 | #endif |