Line data Source code
1 : // Pythia.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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.
5 :
6 : // Function definitions (not found in the header) for the Pythia class.
7 :
8 : #include "Pythia8/Pythia.h"
9 :
10 : // Access time information.
11 : #include <ctime>
12 :
13 : // Allow string and character manipulation.
14 : #include <cctype>
15 :
16 : namespace Pythia8 {
17 :
18 : //==========================================================================
19 :
20 : // The Pythia class.
21 :
22 : //--------------------------------------------------------------------------
23 :
24 : // The current Pythia (sub)version number, to agree with XML version.
25 : const double Pythia::VERSIONNUMBERHEAD = PYTHIA_VERSION;
26 : const double Pythia::VERSIONNUMBERCODE = 8.210;
27 :
28 : //--------------------------------------------------------------------------
29 :
30 : // Constants: could be changed here if desired, but normally should not.
31 : // These are of technical nature, as described for each.
32 :
33 : // Maximum number of tries to produce parton level from given input.
34 : const int Pythia::NTRY = 10;
35 :
36 : // Negative integer to denote that no subrun has been set.
37 : const int Pythia::SUBRUNDEFAULT = -999;
38 :
39 : //--------------------------------------------------------------------------
40 :
41 : // Constructor.
42 :
43 0 : Pythia::Pythia(string xmlDir, bool printBanner) {
44 :
45 : // Initial values for pointers to PDF's.
46 0 : useNewPdfA = false;
47 0 : useNewPdfB = false;
48 0 : useNewPdfHard = false;
49 0 : useNewPdfPomA = false;
50 0 : useNewPdfPomB = false;
51 0 : pdfAPtr = 0;
52 0 : pdfBPtr = 0;
53 0 : pdfHardAPtr = 0;
54 0 : pdfHardBPtr = 0;
55 0 : pdfPomAPtr = 0;
56 0 : pdfPomBPtr = 0;
57 :
58 : // Initial values for pointers to Les Houches Event objects.
59 0 : doLHA = false;
60 0 : useNewLHA = false;
61 0 : lhaUpPtr = 0;
62 :
63 : //Initial value for couplings pointer
64 0 : couplingsPtr = &couplings;
65 :
66 : // Initial value for pointer to external decay handler.
67 0 : decayHandlePtr = 0;
68 :
69 : // Initial value for pointer to user hooks.
70 0 : userHooksPtr = 0;
71 :
72 : // Initial value for pointer to merging hooks.
73 0 : doMerging = false;
74 0 : hasMergingHooks = false;
75 0 : hasOwnMergingHooks = false;
76 0 : mergingHooksPtr = 0;
77 :
78 : // Initial value for pointer to beam shape.
79 0 : useNewBeamShape = false;
80 0 : beamShapePtr = 0;
81 :
82 : // Initial values for pointers to timelike and spacelike showers.
83 0 : useNewTimesDec = false;
84 0 : useNewTimes = false;
85 0 : useNewSpace = false;
86 0 : timesDecPtr = 0;
87 0 : timesPtr = 0;
88 0 : spacePtr = 0;
89 :
90 : // Find path to data files, i.e. xmldoc directory location.
91 : // Environment variable takes precedence, then constructor input,
92 : // and finally the pre-processor constant XMLDIR.
93 0 : xmlPath = "";
94 : const char* PYTHIA8DATA = "PYTHIA8DATA";
95 0 : char* envPath = getenv(PYTHIA8DATA);
96 0 : if (envPath != 0 && *envPath != '\0') {
97 : int i = 0;
98 0 : while (*(envPath+i) != '\0') xmlPath += *(envPath+(i++));
99 0 : } else {
100 0 : if (xmlDir[ xmlDir.length() - 1 ] != '/') xmlDir += "/";
101 0 : xmlPath = xmlDir;
102 0 : ifstream xmlFile((xmlPath + "Index.xml").c_str());
103 : // if (!xmlFile.good()) xmlPath = XMLDIR;
104 0 : xmlFile.close();
105 0 : }
106 0 : if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
107 :
108 : // Read in files with all flags, modes, parms and words.
109 0 : settings.initPtr( &info);
110 0 : string initFile = xmlPath + "Index.xml";
111 0 : isConstructed = settings.init( initFile);
112 0 : if (!isConstructed) {
113 0 : info.errorMsg("Abort from Pythia::Pythia: settings unavailable");
114 0 : return;
115 : }
116 :
117 : // Check that XML version number matches code version number.
118 0 : double versionNumberXML = parm("Pythia:versionNumber");
119 0 : isConstructed = (abs(versionNumberXML - VERSIONNUMBERCODE) < 0.0005);
120 0 : if (!isConstructed) {
121 0 : ostringstream errCode;
122 0 : errCode << fixed << setprecision(3) << ": in code " << VERSIONNUMBERCODE
123 0 : << " but in XML " << versionNumberXML;
124 0 : info.errorMsg("Abort from Pythia::Pythia: unmatched version numbers",
125 0 : errCode.str());
126 : return;
127 0 : }
128 :
129 : // Check that header version number matches code version number.
130 0 : isConstructed = (abs(VERSIONNUMBERHEAD - VERSIONNUMBERCODE) < 0.0005);
131 0 : if (!isConstructed) {
132 0 : ostringstream errCode;
133 0 : errCode << fixed << setprecision(3) << ": in code " << VERSIONNUMBERCODE
134 0 : << " but in header " << VERSIONNUMBERHEAD;
135 0 : info.errorMsg("Abort from Pythia::Pythia: unmatched version numbers",
136 0 : errCode.str());
137 : return;
138 0 : }
139 :
140 : // Read in files with all particle data.
141 0 : particleData.initPtr( &info, &settings, &rndm, couplingsPtr);
142 0 : string dataFile = xmlPath + "ParticleData.xml";
143 0 : isConstructed = particleData.init( dataFile);
144 0 : if (!isConstructed) {
145 0 : info.errorMsg("Abort from Pythia::Pythia: particle data unavailable");
146 0 : return;
147 : }
148 :
149 : // Write the Pythia banner to output.
150 0 : if (printBanner) banner();
151 :
152 : // Not initialized until at the end of the init() call.
153 0 : isInit = false;
154 0 : info.addCounter(0);
155 :
156 0 : }
157 :
158 : //--------------------------------------------------------------------------
159 :
160 : // Destructor.
161 :
162 0 : Pythia::~Pythia() {
163 :
164 : // Delete the PDF's created with new.
165 0 : if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr;
166 0 : if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr;
167 0 : if (useNewPdfA) delete pdfAPtr;
168 0 : if (useNewPdfB) delete pdfBPtr;
169 0 : if (useNewPdfPomA) delete pdfPomAPtr;
170 0 : if (useNewPdfPomB) delete pdfPomBPtr;
171 :
172 : // Delete the Les Houches object created with new.
173 0 : if (useNewLHA) delete lhaUpPtr;
174 :
175 : // Delete the MergingHooks object created with new.
176 0 : if (hasOwnMergingHooks) delete mergingHooksPtr;
177 :
178 : // Delete the BeamShape object created with new.
179 0 : if (useNewBeamShape) delete beamShapePtr;
180 :
181 : // Delete the timelike and spacelike showers created with new.
182 0 : if (useNewTimesDec) delete timesDecPtr;
183 0 : if (useNewTimes && !useNewTimesDec) delete timesPtr;
184 0 : if (useNewSpace) delete spacePtr;
185 :
186 0 : }
187 :
188 : //--------------------------------------------------------------------------
189 :
190 : // Read in one update for a setting or particle data from a single line.
191 :
192 : bool Pythia::readString(string line, bool warn) {
193 :
194 : // Check that constructor worked.
195 0 : if (!isConstructed) return false;
196 :
197 : // If empty line then done.
198 0 : if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return true;
199 :
200 : // If first character is not a letter/digit, then taken to be a comment.
201 0 : int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
202 0 : if (!isalnum(line[firstChar])) return true;
203 :
204 : // Send on particle data to the ParticleData database.
205 0 : if (isdigit(line[firstChar])) {
206 0 : bool passed = particleData.readString(line, warn);
207 0 : if (passed) particleDataBuffer << line << endl;
208 : return passed;
209 : }
210 :
211 : // Everything else sent on to Settings.
212 0 : return settings.readString(line, warn);
213 :
214 0 : }
215 :
216 : //--------------------------------------------------------------------------
217 :
218 : // Read in updates for settings or particle data from user-defined file.
219 :
220 : bool Pythia::readFile(string fileName, bool warn, int subrun) {
221 :
222 : // Check that constructor worked.
223 0 : if (!isConstructed) return false;
224 :
225 : // Open file for reading.
226 0 : const char* cstring = fileName.c_str();
227 0 : ifstream is(cstring);
228 0 : if (!is.good()) {
229 0 : info.errorMsg("Error in Pythia::readFile: did not find file", fileName);
230 0 : return false;
231 : }
232 :
233 : // Hand over real work to next method.
234 0 : return readFile( is, warn, subrun);
235 :
236 0 : }
237 :
238 : //--------------------------------------------------------------------------
239 :
240 : // Read in updates for settings or particle data
241 : // from user-defined stream (or file).
242 :
243 : bool Pythia::readFile(istream& is, bool warn, int subrun) {
244 :
245 : // Check that constructor worked.
246 0 : if (!isConstructed) return false;
247 :
248 : // Read in one line at a time.
249 0 : string line;
250 : bool isCommented = false;
251 : bool accepted = true;
252 : int subrunNow = SUBRUNDEFAULT;
253 0 : while ( getline(is, line) ) {
254 :
255 : // Check whether entering, leaving or inside commented-commands section.
256 0 : int commentLine = readCommented( line);
257 0 : if (commentLine == +1) isCommented = true;
258 0 : else if (commentLine == -1) isCommented = false;
259 0 : else if (isCommented) ;
260 :
261 : else {
262 : // Check whether entered new subrun.
263 0 : int subrunLine = readSubrun( line, warn);
264 0 : if (subrunLine >= 0) subrunNow = subrunLine;
265 :
266 : // Process the line if in correct subrun.
267 0 : if ( (subrunNow == subrun || subrunNow == SUBRUNDEFAULT)
268 0 : && !readString( line, warn) ) accepted = false;
269 : }
270 :
271 : // Reached end of input file.
272 : };
273 0 : return accepted;
274 :
275 0 : }
276 :
277 : //--------------------------------------------------------------------------
278 :
279 : // Routine to pass in pointers to PDF's. Usage optional.
280 :
281 : bool Pythia::setPDFPtr( PDF* pdfAPtrIn, PDF* pdfBPtrIn, PDF* pdfHardAPtrIn,
282 : PDF* pdfHardBPtrIn, PDF* pdfPomAPtrIn, PDF* pdfPomBPtrIn) {
283 :
284 : // Delete any PDF's created in a previous init call.
285 0 : if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr;
286 0 : if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr;
287 0 : if (useNewPdfA) delete pdfAPtr;
288 0 : if (useNewPdfB) delete pdfBPtr;
289 0 : if (useNewPdfPomA) delete pdfPomAPtr;
290 0 : if (useNewPdfPomB) delete pdfPomBPtr;
291 :
292 : // Reset pointers to be empty.
293 0 : useNewPdfA = false;
294 0 : useNewPdfB = false;
295 0 : useNewPdfHard = false;
296 0 : useNewPdfPomA = false;
297 0 : useNewPdfPomB = false;
298 0 : pdfAPtr = 0;
299 0 : pdfBPtr = 0;
300 0 : pdfHardAPtr = 0;
301 0 : pdfHardBPtr = 0;
302 0 : pdfPomAPtr = 0;
303 0 : pdfPomBPtr = 0;
304 :
305 : // Switch off external PDF's by zero as input.
306 0 : if (pdfAPtrIn == 0 && pdfBPtrIn == 0) return true;
307 :
308 : // The two PDF objects cannot be one and the same.
309 0 : if (pdfAPtrIn == pdfBPtrIn) return false;
310 :
311 : // Save pointers.
312 0 : pdfAPtr = pdfAPtrIn;
313 0 : pdfBPtr = pdfBPtrIn;
314 :
315 : // By default same pointers for hard-process PDF's.
316 0 : pdfHardAPtr = pdfAPtrIn;
317 0 : pdfHardBPtr = pdfBPtrIn;
318 :
319 : // Optionally allow separate pointers for hard process.
320 0 : if (pdfHardAPtrIn != 0 && pdfHardBPtrIn != 0) {
321 0 : if (pdfHardAPtrIn == pdfHardBPtrIn) return false;
322 0 : pdfHardAPtr = pdfHardAPtrIn;
323 0 : pdfHardBPtr = pdfHardBPtrIn;
324 0 : }
325 :
326 : // Optionally allow pointers for Pomerons in the proton.
327 0 : if (pdfPomAPtrIn != 0 && pdfPomBPtrIn != 0) {
328 0 : if (pdfPomAPtrIn == pdfPomBPtrIn) return false;
329 0 : pdfPomAPtr = pdfPomAPtrIn;
330 0 : pdfPomBPtr = pdfPomBPtrIn;
331 0 : }
332 :
333 : // Done.
334 0 : return true;
335 0 : }
336 :
337 : //--------------------------------------------------------------------------
338 :
339 : // Routine to initialize with the variable values of the Beams kind.
340 :
341 : bool Pythia::init() {
342 :
343 : // Check that constructor worked.
344 0 : isInit = false;
345 0 : if (!isConstructed) {
346 0 : info.errorMsg("Abort from Pythia::init: constructor "
347 : "initialization failed");
348 0 : return false;
349 : }
350 :
351 : // Early readout, if return false or changed when no beams.
352 0 : doProcessLevel = settings.flag("ProcessLevel:all");
353 :
354 : // Check that changes in Settings and ParticleData have worked.
355 0 : if (settings.readingFailed()) {
356 0 : info.errorMsg("Abort from Pythia::init: some user settings "
357 : "did not make sense");
358 0 : return false;
359 : }
360 0 : if (particleData.readingFailed()) {
361 0 : info.errorMsg("Abort from Pythia::init: some user particle data "
362 : "did not make sense");
363 0 : return false;
364 : }
365 :
366 : // Begin initialization. Find which frame type to use.
367 0 : info.addCounter(1);
368 0 : frameType = mode("Beams:frameType");
369 :
370 : // Initialization with internal processes: read in and set values.
371 0 : if (frameType < 4 ) {
372 0 : doLHA = false;
373 0 : boostType = frameType;
374 0 : idA = mode("Beams:idA");
375 0 : idB = mode("Beams:idB");
376 0 : eCM = parm("Beams:eCM");
377 0 : eA = parm("Beams:eA");
378 0 : eB = parm("Beams:eB");
379 0 : pxA = parm("Beams:pxA");
380 0 : pyA = parm("Beams:pyA");
381 0 : pzA = parm("Beams:pzA");
382 0 : pxB = parm("Beams:pxB");
383 0 : pyB = parm("Beams:pyB");
384 0 : pzB = parm("Beams:pzB");
385 :
386 : // Initialization with a Les Houches Event File or an LHAup object.
387 0 : } else {
388 0 : doLHA = true;
389 0 : boostType = 2;
390 0 : string lhef = word("Beams:LHEF");
391 0 : string lhefHeader = word("Beams:LHEFheader");
392 0 : bool readHeaders = flag("Beams:readLHEFheaders");
393 0 : bool setScales = flag("Beams:setProductionScalesFromLHEF");
394 0 : bool skipInit = flag("Beams:newLHEFsameInit");
395 0 : int nSkipAtInit = mode("Beams:nSkipLHEFatInit");
396 :
397 : // For file input: renew file stream or (re)new Les Houches object.
398 0 : if (frameType == 4) {
399 0 : const char* cstring1 = lhef.c_str();
400 0 : if (useNewLHA && skipInit) lhaUpPtr->newEventFile(cstring1);
401 : else {
402 0 : if (useNewLHA) delete lhaUpPtr;
403 : // Header is optional, so use NULL pointer to indicate no value.
404 0 : const char* cstring2 = (lhefHeader == "void")
405 0 : ? NULL : lhefHeader.c_str();
406 0 : lhaUpPtr = new LHAupLHEF(&info, cstring1, cstring2,
407 0 : readHeaders, setScales);
408 0 : useNewLHA = true;
409 : }
410 :
411 : // Check that file was properly opened.
412 0 : if (!lhaUpPtr->fileFound()) {
413 0 : info.errorMsg("Abort from Pythia::init: "
414 : "Les Houches Event File not found");
415 0 : return false;
416 : }
417 :
418 : // For object input: at least check that not null pointer.
419 0 : } else {
420 0 : if (lhaUpPtr == 0) {
421 0 : info.errorMsg("Abort from Pythia::init: "
422 : "LHAup object not found");
423 0 : return false;
424 : }
425 :
426 : // LHAup object generic abort using fileFound() routine.
427 0 : if (!lhaUpPtr->fileFound()) {
428 0 : info.errorMsg("Abort from Pythia::init: "
429 : "LHAup initialisation error");
430 0 : return false;
431 : }
432 : }
433 :
434 : // Send in pointer to info. Store or replace LHA pointer in other classes.
435 0 : lhaUpPtr->setPtr( &info);
436 0 : processLevel.setLHAPtr( lhaUpPtr);
437 :
438 : // If second time around, only with new file, then simplify.
439 : // Optionally skip ahead a number of events at beginning of file.
440 0 : if (skipInit) {
441 0 : isInit = true;
442 0 : info.addCounter(2);
443 0 : if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit);
444 0 : return true;
445 : }
446 :
447 : // Set up values related to merging hooks.
448 0 : if (frameType == 4 || frameType == 5) {
449 :
450 : // Set up values related to CKKW-L merging.
451 0 : bool doUserMerging = settings.flag("Merging:doUserMerging");
452 0 : bool doMGMerging = settings.flag("Merging:doMGMerging");
453 0 : bool doKTMerging = settings.flag("Merging:doKTMerging");
454 0 : bool doPTLundMerging = settings.flag("Merging:doPTLundMerging");
455 0 : bool doCutBasedMerging = settings.flag("Merging:doCutBasedMerging");
456 : // Set up values related to unitarised CKKW merging
457 0 : bool doUMEPSTree = settings.flag("Merging:doUMEPSTree");
458 0 : bool doUMEPSSubt = settings.flag("Merging:doUMEPSSubt");
459 : // Set up values related to NL3 NLO merging
460 0 : bool doNL3Tree = settings.flag("Merging:doNL3Tree");
461 0 : bool doNL3Loop = settings.flag("Merging:doNL3Loop");
462 0 : bool doNL3Subt = settings.flag("Merging:doNL3Subt");
463 : // Set up values related to unitarised NLO merging
464 0 : bool doUNLOPSTree = settings.flag("Merging:doUNLOPSTree");
465 0 : bool doUNLOPSLoop = settings.flag("Merging:doUNLOPSLoop");
466 0 : bool doUNLOPSSubt = settings.flag("Merging:doUNLOPSSubt");
467 0 : bool doUNLOPSSubtNLO = settings.flag("Merging:doUNLOPSSubtNLO");
468 0 : bool doXSectionEst = settings.flag("Merging:doXSectionEstimate");
469 0 : doMerging = doUserMerging || doMGMerging || doKTMerging
470 0 : || doPTLundMerging || doCutBasedMerging || doUMEPSTree || doUMEPSSubt
471 0 : || doNL3Tree || doNL3Loop || doNL3Subt || doUNLOPSTree
472 0 : || doUNLOPSLoop || doUNLOPSSubt || doUNLOPSSubtNLO || doXSectionEst;
473 :
474 : // Set up MergingHooks object.
475 0 : bool inputMergingHooks = (mergingHooksPtr != 0);
476 0 : if (doMerging && !inputMergingHooks){
477 0 : if (hasOwnMergingHooks && mergingHooksPtr) delete mergingHooksPtr;
478 0 : mergingHooksPtr = new MergingHooks();
479 0 : hasOwnMergingHooks = true;
480 0 : }
481 :
482 0 : hasMergingHooks = (mergingHooksPtr != 0);
483 : // Merging hooks required for merging. If no merging hooks pointer is
484 : // available, exit.
485 0 : if (doMerging && !hasMergingHooks) {
486 0 : info.errorMsg("Abort from Pythia::init: "
487 : "no merging hooks object has been provided");
488 0 : return false;
489 0 : } else if (doMerging) {
490 0 : string lhefIn = (frameType == 4) ? lhef : "";
491 0 : mergingHooksPtr->setLHEInputFile( lhefIn);
492 0 : }
493 : // Initialise counting of Les Houches Events significantly above the
494 : // merging scale.
495 0 : info.setCounter(41,0);
496 0 : }
497 :
498 : // Set LHAinit information (in some external program).
499 0 : if ( !lhaUpPtr->setInit()) {
500 0 : info.errorMsg("Abort from Pythia::init: "
501 : "Les Houches initialization failed");
502 0 : return false;
503 : }
504 :
505 : // Extract beams from values set in an LHAinit object.
506 0 : idA = lhaUpPtr->idBeamA();
507 0 : idB = lhaUpPtr->idBeamB();
508 0 : int idRenameBeams = settings.mode("LesHouches:idRenameBeams");
509 0 : if (abs(idA) == idRenameBeams) idA = 16;
510 0 : if (abs(idB) == idRenameBeams) idB = -16;
511 0 : if (idA == 0 || idB == 0) doProcessLevel = false;
512 0 : eA = lhaUpPtr->eBeamA();
513 0 : eB = lhaUpPtr->eBeamB();
514 :
515 : // Optionally skip ahead a number of events at beginning of file.
516 0 : if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit);
517 0 : }
518 :
519 : // Set up values related to user hooks.
520 0 : hasUserHooks = (userHooksPtr != 0);
521 0 : doVetoProcess = false;
522 0 : doVetoPartons = false;
523 0 : retryPartonLevel = false;
524 0 : if (hasUserHooks) {
525 0 : userHooksPtr->initPtr( &info, &settings, &particleData, &rndm, &beamA,
526 0 : &beamB, &beamPomA, &beamPomB, couplingsPtr, &partonSystems, &sigmaTot);
527 0 : if (!userHooksPtr->initAfterBeams()) {
528 0 : info.errorMsg("Abort from Pythia::init: could not initialise UserHooks");
529 0 : return false;
530 : }
531 0 : doVetoProcess = userHooksPtr->canVetoProcessLevel();
532 0 : doVetoPartons = userHooksPtr->canVetoPartonLevel();
533 0 : retryPartonLevel = userHooksPtr->retryPartonLevel();
534 0 : }
535 :
536 : // Back to common initialization. Reset error counters.
537 0 : nErrEvent = 0;
538 0 : info.setTooLowPTmin(false);
539 0 : info.sigmaReset();
540 :
541 : // Initialize data members extracted from database.
542 0 : doPartonLevel = settings.flag("PartonLevel:all");
543 0 : doHadronLevel = settings.flag("HadronLevel:all");
544 0 : doDiffraction = settings.flag("SoftQCD:all")
545 0 : || settings.flag("SoftQCD:singleDiffractive")
546 0 : || settings.flag("SoftQCD:doubleDiffractive")
547 0 : || settings.flag("SoftQCD:centralDiffractive")
548 0 : || settings.flag("SoftQCD:inelastic");
549 0 : doHardDiff = settings.flag("Diffraction:doHard");
550 0 : doResDec = settings.flag("ProcessLevel:resonanceDecays");
551 0 : doFSRinRes = doPartonLevel && settings.flag("PartonLevel:FSR")
552 0 : && settings.flag("PartonLevel:FSRinResonances");
553 0 : decayRHadrons = settings.flag("RHadrons:allowDecay");
554 0 : doMomentumSpread = settings.flag("Beams:allowMomentumSpread");
555 0 : doVertexSpread = settings.flag("Beams:allowVertexSpread");
556 0 : abortIfVeto = settings.flag("Check:abortIfVeto");
557 0 : checkEvent = settings.flag("Check:event");
558 0 : checkHistory = settings.flag("Check:history");
559 0 : nErrList = settings.mode("Check:nErrList");
560 0 : epTolErr = settings.parm("Check:epTolErr");
561 0 : epTolWarn = settings.parm("Check:epTolWarn");
562 0 : mTolErr = settings.parm("Check:mTolErr");
563 0 : mTolWarn = settings.parm("Check:mTolWarn");
564 :
565 : // Initialise merging hooks.
566 0 : if ( doMerging && (hasMergingHooks || hasOwnMergingHooks) )
567 0 : mergingHooksPtr->init( settings, &info, &particleData, &partonSystems );
568 :
569 : // Initialize the random number generator.
570 0 : if ( settings.flag("Random:setSeed") )
571 0 : rndm.init( settings.mode("Random:seed") );
572 :
573 : // Check that combinations of settings are allowed; change if not.
574 0 : checkSettings();
575 :
576 : // Initialize the SM couplings (needed to initialize resonances).
577 0 : couplingsPtr->init( settings, &rndm );
578 :
579 : // Initialize SLHA interface (including SLHA/BSM couplings).
580 0 : bool useSLHAcouplings = false;
581 0 : slhaInterface.setPtr( &info );
582 :
583 0 : slhaInterface.init( settings, &rndm, couplingsPtr, &particleData,
584 0 : useSLHAcouplings, particleDataBuffer );
585 0 : if (useSLHAcouplings) couplingsPtr = slhaInterface.couplingsPtr;
586 :
587 : // Reset couplingsPtr to the correct memory address.
588 0 : particleData.initPtr( &info, &settings, &rndm, couplingsPtr);
589 0 : if (hasUserHooks) userHooksPtr->initPtr( &info, &settings, &particleData,
590 0 : &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr,
591 0 : &partonSystems, &sigmaTot);
592 :
593 : // Set headers to distinguish the two event listing kinds.
594 0 : int startColTag = settings.mode("Event:startColTag");
595 0 : process.init("(hard process)", &particleData, startColTag);
596 0 : event.init("(complete event)", &particleData, startColTag);
597 :
598 : // Final setup stage of particle data, notably resonance widths.
599 0 : particleData.initWidths( resonancePtrs);
600 :
601 : // Set up R-hadrons particle data, where relevant.
602 0 : rHadrons.init( &info, settings, &particleData, &rndm);
603 :
604 : // Set up objects for timelike and spacelike showers.
605 0 : if (timesDecPtr == 0 || timesPtr == 0) {
606 0 : TimeShower* timesNow = new TimeShower();
607 0 : if (timesDecPtr == 0) {
608 0 : timesDecPtr = timesNow;
609 0 : useNewTimesDec = true;
610 0 : }
611 0 : if (timesPtr == 0) {
612 0 : timesPtr = timesNow;
613 0 : useNewTimes = true;
614 0 : }
615 0 : }
616 0 : if (spacePtr == 0) {
617 0 : spacePtr = new SpaceShower();
618 0 : useNewSpace = true;
619 0 : }
620 :
621 : // Initialize showers, especially for simple showers in decays.
622 0 : timesPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr,
623 0 : &partonSystems, userHooksPtr, mergingHooksPtr);
624 0 : timesDecPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr,
625 0 : &partonSystems, userHooksPtr, mergingHooksPtr);
626 0 : spacePtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr,
627 0 : &partonSystems, userHooksPtr, mergingHooksPtr);
628 0 : timesDecPtr->init( 0, 0);
629 :
630 : // Set up values related to beam shape.
631 0 : if (beamShapePtr == 0) {
632 0 : beamShapePtr = new BeamShape();
633 0 : useNewBeamShape = true;
634 0 : }
635 0 : beamShapePtr->init( settings, &rndm);
636 :
637 : // Check that beams and beam combination can be handled.
638 0 : if (!checkBeams()) {
639 0 : info.errorMsg("Abort from Pythia::init: "
640 : "checkBeams initialization failed");
641 0 : return false;
642 : }
643 :
644 : // Do not set up beam kinematics when no process level.
645 0 : if (!doProcessLevel) boostType = 1;
646 : else {
647 :
648 : // Set up beam kinematics.
649 0 : if (!initKinematics()) {
650 0 : info.errorMsg("Abort from Pythia::init: "
651 : "kinematics initialization failed");
652 0 : return false;
653 : }
654 :
655 : // Set up pointers to PDFs.
656 0 : if (!initPDFs()) {
657 0 : info.errorMsg("Abort from Pythia::init: PDF initialization failed");
658 0 : return false;
659 : }
660 :
661 : // Set up the two beams and the common remnant system.
662 0 : StringFlav* flavSelPtr = hadronLevel.getStringFlavPtr();
663 0 : beamA.init( idA, pzAcm, eA, mA, &info, settings, &particleData, &rndm,
664 0 : pdfAPtr, pdfHardAPtr, isUnresolvedA, flavSelPtr);
665 0 : beamB.init( idB, pzBcm, eB, mB, &info, settings, &particleData, &rndm,
666 0 : pdfBPtr, pdfHardBPtr, isUnresolvedB, flavSelPtr);
667 :
668 : // Optionally set up new alternative beams for these Pomerons.
669 0 : if ( doDiffraction || doHardDiff) {
670 0 : beamPomA.init( 990, 0.5 * eCM, 0.5 * eCM, 0., &info, settings,
671 0 : &particleData, &rndm, pdfPomAPtr, pdfPomAPtr, false, flavSelPtr);
672 0 : beamPomB.init( 990, -0.5 * eCM, 0.5 * eCM, 0., &info, settings,
673 0 : &particleData, &rndm, pdfPomBPtr, pdfPomBPtr, false, flavSelPtr);
674 0 : }
675 : }
676 :
677 : // Send info/pointers to process level for initialization.
678 0 : if ( doProcessLevel && !processLevel.init( &info, settings, &particleData,
679 0 : &rndm, &beamA, &beamB, couplingsPtr, &sigmaTot, doLHA, &slhaInterface,
680 0 : userHooksPtr, sigmaPtrs, phaseSpacePtrs) ) {
681 0 : info.errorMsg("Abort from Pythia::init: "
682 : "processLevel initialization failed");
683 0 : return false;
684 : }
685 :
686 : // Alternatively only initialize resonance decays.
687 0 : if ( !doProcessLevel) processLevel.initDecays( &info, &particleData,
688 0 : &rndm, lhaUpPtr);
689 :
690 : // Send info/pointers to parton level for initialization.
691 0 : if ( doPartonLevel && doProcessLevel && !partonLevel.init( &info, settings,
692 0 : &particleData, &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr,
693 0 : &partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons,
694 0 : userHooksPtr, mergingHooksPtr, false) ) {
695 0 : info.errorMsg("Abort from Pythia::init: "
696 : "partonLevel initialization failed" );
697 0 : return false;
698 : }
699 :
700 : // Alternatively only initialize final-state showers in resonance decays.
701 0 : if ( !doProcessLevel || !doPartonLevel) partonLevel.init( &info, settings,
702 0 : &particleData, &rndm, 0, 0, 0, 0, couplingsPtr, &partonSystems, 0,
703 0 : timesDecPtr, 0, 0, &rHadrons, 0, 0, false);
704 :
705 : // Send info/pointers to parton level for trial shower initialization.
706 0 : if ( doMerging && !trialPartonLevel.init( &info, settings, &particleData,
707 0 : &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr,
708 0 : &partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons,
709 0 : NULL, mergingHooksPtr, true) ) {
710 0 : info.errorMsg("Abort from Pythia::init: "
711 : "trialPartonLevel initialization failed");
712 0 : return false;
713 : }
714 :
715 : // Initialise the merging wrapper class.
716 0 : if (doMerging ) merging.init( &settings, &info, &particleData, &rndm,
717 0 : &beamA, &beamB, mergingHooksPtr, &trialPartonLevel );
718 :
719 : // Send info/pointers to hadron level for initialization.
720 : // Note: forceHadronLevel() can come, so we must always initialize.
721 0 : if ( !hadronLevel.init( &info, settings, &particleData, &rndm,
722 0 : couplingsPtr, timesDecPtr, &rHadrons, decayHandlePtr,
723 0 : handledParticles) ) {
724 0 : info.errorMsg("Abort from Pythia::init: "
725 : "hadronLevel initialization failed");
726 0 : return false;
727 : }
728 :
729 : // Optionally check particle data table for inconsistencies.
730 0 : if ( settings.flag("Check:particleData") )
731 0 : particleData.checkTable( settings.mode("Check:levelParticleData") );
732 :
733 : // Optionally show settings and particle data, changed or all.
734 0 : bool showCS = settings.flag("Init:showChangedSettings");
735 0 : bool showAS = settings.flag("Init:showAllSettings");
736 0 : bool showCPD = settings.flag("Init:showChangedParticleData");
737 0 : bool showCRD = settings.flag("Init:showChangedResonanceData");
738 0 : bool showAPD = settings.flag("Init:showAllParticleData");
739 0 : int show1PD = settings.mode("Init:showOneParticleData");
740 0 : bool showPro = settings.flag("Init:showProcesses");
741 0 : if (showCS) settings.listChanged();
742 0 : if (showAS) settings.listAll();
743 0 : if (show1PD > 0) particleData.list(show1PD);
744 0 : if (showCPD) particleData.listChanged(showCRD);
745 0 : if (showAPD) particleData.listAll();
746 :
747 : // Listing options for the next() routine.
748 0 : nCount = settings.mode("Next:numberCount");
749 0 : nShowLHA = settings.mode("Next:numberShowLHA");
750 0 : nShowInfo = settings.mode("Next:numberShowInfo");
751 0 : nShowProc = settings.mode("Next:numberShowProcess");
752 0 : nShowEvt = settings.mode("Next:numberShowEvent");
753 0 : showSaV = settings.flag("Next:showScaleAndVertex");
754 0 : showMaD = settings.flag("Next:showMothersAndDaughters");
755 :
756 : // Init colour reconnection and junction splitting.
757 0 : colourReconnection.init( &info, settings, &rndm, &particleData,
758 0 : &beamA, &beamB, &partonSystems);
759 0 : junctionSplitting.init(&info, settings, &rndm, &particleData);
760 :
761 : // Flags for colour reconnection.
762 0 : doReconnect = settings.flag("ColourReconnection:reconnect");
763 0 : reconnectMode = settings.mode("ColourReconnection:mode");
764 0 : forceHadronLevelCR = settings.flag("ColourReconnection:forceHadronLevelCR");
765 :
766 : // Succeeded.
767 0 : isInit = true;
768 0 : info.addCounter(2);
769 0 : if (useNewLHA && showPro) lhaUpPtr->listInit();
770 : return true;
771 :
772 0 : }
773 :
774 : //--------------------------------------------------------------------------
775 :
776 : // Check that combinations of settings are allowed; change if not.
777 :
778 : void Pythia::checkSettings() {
779 :
780 : // Double rescattering not allowed if ISR or FSR.
781 0 : if ((settings.flag("PartonLevel:ISR") || settings.flag("PartonLevel:FSR"))
782 0 : && settings.flag("MultipartonInteractions:allowDoubleRescatter")) {
783 0 : info.errorMsg("Warning in Pythia::checkSettings: "
784 : "double rescattering switched off since showering is on");
785 0 : settings.flag("MultipartonInteractions:allowDoubleRescatter", false);
786 0 : }
787 :
788 0 : }
789 :
790 : //--------------------------------------------------------------------------
791 :
792 : // Check that beams and beam combination can be handled. Set up unresolved.
793 :
794 : bool Pythia::checkBeams() {
795 :
796 : // Absolute flavours. If not to do process level then no check needed.
797 0 : int idAabs = abs(idA);
798 0 : int idBabs = abs(idB);
799 0 : if (!doProcessLevel) return true;
800 :
801 : // Neutrino beams always unresolved, charged lepton ones conditionally.
802 0 : bool isLeptonA = (idAabs > 10 && idAabs < 17);
803 0 : bool isLeptonB = (idBabs > 10 && idBabs < 17);
804 0 : bool isUnresLep = !settings.flag("PDF:lepton");
805 0 : isUnresolvedA = isLeptonA && (idAabs%2 == 0 || isUnresLep);
806 0 : isUnresolvedB = isLeptonB && (idBabs%2 == 0 || isUnresLep);
807 :
808 : // Lepton-lepton collisions OK (including neutrinos) if both (un)resolved.
809 0 : if (isLeptonA && isLeptonB && isUnresolvedA == isUnresolvedB) return true;
810 :
811 : // MBR model only implemented for pp/ppbar/pbarp collisions.
812 0 : int PomFlux = settings.mode("Diffraction:PomFlux");
813 0 : if (PomFlux == 5) {
814 0 : bool ispp = (idAabs == 2212 && idBabs == 2212);
815 0 : bool ispbarpbar = (idA == -2212 && idB == -2212);
816 0 : if (ispp && !ispbarpbar) return true;
817 0 : info.errorMsg("Error in Pythia::init: cannot handle this beam combination"
818 : " with PomFlux == 5");
819 0 : return false;
820 : }
821 :
822 : // Hadron-hadron collisions OK, with Pomeron counted as hadron.
823 0 : bool isHadronA = (idAabs == 2212) || (idAabs == 2112) || (idA == 111)
824 0 : || (idAabs == 211) || (idA == 990);
825 0 : bool isHadronB = (idBabs == 2212) || (idBabs == 2112) || (idB == 111)
826 0 : || (idBabs == 211) || (idB == 990);
827 0 : if (isHadronA && isHadronB) return true;
828 :
829 : // Lepton-hadron collisions OK for DIS processes, although still primitive.
830 0 : if ( (isLeptonA && isHadronB) || (isHadronA && isLeptonB) ) {
831 0 : bool doDIS = settings.flag("WeakBosonExchange:all")
832 0 : || settings.flag("WeakBosonExchange:ff2ff(t:gmZ)")
833 0 : || settings.flag("WeakBosonExchange:ff2ff(t:W)");
834 0 : if (doDIS) return true;
835 0 : }
836 :
837 : // If no case above then failed.
838 0 : info.errorMsg("Error in Pythia::init: cannot handle this beam combination");
839 0 : return false;
840 :
841 0 : }
842 :
843 : //--------------------------------------------------------------------------
844 :
845 : // Calculate kinematics at initialization. Store beam four-momenta.
846 :
847 : bool Pythia::initKinematics() {
848 :
849 : // Find masses. Initial guess that we are in CM frame.
850 0 : mA = particleData.m0(idA);
851 0 : mB = particleData.m0(idB);
852 0 : betaZ = 0.;
853 0 : gammaZ = 1.;
854 :
855 : // Collinear beams not in CM frame: find CM energy.
856 0 : if (boostType == 2) {
857 0 : eA = max(eA, mA);
858 0 : eB = max(eB, mB);
859 0 : pzA = sqrt(eA*eA - mA*mA);
860 0 : pzB = -sqrt(eB*eB - mB*mB);
861 0 : pAinit = Vec4( 0., 0., pzA, eA);
862 0 : pBinit = Vec4( 0., 0., pzB, eB);
863 0 : eCM = sqrt( pow2(eA + eB) - pow2(pzA + pzB) );
864 :
865 : // Find boost to rest frame.
866 0 : betaZ = (pzA + pzB) / (eA + eB);
867 0 : gammaZ = (eA + eB) / eCM;
868 0 : if (abs(betaZ) < 1e-10) boostType = 1;
869 : }
870 :
871 : // Completely general beam directions: find CM energy.
872 0 : else if (boostType == 3) {
873 0 : eA = sqrt( pxA*pxA + pyA*pyA + pzA*pzA + mA*mA);
874 0 : eB = sqrt( pxB*pxB + pyB*pyB + pzB*pzB + mB*mB);
875 0 : pAinit = Vec4( pxA, pyA, pzA, eA);
876 0 : pBinit = Vec4( pxB, pyB, pzB, eB);
877 0 : eCM = (pAinit + pBinit).mCalc();
878 :
879 : // Find boost+rotation needed to move from/to CM frame.
880 0 : MfromCM.reset();
881 0 : MfromCM.fromCMframe( pAinit, pBinit);
882 0 : MtoCM = MfromCM;
883 0 : MtoCM.invert();
884 0 : }
885 :
886 : // Fail if CM energy below beam masses.
887 0 : if (eCM < mA + mB) {
888 0 : info.errorMsg("Error in Pythia::initKinematics: too low energy");
889 0 : return false;
890 : }
891 :
892 : // Set up CM-frame kinematics with beams along +-z axis.
893 0 : pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
894 0 : * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
895 0 : pzBcm = -pzAcm;
896 0 : eA = sqrt(mA*mA + pzAcm*pzAcm);
897 0 : eB = sqrt(mB*mB + pzBcm*pzBcm);
898 :
899 : // If in CM frame then store beam four-vectors (else already done above).
900 0 : if (boostType != 2 && boostType != 3) {
901 0 : pAinit = Vec4( 0., 0., pzAcm, eA);
902 0 : pBinit = Vec4( 0., 0., pzBcm, eB);
903 0 : }
904 :
905 : // Store main info for access in process generation.
906 0 : info.setBeamA( idA, pzAcm, eA, mA);
907 0 : info.setBeamB( idB, pzBcm, eB, mB);
908 0 : info.setECM( eCM);
909 :
910 : // Must allow for generic boost+rotation when beam momentum spread.
911 0 : if (doMomentumSpread) boostType = 3;
912 :
913 : // Done.
914 0 : return true;
915 :
916 0 : }
917 :
918 : //--------------------------------------------------------------------------
919 :
920 : // Set up pointers to PDFs.
921 :
922 : bool Pythia::initPDFs() {
923 :
924 : // Delete any PDF's created in a previous init call.
925 0 : if (useNewPdfHard) {
926 0 : if (pdfHardAPtr != pdfAPtr) {
927 0 : delete pdfHardAPtr;
928 0 : pdfHardAPtr = 0;
929 0 : }
930 0 : if (pdfHardBPtr != pdfBPtr) {
931 0 : delete pdfHardBPtr;
932 0 : pdfHardBPtr = 0;
933 0 : }
934 0 : useNewPdfHard = false;
935 0 : }
936 0 : if (useNewPdfA) {
937 0 : delete pdfAPtr;
938 0 : useNewPdfA = false;
939 0 : pdfAPtr = 0;
940 0 : }
941 0 : if (useNewPdfB) {
942 0 : delete pdfBPtr;
943 0 : useNewPdfB = false;
944 0 : pdfBPtr = 0;
945 0 : }
946 0 : if (useNewPdfPomA) {
947 0 : delete pdfPomAPtr;
948 0 : useNewPdfPomA = false;
949 0 : pdfPomAPtr = 0;
950 0 : }
951 0 : if (useNewPdfPomB) {
952 0 : delete pdfPomBPtr;
953 0 : useNewPdfPomB = false;
954 0 : pdfPomBPtr = 0;
955 0 : }
956 :
957 : // Set up the PDF's, if not already done.
958 0 : if (pdfAPtr == 0) {
959 0 : pdfAPtr = getPDFPtr(idA);
960 0 : if (pdfAPtr == 0 || !pdfAPtr->isSetup()) {
961 0 : info.errorMsg("Error in Pythia::init: "
962 : "could not set up PDF for beam A");
963 0 : return false;
964 : }
965 0 : pdfHardAPtr = pdfAPtr;
966 0 : useNewPdfA = true;
967 0 : }
968 0 : if (pdfBPtr == 0) {
969 0 : pdfBPtr = getPDFPtr(idB, 1, "B");
970 0 : if (pdfBPtr == 0 || !pdfBPtr->isSetup()) {
971 0 : info.errorMsg("Error in Pythia::init: "
972 : "could not set up PDF for beam B");
973 0 : return false;
974 : }
975 0 : pdfHardBPtr = pdfBPtr;
976 0 : useNewPdfB = true;
977 0 : }
978 :
979 : // Optionally set up separate PDF's for hard process.
980 0 : if (settings.flag("PDF:useHard") && useNewPdfA && useNewPdfB) {
981 0 : pdfHardAPtr = getPDFPtr(idA, 2);
982 0 : if (!pdfHardAPtr->isSetup()) return false;
983 0 : pdfHardBPtr = getPDFPtr(idB, 2, "B");
984 0 : if (!pdfHardBPtr->isSetup()) return false;
985 0 : useNewPdfHard = true;
986 0 : }
987 :
988 : // Optionally set up Pomeron PDF's for diffractive physics.
989 0 : if ( doDiffraction || doHardDiff) {
990 0 : if (pdfPomAPtr == 0) {
991 0 : pdfPomAPtr = getPDFPtr(990);
992 0 : useNewPdfPomA = true;
993 0 : }
994 0 : if (pdfPomBPtr == 0) {
995 0 : pdfPomBPtr = getPDFPtr(990);
996 0 : useNewPdfPomB = true;
997 0 : }
998 : }
999 :
1000 : // Done.
1001 0 : return true;
1002 :
1003 0 : }
1004 :
1005 : //--------------------------------------------------------------------------
1006 :
1007 : // Main routine to generate the next event, using internal machinery.
1008 :
1009 : bool Pythia::next() {
1010 :
1011 : // Check that constructor worked.
1012 0 : if (!isConstructed) return false;
1013 :
1014 : // Regularly print how many events have been generated.
1015 0 : int nPrevious = info.getCounter(3);
1016 0 : if (nCount > 0 && nPrevious > 0 && nPrevious%nCount == 0)
1017 0 : cout << "\n Pythia::next(): " << nPrevious
1018 0 : << " events have been generated " << endl;
1019 :
1020 : // Set/reset info counters specific to each event.
1021 0 : info.addCounter(3);
1022 0 : for (int i = 10; i < 13; ++i) info.setCounter(i);
1023 :
1024 : // Simpler option when no hard process, i.e. mainly hadron level.
1025 0 : if (!doProcessLevel) {
1026 :
1027 : // Optionally fetch in resonance decays from LHA interface.
1028 0 : if (doLHA && !processLevel.nextLHAdec( event)) {
1029 0 : if (info.atEndOfFile()) info.errorMsg("Abort from Pythia::next:"
1030 : " reached end of Les Houches Events File");
1031 0 : return false;
1032 : }
1033 :
1034 : // Reset info array (while event record contains data).
1035 0 : info.clear();
1036 :
1037 : // Set correct energy for system.
1038 0 : Vec4 pSum = 0.;
1039 0 : for (int i = 1; i < event.size(); ++i)
1040 0 : if (event[i].isFinal()) pSum += event[i].p();
1041 0 : event[0].p( pSum );
1042 0 : event[0].m( pSum.mCalc() );
1043 :
1044 : // Generate hadronization and decays.
1045 0 : bool status = forceHadronLevel();
1046 0 : if (status) info.addCounter(4);
1047 0 : if (status && nPrevious < nShowEvt) event.list(showSaV, showMaD);
1048 : return status;
1049 0 : }
1050 :
1051 : // Reset arrays.
1052 0 : info.clear();
1053 0 : process.clear();
1054 0 : event.clear();
1055 0 : partonSystems.clear();
1056 0 : beamA.clear();
1057 0 : beamB.clear();
1058 0 : beamPomA.clear();
1059 0 : beamPomB.clear();
1060 :
1061 : // Pick current beam valence flavours (for pi0, K0S, K0L, Pomeron).
1062 0 : beamA.newValenceContent();
1063 0 : beamB.newValenceContent();
1064 0 : if ( doDiffraction || doHardDiff) {
1065 0 : beamPomA.newValenceContent();
1066 0 : beamPomB.newValenceContent();
1067 0 : }
1068 :
1069 : // Can only generate event if initialization worked.
1070 0 : if (!isInit) {
1071 0 : info.errorMsg("Abort from Pythia::next: "
1072 : "not properly initialized so cannot generate events");
1073 0 : return false;
1074 : }
1075 :
1076 : // Pick beam momentum spread and beam vertex.
1077 0 : if (doMomentumSpread || doVertexSpread) beamShapePtr->pick();
1078 :
1079 : // Recalculate kinematics when beam momentum spread.
1080 0 : if (doMomentumSpread) nextKinematics();
1081 :
1082 : // Outer loop over hard processes; only relevant for user-set vetoes.
1083 : for ( ; ; ) {
1084 :
1085 0 : info.addCounter(10);
1086 : bool hasVetoed = false;
1087 : bool hasVetoedDiff = false;
1088 :
1089 : // Provide the hard process that starts it off. Only one try.
1090 0 : info.clear();
1091 0 : process.clear();
1092 :
1093 0 : if ( !processLevel.next( process) ) {
1094 0 : if (doLHA && info.atEndOfFile()) info.errorMsg("Abort from "
1095 : "Pythia::next: reached end of Les Houches Events File");
1096 0 : else info.errorMsg("Abort from Pythia::next: "
1097 : "processLevel failed; giving up");
1098 0 : return false;
1099 : }
1100 :
1101 0 : info.addCounter(11);
1102 :
1103 : // Possibility for a user veto of the process-level event.
1104 0 : if (doVetoProcess) {
1105 0 : hasVetoed = userHooksPtr->doVetoProcessLevel( process);
1106 0 : if (hasVetoed) {
1107 0 : if (abortIfVeto) return false;
1108 0 : continue;
1109 : }
1110 : }
1111 :
1112 : // Possibility to perform matrix element merging for this event.
1113 0 : if (doMerging) {
1114 0 : int veto = merging.mergeProcess( process );
1115 : // Apply possible merging scale cut.
1116 0 : if ( veto == -1 ) {
1117 : hasVetoed = true;
1118 0 : if (abortIfVeto) return false;
1119 0 : continue;
1120 : // Exit because of vanishing no-emission probability.
1121 0 : } else if ( veto == 0 ) break;
1122 :
1123 : // Redo resonance decays after the merging, in case the resonance
1124 : // structure has been changed because of reclusterings.
1125 0 : if (veto == 2 && doResDec) processLevel.nextDecays( process);
1126 0 : }
1127 :
1128 : // Possibility to stop the generation at this stage.
1129 0 : if (!doPartonLevel) {
1130 0 : boostAndVertex( true, true);
1131 0 : processLevel.accumulate();
1132 0 : info.addCounter(4);
1133 0 : if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1134 0 : if (nPrevious < nShowInfo) info.list();
1135 0 : if (nPrevious < nShowProc) process.list(showSaV, showMaD);
1136 0 : return true;
1137 : }
1138 :
1139 : // Save spare copy of process record in case of problems.
1140 0 : Event processSave = process;
1141 0 : int sizeMPI = info.sizeMPIarrays();
1142 0 : info.addCounter(12);
1143 0 : for (int i = 14; i < 19; ++i) info.setCounter(i);
1144 :
1145 : // Allow up to ten tries for parton- and hadron-level processing.
1146 : bool physical = true;
1147 0 : for (int iTry = 0; iTry < NTRY; ++iTry) {
1148 :
1149 0 : info.addCounter(14);
1150 : physical = true;
1151 : hasVetoed = false;
1152 :
1153 : // Restore original process record if problems.
1154 0 : if (iTry > 0) process = processSave;
1155 0 : if (iTry > 0) info.resizeMPIarrays( sizeMPI);
1156 :
1157 : // Reset event record and (extracted partons from) beam remnants.
1158 0 : event.clear();
1159 0 : beamA.clear();
1160 0 : beamB.clear();
1161 0 : beamPomA.clear();
1162 0 : beamPomB.clear();
1163 0 : partonSystems.clear();
1164 :
1165 : // Parton-level evolution: ISR, FSR, MPI.
1166 0 : if ( !partonLevel.next( process, event) ) {
1167 :
1168 : // Abort event generation if parton level is set to abort.
1169 0 : if (info.getAbortPartonLevel()) return false;
1170 :
1171 : // Skip to next hard process for failure owing to deliberate veto,
1172 : // or alternatively retry for the same hard process.
1173 0 : hasVetoed = partonLevel.hasVetoed();
1174 0 : if (hasVetoed) {
1175 0 : if (retryPartonLevel) {
1176 0 : --iTry;
1177 0 : continue;
1178 : }
1179 0 : if (abortIfVeto) return false;
1180 0 : break;
1181 : }
1182 :
1183 : // If hard diffractive event has been discarded retry partonLevel.
1184 0 : hasVetoedDiff = partonLevel.hasVetoedDiff();
1185 0 : if (hasVetoedDiff) {
1186 0 : info.errorMsg("Warning in Pythia::next: "
1187 : "discarding hard diffractive event from partonLevel; try again");
1188 0 : break;
1189 : }
1190 :
1191 : // Else make a new try for other failures.
1192 0 : info.errorMsg("Error in Pythia::next: "
1193 : "partonLevel failed; try again");
1194 : physical = false;
1195 0 : continue;
1196 : }
1197 0 : info.addCounter(15);
1198 :
1199 : // Possibility for a user veto of the parton-level event.
1200 0 : if (doVetoPartons) {
1201 0 : hasVetoed = userHooksPtr->doVetoPartonLevel( event);
1202 0 : if (hasVetoed) {
1203 0 : if (abortIfVeto) return false;
1204 0 : break;
1205 : }
1206 : }
1207 :
1208 : // Boost to lab frame (before decays, for vertices).
1209 0 : boostAndVertex( true, true);
1210 :
1211 : // Possibility to stop the generation at this stage.
1212 0 : if (!doHadronLevel) {
1213 0 : processLevel.accumulate();
1214 0 : partonLevel.accumulate();
1215 : // Optionally check final event for problems.
1216 0 : if (checkEvent && !check()) {
1217 0 : info.errorMsg("Abort from Pythia::next: "
1218 : "check of event revealed problems");
1219 0 : return false;
1220 : }
1221 0 : info.addCounter(4);
1222 0 : if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1223 0 : if (nPrevious < nShowInfo) info.list();
1224 0 : if (nPrevious < nShowProc) process.list(showSaV, showMaD);
1225 0 : if (nPrevious < nShowEvt) event.list(showSaV, showMaD);
1226 0 : return true;
1227 : }
1228 :
1229 : // Hadron-level: hadronization, decays.
1230 0 : info.addCounter(16);
1231 0 : if ( !hadronLevel.next( event) ) {
1232 0 : info.errorMsg("Error in Pythia::next: "
1233 : "hadronLevel failed; try again");
1234 : physical = false;
1235 0 : continue;
1236 : }
1237 :
1238 : // If R-hadrons have been formed, then (optionally) let them decay.
1239 0 : if (decayRHadrons && rHadrons.exist() && !doRHadronDecays()) {
1240 0 : info.errorMsg("Error in Pythia::next: "
1241 : "decayRHadrons failed; try again");
1242 : physical = false;
1243 0 : continue;
1244 : }
1245 0 : info.addCounter(17);
1246 :
1247 : // Optionally check final event for problems.
1248 0 : if (checkEvent && !check()) {
1249 0 : info.errorMsg("Error in Pythia::next: "
1250 : "check of event revealed problems");
1251 : physical = false;
1252 0 : continue;
1253 : }
1254 :
1255 : // Stop parton- and hadron-level looping if you got this far.
1256 0 : info.addCounter(18);
1257 0 : break;
1258 : }
1259 :
1260 : // If event vetoed then to make a new try.
1261 0 : if (hasVetoed || hasVetoedDiff) {
1262 0 : if (abortIfVeto) return false;
1263 0 : continue;
1264 : }
1265 :
1266 : // If event failed any other way (after ten tries) then give up.
1267 0 : if (!physical) {
1268 0 : info.errorMsg("Abort from Pythia::next: "
1269 : "parton+hadronLevel failed; giving up");
1270 0 : return false;
1271 : }
1272 :
1273 : // Process- and parton-level statistics. Event scale.
1274 0 : processLevel.accumulate();
1275 0 : partonLevel.accumulate();
1276 0 : event.scale( process.scale() );
1277 :
1278 : // End of outer loop over hard processes. Done with normal option.
1279 0 : info.addCounter(13);
1280 0 : break;
1281 0 : }
1282 :
1283 : // List events.
1284 0 : if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1285 0 : if (nPrevious < nShowInfo) info.list();
1286 0 : if (nPrevious < nShowProc) process.list(showSaV,showMaD);
1287 0 : if (nPrevious < nShowEvt) event.list(showSaV, showMaD);
1288 :
1289 : // Done.
1290 0 : info.addCounter(4);
1291 0 : return true;
1292 :
1293 0 : }
1294 :
1295 : //--------------------------------------------------------------------------
1296 :
1297 : // Generate only the hadronization/decay stage, using internal machinery.
1298 : // The "event" instance should already contain a parton-level configuration.
1299 :
1300 : bool Pythia::forceHadronLevel(bool findJunctions) {
1301 :
1302 : // Can only generate event if initialization worked.
1303 0 : if (!isInit) {
1304 0 : info.errorMsg("Abort from Pythia::forceHadronLevel: "
1305 : "not properly initialized so cannot generate events");
1306 0 : return false;
1307 : }
1308 :
1309 : // Check whether any junctions in system. (Normally done in ProcessLevel.)
1310 : // Avoid it if there are no final-state coloured partons.
1311 0 : if (findJunctions) {
1312 0 : event.clearJunctions();
1313 0 : for (int i = 0; i < event.size(); ++i)
1314 0 : if (event[i].isFinal()
1315 0 : && (event[i].col() != 0 || event[i].acol() != 0)) {
1316 0 : processLevel.findJunctions( event);
1317 0 : break;
1318 : }
1319 0 : }
1320 :
1321 : // Allow for CR before the hadronization.
1322 0 : if (forceHadronLevelCR) {
1323 :
1324 : // Setup parton system for SK-I and SK-II colour reconnection.
1325 : // Require all final state particles to have the Ws as mothers.
1326 0 : if (reconnectMode == 3 || reconnectMode == 4) {
1327 0 : partonSystems.clear();
1328 0 : partonSystems.addSys();
1329 0 : partonSystems.addSys();
1330 0 : for (int i = 5;i < event.size();++i) {
1331 0 : if (event[i].mother1() - 3 < 0 || event[i].mother1() - 3 > 1) {
1332 0 : info.errorMsg("Error from Pythia::forceHadronLevel: "
1333 : " Event is not setup correctly for SK-I or SK-II CR");
1334 0 : return false;
1335 : }
1336 0 : partonSystems.addOut(event[i].mother1() - 3,i);
1337 : }
1338 : }
1339 :
1340 : // save spare copy of event in case of failure.
1341 0 : Event spareEvent = event;
1342 : bool colCorrect = false;
1343 :
1344 : // Allow up to ten tries for CR.
1345 0 : for (int iTry = 0; iTry < NTRY; ++ iTry) {
1346 0 : colourReconnection.next(event, 0);
1347 0 : if (junctionSplitting.checkColours(event)) {
1348 : colCorrect = true;
1349 0 : break;
1350 : }
1351 0 : else event = spareEvent;
1352 : }
1353 :
1354 0 : if (!colCorrect) {
1355 0 : info.errorMsg("Error in Pythia::forceHadronLevel: "
1356 : "Colour reconnection failed.");
1357 0 : return false;
1358 : }
1359 0 : }
1360 :
1361 : // Save spare copy of event in case of failure.
1362 0 : Event spareEvent = event;
1363 :
1364 : // Allow up to ten tries for hadron-level processing.
1365 : bool physical = true;
1366 0 : for (int iTry = 0; iTry < NTRY; ++ iTry) {
1367 : physical = true;
1368 :
1369 : // Check whether any resonances need to be handled at process level.
1370 0 : if (doResDec) {
1371 0 : process = event;
1372 0 : processLevel.nextDecays( process);
1373 :
1374 : // Allow for showers if decays happened at process level.
1375 0 : if (process.size() > event.size()) {
1376 0 : if (doFSRinRes) {
1377 0 : partonLevel.setupShowerSys( process, event);
1378 0 : partonLevel.resonanceShowers( process, event, false);
1379 0 : } else event = process;
1380 : }
1381 : }
1382 :
1383 : // Hadron-level: hadronization, decays.
1384 0 : if (hadronLevel.next( event)) break;
1385 :
1386 : // If failure then warn, restore original configuration and try again.
1387 0 : info.errorMsg("Error in Pythia::forceHadronLevel: "
1388 : "hadronLevel failed; try again");
1389 : physical = false;
1390 0 : event = spareEvent;
1391 : }
1392 :
1393 : // Done for simpler option.
1394 0 : if (!physical) {
1395 0 : info.errorMsg("Abort from Pythia::forceHadronLevel: "
1396 : "hadronLevel failed; giving up");
1397 0 : return false;
1398 : }
1399 :
1400 : // Optionally check final event for problems.
1401 0 : if (checkEvent && !check()) {
1402 0 : info.errorMsg("Abort from Pythia::forceHadronLevel: "
1403 : "check of event revealed problems");
1404 0 : return false;
1405 : }
1406 :
1407 : // Done.
1408 0 : return true;
1409 :
1410 0 : }
1411 :
1412 : //--------------------------------------------------------------------------
1413 :
1414 : // Recalculate kinematics for each event when beam momentum has a spread.
1415 :
1416 : void Pythia::nextKinematics() {
1417 :
1418 : // Read out momentum shift to give current beam momenta.
1419 0 : pAnow = pAinit + beamShapePtr->deltaPA();
1420 0 : pAnow.e( sqrt(pAnow.pAbs2() + mA * mA) );
1421 0 : pBnow = pBinit + beamShapePtr->deltaPB();
1422 0 : pBnow.e( sqrt(pBnow.pAbs2() + mB * mB) );
1423 :
1424 : // Construct CM frame kinematics.
1425 0 : eCM = (pAnow + pBnow).mCalc();
1426 0 : pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
1427 0 : * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
1428 0 : pzBcm = -pzAcm;
1429 0 : eA = sqrt(mA*mA + pzAcm*pzAcm);
1430 0 : eB = sqrt(mB*mB + pzBcm*pzBcm);
1431 :
1432 : // Set relevant info for other classes to use.
1433 0 : info.setBeamA( idA, pzAcm, eA, mA);
1434 0 : info.setBeamB( idB, pzBcm, eB, mB);
1435 0 : info.setECM( eCM);
1436 0 : beamA.newPzE( pzAcm, eA);
1437 0 : beamB.newPzE( pzBcm, eB);
1438 :
1439 : // Set boost/rotation matrices from/to CM frame.
1440 0 : MfromCM.reset();
1441 0 : MfromCM.fromCMframe( pAnow, pBnow);
1442 0 : MtoCM = MfromCM;
1443 0 : MtoCM.invert();
1444 :
1445 0 : }
1446 :
1447 : //--------------------------------------------------------------------------
1448 :
1449 : // Boost from CM frame to lab frame, or inverse. Set production vertex.
1450 :
1451 : void Pythia::boostAndVertex( bool toLab, bool setVertex) {
1452 :
1453 : // Boost process from CM frame to lab frame.
1454 0 : if (toLab) {
1455 0 : if (boostType == 2) process.bst(0., 0., betaZ, gammaZ);
1456 0 : else if (boostType == 3) process.rotbst(MfromCM);
1457 :
1458 : // Boost nonempty event from CM frame to lab frame.
1459 0 : if (event.size() > 0) {
1460 0 : if (boostType == 2) event.bst(0., 0., betaZ, gammaZ);
1461 0 : else if (boostType == 3) event.rotbst(MfromCM);
1462 : }
1463 :
1464 : // Boost process from lab frame to CM frame.
1465 : } else {
1466 0 : if (boostType == 2) process.bst(0., 0., -betaZ, gammaZ);
1467 0 : else if (boostType == 3) process.rotbst(MtoCM);
1468 :
1469 : // Boost nonempty event from lab frame to CM frame.
1470 0 : if (event.size() > 0) {
1471 0 : if (boostType == 2) event.bst(0., 0., -betaZ, gammaZ);
1472 0 : else if (boostType == 3) event.rotbst(MtoCM);
1473 : }
1474 : }
1475 :
1476 : // Set production vertex; assumes particles are in lab frame and at origin.
1477 0 : if (setVertex && doVertexSpread) {
1478 0 : Vec4 vertex = beamShapePtr->vertex();
1479 0 : for (int i = 0; i < process.size(); ++i) process[i].vProd( vertex);
1480 0 : for (int i = 0; i < event.size(); ++i) event[i].vProd( vertex);
1481 0 : }
1482 :
1483 0 : }
1484 :
1485 : //--------------------------------------------------------------------------
1486 :
1487 : // Perform R-hadron decays, either as part of normal evolution or forced.
1488 :
1489 : bool Pythia::doRHadronDecays( ) {
1490 :
1491 : // Check if R-hadrons exist to be processed.
1492 0 : if ( !rHadrons.exist() ) return true;
1493 :
1494 : // Do the R-hadron decay itself.
1495 0 : if ( !rHadrons.decay( event) ) return false;
1496 :
1497 : // Perform showers in resonance decay chains.
1498 0 : if ( !partonLevel.resonanceShowers( process, event, false) ) return false;
1499 :
1500 : // Subsequent hadronization and decays.
1501 0 : if ( !hadronLevel.next( event) ) return false;
1502 :
1503 : // Done.
1504 0 : return true;
1505 :
1506 0 : }
1507 :
1508 : //--------------------------------------------------------------------------
1509 :
1510 : // Print statistics on event generation.
1511 :
1512 : void Pythia::stat() {
1513 :
1514 : // Read out settings for what to include.
1515 0 : bool showPrL = settings.flag("Stat:showProcessLevel");
1516 0 : bool showPaL = settings.flag("Stat:showPartonLevel");
1517 0 : bool showErr = settings.flag("Stat:showErrors");
1518 0 : bool reset = settings.flag("Stat:reset");
1519 :
1520 : // Statistics on cross section and number of events.
1521 0 : if (doProcessLevel) {
1522 0 : if (showPrL) processLevel.statistics(false);
1523 0 : if (reset) processLevel.resetStatistics();
1524 : }
1525 :
1526 : // Statistics from other classes, currently multiparton interactions.
1527 0 : if (showPaL) partonLevel.statistics(false);
1528 0 : if (reset) partonLevel.resetStatistics();
1529 :
1530 : // Merging statistics.
1531 0 : if (doMerging) merging.statistics();
1532 :
1533 : // Summary of which and how many warnings/errors encountered.
1534 0 : if (showErr) info.errorStatistics();
1535 0 : if (reset) info.errorReset();
1536 :
1537 0 : }
1538 :
1539 : //--------------------------------------------------------------------------
1540 :
1541 : // Write the Pythia banner, with symbol and version information.
1542 :
1543 : void Pythia::banner(ostream& os) {
1544 :
1545 : // Read in version number and last date of change.
1546 0 : double versionNumber = settings.parm("Pythia:versionNumber");
1547 0 : int versionDate = settings.mode("Pythia:versionDate");
1548 0 : string month[12] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun",
1549 0 : "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
1550 :
1551 : // Get date and time.
1552 0 : time_t t = time(0);
1553 0 : char dateNow[12];
1554 0 : strftime(dateNow,12,"%d %b %Y",localtime(&t));
1555 0 : char timeNow[9];
1556 0 : strftime(timeNow,9,"%H:%M:%S",localtime(&t));
1557 :
1558 0 : os << "\n"
1559 0 : << " *-------------------------------------------"
1560 0 : << "-----------------------------------------* \n"
1561 0 : << " | "
1562 0 : << " | \n"
1563 0 : << " | *----------------------------------------"
1564 0 : << "--------------------------------------* | \n"
1565 0 : << " | | "
1566 0 : << " | | \n"
1567 0 : << " | | "
1568 0 : << " | | \n"
1569 0 : << " | | PPP Y Y TTTTT H H III A "
1570 0 : << " Welcome to the Lund Monte Carlo! | | \n"
1571 0 : << " | | P P Y Y T H H I A A "
1572 0 : << " This is PYTHIA version " << fixed << setprecision(3)
1573 0 : << setw(5) << versionNumber << " | | \n"
1574 0 : << " | | PPP Y T HHHHH I AAAAA"
1575 0 : << " Last date of change: " << setw(2) << versionDate%100
1576 0 : << " " << month[ (versionDate/100)%100 - 1 ]
1577 0 : << " " << setw(4) << versionDate/10000 << " | | \n"
1578 0 : << " | | P Y T H H I A A"
1579 0 : << " | | \n"
1580 0 : << " | | P Y T H H III A A"
1581 0 : << " Now is " << dateNow << " at " << timeNow << " | | \n"
1582 0 : << " | | "
1583 0 : << " | | \n"
1584 0 : << " | | Torbjorn Sjostrand; Department of As"
1585 0 : << "tronomy and Theoretical Physics, | | \n"
1586 0 : << " | | Lund University, Solvegatan 14A, S"
1587 0 : << "E-223 62 Lund, Sweden; | | \n"
1588 0 : << " | | e-mail: torbjorn@thep.lu.se "
1589 0 : << " | | \n"
1590 0 : << " | | Jesper Roy Christiansen; Department "
1591 0 : << "of Astronomy and Theoretical Physics, | | \n"
1592 0 : << " | | Lund University, Solvegatan 14A, S"
1593 0 : << "E-223 62 Lund, Sweden; | | \n"
1594 0 : << " | | e-mail: Jesper.Roy.Christiansen@th"
1595 0 : << "ep.lu.se | | \n"
1596 0 : << " | | Nishita Desai; Institut fuer Theoret"
1597 0 : << "ische Physik, | | \n"
1598 0 : << " | | Universitaet Heidelberg, Philosophe"
1599 0 : << "nweg 16, D-69120 Heidelberg, Germany; | | \n"
1600 0 : << " | | e-mail: n.desai@thphys.uni-heidelb"
1601 0 : << "erg.de | | \n"
1602 0 : << " | | Philip Ilten; Massachusetts Institut"
1603 0 : << "e of Technology, | | \n"
1604 0 : << " | | stationed at CERN, CH-1211 Geneva "
1605 0 : << "23, Switzerland; | | \n"
1606 0 : << " | | e-mail: philten@cern.ch "
1607 0 : << " | | \n"
1608 0 : << " | | Stephen Mrenna; Computing Division, "
1609 0 : << "Simulations Group, | | \n"
1610 0 : << " | | Fermi National Accelerator Laborat"
1611 0 : << "ory, MS 234, Batavia, IL 60510, USA; | | \n"
1612 0 : << " | | e-mail: mrenna@fnal.gov "
1613 0 : << " | | \n"
1614 0 : << " | | Stefan Prestel; Theoretical Physics "
1615 0 : << "Group, | | \n"
1616 0 : << " | | SLAC National Accelerator Laborato"
1617 0 : << "ry, Menlo Park, CA 94025, USA; | | \n"
1618 0 : << " | | e-mail: prestel@slac.stanford.edu "
1619 0 : << " | | \n"
1620 0 : << " | | Christine O. Rasmussen; Department o"
1621 0 : << "f Astronomy and Theoretical Physics, | | \n"
1622 0 : << " | | Lund University, Solvegatan 14A, S"
1623 0 : << "E-223 62 Lund, Sweden; | | \n"
1624 0 : << " | | e-mail: christine.rasmussen@thep.l"
1625 0 : << "u.se | | \n"
1626 0 : << " | | Peter Skands; School of Physics, "
1627 0 : << " | | \n"
1628 0 : << " | | Monash University, PO Box 27, 3800"
1629 0 : << " Melbourne, Australia; | | \n"
1630 0 : << " | | e-mail: peter.skands@monash.edu "
1631 0 : << " | | \n"
1632 0 : << " | | "
1633 0 : << " | | \n"
1634 0 : << " | | The main program reference is 'An Int"
1635 0 : << "roduction to PYTHIA 8.2', | | \n"
1636 0 : << " | | T. Sjostrand et al, Comput. Phys. Com"
1637 0 : << "mun. 191 (2005) 159 | | \n"
1638 0 : << " | | [arXiv:1410.3012 [hep-ph]] "
1639 0 : << " | | \n"
1640 0 : << " | | "
1641 0 : << " | | \n"
1642 0 : << " | | The main physics reference is the 'PY"
1643 0 : << "THIA 6.4 Physics and Manual', | | \n"
1644 0 : << " | | T. Sjostrand, S. Mrenna and P. Skands"
1645 0 : << ", JHEP05 (2006) 026 [hep-ph/0603175] | | \n"
1646 0 : << " | | "
1647 0 : << " | | \n"
1648 0 : << " | | An archive of program versions and do"
1649 0 : << "cumentation is found on the web: | | \n"
1650 0 : << " | | http://www.thep.lu.se/Pythia "
1651 0 : << " | | \n"
1652 0 : << " | | "
1653 0 : << " | | \n"
1654 0 : << " | | This program is released under the GN"
1655 0 : << "U General Public Licence version 2. | | \n"
1656 0 : << " | | Please respect the MCnet Guidelines f"
1657 0 : << "or Event Generator Authors and Users. | | \n"
1658 0 : << " | | "
1659 0 : << " | | \n"
1660 0 : << " | | Disclaimer: this program comes withou"
1661 0 : << "t any guarantees. | | \n"
1662 0 : << " | | Beware of errors and use common sense"
1663 0 : << " when interpreting results. | | \n"
1664 0 : << " | | "
1665 0 : << " | | \n"
1666 0 : << " | | Copyright (C) 2015 Torbjorn Sjostrand"
1667 0 : << " | | \n"
1668 0 : << " | | "
1669 0 : << " | | \n"
1670 0 : << " | | "
1671 0 : << " | | \n"
1672 0 : << " | *----------------------------------------"
1673 0 : << "--------------------------------------* | \n"
1674 0 : << " | "
1675 0 : << " | \n"
1676 0 : << " *-------------------------------------------"
1677 0 : << "-----------------------------------------* \n" << endl;
1678 :
1679 0 : }
1680 :
1681 : //--------------------------------------------------------------------------
1682 :
1683 : // Check for lines in file that mark the beginning of new subrun.
1684 :
1685 : int Pythia::readSubrun(string line, bool warn, ostream& os) {
1686 :
1687 : // If empty line then done.
1688 0 : int subrunLine = SUBRUNDEFAULT;
1689 0 : if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos)
1690 0 : return subrunLine;
1691 :
1692 : // If first character is not a letter, then done.
1693 0 : string lineNow = line;
1694 0 : int firstChar = lineNow.find_first_not_of(" \n\t\v\b\r\f\a");
1695 0 : if (!isalpha(lineNow[firstChar])) return subrunLine;
1696 :
1697 : // Replace an equal sign by a blank to make parsing simpler.
1698 0 : while (lineNow.find("=") != string::npos) {
1699 0 : int firstEqual = lineNow.find_first_of("=");
1700 0 : lineNow.replace(firstEqual, 1, " ");
1701 : }
1702 :
1703 : // Get first word of a line.
1704 0 : istringstream splitLine(lineNow);
1705 0 : string name;
1706 0 : splitLine >> name;
1707 :
1708 : // Replace two colons by one (:: -> :) to allow for such mistakes.
1709 0 : while (name.find("::") != string::npos) {
1710 0 : int firstColonColon = name.find_first_of("::");
1711 0 : name.replace(firstColonColon, 2, ":");
1712 : }
1713 :
1714 :
1715 : // Convert to lowercase.
1716 0 : for (int i = 0; i < int(name.length()); ++i) name[i] = tolower(name[i]);
1717 :
1718 : // If no match then done.
1719 0 : if (name != "main:subrun") return subrunLine;
1720 :
1721 : // Else find new subrun number and return it.
1722 0 : splitLine >> subrunLine;
1723 0 : if (!splitLine) {
1724 0 : if (warn) os << "\n PYTHIA Warning: Main:subrun number not"
1725 0 : << " recognized; skip:\n " << line << endl;
1726 0 : subrunLine = SUBRUNDEFAULT;
1727 0 : }
1728 0 : return subrunLine;
1729 :
1730 0 : }
1731 :
1732 : //--------------------------------------------------------------------------
1733 :
1734 : // Check for lines in file that mark the beginning or end of commented section.
1735 : // Return +1 for beginning, -1 for end, 0 else.
1736 :
1737 : int Pythia::readCommented(string line) {
1738 :
1739 : // If less than two nontrivial characters on line then done.
1740 0 : if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return 0;
1741 0 : int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
1742 0 : if (int(line.size()) < firstChar + 2) return 0;
1743 :
1744 : // If first two nontrivial characters are /* or */ then done.
1745 0 : if (line.substr(firstChar, 2) == "/*") return +1;
1746 0 : if (line.substr(firstChar, 2) == "*/") return -1;
1747 :
1748 : // Else done.
1749 0 : return 0;
1750 :
1751 0 : }
1752 :
1753 : //--------------------------------------------------------------------------
1754 :
1755 : // Check that the final event makes sense: no unknown id codes;
1756 : // charge and energy-momentum conserved.
1757 :
1758 : bool Pythia::check(ostream& os) {
1759 :
1760 : // Reset.
1761 : bool physical = true;
1762 : bool listVertices = false;
1763 : bool listHistory = false;
1764 : bool listSystems = false;
1765 : bool listBeams = false;
1766 0 : iErrId.resize(0);
1767 0 : iErrCol.resize(0);
1768 0 : iErrEpm.resize(0);
1769 0 : iErrNan.resize(0);
1770 0 : iErrNanVtx.resize(0);
1771 0 : Vec4 pSum;
1772 : double chargeSum = 0.;
1773 :
1774 : // Incoming beams counted with negative momentum and charge.
1775 0 : if (doProcessLevel) {
1776 0 : pSum = - (event[1].p() + event[2].p());
1777 0 : chargeSum = - (event[1].charge() + event[2].charge());
1778 :
1779 : // If no ProcessLevel then sum final state of process record.
1780 0 : } else if (process.size() > 0) {
1781 0 : pSum = - process[0].p();
1782 0 : for (int i = 0; i < process.size(); ++i)
1783 0 : if (process[i].isFinal()) chargeSum -= process[i].charge();
1784 :
1785 : // If process not filled, then use outgoing primary in event.
1786 0 : } else {
1787 0 : pSum = - event[0].p();
1788 0 : for (int i = 1; i < event.size(); ++i)
1789 0 : if (event[i].statusAbs() < 10 || event[i].statusAbs() == 23)
1790 0 : chargeSum -= event[i].charge();
1791 : }
1792 0 : double eLab = abs(pSum.e());
1793 :
1794 : // Loop over particles in the event.
1795 0 : for (int i = 0; i < event.size(); ++i) {
1796 :
1797 : // Look for any unrecognized particle codes.
1798 0 : int id = event[i].id();
1799 0 : if (id == 0 || !particleData.isParticle(id)) {
1800 0 : ostringstream errCode;
1801 0 : errCode << ", i = " << i << ", id = " << id;
1802 0 : info.errorMsg("Error in Pythia::check: "
1803 0 : "unknown particle code", errCode.str());
1804 : physical = false;
1805 0 : iErrId.push_back(i);
1806 :
1807 : // Check that colour assignments are the expected ones.
1808 0 : } else {
1809 0 : int colType = event[i].colType();
1810 0 : int col = event[i].col();
1811 0 : int acol = event[i].acol();
1812 0 : if ( (colType == 0 && (col > 0 || acol > 0))
1813 0 : || (colType == 1 && (col <= 0 || acol > 0))
1814 0 : || (colType == -1 && (col > 0 || acol <= 0))
1815 0 : || (colType == 2 && (col <= 0 || acol <= 0)) ) {
1816 0 : ostringstream errCode;
1817 0 : errCode << ", i = " << i << ", id = " << id << " cols = " << col
1818 0 : << " " << acol;
1819 0 : info.errorMsg("Error in Pythia::check: "
1820 0 : "incorrect colours", errCode.str());
1821 : physical = false;
1822 0 : iErrCol.push_back(i);
1823 0 : }
1824 : }
1825 :
1826 : // Look for particles with mismatched or not-a-number energy/momentum/mass.
1827 0 : if (abs(event[i].px()) >= 0. && abs(event[i].py()) >= 0.
1828 0 : && abs(event[i].pz()) >= 0. && abs(event[i].e()) >= 0.
1829 0 : && abs(event[i].m()) >= 0.) {
1830 0 : double errMass = abs(event[i].mCalc() - event[i].m())
1831 0 : / max( 1.0, event[i].e());
1832 0 : if (errMass > mTolErr) {
1833 0 : info.errorMsg("Error in Pythia::check: "
1834 : "unmatched particle energy/momentum/mass");
1835 : physical = false;
1836 0 : iErrEpm.push_back(i);
1837 0 : } else if (errMass > mTolWarn) {
1838 0 : info.errorMsg("Warning in Pythia::check: "
1839 : "not quite matched particle energy/momentum/mass");
1840 0 : }
1841 0 : } else {
1842 0 : info.errorMsg("Error in Pythia::check: "
1843 : "not-a-number energy/momentum/mass");
1844 : physical = false;
1845 0 : iErrNan.push_back(i);
1846 : }
1847 :
1848 : // Look for particles with not-a-number vertex/lifetime.
1849 0 : if (abs(event[i].xProd()) >= 0. && abs(event[i].yProd()) >= 0.
1850 0 : && abs(event[i].zProd()) >= 0. && abs(event[i].tProd()) >= 0.
1851 0 : && abs(event[i].tau()) >= 0.) ;
1852 : else {
1853 0 : info.errorMsg("Error in Pythia::check: "
1854 : "not-a-number vertex/lifetime");
1855 : physical = false;
1856 : listVertices = true;
1857 0 : iErrNanVtx.push_back(i);
1858 : }
1859 :
1860 : // Add final-state four-momentum and charge.
1861 0 : if (event[i].isFinal()) {
1862 0 : pSum += event[i].p();
1863 0 : chargeSum += event[i].charge();
1864 0 : }
1865 :
1866 : // End of particle loop.
1867 : }
1868 :
1869 : // Check energy-momentum/charge conservation.
1870 0 : double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
1871 0 : + abs(pSum.pz());
1872 0 : if (epDev > epTolErr * eLab) {
1873 0 : info.errorMsg("Error in Pythia::check: energy-momentum not conserved");
1874 : physical = false;
1875 0 : } else if (epDev > epTolWarn * eLab) {
1876 0 : info.errorMsg("Warning in Pythia::check: "
1877 : "energy-momentum not quite conserved");
1878 0 : }
1879 0 : if (abs(chargeSum) > 0.1) {
1880 0 : info.errorMsg("Error in Pythia::check: charge not conserved");
1881 : physical = false;
1882 0 : }
1883 :
1884 : // Check that beams and event records agree on incoming partons.
1885 : // Only meaningful for resolved beams.
1886 0 : if (info.isResolved() && !info.hasUnresolvedBeams())
1887 0 : for (int iSys = 0; iSys < beamA.sizeInit(); ++iSys) {
1888 0 : int eventANw = partonSystems.getInA(iSys);
1889 0 : int eventBNw = partonSystems.getInB(iSys);
1890 0 : int beamANw = beamA[iSys].iPos();
1891 0 : int beamBNw = beamB[iSys].iPos();
1892 0 : if (eventANw != beamANw || eventBNw != beamBNw) {
1893 0 : info.errorMsg("Error in Pythia::check: "
1894 : "event and beams records disagree");
1895 : physical = false;
1896 : listSystems = true;
1897 : listBeams = true;
1898 0 : }
1899 0 : }
1900 :
1901 : // Check that mother and daughter information match for each particle.
1902 0 : vector<int> noMot;
1903 0 : vector<int> noDau;
1904 0 : vector< pair<int,int> > noMotDau;
1905 0 : if (checkHistory) {
1906 :
1907 : // Loop through the event and check that there are beam particles.
1908 : bool hasBeams = false;
1909 0 : for (int i = 0; i < event.size(); ++i) {
1910 0 : int status = event[i].status();
1911 0 : if (abs(status) == 12) hasBeams = true;
1912 :
1913 : // Check that mother and daughter lists not empty where not expected to.
1914 0 : vector<int> mList = event[i].motherList();
1915 0 : vector<int> dList = event[i].daughterList();
1916 0 : if (mList.size() == 0 && abs(status) != 11 && abs(status) != 12)
1917 0 : noMot.push_back(i);
1918 0 : if (dList.size() == 0 && status < 0 && status != -11)
1919 0 : noDau.push_back(i);
1920 :
1921 : // Check that the particle appears in the daughters list of each mother.
1922 0 : for (int j = 0; j < int(mList.size()); ++j) {
1923 0 : if ( event[mList[j]].daughter1() <= i
1924 0 : && event[mList[j]].daughter2() >= i ) continue;
1925 0 : vector<int> dmList = event[mList[j]].daughterList();
1926 : bool foundMatch = false;
1927 0 : for (int k = 0; k < int(dmList.size()); ++k)
1928 0 : if (dmList[k] == i) {
1929 : foundMatch = true;
1930 0 : break;
1931 : }
1932 0 : if (!hasBeams && mList.size() == 1 && mList[0] == 0) foundMatch = true;
1933 0 : if (!foundMatch) {
1934 : bool oldPair = false;
1935 0 : for (int k = 0; k < int(noMotDau.size()); ++k)
1936 0 : if (noMotDau[k].first == mList[j] && noMotDau[k].second == i) {
1937 : oldPair = true;
1938 0 : break;
1939 : }
1940 0 : if (!oldPair) noMotDau.push_back( make_pair( mList[j], i) );
1941 0 : }
1942 0 : }
1943 :
1944 : // Check that the particle appears in the mothers list of each daughter.
1945 0 : for (int j = 0; j < int(dList.size()); ++j) {
1946 0 : if ( event[dList[j]].statusAbs() > 80
1947 0 : && event[dList[j]].statusAbs() < 90
1948 0 : && event[dList[j]].mother1() <= i
1949 0 : && event[dList[j]].mother2() >= i) continue;
1950 0 : vector<int> mdList = event[dList[j]].motherList();
1951 : bool foundMatch = false;
1952 0 : for (int k = 0; k < int(mdList.size()); ++k)
1953 0 : if (mdList[k] == i) {
1954 : foundMatch = true;
1955 0 : break;
1956 : }
1957 0 : if (!foundMatch) {
1958 : bool oldPair = false;
1959 0 : for (int k = 0; k < int(noMotDau.size()); ++k)
1960 0 : if (noMotDau[k].first == i && noMotDau[k].second == dList[j]) {
1961 : oldPair = true;
1962 0 : break;
1963 : }
1964 0 : if (!oldPair) noMotDau.push_back( make_pair( i, dList[j]) );
1965 0 : }
1966 0 : }
1967 0 : }
1968 :
1969 : // Warn if any errors were found.
1970 0 : if (noMot.size() > 0 || noDau.size() > 0 || noMotDau.size() > 0) {
1971 0 : info.errorMsg("Error in Pythia::check: "
1972 : "mismatch in daughter and mother lists");
1973 : physical = false;
1974 : listHistory = true;
1975 0 : }
1976 0 : }
1977 :
1978 : // Done for sensible events.
1979 0 : if (physical) return true;
1980 :
1981 : // Print (the first few) flawed events: local info.
1982 0 : if (nErrEvent < nErrList) {
1983 0 : os << "\n PYTHIA erroneous event info: \n";
1984 0 : if (iErrId.size() > 0) {
1985 0 : os << " unknown particle codes in lines ";
1986 0 : for (int i = 0; i < int(iErrId.size()); ++i)
1987 0 : os << iErrId[i] << " ";
1988 0 : os << "\n";
1989 : }
1990 0 : if (iErrCol.size() > 0) {
1991 0 : os << " incorrect colour assignments in lines ";
1992 0 : for (int i = 0; i < int(iErrCol.size()); ++i)
1993 0 : os << iErrCol[i] << " ";
1994 0 : os << "\n";
1995 : }
1996 0 : if (iErrEpm.size() > 0) {
1997 0 : os << " mismatch between energy/momentum/mass in lines ";
1998 0 : for (int i = 0; i < int(iErrEpm.size()); ++i)
1999 0 : os << iErrEpm[i] << " ";
2000 0 : os << "\n";
2001 : }
2002 0 : if (iErrNan.size() > 0) {
2003 0 : os << " not-a-number energy/momentum/mass in lines ";
2004 0 : for (int i = 0; i < int(iErrNan.size()); ++i)
2005 0 : os << iErrNan[i] << " ";
2006 0 : os << "\n";
2007 : }
2008 0 : if (iErrNanVtx.size() > 0) {
2009 0 : os << " not-a-number vertex/lifetime in lines ";
2010 0 : for (int i = 0; i < int(iErrNanVtx.size()); ++i)
2011 0 : os << iErrNanVtx[i] << " ";
2012 0 : os << "\n";
2013 : }
2014 0 : if (epDev > epTolErr * eLab) os << scientific << setprecision(3)
2015 0 : << " total energy-momentum non-conservation = " << epDev << "\n";
2016 0 : if (abs(chargeSum) > 0.1) os << fixed << setprecision(2)
2017 0 : << " total charge non-conservation = " << chargeSum << "\n";
2018 0 : if (noMot.size() > 0) {
2019 0 : os << " missing mothers for particles ";
2020 0 : for (int i = 0; i < int(noMot.size()); ++i) os << noMot[i] << " ";
2021 0 : os << "\n";
2022 : }
2023 0 : if (noDau.size() > 0) {
2024 0 : os << " missing daughters for particles ";
2025 0 : for (int i = 0; i < int(noDau.size()); ++i) os << noDau[i] << " ";
2026 0 : os << "\n";
2027 : }
2028 0 : if (noMotDau.size() > 0) {
2029 0 : os << " inconsistent history for (mother,daughter) pairs ";
2030 0 : for (int i = 0; i < int(noMotDau.size()); ++i)
2031 0 : os << "(" << noMotDau[i].first << "," << noMotDau[i].second << ") ";
2032 0 : os << "\n";
2033 : }
2034 :
2035 : // Print (the first few) flawed events: standard listings.
2036 0 : info.list();
2037 0 : event.list(listVertices, listHistory);
2038 0 : if (listSystems) partonSystems.list();
2039 0 : if (listBeams) beamA.list();
2040 0 : if (listBeams) beamB.list();
2041 : }
2042 :
2043 : // Update error counter. Done also for flawed event.
2044 0 : ++nErrEvent;
2045 0 : return false;
2046 :
2047 0 : }
2048 :
2049 : //--------------------------------------------------------------------------
2050 :
2051 : // Routine to set up a PDF pointer.
2052 :
2053 : PDF* Pythia::getPDFPtr(int idIn, int sequence, string beam) {
2054 :
2055 : // Temporary pointer to be returned.
2056 : PDF* tempPDFPtr = 0;
2057 :
2058 : // One option is to treat a Pomeron like a pi0.
2059 0 : if (idIn == 990 && settings.mode("PDF:PomSet") == 2) idIn = 111;
2060 :
2061 : // Proton beam, normal or hard choice. Also used for neutron.
2062 0 : if (abs(idIn) == 2212 || abs(idIn) == 2112) {
2063 0 : string pSet = settings.word("PDF:p"
2064 0 : + string(sequence == 1 ? "" : "Hard") + "Set" + beam);
2065 0 : if (pSet == "void" && sequence != 1 && beam == "B")
2066 0 : pSet = settings.word("PDF:pHardSet");
2067 0 : if (pSet == "void") pSet = settings.word("PDF:pSet");
2068 0 : istringstream pSetStream(pSet);
2069 0 : int pSetInt(0);
2070 0 : pSetStream >> pSetInt;
2071 :
2072 : // Use sets from LHAPDF.
2073 0 : if (pSetInt == 0)
2074 0 : tempPDFPtr = new LHAPDF(idIn, pSet, &info);
2075 :
2076 : // Use internal sets.
2077 0 : else if (pSetInt == 1) tempPDFPtr = new GRV94L(idIn);
2078 0 : else if (pSetInt == 2) tempPDFPtr = new CTEQ5L(idIn);
2079 0 : else if (pSetInt <= 6)
2080 0 : tempPDFPtr = new MSTWpdf(idIn, pSetInt - 2, xmlPath, &info);
2081 0 : else if (pSetInt <= 12)
2082 0 : tempPDFPtr = new CTEQ6pdf(idIn, pSetInt - 6, xmlPath, &info);
2083 0 : else if (pSetInt <= 16)
2084 0 : tempPDFPtr = new NNPDF(idIn, pSetInt - 12, xmlPath, &info);
2085 : else tempPDFPtr = 0;
2086 0 : }
2087 :
2088 : // Pion beam (or, in one option, Pomeron beam).
2089 0 : else if (abs(idIn) == 211 || idIn == 111) {
2090 0 : string pSet = settings.word("PDF:piSet" + beam);
2091 0 : istringstream pSetStream(pSet);
2092 0 : int pSetInt(0);
2093 0 : pSetStream >> pSetInt;
2094 :
2095 : // Use sets from LHAPDF.
2096 0 : if (pSetInt == 0)
2097 0 : tempPDFPtr = new LHAPDF(idIn, pSet, &info);
2098 :
2099 : // Use internal set.
2100 0 : else if (pSetInt == 1) tempPDFPtr = new GRVpiL(idIn);
2101 : else tempPDFPtr = 0;
2102 0 : }
2103 :
2104 : // Pomeron beam, if not treated like a pi0 beam.
2105 0 : else if (idIn == 990) {
2106 0 : int pomSet = settings.mode("PDF:PomSet");
2107 0 : double rescale = settings.parm("PDF:PomRescale");
2108 :
2109 : // A generic Q2-independent parametrization.
2110 0 : if (pomSet == 1) {
2111 0 : double gluonA = settings.parm("PDF:PomGluonA");
2112 0 : double gluonB = settings.parm("PDF:PomGluonB");
2113 0 : double quarkA = settings.parm("PDF:PomQuarkA");
2114 0 : double quarkB = settings.parm("PDF:PomQuarkB");
2115 0 : double quarkFrac = settings.parm("PDF:PomQuarkFrac");
2116 0 : double strangeSupp = settings.parm("PDF:PomStrangeSupp");
2117 0 : tempPDFPtr = new PomFix( 990, gluonA, gluonB, quarkA, quarkB,
2118 : quarkFrac, strangeSupp);
2119 0 : }
2120 :
2121 : // The H1 Q2-dependent parametrizations. Initialization requires files.
2122 0 : else if (pomSet == 3 || pomSet == 4)
2123 0 : tempPDFPtr = new PomH1FitAB( 990, pomSet - 2, rescale, xmlPath, &info);
2124 0 : else if (pomSet == 5)
2125 0 : tempPDFPtr = new PomH1Jets( 990, rescale, xmlPath, &info);
2126 0 : else if (pomSet == 6)
2127 0 : tempPDFPtr = new PomH1FitAB( 990, 3, rescale, xmlPath, &info);
2128 0 : }
2129 :
2130 : // Lepton beam: neutrino, resolved charged lepton or unresolved ditto.
2131 0 : else if (abs(idIn) > 10 && abs(idIn) < 17) {
2132 0 : if (abs(idIn)%2 == 0) tempPDFPtr = new NeutrinoPoint(idIn);
2133 0 : else if (settings.flag("PDF:lepton")) tempPDFPtr = new Lepton(idIn);
2134 0 : else tempPDFPtr = new LeptonPoint(idIn);
2135 : }
2136 :
2137 : // Optionally allow extrapolation beyond x and Q2 limits.
2138 0 : if (tempPDFPtr)
2139 0 : tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolate") );
2140 :
2141 : // Done.
2142 0 : return tempPDFPtr;
2143 0 : }
2144 :
2145 : //==========================================================================
2146 :
2147 : } // end namespace Pythia8
|
