1 // Pythia.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2008 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 // Function definitions (not found in the header) for the Pythia class.
10 // Access time information.
13 // Allow string and character manipulation.
18 //**************************************************************************
24 // Constants: could be changed here if desired, but normally should not.
25 // These are of technical nature, as described for each.
27 // Maximum number of tries to produce parton level from given input.
28 const int Pythia::NTRY = 10;
30 // Negative integer to denote that no subrun has been set.
31 const int Pythia::SUBRUNDEFAULT = -999;
37 Pythia::Pythia(string xmlDir) {
39 // Initial values for pointers to PDF's.
42 useNewPdfHard = false;
48 // Initial values for pointers to Les Houches Event objects.
53 // Initial value for pointer to external decay handler.
56 // Initial value for pointer to user hooks.
59 // Initial value for pointer to beam shape.
60 useNewBeamShape = false;
63 // Initial values for pointers to timelike and spacelike showers.
70 // Send Info pointer to handle error printout/counting correctly.
71 settings.initPtr( &info);
72 ParticleDataEntry::initPtr( &info);
73 particleData.initPtr( &info);
75 // Find path to data files, i.e. xmldoc directory location.
76 // Environment variable takes precedence, else use constructor input.
78 const char* PYTHIA8DATA = "PYTHIA8DATA";
79 char* envPath = getenv(PYTHIA8DATA);
80 if (envPath != 0 && *envPath != '\0') {
82 while (*(envPath+i) != '\0') path += *(envPath+(i++));
85 if (path[ path.length() - 1 ] != '/') path += "/";
87 // Read in files with all flags, modes, parms and words.
88 string initFile = path + "Index.xml";
89 isConstructed = settings.init( initFile);
91 info.errorMsg("Abort from Pythia::Pythia: settings unavailable");
95 // Read in files with all particle data.
96 string dataFile = path + "ParticleData.xml";
97 isConstructed = particleData.init( dataFile);
99 info.errorMsg("Abort from Pythia::Pythia: particle data unavailable");
103 // Write the Pythia banner to output.
106 // Set headers to distinguish the two event listing kinds.
107 process.init("(hard process)");
108 event.init("(complete event)");
110 // Not initialized until at the end of initInternal.
121 // Delete the PDF's created with new.
122 if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr;
123 if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr;
124 if (useNewPdfA) delete pdfAPtr;
125 if (useNewPdfB) delete pdfBPtr;
127 // Delete the Les Houches object created with new.
128 if (useNewLHA) delete lhaUpPtr;
130 // Delete the BeamShape object created with new.
131 if (useNewBeamShape) delete beamShapePtr;
133 // Delete the timelike and spacelike showers created with new.
134 if (useNewTimes) delete timesPtr;
135 if (useNewSpace) delete spacePtr;
141 // Read in one update for a setting or particle data from a single line.
143 bool Pythia::readString(string line, bool warn) {
145 // Check that constructor worked.
146 if (!isConstructed) return false;
148 // If empty line then done.
149 if (line.find_first_not_of(" ") == string::npos) return true;
151 // If first character is not a letter/digit, then taken to be a comment.
152 int firstChar = line.find_first_not_of(" ");
153 if (!isalnum(line[firstChar])) return true;
155 // Send on particle data to the ParticleData database.
156 if (isdigit(line[firstChar]))
157 return particleData.readString(line, warn);
159 // Everything else sent on to Settings.
160 return settings.readString(line, warn);
166 // Read in updates for settings or particle data from user-defined file.
168 bool Pythia::readFile(string fileName, bool warn, int subrun) {
170 // Check that constructor worked.
171 if (!isConstructed) return false;
173 // Open file with updates.
174 const char* cstring = fileName.c_str();
175 ifstream is(cstring);
177 info.errorMsg("Error in Pythia::readFile: did not find file", fileName);
181 // Read in one line at a time.
183 bool accepted = true;
184 int subrunNow = SUBRUNDEFAULT;
185 while ( getline(is, line) ) {
187 // Check whether entered new subrun.
188 int subrunLine = readSubrun( line, warn);
189 if (subrunLine >= 0) subrunNow = subrunLine;
191 // Process the line if in correct subrun.
192 if ( (subrunNow == subrun || subrunNow == SUBRUNDEFAULT)
193 && !readString( line, warn) ) accepted = false;
195 // Reached end of input file.
202 // Routine to pass in pointers to PDF's. Usage optional.
204 bool Pythia::setPDFPtr( PDF* pdfAPtrIn, PDF* pdfBPtrIn, PDF* pdfHardAPtrIn,
205 PDF* pdfHardBPtrIn) {
207 // The routine can have no effect if PDF's already assigned.
208 if (pdfAPtr != 0 || pdfBPtr != 0) return false;
210 // The two PDF objects cannot be one and the same, or unassigned.
211 if (pdfAPtrIn == pdfBPtrIn || pdfAPtrIn == 0 || pdfBPtrIn == 0) return false;
217 // By default same pointers for hard-process PDF's.
218 pdfHardAPtr = pdfAPtrIn;
219 pdfHardBPtr = pdfBPtrIn;
221 // Optionally allow separate pointers for hard process.
222 if (pdfHardAPtrIn == 0 || pdfHardBPtrIn == 0) return true;
223 if (pdfHardAPtrIn == pdfHardBPtrIn) return false;
224 pdfHardAPtr = pdfHardAPtrIn;
225 pdfHardBPtr = pdfHardBPtrIn;
233 // Routine to initialize with CM energy.
235 bool Pythia::init( int idAin, int idBin, double eCMin) {
237 // Read in and set values.
244 // Send on to common initialization.
245 bool status = initInternal();
246 if (!status) info.errorMsg("Abort from Pythia::init: initialization failed");
253 // Routine to initialize with two collinear beams, energies specified.
255 bool Pythia::init( int idAin, int idBin, double eAin, double eBin) {
257 // Read in and set values.
265 // Send on to common initialization.
266 bool status = initInternal();
267 if (!status) info.errorMsg("Abort from Pythia::init: initialization failed");
274 // Routine to initialize with two beams specified by three-momenta.
276 bool Pythia::init( int idAin, int idBin, double pxAin, double pyAin,
277 double pzAin, double pxBin, double pyBin, double pzBin) {
279 // Read in and set values.
291 // Send on to common initialization.
292 bool status = initInternal();
293 if (!status) info.errorMsg("Abort from Pythia::init: initialization failed");
300 // Routine to initialize when all info is given in a Les Houches Event File.
302 bool Pythia::init( string LesHouchesEventFile, bool skipInit) {
304 // Destroy any previous LHAup object.
305 if (useNewLHA) delete lhaUpPtr;
307 // Create LHAup object. Send in pointer to info.
308 const char* cstring = LesHouchesEventFile.c_str();
309 lhaUpPtr = new LHAupLHEF(cstring);
310 lhaUpPtr->setPtr( &info);
314 // Store or replace LHA pointer in other classes.
315 processLevel.setLHAPtr( lhaUpPtr);
317 // If second time around, only with new file, then simplify.
318 if (skipInit) return true;
320 // Set LHAinit information (in some external program).
321 if (!lhaUpPtr->setInit()) {
322 info.errorMsg("Abort from Pythia::init: "
323 "Les Houches initialization failed");
327 // Extract beams from values set in an LHAinit object.
328 idA = lhaUpPtr->idBeamA();
329 idB = lhaUpPtr->idBeamB();
330 eA = lhaUpPtr->eBeamA();
331 eB = lhaUpPtr->eBeamB();
334 // Now do normal initialization. List info if there.
335 bool status = initInternal();
336 lhaUpPtr->listInit();
337 if (!status) info.errorMsg("Abort from Pythia::init: initialization failed");
344 // Routine to initialize with the variable values of the Beams kind.
346 bool Pythia::init() {
348 // Check if to read from Les Houches Event File set, and is so send on.
349 if (mode("Beams:frameType") == 4) {
350 string lhef = word("Beams:LHEF");
351 bool skipInit = flag("Main:LHEFskipInit");
352 return init( lhef, skipInit);
355 // Read in and set values.
356 idA = mode("Beams:idA");
357 idB = mode("Beams:idB");
358 frameType = mode("Beams:frameType");
359 eCM = parm("Beams:eCM");
360 eA = parm("Beams:eA");
361 eB = parm("Beams:eB");
362 pxA = parm("Beams:pxA");
363 pyA = parm("Beams:pyA");
364 pzA = parm("Beams:pzA");
365 pxB = parm("Beams:pxB");
366 pyB = parm("Beams:pyB");
367 pzB = parm("Beams:pzB");
370 // Send on to common initialization.
371 bool status = initInternal();
372 if (!status) info.errorMsg("Abort from Pythia::init: initialization failed");
379 // Routine to initialize when beam info is given in an LHAup object.
381 bool Pythia::init( LHAup* lhaUpPtrIn) {
383 // Save and set flag for subsequent usage of LHAup object.
384 lhaUpPtr = lhaUpPtrIn;
387 // Store LHA pointer in other classes.
388 processLevel.setLHAPtr( lhaUpPtr);
390 // Set LHAinit information (in some external program).
391 if (!lhaUpPtr->setInit()) {
392 info.errorMsg("Abort from Pythia::init: "
393 "Les Houches initialization failed");
397 // Extract beams from values set in an LHAinit object.
398 idA = lhaUpPtr->idBeamA();
399 idB = lhaUpPtr->idBeamB();
400 eA = lhaUpPtr->eBeamA();
401 eB = lhaUpPtr->eBeamB();
404 // Now do normal initialization. List info if there.
405 bool status = initInternal();
406 lhaUpPtr->listInit();
407 if (!status) info.errorMsg("Abort from Pythia::init: initialization failed");
414 // Main routine to initialize the generation process.
415 // (The alternative init forms end up in this one.)
417 bool Pythia::initInternal() {
419 // Check that constructor worked.
421 if (!isConstructed) return false;
423 // Reset error counters.
427 // Initialize data members extracted from database.
428 doProcessLevel = settings.flag("ProcessLevel:all");
429 doPartonLevel = settings.flag("PartonLevel:all") && doProcessLevel;
430 doHadronLevel = settings.flag("HadronLevel:all");
431 doMomentumSpread = settings.flag("Beams:allowMomentumSpread");
432 doVertexSpread = settings.flag("Beams:allowVertexSpread");
433 checkEvent = settings.flag("Check:event");
434 nErrList = settings.mode("Check:nErrList");
435 epTolErr = settings.parm("Check:epTolErr");
436 epTolWarn = settings.parm("Check:epTolWarn");
438 // Initialize the random number generator.
439 if ( settings.flag("Random:setSeed") )
440 Rndm::init( settings.mode("Random:seed") );
442 // Initialize tunes to e+e- and pp/ppbar data.
445 // Initialize couplings (needed to initialize resonances).
446 AlphaEM::initStatic();
447 CoupEW::initStatic();
450 // Initialize some aspects of particle data, including resonances.
451 ParticleDataEntry::initStatic();
452 particleData.initBWmass();
453 particleData.initResonances(resonancePtrs);
455 // Set up values related to user hooks.
456 hasUserHooks = (userHooksPtr > 0);
457 doVetoProcess = (hasUserHooks)
458 ? userHooksPtr->canVetoProcessLevel() : false;
459 doVetoPartons = (hasUserHooks)
460 ? userHooksPtr->canVetoPartonLevel() : false;
462 // Set up values related to beam shape.
463 if (beamShapePtr == 0) {
464 beamShapePtr = new BeamShape();
465 useNewBeamShape = true;
467 beamShapePtr->init();
469 // Set up objects for timelike and spacelike showers.
470 if (timesDecPtr == 0 || timesPtr == 0) {
471 TimeShower* timesNow = new TimeShower();
472 if (timesDecPtr == 0) timesDecPtr = timesNow;
473 if (timesPtr == 0) timesPtr = timesNow;
477 spacePtr = new SpaceShower();
481 // Initialize showers, especially for simple showers in decays.
482 timesPtr->initPtr( &info);
483 timesDecPtr->initPtr( &info);
484 spacePtr->initPtr( &info);
485 timesDecPtr->init( 0, 0);
487 // Check that beams and beam combination can be handled.
488 // Only allow neutrinos as beams when leptons unresolved.
489 bool canHandleBeams = false;
490 int idAabs = abs(idA);
491 int idBabs = abs(idB);
492 if (doProcessLevel) {
493 if (idAabs == 2212 && idBabs == 2212) canHandleBeams = true;
494 else if ( idAabs == idBabs && (idAabs == 11 || idAabs == 13
495 || idAabs == 15) ) canHandleBeams = true;
496 else if ( idAabs > 10 && idAabs < 17 && idA * idB < 0
497 && !settings.flag("PDF:lepton") ) {
498 if (idAabs == idBabs) canHandleBeams = true;
499 int idMax = max(idAabs, idBabs);
500 int idMin = min(idAabs, idBabs);
501 if (idMax - idMin == 1 && idMax%2 == 0) canHandleBeams = true;
503 if (!canHandleBeams) {
504 info.errorMsg("Error in Pythia::init: "
505 "cannot handle this beam combination");
510 // Do not set up beam kinematics when no process level.
511 if (!doProcessLevel) frameType = 1;
514 // Set up beam kinematics.
515 if (!initKinematics()) return false;
517 // Set up the PDF's, if not already done.
519 pdfAPtr = getPDFPtr(idA);
520 if (!pdfAPtr->isSetup()) return false;
521 pdfHardAPtr = pdfAPtr;
525 pdfBPtr = getPDFPtr(idB);
526 if (!pdfBPtr->isSetup()) return false;
527 pdfHardBPtr = pdfBPtr;
531 // Optionally set up separate PDF's for hard process.
532 if (settings.flag("PDF:useHard")) {
533 pdfHardAPtr = getPDFPtr(idA, 2);
534 if (!pdfHardAPtr->isSetup()) return false;
535 pdfHardBPtr = getPDFPtr(idB, 2);
536 if (!pdfHardBPtr->isSetup()) return false;
537 useNewPdfHard = true;
540 // Set up the two beams and the common remnant system.
541 StringFlav* flavSel = hadronLevel.getStringFlavPtr();
542 bool isUnresolvedA = ( ParticleDataTable::isLepton(idA)
543 && !settings.flag("PDF:lepton") );
544 bool isUnresolvedB = ( ParticleDataTable::isLepton(idB)
545 && !settings.flag("PDF:lepton") );
546 beamA.init( idA, pzAcm, eA, mA, &info, pdfAPtr, pdfHardAPtr,
547 isUnresolvedA, flavSel);
548 beamB.init( idB, pzBcm, eB, mB, &info, pdfBPtr, pdfHardBPtr,
549 isUnresolvedB, flavSel);
552 // Send info/pointers to process level for initialization.
553 if ( doProcessLevel && !processLevel.init( &info, &beamA, &beamB,
554 &sigmaTot, doLHA, &slha, userHooksPtr, sigmaPtrs) ) return false;
556 // Send info/pointers to parton level for initialization.
557 if ( doPartonLevel && !partonLevel.init( &info, &beamA, &beamB,
558 &sigmaTot, timesDecPtr, timesPtr, spacePtr, userHooksPtr) )
561 // Send info/pointers to hadron level for initialization.
562 // Note: forceHadronLevel() can come, so we must always initialize.
563 if ( !hadronLevel.init( &info, timesDecPtr, decayHandlePtr,
564 handledParticles) ) return false;
566 // Optionally check particle data table for inconsistencies.
567 if ( settings.flag("Check:particleData") )
568 particleData.checkTable( settings.mode("Check:levelParticleData") );
577 // Calculate kinematics at initialization. Store beam four-momenta.
579 bool Pythia::initKinematics() {
581 // Find masses. Initial guess that we are in CM frame.
582 mA = ParticleDataTable::m0(idA);
583 mB = ParticleDataTable::m0(idB);
587 // Collinear beams not in CM frame: find CM energy.
588 if (frameType == 2) {
591 pzA = sqrt(eA*eA - mA*mA);
592 pzB = -sqrt(eB*eB - mB*mB);
593 pAinit = Vec4( 0., 0., pzA, eA);
594 pBinit = Vec4( 0., 0., pzB, eB);
595 eCM = sqrt( pow2(eA + eB) - pow2(pzA + pzB) );
597 // Find boost to rest frame.
598 betaZ = (pzA + pzB) / (eA + eB);
599 gammaZ = (eA + eB) / eCM;
600 if (abs(betaZ) < 1e-10) frameType = 1;
603 // Completely general beam directions: find CM energy.
604 else if (frameType == 3) {
605 eA = sqrt( pxA*pxA + pyA*pyA + pzA*pzA + mA*mA);
606 eB = sqrt( pxB*pxB + pyB*pyB + pzB*pzB + mB*mB);
607 pAinit = Vec4( pxA, pyA, pzA, eA);
608 pBinit = Vec4( pxB, pyB, pzB, eB);
609 eCM = (pAinit + pBinit).mCalc();
611 // Find boost+rotation needed to move from/to CM frame.
613 MfromCM.fromCMframe( pAinit, pBinit);
618 // Fail if CM energy below beam masses.
619 if (eCM < mA + mB) return false;
621 // Set up CM-frame kinematics with beams along +-z axis.
622 pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
623 * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
625 eA = sqrt(mA*mA + pzAcm*pzAcm);
626 eB = sqrt(mB*mB + pzBcm*pzBcm);
628 // If in CM frame then store beam four-vectors (else already done above).
629 if (frameType != 2 && frameType != 3) {
630 pAinit = Vec4( 0., 0., pzAcm, eA);
631 pBinit = Vec4( 0., 0., pzBcm, eB);
634 // Store main info for access in process generation.
635 info.setBeamA( idA, pzAcm, eA, mA);
636 info.setBeamB( idB, pzBcm, eB, mB);
639 // Must allow for generic boost+rotation when beam momentum spread.
640 if (doMomentumSpread) frameType = 3;
649 // Initialize tunes to e+e- and pp/ppbar data.
651 void Pythia::initTunes() {
653 // Modes to use. Fast return if all is default.
654 int eeTune = settings.mode("Tune:ee");
655 int ppTune = settings.mode("Tune:pp");
656 if (eeTune == 0 && ppTune == 0) return;
658 // Marc Montull's tune to particle composition at LEP1.
660 settings.parm("StringZ:aLund", 0.76 );
661 settings.parm("StringZ:bLund", 0.58 ); // default
662 settings.parm("StringFlav:probStoUD", 0.22 );
663 settings.parm("StringFlav:probQQtoQ", 0.08 );
664 settings.parm("StringFlav:probSQtoQQ", 0.75 );
665 settings.parm("StringFlav:probQQ1toQQ0", 0.025 );
666 settings.parm("StringFlav:mesonUDvector", 0.5 );
667 settings.parm("StringFlav:mesonSvector", 0.6 );
668 settings.parm("StringFlav:mesonCvector", 1.5 );
669 settings.parm("StringFlav:mesonBvector", 2.5 );
670 settings.parm("StringFlav:etaSup", 0.60 );
671 settings.parm("StringFlav:etaPrimeSup", 0.15 );
672 settings.parm("StringFlav:popcornSpair", 1.0 );
673 settings.parm("StringFlav:popcornSmeson", 1.0 );
680 // Main routine to generate the next event, using internal machinery.
682 bool Pythia::next() {
684 // Simpler option when only HadronLevel to be generated.
685 if (!doProcessLevel) {
687 // Set correct energy for system.
689 for (int i = 1; i < event.size(); ++i)
690 if (event[i].isFinal()) pSum += event[i].p();
692 event[0].m( pSum.mCalc() );
694 // Generate hadronization and decays.
695 return forceHadronLevel();
703 // Can only generate event if initialization worked.
705 info.errorMsg("Abort from Pythia::next: "
706 "not properly initialized so cannot generate events");
710 // Pick beam momentum spread and beam vertex.
711 if (doMomentumSpread || doVertexSpread) beamShapePtr->pick();
713 // Recalculate kinematics when beam momentum spread.
714 if (doMomentumSpread) nextKinematics();
716 // Outer loop over hard processes; only relevant for user-set vetoes.
718 bool hasVetoed = false;
720 // Provide the hard process that starts it off. Only one try.
723 if ( !processLevel.next( process) ) {
724 info.errorMsg("Abort from Pythia::next: "
725 "processLevel failed; giving up");
729 // Possibility for a user veto of the process-level event.
731 hasVetoed = userHooksPtr->doVetoProcessLevel( process);
732 if (hasVetoed) continue;
735 // Possibility to stop the generation at this stage.
736 if (!doPartonLevel) {
737 boostAndVertex( true, true);
738 processLevel.accumulate();
742 // Allow up to ten tries for parton- and hadron-level processing.
743 bool physical = true;
744 bool hasBoosted = false;
745 for (int iTry = 0; iTry < NTRY; ++ iTry) {
749 // Restore process record to CM frame if it was boosted.
751 boostAndVertex( false, false);
755 // Reset event record and (extracted partons from) beam remnants.
760 // Parton-level evolution: ISR, FSR, MI.
761 if ( !partonLevel.next( process, event) ) {
762 // Skip to next hard process for failure owing to deliberate veto.
763 hasVetoed = partonLevel.hasVetoed();
764 if (hasVetoed) break;
765 // Else make a new try for other failures.
766 info.errorMsg("Error in Pythia::next: "
767 "partonLevel failed; try again");
772 // Possibility for a user veto of the parton-level event.
774 hasVetoed = userHooksPtr->doVetoPartonLevel( event);
775 if (hasVetoed) break;
778 // Boost to lab frame (before decays, for vertices).
779 boostAndVertex( true, true);
780 if (frameType == 2 || frameType == 3) hasBoosted = true;
782 // Possibility to stop the generation at this stage.
783 if (!doHadronLevel) {
784 processLevel.accumulate();
785 partonLevel.accumulate();
787 // Optionally check final event for problems.
788 if (checkEvent && !check()) {
789 info.errorMsg("Abort from Pythia::next: "
790 "check of event revealed problems");
796 // Hadron-level: hadronization, decays.
797 if ( !hadronLevel.next( event) ) {
798 info.errorMsg("Error in Pythia::next: "
799 "hadronLevel failed; try again");
804 // Stop parton- and hadron-level looping if you got this far.
808 // If event vetoed then to make a new try.
809 if (hasVetoed) continue;
811 // If event failed any other way (after ten tries) then give up.
813 info.errorMsg("Abort from Pythia::next: "
814 "parton+hadronLevel failed; giving up");
818 // Process- and parton-level statistics.
819 processLevel.accumulate();
820 partonLevel.accumulate();
822 // End of outer loop over hard processes. Done with normal option.
826 // Optionally check final event for problems.
827 if (checkEvent && !check()) {
828 info.errorMsg("Abort from Pythia::next: "
829 "check of event revealed problems");
840 // Generate only the hadronization/decay stage, using internal machinery.
841 // The "event" instance should already contain a parton-level configuration.
843 bool Pythia::forceHadronLevel() {
845 // Can only generate event if initialization worked.
847 info.errorMsg("Abort from Pythia::forceHadronLevel: "
848 "not properly initialized so cannot generate events");
852 // Check whether any junctions in system. (Normally done in ProcessLevel.)
853 event.clearJunctions();
854 processLevel.findJunctions( event);
856 // Save spare copy of event in case of failure.
857 Event spareEvent = event;
859 // Allow up to ten tries for hadron-level processing.
860 bool physical = true;
861 for (int iTry = 0; iTry < NTRY; ++ iTry) {
864 // Hadron-level: hadronization, decays.
865 if (hadronLevel.next( event)) break;
867 // If failure then warn, restore original configuration and try again.
868 info.errorMsg("Error in Pythia::forceHadronLevel: "
869 "hadronLevel failed; try again");
874 // Done for simpler option.
876 info.errorMsg("Abort from Pythia::forceHadronLevel: "
877 "hadronLevel failed; giving up");
881 // Optionally check final event for problems.
882 if (checkEvent && !check()) {
883 info.errorMsg("Abort from Pythia::forceHadronLevel: "
884 "check of event revealed problems");
895 // Recalculate kinematics for each event when beam momentum has a spread.
897 void Pythia::nextKinematics() {
899 // Read out momentum shift to give current beam momenta.
900 pAnow = pAinit + beamShapePtr->deltaPA();
901 pAnow.e( sqrt(pAnow.pAbs2() + mA * mA) );
902 pBnow = pBinit + beamShapePtr->deltaPB();
903 pBnow.e( sqrt(pBnow.pAbs2() + mB * mB) );
905 // Construct CM frame kinematics.
906 eCM = (pAnow + pBnow).mCalc();
907 pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
908 * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
910 eA = sqrt(mA*mA + pzAcm*pzAcm);
911 eB = sqrt(mB*mB + pzBcm*pzBcm);
913 // Set relevant info for other classes to use.
914 info.setBeamA( idA, pzAcm, eA, mA);
915 info.setBeamB( idB, pzBcm, eB, mB);
917 beamA.newPzE( pzAcm, eA);
918 beamB.newPzE( pzBcm, eB);
920 // Set boost/rotation matrices from/to CM frame.
922 MfromCM.fromCMframe( pAnow, pBnow);
930 // Boost from CM frame to lab frame, or inverse. Set production vertex.
932 void Pythia::boostAndVertex( bool toLab, bool setVertex) {
934 // Boost process from CM frame to lab frame.
936 if (frameType == 2) process.bst(0., 0., betaZ, gammaZ);
937 else if (frameType == 3) process.rotbst(MfromCM);
939 // Boost nonempty event from CM frame to lab frame.
940 if (event.size() > 0) {
941 if (frameType == 2) event.bst(0., 0., betaZ, gammaZ);
942 else if (frameType == 3) event.rotbst(MfromCM);
945 // Boost process from lab frame to CM frame.
947 if (frameType == 2) process.bst(0., 0., -betaZ, gammaZ);
948 else if (frameType == 3) process.rotbst(MtoCM);
950 // Boost nonempty event from lab frame to CM frame.
951 if (event.size() > 0) {
952 if (frameType == 2) event.bst(0., 0., -betaZ, gammaZ);
953 else if (frameType == 3) event.rotbst(MtoCM);
957 // Set production vertex; assumes particles are in lab frame and at origin
958 if (setVertex && doVertexSpread) {
959 Vec4 vertex = beamShapePtr->vertex();
960 for (int i = 0; i < process.size(); ++i) process[i].vProd( vertex);
961 for (int i = 0; i < event.size(); ++i) event[i].vProd( vertex);
968 // Print statistics on event generation.
970 void Pythia::statistics(bool all, bool reset) {
972 // Statistics on cross section and number of events.
973 if (doProcessLevel) processLevel.statistics(reset);
975 // Statistics from other classes, e.g. multiple interactions.
976 if (all) partonLevel.statistics(reset);
978 // Summary of which and how many warnings/errors encountered.
979 info.errorStatistics();
980 if (reset) info.errorReset();
986 // Write the Pythia banner, with symbol and version information.
988 void Pythia::banner(ostream& os) {
990 // Read in version number and last date of change.
991 double versionNumber = Settings::parm("Pythia:versionNumber");
992 int versionDate = Settings::mode("Pythia:versionDate");
993 string month[12] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun",
994 "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
996 // Get date and time.
999 strftime(dateNow,12,"%d %b %Y",localtime(&t));
1001 strftime(timeNow,9,"%H:%M:%S",localtime(&t));
1005 << " *-------------------------------------------"
1006 << "-----------------------------------------* \n"
1009 << " | *----------------------------------------"
1010 << "--------------------------------------* | \n"
1015 << " | | PPP Y Y TTTTT H H III A "
1016 << " Welcome to the Lund Monte Carlo! | | \n"
1017 << " | | P P Y Y T H H I A A "
1018 << " This is PYTHIA version " << fixed << setprecision(3)
1019 << setw(5) << versionNumber << " | | \n"
1020 << " | | PPP Y T HHHHH I AAAAA"
1021 << " Last date of change: " << setw(2) << versionDate%100
1022 << " " << month[ (versionDate/100)%100 - 1 ]
1023 << " " << setw(4) << versionDate/10000 << " | | \n"
1024 << " | | P Y T H H I A A"
1026 << " | | P Y T H H III A A"
1027 << " Now is " << dateNow << " at " << timeNow << " | | \n"
1030 << " | | Main author: Torbjorn Sjostrand; CERN"
1031 << "/PH, CH-1211 Geneva, Switzerland, | | \n"
1032 << " | | and Department of Theoretical Physi"
1033 << "cs, Lund University, Lund, Sweden; | | \n"
1034 << " | | phone: + 41 - 22 - 767 82 27; e-mai"
1035 << "l: torbjorn@thep.lu.se | | \n"
1036 << " | | Author: Stephen Mrenna; Computing Div"
1037 << "ision, Simulations Group, | | \n"
1038 << " | | Fermi National Accelerator Laborato"
1039 << "ry, MS 234, Batavia, IL 60510, USA; | | \n"
1040 << " | | phone: + 1 - 630 - 840 - 2556; e-ma"
1041 << "il: mrenna@fnal.gov | | \n"
1042 << " | | Author: Peter Skands; CERN/PH, CH-121"
1043 << "1 Geneva, Switzerland, | | \n"
1044 << " | | and Theoretical Physics Department,"
1046 << " | | Fermi National Accelerator Laborato"
1047 << "ry, MS 106, Batavia, IL 60510, USA; | | \n"
1048 << " | | phone: + 41 - 22 - 767 24 59; e-mai"
1049 << "l: skands@fnal.gov | | \n"
1052 << " | | The main program reference is the 'Br"
1053 << "ief Introduction to PYTHIA 8.1', | | \n"
1054 << " | | T. Sjostrand, S. Mrenna and P. Skands"
1055 << ", arXiv:0710.3820 | | \n"
1058 << " | | The main physics reference is the 'PY"
1059 << "THIA 6.4 Physics and Manual', | | \n"
1060 << " | | T. Sjostrand, S. Mrenna and P. Skands"
1061 << ", JHEP05 (2006) 026 [hep-ph/0603175]. | | \n"
1064 << " | | An archive of program versions and do"
1065 << "cumentation is found on the web: | | \n"
1066 << " | | http://www.thep.lu.se/~torbjorn/Pythi"
1070 << " | | This program is released under the GN"
1071 << "U General Public Licence version 2. | | \n"
1072 << " | | Please respect the MCnet Guidelines f"
1073 << "or Event Generator Authors and Users. | | \n"
1076 << " | | Disclaimer: this program comes withou"
1077 << "t any guarantees. | | \n"
1078 << " | | Beware of errors and use common sense"
1079 << " when interpreting results. | | \n"
1082 << " | | Copyright (C) 2008 Torbjorn Sjostrand"
1088 << " | *----------------------------------------"
1089 << "--------------------------------------* | \n"
1092 << " *-------------------------------------------"
1093 << "-----------------------------------------* \n" << endl;
1099 // Check for lines in file that mark the beginning of new subrun.
1101 int Pythia::readSubrun(string line, bool warn, ostream& os) {
1103 // If empty line then done.
1104 int subrunLine = SUBRUNDEFAULT;
1105 if (line.find_first_not_of(" ") == string::npos) return subrunLine;
1107 // If first character is not a letter, then done.
1108 string lineNow = line;
1109 int firstChar = lineNow.find_first_not_of(" ");
1110 if (!isalpha(lineNow[firstChar])) return subrunLine;
1112 // Replace an equal sign by a blank to make parsing simpler.
1113 while (lineNow.find("=") != string::npos) {
1114 int firstEqual = lineNow.find_first_of("=");
1115 lineNow.replace(firstEqual, 1, " ");
1118 // Get first word of a line.
1119 istringstream splitLine(lineNow);
1123 // Replace two colons by one (:: -> :) to allow for such mistakes.
1124 while (name.find("::") != string::npos) {
1125 int firstColonColon = name.find_first_of("::");
1126 name.replace(firstColonColon, 2, ":");
1129 // Convert to lowercase.
1130 for (int i = 0; i < int(name.length()); ++i)
1131 name[i] = std::tolower(name[i]);
1133 // If no match then done.
1134 if (name != "main:subrun") return subrunLine;
1136 // Else find new subrun number and return it.
1137 splitLine >> subrunLine;
1139 if (warn) os << "\n PYTHIA Warning: Main:subrun number not"
1140 << " recognized; skip:\n " << line << endl;
1141 subrunLine = SUBRUNDEFAULT;
1149 // Check that the final event makes sense: no unknown id codes;
1150 // charge and energy-momentum conserved.
1152 bool Pythia::check(ostream& os) {
1155 bool physical = true;
1160 double chargeSum = 0.;
1162 // Incoming beams counted with negative momentum and charge.
1163 if (doProcessLevel) {
1164 pSum = - (event[1].p() + event[2].p());
1165 chargeSum = - (event[1].charge() + event[2].charge());
1167 // If no ProcessLevel then sum momentum and charge in initial state.
1169 pSum = - event[0].p();
1170 for (int i = 0; i < process.size(); ++i)
1171 if (process[i].isFinal()) chargeSum -= process[i].charge();
1173 double eLab = abs(pSum.e());
1175 // Loop over particles in the event.
1176 for (int i = 0; i < event.size(); ++i) {
1178 // Look for any unrecognized particle codes.
1179 int id = event[i].id();
1180 if (id == 0 || !ParticleDataTable::isParticle(id)) {
1181 ostringstream errCode;
1182 errCode << ", id = " << id;
1183 info.errorMsg("Error in Pythia::check: "
1184 "unknown particle code", errCode.str());
1186 iErrId.push_back(i);
1188 // Check that colour assignments are the expected ones.
1190 int colType = event[i].colType();
1191 int col = event[i].col();
1192 int acol = event[i].acol();
1193 if ( (colType == 0 && (col > 0 || acol > 0))
1194 || (colType == 1 && (col <= 0 || acol > 0))
1195 || (colType == -1 && (col > 0 || acol <= 0))
1196 || (colType == 2 && (col <= 0 || acol <= 0)) ) {
1197 ostringstream errCode;
1198 errCode << ", id = " << id << " cols = " << col << " " << acol;
1199 info.errorMsg("Error in Pythia::check: "
1200 "incorrect colours", errCode.str());
1202 iErrCol.push_back(i);
1206 // Look for particles with not-a-number energy/momentum/mass.
1207 if (abs(event[i].px()) >= 0. && abs(event[i].py()) >= 0.
1208 && abs(event[i].pz()) >= 0. && abs(event[i].e()) >= 0.
1209 && abs(event[i].m()) >= 0.) ;
1211 info.errorMsg("Error in Pythia::check: "
1212 "not-a-number energy/momentum/mass");
1214 iErrNan.push_back(i);
1217 // Add final-state four-momentum and charge.
1218 if (event[i].isFinal()) {
1219 pSum += event[i].p();
1220 chargeSum += event[i].charge();
1223 // End of particle loop.
1226 // Check energy-momentum/charge conservation.
1227 double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
1229 if (epDev > epTolErr * eLab) {
1230 info.errorMsg("Error in Pythia::check: energy-momentum not conserved");
1232 } else if (epDev > epTolWarn * eLab) {
1233 info.errorMsg("Warning in Pythia::check: "
1234 "energy-momentum not quite conserved");
1236 if (abs(chargeSum) > 0.1) {
1237 info.errorMsg("Error in Pythia::check: charge not conserved");
1241 // Done for sensible events.
1242 if (physical) return true;
1244 // Print (the first few) flawed events.
1245 if (nErrEvent < nErrList) {
1246 os << " PYTHIA erroneous event info: \n";
1247 if (iErrId.size() > 0) {
1248 os << " unknown particle codes in lines ";
1249 for (int i = 0; i < int(iErrId.size()); ++i)
1250 os << iErrId[i] << " ";
1253 if (iErrCol.size() > 0) {
1254 os << " incorrect colour assignments in lines ";
1255 for (int i = 0; i < int(iErrCol.size()); ++i)
1256 os << iErrCol[i] << " ";
1259 if (iErrNan.size() > 0) {
1260 os << " not-a-number energy/momentum/mass in lines ";
1261 for (int i = 0; i < int(iErrNan.size()); ++i)
1262 os << iErrNan[i] << " ";
1265 if (epDev > epTolErr * eLab) os << scientific << setprecision(3)
1266 << " total energy-momentum non-conservation = " << epDev << "\n";
1267 if (abs(chargeSum) > 0.1) os << fixed << setprecision(2)
1268 << " total charge non-conservation = " << chargeSum << "\n";
1273 // Update error counter. Done also for flawed event.
1281 // Routine to set up a PDF pointer.
1283 PDF* Pythia::getPDFPtr(int idIn, int sequence) {
1287 // Proton beam, normal choice.
1288 if (abs(idIn) == 2212 && sequence == 1) {
1289 int pSet = settings.mode("PDF:pSet");
1290 bool useLHAPDF = settings.flag("PDF:useLHAPDF");
1292 // Use internal sets.
1294 if (pSet == 1) tempPDFPtr = new GRV94L(idIn);
1295 else tempPDFPtr = new CTEQ5L(idIn);
1298 // Use sets from LHAPDF.
1300 string LHAPDFset = settings.word("PDF:LHAPDFset");
1301 int LHAPDFmember = settings.mode("PDF:LHAPDFmember");
1302 tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 1, &info);
1304 // Optionally allow extrapolation beyond x and Q2 limits.
1305 tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
1309 // Proton beam, special choice for the hard process..
1310 else if (abs(idIn) == 2212) {
1311 int pSet = settings.mode("PDF:pHardSet");
1312 bool useLHAPDF = settings.flag("PDF:useHardLHAPDF");
1314 // Use internal sets.
1316 if (pSet == 1) tempPDFPtr = new GRV94L(idIn);
1317 else tempPDFPtr = new CTEQ5L(idIn);
1320 // Use sets from LHAPDF.
1322 string LHAPDFset = settings.word("PDF:hardLHAPDFset");
1323 int LHAPDFmember = settings.mode("PDF:hardLHAPDFmember");
1324 tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 2, &info);
1326 // Optionally allow extrapolation beyond x and Q2 limits.
1327 tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
1331 // Lepton beam; resolved or not.
1333 if (settings.flag("PDF:lepton") && abs(idIn)%2 == 1)
1334 tempPDFPtr = new Lepton(idIn);
1335 else tempPDFPtr = new LeptonPoint(idIn);
1344 } // end namespace Pythia8