Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

Pythia.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago
first import of structure, PYTHIA8 and DELPHES
File size: 81.0 KB

Line
1	// Pythia.cc is a part of the PYTHIA event generator.
2	// Copyright (C) 2012 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 "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::VERSIONNUMBERCODE = 8.170;
26
27	//--------------------------------------------------------------------------
28
29	// Constants: could be changed here if desired, but normally should not.
30	// These are of technical nature, as described for each.
31
32	// Maximum number of tries to produce parton level from given input.
33	const int Pythia::NTRY = 10;
34
35	// Negative integer to denote that no subrun has been set.
36	const int Pythia::SUBRUNDEFAULT = -999;
37
38	//--------------------------------------------------------------------------
39
40	// Constructor.
41
42	Pythia::Pythia(string xmlDir) {
43
44	// Initial values for pointers to PDF's.
45	useNewPdfA = false;
46	useNewPdfB = false;
47	useNewPdfHard = false;
48	useNewPdfPomA = false;
49	useNewPdfPomB = false;
50	pdfAPtr = 0;
51	pdfBPtr = 0;
52	pdfHardAPtr = 0;
53	pdfHardBPtr = 0;
54	pdfPomAPtr = 0;
55	pdfPomBPtr = 0;
56
57	// Initial values for pointers to Les Houches Event objects.
58	doLHA = false;
59	useNewLHA = false;
60	lhaUpPtr = 0;
61
62	//Initial value for couplings pointer
63	couplingsPtr = &couplings;
64
65	// Initial value for pointer to external decay handler.
66	decayHandlePtr = 0;
67
68	// Initial value for pointer to user hooks.
69	userHooksPtr = 0;
70
71	// Initial value for pointer to merging hooks.
72	doMerging = false;
73	hasMergingHooks = false;
74	hasOwnMergingHooks = false;
75	mergingHooksPtr = 0;
76
77	// Initial value for pointer to beam shape.
78	useNewBeamShape = false;
79	beamShapePtr = 0;
80
81	// Initial values for pointers to timelike and spacelike showers.
82	useNewTimes = false;
83	useNewSpace = false;
84	timesDecPtr = 0;
85	timesPtr = 0;
86	spacePtr = 0;
87
88	// Find path to data files, i.e. xmldoc directory location.
89	// Environment variable takes precedence, else use constructor input.
90	xmlPath = "";
91	const char* PYTHIA8DATA = "PYTHIA8DATA";
92	char* envPath = getenv(PYTHIA8DATA);
93	if (envPath != 0 && *envPath != '\0') {
94	int i = 0;
95	while ((envPath+i) != '\0') xmlPath += (envPath+(i++));
96	}
97	else xmlPath = xmlDir;
98	if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
99
100	// Read in files with all flags, modes, parms and words.
101	settings.initPtr( &info);
102	string initFile = xmlPath + "Index.xml";
103	isConstructed = settings.init( initFile);
104	if (!isConstructed) {
105	info.errorMsg("Abort from Pythia::Pythia: settings unavailable");
106	return;
107	}
108
109	// Check that XML version number matches code version number.
110	double versionNumberXML = parm("Pythia:versionNumber");
111	isConstructed = (abs(versionNumberXML - VERSIONNUMBERCODE) < 0.0005);
112	if (!isConstructed) {
113	ostringstream errCode;
114	errCode << fixed << setprecision(3) << ": in code " << VERSIONNUMBERCODE
115	<< " but in XML " << versionNumberXML;
116	info.errorMsg("Abort from Pythia::Pythia: unmatched version numbers",
117	errCode.str());
118	return;
119	}
120
121	// Read in files with all particle data.
122	particleData.initPtr( &info, &settings, &rndm, couplingsPtr);
123	string dataFile = xmlPath + "ParticleData.xml";
124	isConstructed = particleData.init( dataFile);
125	if (!isConstructed) {
126	info.errorMsg("Abort from Pythia::Pythia: particle data unavailable");
127	return;
128	}
129
130	// Write the Pythia banner to output.
131	banner();
132
133	// Not initialized until at the end of the init() call.
134	isInit = false;
135	info.addCounter(0);
136
137	}
138
139	//--------------------------------------------------------------------------
140
141	// Destructor.
142
143	Pythia::~Pythia() {
144
145	// Delete the PDF's created with new.
146	if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr;
147	if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr;
148	if (useNewPdfA) delete pdfAPtr;
149	if (useNewPdfB) delete pdfBPtr;
150	if (useNewPdfPomA) delete pdfPomAPtr;
151	if (useNewPdfPomB) delete pdfPomBPtr;
152
153	// Delete the Les Houches object created with new.
154	if (useNewLHA) delete lhaUpPtr;
155
156	// Delete the MergingHooks object created with new.
157	if (hasOwnMergingHooks) delete mergingHooksPtr;
158
159	// Delete the BeamShape object created with new.
160	if (useNewBeamShape) delete beamShapePtr;
161
162	// Delete the timelike and spacelike showers created with new.
163	if (useNewTimes) delete timesPtr;
164	if (useNewSpace) delete spacePtr;
165
166	}
167
168	//--------------------------------------------------------------------------
169
170	// Read in one update for a setting or particle data from a single line.
171
172	bool Pythia::readString(string line, bool warn) {
173
174	// Check that constructor worked.
175	if (!isConstructed) return false;
176
177	// If empty line then done.
178	if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return true;
179
180	// If first character is not a letter/digit, then taken to be a comment.
181	int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
182	if (!isalnum(line[firstChar])) return true;
183
184	// Send on particle data to the ParticleData database.
185	if (isdigit(line[firstChar]))
186	return particleData.readString(line, warn);
187
188	// Everything else sent on to Settings.
189	return settings.readString(line, warn);
190
191	}
192
193	//--------------------------------------------------------------------------
194
195	// Read in updates for settings or particle data from user-defined file.
196
197	bool Pythia::readFile(string fileName, bool warn, int subrun) {
198
199	// Check that constructor worked.
200	if (!isConstructed) return false;
201
202	// Open file for reading.
203	const char* cstring = fileName.c_str();
204	ifstream is(cstring);
205	if (!is.good()) {
206	info.errorMsg("Error in Pythia::readFile: did not find file", fileName);
207	return false;
208	}
209
210	// Hand over real work to next method.
211	return readFile( is, warn, subrun);
212
213	}
214
215	//--------------------------------------------------------------------------
216
217	// Read in updates for settings or particle data
218	// from user-defined stream (or file).
219
220	bool Pythia::readFile(istream& is, bool warn, int subrun) {
221
222	// Check that constructor worked.
223	if (!isConstructed) return false;
224
225	// Read in one line at a time.
226	string line;
227	bool accepted = true;
228	int subrunNow = SUBRUNDEFAULT;
229	while ( getline(is, line) ) {
230
231	// Check whether entered new subrun.
232	int subrunLine = readSubrun( line, warn);
233	if (subrunLine >= 0) subrunNow = subrunLine;
234
235	// Process the line if in correct subrun.
236	if ( (subrunNow == subrun \|\| subrunNow == SUBRUNDEFAULT)
237	&& !readString( line, warn) ) accepted = false;
238
239	// Reached end of input file.
240	};
241	return accepted;
242
243	}
244
245	//--------------------------------------------------------------------------
246
247	// Routine to pass in pointers to PDF's. Usage optional.
248
249	bool Pythia::setPDFPtr( PDF* pdfAPtrIn, PDF* pdfBPtrIn, PDF* pdfHardAPtrIn,
250	PDF* pdfHardBPtrIn, PDF* pdfPomAPtrIn, PDF* pdfPomBPtrIn) {
251
252	// Delete any PDF's created in a previous init call.
253	if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr;
254	if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr;
255	if (useNewPdfA) delete pdfAPtr;
256	if (useNewPdfB) delete pdfBPtr;
257	if (useNewPdfPomA) delete pdfPomAPtr;
258	if (useNewPdfPomB) delete pdfPomBPtr;
259
260	// Reset pointers to be empty.
261	useNewPdfA = false;
262	useNewPdfB = false;
263	useNewPdfHard = false;
264	useNewPdfPomA = false;
265	useNewPdfPomB = false;
266	pdfAPtr = 0;
267	pdfBPtr = 0;
268	pdfHardAPtr = 0;
269	pdfHardBPtr = 0;
270	pdfPomAPtr = 0;
271	pdfPomBPtr = 0;
272
273	// Switch off external PDF's by zero as input.
274	if (pdfAPtrIn == 0 && pdfBPtrIn == 0) return true;
275
276	// The two PDF objects cannot be one and the same.
277	if (pdfAPtrIn == pdfBPtrIn) return false;
278
279	// Save pointers.
280	pdfAPtr = pdfAPtrIn;
281	pdfBPtr = pdfBPtrIn;
282
283	// By default same pointers for hard-process PDF's.
284	pdfHardAPtr = pdfAPtrIn;
285	pdfHardBPtr = pdfBPtrIn;
286
287	// Optionally allow separate pointers for hard process.
288	if (pdfHardAPtrIn != 0 && pdfHardBPtrIn != 0) {
289	if (pdfHardAPtrIn == pdfHardBPtrIn) return false;
290	pdfHardAPtr = pdfHardAPtrIn;
291	pdfHardBPtr = pdfHardBPtrIn;
292	}
293
294	// Optionally allow pointers for Pomerons in the proton.
295	if (pdfPomAPtrIn != 0 && pdfPomBPtrIn != 0) {
296	if (pdfPomAPtrIn == pdfPomBPtrIn) return false;
297	pdfPomAPtr = pdfPomAPtrIn;
298	pdfPomBPtr = pdfPomBPtrIn;
299	}
300
301	// Done.
302	return true;
303	}
304
305	//--------------------------------------------------------------------------
306
307	// Routine to initialize with the variable values of the Beams kind.
308
309	bool Pythia::init() {
310
311	// Check that constructor worked.
312	isInit = false;
313	if (!isConstructed) {
314	info.errorMsg("Abort from Pythia::init: constructur "
315	"initialization failed");
316	return false;
317	}
318
319	// Begin initialization. Find which frame type to use.
320	info.addCounter(1);
321	frameType = mode("Beams:frameType");
322
323	// Initialization with internal processes: read in and set values.
324	if (frameType < 4 ) {
325	doLHA = false;
326	boostType = frameType;
327	idA = mode("Beams:idA");
328	idB = mode("Beams:idB");
329	eCM = parm("Beams:eCM");
330	eA = parm("Beams:eA");
331	eB = parm("Beams:eB");
332	pxA = parm("Beams:pxA");
333	pyA = parm("Beams:pyA");
334	pzA = parm("Beams:pzA");
335	pxB = parm("Beams:pxB");
336	pyB = parm("Beams:pyB");
337	pzB = parm("Beams:pzB");
338
339	// Initialization with a Les Houches Event File or an LHAup object.
340	} else {
341	doLHA = true;
342	boostType = 2;
343	string lhef = word("Beams:LHEF");
344	string lhefHeader = word("Beams:LHEFheader");
345	bool readHeaders = flag("Beams:readLHEFheaders");
346	bool skipInit = flag("Beams:newLHEFsameInit");
347	int nSkipAtInit = mode("Beams:nSkipLHEFatInit");
348
349	// For file input: renew file stream or (re)new Les Houches object.
350	if (frameType == 4) {
351	const char* cstring1 = lhef.c_str();
352	if (useNewLHA && skipInit) lhaUpPtr->newEventFile(cstring1);
353	else {
354	if (useNewLHA) delete lhaUpPtr;
355	// Header is optional, so use NULL pointer to indicate no value.
356	const char* cstring2 = (lhefHeader == "void") ?
357	NULL : lhefHeader.c_str();
358	lhaUpPtr = new LHAupLHEF(cstring1, cstring2, readHeaders);
359	useNewLHA = true;
360	}
361
362	// Check that file was properly opened.
363	if (!lhaUpPtr->fileFound()) {
364	info.errorMsg("Abort from Pythia::init: "
365	"Les Houches Event File not found");
366	return false;
367	}
368
369	// For object input: at least check that not null pointer.
370	} else {
371	if (lhaUpPtr == 0) {
372	info.errorMsg("Abort from Pythia::init: "
373	"LHAup object not found");
374	return false;
375	}
376
377	// LHAup object generic abort using fileFound() routine
378	if (!lhaUpPtr->fileFound()) {
379	info.errorMsg("Abort from Pythia::init: "
380	"LHAup initialisation error");
381	return false;
382	}
383	}
384
385	// Send in pointer to info. Store or replace LHA pointer in other classes.
386	lhaUpPtr->setPtr( &info);
387	processLevel.setLHAPtr( lhaUpPtr);
388
389	// If second time around, only with new file, then simplify.
390	// Optionally skip ahead a number of events at beginning of file.
391	if (skipInit) {
392	isInit = true;
393	info.addCounter(2);
394	if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit);
395	return true;
396	}
397
398	// Set up values related to merging hooks.
399	if (frameType == 4) {
400	// Set up values related to CKKW-L merging.
401	doUserMerging = settings.flag("Merging:doUserMerging");
402	doMGMerging = settings.flag("Merging:doMGMerging");
403	doKTMerging = settings.flag("Merging:doKTMerging");
404	doPTLundMerging = settings.flag("Merging:doPTLundMerging");
405	doCutBasedMerging = settings.flag("Merging:doCutBasedMerging");
406	doMerging = doUserMerging \|\| doMGMerging \|\| doKTMerging
407	\|\| doPTLundMerging \|\| doCutBasedMerging;
408	// Set up MergingHooks object
409	if (!doUserMerging){
410	if (hasOwnMergingHooks && mergingHooksPtr) delete mergingHooksPtr;
411	mergingHooksPtr = new MergingHooks();
412	hasOwnMergingHooks = true;
413	}
414	hasMergingHooks = (mergingHooksPtr != 0);
415	// Merging hooks required for merging. If no merging hooks pointer is
416	// available, exit.
417	if (!hasMergingHooks) {
418	info.errorMsg("Abort from Pythia::init: "
419	"no merging hooks object has been provided");
420	return false;
421	} else mergingHooksPtr->setLHEInputFile( lhef);
422	// Initialise counting of Les Houches Events significantly above the
423	// merging scale.
424	info.setCounter(41,0);
425	}
426
427	// Set LHAinit information (in some external program).
428	if (settings.flag("ProcessLevel:all") && !lhaUpPtr->setInit()) {
429	info.errorMsg("Abort from Pythia::init: "
430	"Les Houches initialization failed");
431	return false;
432	}
433
434	// Extract beams from values set in an LHAinit object.
435	idA = lhaUpPtr->idBeamA();
436	idB = lhaUpPtr->idBeamB();
437	eA = lhaUpPtr->eBeamA();
438	eB = lhaUpPtr->eBeamB();
439	// Optionally skip ahead a number of events at beginning of file.
440	if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit);
441	}
442
443	// Set up values related to user hooks.
444	hasUserHooks = (userHooksPtr != 0);
445	doVetoProcess = false;
446	doVetoPartons = false;
447	if (hasUserHooks) {
448	userHooksPtr->initPtr( &info, &settings, &particleData, &rndm,
449	&beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr, &partonSystems,
450	&sigmaTot);
451	if (!userHooksPtr->initAfterBeams()) {
452	info.errorMsg("Abort from Pythia::init: could not initialise UserHooks");
453	return false;
454	}
455	doVetoProcess = userHooksPtr->canVetoProcessLevel();
456	doVetoPartons = userHooksPtr->canVetoPartonLevel();
457	}
458
459	// Back to common initialization. Reset error counters.
460	nErrEvent = 0;
461	info.errorReset();
462	info.setTooLowPTmin(false);
463	info.sigmaReset();
464
465	// Initialize data members extracted from database.
466	doProcessLevel = settings.flag("ProcessLevel:all");
467	doPartonLevel = settings.flag("PartonLevel:all");
468	doHadronLevel = settings.flag("HadronLevel:all");
469	doDiffraction = settings.flag("SoftQCD:all")
470	\|\| settings.flag("SoftQCD:singleDiffractive")
471	\|\| settings.flag("SoftQCD:doubleDiffractive")
472	\|\| settings.flag("SoftQCD:centralDiffractive");
473	doResDec = settings.flag("Standalone:allowResDec");
474	doFSRinRes = doPartonLevel && settings.flag("PartonLevel:FSR")
475	&& settings.flag("PartonLevel:FSRinResonances");
476	decayRHadrons = settings.flag("RHadrons:allowDecay");
477	doMomentumSpread = settings.flag("Beams:allowMomentumSpread");
478	doVertexSpread = settings.flag("Beams:allowVertexSpread");
479	abortIfVeto = settings.flag("Check:abortIfVeto");
480	checkEvent = settings.flag("Check:event");
481	checkHistory = settings.flag("Check:history");
482	nErrList = settings.mode("Check:nErrList");
483	epTolErr = settings.parm("Check:epTolErr");
484	epTolWarn = settings.parm("Check:epTolWarn");
485
486
487	// Initialise merging hooks.
488	if (hasMergingHooks && doMerging )
489	mergingHooksPtr->init( settings, &info, &particleData );
490
491	// Initialize the random number generator.
492	if ( settings.flag("Random:setSeed") )
493	rndm.init( settings.mode("Random:seed") );
494
495	// Check that combinations of settings are allowed; change if not.
496	checkSettings();
497
498	// Initialize couplings (needed to initialize resonances).
499	// Check if SUSY couplings need to be read in
500	if( !initSLHA()) info.errorMsg("Error in Pythia::init: "
501	"Could not read SLHA file");
502	if (couplingsPtr->isSUSY) {
503	// Initialize the SM and SUSY.
504	coupSUSY.init( settings, &rndm);
505	coupSUSY.initSUSY(&slha, &settings, &particleData);
506	couplingsPtr = (Couplings *) &coupSUSY;
507	} else {
508	// Initialize the SM couplings.
509	couplingsPtr->init( settings, &rndm);
510	}
511
512	// Reset couplingsPtr to the correct place.
513	particleData.initPtr( &info, &settings, &rndm, couplingsPtr);
514
515	// Set headers to distinguish the two event listing kinds.
516	int startColTag = settings.mode("Event:startColTag");
517	process.init("(hard process)", &particleData, startColTag);
518	event.init("(complete event)", &particleData, startColTag);
519
520	// Final setup stage of particle data, notably resonance widths.
521	particleData.initWidths( resonancePtrs);
522
523	// Set up R-hadrons particle data, where relevant.
524	rHadrons.init( &info, settings, &particleData, &rndm);
525
526	// Set up objects for timelike and spacelike showers.
527	if (timesDecPtr == 0 \|\| timesPtr == 0) {
528	TimeShower* timesNow = new TimeShower();
529	if (timesDecPtr == 0) timesDecPtr = timesNow;
530	if (timesPtr == 0) timesPtr = timesNow;
531	useNewTimes = true;
532	}
533	if (spacePtr == 0) {
534	spacePtr = new SpaceShower();
535	useNewSpace = true;
536	}
537
538	// Initialize showers, especially for simple showers in decays.
539	timesPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr,
540	&partonSystems, userHooksPtr);
541	timesDecPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr,
542	&partonSystems, userHooksPtr);
543	spacePtr->initPtr( &info, &settings, &particleData, &rndm,
544	&partonSystems, userHooksPtr);
545	timesDecPtr->init( 0, 0);
546
547	// Set up values related to beam shape.
548	if (beamShapePtr == 0) {
549	beamShapePtr = new BeamShape();
550	useNewBeamShape = true;
551	}
552	beamShapePtr->init( settings, &rndm);
553
554	// Check that beams and beam combination can be handled.
555	if (!checkBeams()) {
556	info.errorMsg("Abort from Pythia::init: "
557	"checkBeams initialization failed");
558	return false;
559	}
560
561	// Do not set up beam kinematics when no process level.
562	if (!doProcessLevel) boostType = 1;
563	else {
564
565	// Set up beam kinematics.
566	if (!initKinematics()) {
567	info.errorMsg("Abort from Pythia::init: "
568	"kinematics initialization failed");
569	return false;
570	}
571
572	// Set up pointers to PDFs.
573	if (!initPDFs()) {
574	info.errorMsg("Abort from Pythia::init: PDF initialization failed");
575	return false;
576	}
577
578	// Set up the two beams and the common remnant system.
579	StringFlav* flavSelPtr = hadronLevel.getStringFlavPtr();
580	beamA.init( idA, pzAcm, eA, mA, &info, settings, &particleData, &rndm,
581	pdfAPtr, pdfHardAPtr, isUnresolvedA, flavSelPtr);
582	beamB.init( idB, pzBcm, eB, mB, &info, settings, &particleData, &rndm,
583	pdfBPtr, pdfHardBPtr, isUnresolvedB, flavSelPtr);
584
585	// Optionally set up new alternative beams for these Pomerons.
586	if ( doDiffraction) {
587	beamPomA.init( 990, 0.5 * eCM, 0.5 * eCM, 0., &info, settings,
588	&particleData, &rndm, pdfPomAPtr, pdfPomAPtr, false, flavSelPtr);
589	beamPomB.init( 990, -0.5 * eCM, 0.5 * eCM, 0., &info, settings,
590	&particleData, &rndm, pdfPomBPtr, pdfPomBPtr, false, flavSelPtr);
591	}
592	}
593
594	// Send info/pointers to process level for initialization.
595	if ( doProcessLevel && !processLevel.init( &info, settings, &particleData,
596	&rndm, &beamA, &beamB, couplingsPtr, &sigmaTot, doLHA, &slha, userHooksPtr,
597	sigmaPtrs) ) {
598	info.errorMsg("Abort from Pythia::init: "
599	"processLevel initialization failed");
600	return false;
601	}
602
603	// Optionally only initialize resonance decays.
604	if ( !doProcessLevel && doResDec) processLevel.initDecays( &info,
605	&particleData, &rndm, lhaUpPtr);
606
607	// Send info/pointers to parton level for initialization.
608	if ( doPartonLevel && doProcessLevel && !partonLevel.init( &info, settings,
609	&particleData, &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr,
610	&partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons,
611	userHooksPtr, mergingHooksPtr, false) ) {
612	info.errorMsg("Abort from Pythia::init: "
613	"partonLevel initialization failed" );
614	return false;
615	}
616
617	// Optionally only initialize final-state showers in resonance decays.
618	if ( (!doProcessLevel \|\| !doPartonLevel) && doFSRinRes) partonLevel.init(
619	&info, settings, &particleData, &rndm, 0, 0, 0, 0, couplingsPtr,
620	&partonSystems, 0, timesDecPtr, 0, 0, &rHadrons, 0, 0, false);
621
622	// Send info/pointers to parton level for trial shower initialization.
623	if ( doMerging
624	&& !trialPartonLevel.init( &info, settings, &particleData,
625	&rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr,
626	&partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons,
627	NULL, mergingHooksPtr, true) ) {
628	info.errorMsg("Abort from Pythia::init: "
629	"trialPartonLevel initialization failed");
630	return false;
631	}
632
633	// Send info/pointers to hadron level for initialization.
634	// Note: forceHadronLevel() can come, so we must always initialize.
635	if ( !hadronLevel.init( &info, settings, &particleData, &rndm,
636	couplingsPtr, timesDecPtr, &rHadrons, decayHandlePtr,
637	handledParticles) ) {
638	info.errorMsg("Abort from Pythia::init: "
639	"hadronLevel initialization failed");
640	return false;
641	}
642
643	// Optionally check particle data table for inconsistencies.
644	if ( settings.flag("Check:particleData") )
645	particleData.checkTable( settings.mode("Check:levelParticleData") );
646
647	// Optionally show settings and particle data, changed or all.
648	bool showCS = settings.flag("Init:showChangedSettings");
649	bool showAS = settings.flag("Init:showAllSettings");
650	bool showCPD = settings.flag("Init:showChangedParticleData");
651	bool showCRD = settings.flag("Init:showChangedResonanceData");
652	bool showAPD = settings.flag("Init:showAllParticleData");
653	int show1PD = settings.mode("Init:showOneParticleData");
654	bool showPro = settings.flag("Init:showProcesses");
655	if (showCS) settings.listChanged();
656	if (showAS) settings.listAll();
657	if (show1PD > 0) particleData.list(show1PD);
658	if (showCPD) particleData.listChanged(showCRD);
659	if (showAPD) particleData.listAll();
660
661	// Listing options for the next() routine.
662	nCount = settings.mode("Next:numberCount");
663	nShowLHA = settings.mode("Next:numberShowLHA");
664	nShowInfo = settings.mode("Next:numberShowInfo");
665	nShowProc = settings.mode("Next:numberShowProcess");
666	nShowEvt = settings.mode("Next:numberShowEvent");
667	showSaV = settings.flag("Next:showScaleAndVertex");
668	showMaD = settings.flag("Next:showMothersAndDaughters");
669
670	// Succeeded.
671	isInit = true;
672	info.addCounter(2);
673	if (useNewLHA && showPro) lhaUpPtr->listInit();
674	return true;
675
676	}
677
678	//--------------------------------------------------------------------------
679
680	// Routine to initialize with CM energy.
681
682	bool Pythia::init( int idAin, int idBin, double eCMin) {
683
684	// Set arguments in Settings database.
685	settings.mode("Beams:idA", idAin);
686	settings.mode("Beams:idB", idBin);
687	settings.mode("Beams:frameType", 1);
688	settings.parm("Beams:eCM", eCMin);
689
690	// Send on to common initialization.
691	return init();
692
693	}
694
695	//--------------------------------------------------------------------------
696
697	// Routine to initialize with two collinear beams, energies specified.
698
699	bool Pythia::init( int idAin, int idBin, double eAin, double eBin) {
700
701	// Set arguments in Settings database.
702	settings.mode("Beams:idA", idAin);
703	settings.mode("Beams:idB", idBin);
704	settings.mode("Beams:frameType", 2);
705	settings.parm("Beams:eA", eAin);
706	settings.parm("Beams:eB", eBin);
707
708	// Send on to common initialization.
709	return init();
710
711	}
712
713	//--------------------------------------------------------------------------
714
715	// Routine to initialize with two beams specified by three-momenta.
716
717	bool Pythia::init( int idAin, int idBin, double pxAin, double pyAin,
718	double pzAin, double pxBin, double pyBin, double pzBin) {
719
720	// Set arguments in Settings database.
721	settings.mode("Beams:idA", idAin);
722	settings.mode("Beams:idB", idBin);
723	settings.mode("Beams:frameType", 3);
724	settings.parm("Beams:pxA", pxAin);
725	settings.parm("Beams:pyA", pyAin);
726	settings.parm("Beams:pzA", pzAin);
727	settings.parm("Beams:pxB", pxBin);
728	settings.parm("Beams:pyB", pyBin);
729	settings.parm("Beams:pzB", pzBin);
730
731	// Send on to common initialization.
732	return init();
733
734	}
735
736	//--------------------------------------------------------------------------
737
738	// Routine to initialize when all info is given in a Les Houches Event File.
739
740	bool Pythia::init( string LesHouchesEventFile, bool skipInit) {
741
742	// Set arguments in Settings database.
743	settings.mode("Beams:frameType", 4);
744	settings.word("Beams:LHEF", LesHouchesEventFile);
745	settings.flag("Beams:newLHEFsameInit", skipInit);
746
747	// Send on to common initialization.
748	return init();
749
750	}
751
752	//--------------------------------------------------------------------------
753
754	// Routine to initialize when beam info is given in an LHAup object.
755
756	bool Pythia::init( LHAup* lhaUpPtrIn) {
757
758	// Store LHAup object and set arguments in Settings database.
759	setLHAupPtr( lhaUpPtrIn);
760	settings.mode("Beams:frameType", 5);
761
762	// Send on to common initialization.
763	return init();
764
765	}
766
767	//--------------------------------------------------------------------------
768
769	// Check that combinations of settings are allowed; change if not.
770
771	void Pythia::checkSettings() {
772
773	// Double rescattering not allowed if ISR or FSR.
774	if ((settings.flag("PartonLevel:ISR") \|\| settings.flag("PartonLevel:FSR"))
775	&& settings.flag("MultipartonInteractions:allowDoubleRescatter")) {
776	info.errorMsg("Warning in Pythia::checkSettings: "
777	"double rescattering switched off since showering is on");
778	settings.flag("MultipartonInteractions:allowDoubleRescatter", false);
779	}
780
781	}
782
783	//--------------------------------------------------------------------------
784
785	// Check that beams and beam combination can be handled. Set up unresolved.
786
787	bool Pythia::checkBeams() {
788
789	// Absolute flavours. If not to do process level then no check needed.
790	int idAabs = abs(idA);
791	int idBabs = abs(idB);
792	if (!doProcessLevel) return true;
793
794	// Neutrino beams always unresolved, charged lepton ones conditionally.
795	bool isLeptonA = (idAabs > 10 && idAabs < 17);
796	bool isLeptonB = (idBabs > 10 && idBabs < 17);
797	bool isUnresLep = !settings.flag("PDF:lepton");
798	isUnresolvedA = isLeptonA && (idAabs%2 == 0 \|\| isUnresLep);
799	isUnresolvedB = isLeptonB && (idBabs%2 == 0 \|\| isUnresLep);
800
801	// Lepton-lepton collisions OK (including neutrinos) if both (un)resolved.
802	if (isLeptonA && isLeptonB && isUnresolvedA == isUnresolvedB) return true;
803
804	// MBR model only implemented for pp/ppbar/pbarp collisions.
805	int PomFlux = settings.mode("Diffraction:PomFlux");
806	if (PomFlux == 5) {
807	bool ispp = (idAabs == 2212 && idBabs == 2212);
808	bool ispbarpbar = (idA == -2212 && idB == -2212);
809	if (ispp && !ispbarpbar) return true;
810	info.errorMsg("Error in Pythia::init: cannot handle this beam combination"
811	" with PomFlux == 5");
812	return false;
813	}
814
815	// Hadron-hadron collisions OK, with Pomeron counted as hadron.
816	bool isHadronA = (idAabs == 2212) \|\| (idA == 111) \|\| (idAabs == 211)
817	\|\| (idA == 990);
818	bool isHadronB = (idBabs == 2212) \|\| (idA == 111)\|\| (idBabs == 211)
819	\|\| (idB == 990);
820	if (isHadronA && isHadronB) return true;
821
822	// If no case above then failed.
823	info.errorMsg("Error in Pythia::init: cannot handle this beam combination");
824	return false;
825
826	}
827
828	//--------------------------------------------------------------------------
829
830	// Calculate kinematics at initialization. Store beam four-momenta.
831
832	bool Pythia::initKinematics() {
833
834	// Find masses. Initial guess that we are in CM frame.
835	mA = particleData.m0(idA);
836	mB = particleData.m0(idB);
837	betaZ = 0.;
838	gammaZ = 1.;
839
840	// Collinear beams not in CM frame: find CM energy.
841	if (boostType == 2) {
842	eA = max(eA, mA);
843	eB = max(eB, mB);
844	pzA = sqrt(eAeA - mAmA);
845	pzB = -sqrt(eBeB - mBmB);
846	pAinit = Vec4( 0., 0., pzA, eA);
847	pBinit = Vec4( 0., 0., pzB, eB);
848	eCM = sqrt( pow2(eA + eB) - pow2(pzA + pzB) );
849
850	// Find boost to rest frame.
851	betaZ = (pzA + pzB) / (eA + eB);
852	gammaZ = (eA + eB) / eCM;
853	if (abs(betaZ) < 1e-10) boostType = 1;
854	}
855
856	// Completely general beam directions: find CM energy.
857	else if (boostType == 3) {
858	eA = sqrt( pxApxA + pyApyA + pzApzA + mAmA);
859	eB = sqrt( pxBpxB + pyBpyB + pzBpzB + mBmB);
860	pAinit = Vec4( pxA, pyA, pzA, eA);
861	pBinit = Vec4( pxB, pyB, pzB, eB);
862	eCM = (pAinit + pBinit).mCalc();
863
864	// Find boost+rotation needed to move from/to CM frame.
865	MfromCM.reset();
866	MfromCM.fromCMframe( pAinit, pBinit);
867	MtoCM = MfromCM;
868	MtoCM.invert();
869	}
870
871	// Fail if CM energy below beam masses.
872	if (eCM < mA + mB) {
873	info.errorMsg("Error in Pythia::initKinematics: too low energy");
874	return false;
875	}
876
877	// Set up CM-frame kinematics with beams along +-z axis.
878	pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
879	* (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
880	pzBcm = -pzAcm;
881	eA = sqrt(mAmA + pzAcmpzAcm);
882	eB = sqrt(mBmB + pzBcmpzBcm);
883
884	// If in CM frame then store beam four-vectors (else already done above).
885	if (boostType != 2 && boostType != 3) {
886	pAinit = Vec4( 0., 0., pzAcm, eA);
887	pBinit = Vec4( 0., 0., pzBcm, eB);
888	}
889
890	// Store main info for access in process generation.
891	info.setBeamA( idA, pzAcm, eA, mA);
892	info.setBeamB( idB, pzBcm, eB, mB);
893	info.setECM( eCM);
894
895	// Must allow for generic boost+rotation when beam momentum spread.
896	if (doMomentumSpread) boostType = 3;
897
898	// Done.
899	return true;
900
901	}
902
903	//--------------------------------------------------------------------------
904
905	// Set up pointers to PDFs.
906
907	bool Pythia::initPDFs() {
908
909	// Delete any PDF's created in a previous init call.
910	if (useNewPdfHard) {
911	if (pdfHardAPtr != pdfAPtr) {
912	delete pdfHardAPtr;
913	pdfHardAPtr = 0;
914	}
915	if (pdfHardBPtr != pdfBPtr) {
916	delete pdfHardBPtr;
917	pdfHardBPtr = 0;
918	}
919	useNewPdfHard = false;
920	}
921	if (useNewPdfA) {
922	delete pdfAPtr;
923	useNewPdfA = false;
924	pdfAPtr = 0;
925	}
926	if (useNewPdfB) {
927	delete pdfBPtr;
928	useNewPdfB = false;
929	pdfBPtr = 0;
930	}
931	if (useNewPdfPomA) {
932	delete pdfPomAPtr;
933	useNewPdfPomA = false;
934	pdfPomAPtr = 0;
935	}
936	if (useNewPdfPomB) {
937	delete pdfPomBPtr;
938	useNewPdfPomB = false;
939	pdfPomBPtr = 0;
940	}
941
942	// Set up the PDF's, if not already done.
943	if (pdfAPtr == 0) {
944	pdfAPtr = getPDFPtr(idA);
945	if (pdfAPtr == 0 \|\| !pdfAPtr->isSetup()) {
946	info.errorMsg("Error in Pythia::init: "
947	"could not set up PDF for beam A");
948	return false;
949	}
950	pdfHardAPtr = pdfAPtr;
951	useNewPdfA = true;
952	}
953	if (pdfBPtr == 0) {
954	pdfBPtr = getPDFPtr(idB);
955	if (pdfBPtr == 0 \|\| !pdfBPtr->isSetup()) {
956	info.errorMsg("Error in Pythia::init: "
957	"could not set up PDF for beam B");
958	return false;
959	}
960	pdfHardBPtr = pdfBPtr;
961	useNewPdfB = true;
962	}
963
964	// Optionally set up separate PDF's for hard process.
965	if (settings.flag("PDF:useHard") && useNewPdfA && useNewPdfB) {
966	pdfHardAPtr = getPDFPtr(idA, 2);
967	if (!pdfHardAPtr->isSetup()) return false;
968	pdfHardBPtr = getPDFPtr(idB, 2);
969	if (!pdfHardBPtr->isSetup()) return false;
970	useNewPdfHard = true;
971	}
972
973	// Optionally set up Pomeron PDF's for diffractive physics.
974	if ( doDiffraction) {
975	if (pdfPomAPtr == 0) {
976	pdfPomAPtr = getPDFPtr(990);
977	useNewPdfPomA = true;
978	}
979	if (pdfPomBPtr == 0) {
980	pdfPomBPtr = getPDFPtr(990);
981	useNewPdfPomB = true;
982	}
983	}
984
985	// Done.
986	return true;
987
988	}
989
990	//--------------------------------------------------------------------------
991
992	// Main routine to generate the next event, using internal machinery.
993
994	bool Pythia::next() {
995
996	// Regularly print how many events have been generated.
997	int nPrevious = info.getCounter(3);
998	if (nCount > 0 && nPrevious > 0 && nPrevious%nCount == 0)
999	cout << "\n Pythia::next(): " << nPrevious
1000	<< " events have been generated " << endl;
1001
1002	// Set/reset info counters specific to each event.
1003	info.addCounter(3);
1004	for (int i = 10; i < 13; ++i) info.setCounter(i);
1005
1006	// Simpler option when no hard process, i.e. mainly hadron level.
1007	if (!doProcessLevel) {
1008
1009	// Optionally fetch in resonance decays from LHA interface.
1010	if (doLHA && !processLevel.nextLHAdec( event)) {
1011	if (info.atEndOfFile()) info.errorMsg("Abort from Pythia::next:"
1012	" reached end of Les Houches Events File");
1013	return false;
1014	}
1015
1016	// Reset info array (while event record contains data).
1017	info.clear();
1018
1019	// Set correct energy for system.
1020	Vec4 pSum = 0.;
1021	for (int i = 1; i < event.size(); ++i)
1022	if (event[i].isFinal()) pSum += event[i].p();
1023	event[0].p( pSum );
1024	event[0].m( pSum.mCalc() );
1025
1026	// Generate hadronization and decays.
1027	bool status = forceHadronLevel();
1028	if (status) info.addCounter(4);
1029	if (status && nPrevious < nShowEvt) event.list(showSaV, showMaD);
1030	return status;
1031	}
1032
1033	// Reset arrays.
1034	info.clear();
1035	process.clear();
1036	event.clear();
1037	partonSystems.clear();
1038	beamA.clear();
1039	beamB.clear();
1040	beamPomA.clear();
1041	beamPomB.clear();
1042
1043	// Can only generate event if initialization worked.
1044	if (!isInit) {
1045	info.errorMsg("Abort from Pythia::next: "
1046	"not properly initialized so cannot generate events");
1047	return false;
1048	}
1049
1050	// Pick beam momentum spread and beam vertex.
1051	if (doMomentumSpread \|\| doVertexSpread) beamShapePtr->pick();
1052
1053	// Recalculate kinematics when beam momentum spread.
1054	if (doMomentumSpread) nextKinematics();
1055
1056	// Outer loop over hard processes; only relevant for user-set vetoes.
1057	for ( ; ; ) {
1058	info.addCounter(10);
1059	bool hasVetoed = false;
1060
1061	// Provide the hard process that starts it off. Only one try.
1062	info.clear();
1063	process.clear();
1064
1065	if ( !processLevel.next( process) ) {
1066	if (doLHA && info.atEndOfFile()) info.errorMsg("Abort from "
1067	"Pythia::next: reached end of Les Houches Events File");
1068	else info.errorMsg("Abort from Pythia::next: "
1069	"processLevel failed; giving up");
1070	return false;
1071	}
1072
1073	info.addCounter(11);
1074
1075	// Possibility for a user veto of the process-level event.
1076	if (doVetoProcess) {
1077	hasVetoed = userHooksPtr->doVetoProcessLevel( process);
1078	if (hasVetoed) {
1079	if (abortIfVeto) return false;
1080	continue;
1081	}
1082	}
1083
1084	// Possibility to perform CKKW-L merging on this event.
1085	if(doMerging) {
1086	hasVetoed = mergeProcess();
1087	if (hasVetoed) {
1088	if (abortIfVeto) return false;
1089	continue;
1090	}
1091	}
1092
1093	// Possibility to stop the generation at this stage.
1094	if (!doPartonLevel) {
1095	boostAndVertex( true, true);
1096	processLevel.accumulate();
1097	info.addCounter(4);
1098	if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1099	if (nPrevious < nShowInfo) info.list();
1100	if (nPrevious < nShowProc) process.list(showSaV, showMaD);
1101	return true;
1102	}
1103
1104	// Save spare copy of process record in case of problems.
1105	Event processSave = process;
1106	int sizeMPI = info.sizeMPIarrays();
1107	info.addCounter(12);
1108	for (int i = 14; i < 19; ++i) info.setCounter(i);
1109
1110	// Allow up to ten tries for parton- and hadron-level processing.
1111	bool physical = true;
1112	for (int iTry = 0; iTry < NTRY; ++ iTry) {
1113	info.addCounter(14);
1114	physical = true;
1115
1116	// Restore original process record if problems.
1117	if (iTry > 0) process = processSave;
1118	if (iTry > 0) info.resizeMPIarrays( sizeMPI);
1119
1120	// Reset event record and (extracted partons from) beam remnants.
1121	event.clear();
1122	beamA.clear();
1123	beamB.clear();
1124	beamPomA.clear();
1125	beamPomB.clear();
1126	partonSystems.clear();
1127
1128	// Parton-level evolution: ISR, FSR, MPI.
1129	if ( !partonLevel.next( process, event) ) {
1130	// Skip to next hard process for failure owing to deliberate veto.
1131	hasVetoed = partonLevel.hasVetoed();
1132	if (hasVetoed) {
1133	if (abortIfVeto) return false;
1134	break;
1135	}
1136	// Else make a new try for other failures.
1137	info.errorMsg("Error in Pythia::next: "
1138	"partonLevel failed; try again");
1139	physical = false;
1140	continue;
1141	}
1142	info.addCounter(15);
1143
1144	// Possibility for a user veto of the parton-level event.
1145	if (doVetoPartons) {
1146	hasVetoed = userHooksPtr->doVetoPartonLevel( event);
1147	if (hasVetoed) {
1148	if (abortIfVeto) return false;
1149	break;
1150	}
1151	}
1152
1153	// Boost to lab frame (before decays, for vertices).
1154	boostAndVertex( true, true);
1155
1156	// Possibility to stop the generation at this stage.
1157	if (!doHadronLevel) {
1158	processLevel.accumulate();
1159	partonLevel.accumulate();
1160
1161	// Optionally check final event for problems.
1162	if (checkEvent && !check()) {
1163	info.errorMsg("Abort from Pythia::next: "
1164	"check of event revealed problems");
1165	return false;
1166	}
1167	info.addCounter(4);
1168	if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1169	if (nPrevious < nShowInfo) info.list();
1170	if (nPrevious < nShowProc) process.list(showSaV, showMaD);
1171	if (nPrevious < nShowEvt) event.list(showSaV, showMaD);
1172	return true;
1173	}
1174
1175	// Hadron-level: hadronization, decays.
1176	info.addCounter(16);
1177	if ( !hadronLevel.next( event) ) {
1178	info.errorMsg("Error in Pythia::next: "
1179	"hadronLevel failed; try again");
1180	physical = false;
1181	continue;
1182	}
1183
1184	// If R-hadrons have been formed, then (optionally) let them decay.
1185	if (decayRHadrons && rHadrons.exist() && !doRHadronDecays()) {
1186	info.errorMsg("Error in Pythia::next: "
1187	"decayRHadrons failed; try again");
1188	physical = false;
1189	continue;
1190	}
1191	info.addCounter(17);
1192
1193	// Optionally check final event for problems.
1194	if (checkEvent && !check()) {
1195	info.errorMsg("Error in Pythia::next: "
1196	"check of event revealed problems");
1197	physical = false;
1198	continue;
1199	}
1200
1201	// Stop parton- and hadron-level looping if you got this far.
1202	info.addCounter(18);
1203	break;
1204	}
1205
1206	// If event vetoed then to make a new try.
1207	if (hasVetoed) {
1208	if (abortIfVeto) return false;
1209	continue;
1210	}
1211
1212	// If event failed any other way (after ten tries) then give up.
1213	if (!physical) {
1214	info.errorMsg("Abort from Pythia::next: "
1215	"parton+hadronLevel failed; giving up");
1216	return false;
1217	}
1218
1219	// Process- and parton-level statistics. Event scale.
1220	processLevel.accumulate();
1221	partonLevel.accumulate();
1222	event.scale( process.scale() );
1223
1224	// End of outer loop over hard processes. Done with normal option.
1225	info.addCounter(13);
1226	break;
1227	}
1228
1229	// List events.
1230	if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1231	if (nPrevious < nShowInfo) info.list();
1232	if (nPrevious < nShowProc) process.list(showSaV,showMaD);
1233	if (nPrevious < nShowEvt) event.list(showSaV, showMaD);
1234
1235	// Done.
1236	info.addCounter(4);
1237	return true;
1238
1239	}
1240
1241	//--------------------------------------------------------------------------
1242
1243	// Generate only the hadronization/decay stage, using internal machinery.
1244	// The "event" instance should already contain a parton-level configuration.
1245
1246	bool Pythia::forceHadronLevel(bool findJunctions) {
1247
1248	// Can only generate event if initialization worked.
1249	if (!isInit) {
1250	info.errorMsg("Abort from Pythia::forceHadronLevel: "
1251	"not properly initialized so cannot generate events");
1252	return false;
1253	}
1254
1255	// Check whether any junctions in system. (Normally done in ProcessLevel.)
1256	if (findJunctions) {
1257	event.clearJunctions();
1258	processLevel.findJunctions( event);
1259	}
1260
1261	// Save spare copy of event in case of failure.
1262	Event spareEvent = event;
1263
1264	// Allow up to ten tries for hadron-level processing.
1265	bool physical = true;
1266	for (int iTry = 0; iTry < NTRY; ++ iTry) {
1267	physical = true;
1268
1269	// Check whether any resonances need to be handled at process level.
1270	if (doResDec) {
1271	process = event;
1272	processLevel.nextDecays( process);
1273	bool hasDecays = false;
1274	for (int i = 1; i < process.size(); ++i)
1275	if (process[i].status() < 0) hasDecays = true;
1276
1277	// Allow for showers if decays happened at process level.
1278	if (hasDecays) {
1279	if (doFSRinRes) {
1280	partonLevel.setupShowerSys( process, event);
1281	partonLevel.resonanceShowers( process, event, false);
1282	} else event = process;
1283	}
1284	}
1285
1286	// Hadron-level: hadronization, decays.
1287	if (hadronLevel.next( event)) break;
1288
1289	// If failure then warn, restore original configuration and try again.
1290	info.errorMsg("Error in Pythia::forceHadronLevel: "
1291	"hadronLevel failed; try again");
1292	physical = false;
1293	event = spareEvent;
1294	}
1295
1296	// Done for simpler option.
1297	if (!physical) {
1298	info.errorMsg("Abort from Pythia::forceHadronLevel: "
1299	"hadronLevel failed; giving up");
1300	return false;
1301	}
1302
1303	// Optionally check final event for problems.
1304	if (checkEvent && !check()) {
1305	info.errorMsg("Abort from Pythia::forceHadronLevel: "
1306	"check of event revealed problems");
1307	return false;
1308	}
1309
1310	// Done.
1311	return true;
1312
1313	}
1314
1315	//--------------------------------------------------------------------------
1316
1317	// Recalculate kinematics for each event when beam momentum has a spread.
1318
1319	void Pythia::nextKinematics() {
1320
1321	// Read out momentum shift to give current beam momenta.
1322	pAnow = pAinit + beamShapePtr->deltaPA();
1323	pAnow.e( sqrt(pAnow.pAbs2() + mA * mA) );
1324	pBnow = pBinit + beamShapePtr->deltaPB();
1325	pBnow.e( sqrt(pBnow.pAbs2() + mB * mB) );
1326
1327	// Construct CM frame kinematics.
1328	eCM = (pAnow + pBnow).mCalc();
1329	pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
1330	* (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
1331	pzBcm = -pzAcm;
1332	eA = sqrt(mAmA + pzAcmpzAcm);
1333	eB = sqrt(mBmB + pzBcmpzBcm);
1334
1335	// Set relevant info for other classes to use.
1336	info.setBeamA( idA, pzAcm, eA, mA);
1337	info.setBeamB( idB, pzBcm, eB, mB);
1338	info.setECM( eCM);
1339	beamA.newPzE( pzAcm, eA);
1340	beamB.newPzE( pzBcm, eB);
1341
1342	// Set boost/rotation matrices from/to CM frame.
1343	MfromCM.reset();
1344	MfromCM.fromCMframe( pAnow, pBnow);
1345	MtoCM = MfromCM;
1346	MtoCM.invert();
1347
1348	}
1349
1350	//--------------------------------------------------------------------------
1351
1352	// Boost from CM frame to lab frame, or inverse. Set production vertex.
1353
1354	void Pythia::boostAndVertex( bool toLab, bool setVertex) {
1355
1356	// Boost process from CM frame to lab frame.
1357	if (toLab) {
1358	if (boostType == 2) process.bst(0., 0., betaZ, gammaZ);
1359	else if (boostType == 3) process.rotbst(MfromCM);
1360
1361	// Boost nonempty event from CM frame to lab frame.
1362	if (event.size() > 0) {
1363	if (boostType == 2) event.bst(0., 0., betaZ, gammaZ);
1364	else if (boostType == 3) event.rotbst(MfromCM);
1365	}
1366
1367	// Boost process from lab frame to CM frame.
1368	} else {
1369	if (boostType == 2) process.bst(0., 0., -betaZ, gammaZ);
1370	else if (boostType == 3) process.rotbst(MtoCM);
1371
1372	// Boost nonempty event from lab frame to CM frame.
1373	if (event.size() > 0) {
1374	if (boostType == 2) event.bst(0., 0., -betaZ, gammaZ);
1375	else if (boostType == 3) event.rotbst(MtoCM);
1376	}
1377	}
1378
1379	// Set production vertex; assumes particles are in lab frame and at origin.
1380	if (setVertex && doVertexSpread) {
1381	Vec4 vertex = beamShapePtr->vertex();
1382	for (int i = 0; i < process.size(); ++i) process[i].vProd( vertex);
1383	for (int i = 0; i < event.size(); ++i) event[i].vProd( vertex);
1384	}
1385
1386	}
1387
1388	//--------------------------------------------------------------------------
1389
1390	// Perform R-hadron decays, either as part of normal evolution or forced.
1391
1392	bool Pythia::doRHadronDecays( ) {
1393
1394	// Check if R-hadrons exist to be processed.
1395	if ( !rHadrons.exist() ) return true;
1396
1397	// Do the R-hadron decay itself.
1398	if ( !rHadrons.decay( event) ) return false;
1399
1400	// Perform showers in resonance decay chains.
1401	if ( !partonLevel.resonanceShowers( process, event, false) ) return false;
1402
1403	// Subsequent hadronization and decays.
1404	if ( !hadronLevel.next( event) ) return false;
1405
1406	// Done.
1407	return true;
1408
1409	}
1410
1411	//--------------------------------------------------------------------------
1412
1413	// Print statistics on event generation.
1414
1415	void Pythia::stat() {
1416
1417	// Read out settings for what to include.
1418	bool showPrL = settings.flag("Stat:showProcessLevel");
1419	bool showPaL = settings.flag("Stat:showPartonLevel");
1420	bool showErr = settings.flag("Stat:showErrors");
1421	bool reset = settings.flag("Stat:reset");
1422
1423	// Statistics on cross section and number of events.
1424	if (doProcessLevel) {
1425	if (showPrL) processLevel.statistics(false);
1426	if (reset) processLevel.resetStatistics();
1427	}
1428
1429	// Statistics from other classes, currently multiparton interactions.
1430	if (showPaL) partonLevel.statistics(false);
1431	if (reset) partonLevel.resetStatistics();
1432
1433	// Warning if significant part of events considerably above the
1434	// merging scale value.
1435	if (hasMergingHooks && doMerging && mergingHooksPtr->stats() ) {
1436	string message="Warning in MergingHooks::stats: Most LH";
1437	message+=" Events significantly above Merging:TMS cut. Please check.";
1438	info.errorMsg(message);
1439	}
1440
1441	// Summary of which and how many warnings/errors encountered.
1442	if (showErr) info.errorStatistics();
1443	if (reset) info.errorReset();
1444
1445	}
1446
1447	//--------------------------------------------------------------------------
1448
1449	// Print statistics on event generation - deprecated version.
1450
1451	void Pythia::statistics(bool all, bool reset) {
1452
1453	// Statistics on cross section and number of events.
1454	if (doProcessLevel) processLevel.statistics(reset);
1455
1456	// Statistics from other classes, e.g. multiparton interactions.
1457	if (all) partonLevel.statistics(reset);
1458
1459	// Warning if significant part of events considerably above the
1460	// merging scale value.
1461	if (hasMergingHooks && doMerging && mergingHooksPtr->stats() ) {
1462	string message="Warning in MergingHooks::stats: Most LH";
1463	message+=" Events significantly above Merging:TMS cut. Please check.";
1464	info.errorMsg(message);
1465	}
1466
1467	// Summary of which and how many warnings/errors encountered.
1468	info.errorStatistics();
1469	if (reset) info.errorReset();
1470
1471	}
1472
1473	//--------------------------------------------------------------------------
1474
1475	// Write the Pythia banner, with symbol and version information.
1476
1477	void Pythia::banner(ostream& os) {
1478
1479	// Read in version number and last date of change.
1480	double versionNumber = settings.parm("Pythia:versionNumber");
1481	int versionDate = settings.mode("Pythia:versionDate");
1482	string month[12] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun",
1483	"Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
1484
1485	// Get date and time.
1486	time_t t = time(0);
1487	char dateNow[12];
1488	strftime(dateNow,12,"%d %b %Y",localtime(&t));
1489	char timeNow[9];
1490	strftime(timeNow,9,"%H:%M:%S",localtime(&t));
1491
1492	os << "\n"
1493	<< " *-------------------------------------------"
1494	<< "-----------------------------------------* \n"
1495	<< " \| "
1496	<< " \| \n"
1497	<< " \| *----------------------------------------"
1498	<< "--------------------------------------* \| \n"
1499	<< " \| \| "
1500	<< " \| \| \n"
1501	<< " \| \| "
1502	<< " \| \| \n"
1503	<< " \| \| PPP Y Y TTTTT H H III A "
1504	<< " Welcome to the Lund Monte Carlo! \| \| \n"
1505	<< " \| \| P P Y Y T H H I A A "
1506	<< " This is PYTHIA version " << fixed << setprecision(3)
1507	<< setw(5) << versionNumber << " \| \| \n"
1508	<< " \| \| PPP Y T HHHHH I AAAAA"
1509	<< " Last date of change: " << setw(2) << versionDate%100
1510	<< " " << month[ (versionDate/100)%100 - 1 ]
1511	<< " " << setw(4) << versionDate/10000 << " \| \| \n"
1512	<< " \| \| P Y T H H I A A"
1513	<< " \| \| \n"
1514	<< " \| \| P Y T H H III A A"
1515	<< " Now is " << dateNow << " at " << timeNow << " \| \| \n"
1516	<< " \| \| "
1517	<< " \| \| \n"
1518	<< " \| \| Torbjorn Sjostrand; Department of As"
1519	<< "tronomy and Theoretical Physics, \| \| \n"
1520	<< " \| \| Lund University, Solvegatan 14A, S"
1521	<< "E-223 62 Lund, Sweden; \| \| \n"
1522	<< " \| \| phone: + 46 - 46 - 222 48 16; e-ma"
1523	<< "il: torbjorn@thep.lu.se \| \| \n"
1524	<< " \| \| Stefan Ask; Department of Physics, U"
1525	<< "niversity of Cambridge, \| \| \n"
1526	<< " \| \| Cavendish Laboratory, JJ Thomson A"
1527	<< "ve., Cambridge CB3 0HE, UK; \| \| \n"
1528	<< " \| \| phone: + 41 - 22 - 767 6707; e-mai"
1529	<< "l: Stefan.Ask@cern.ch \| \| \n"
1530	<< " \| \| Richard Corke; Niels Bohr Institute,"
1531	<< " University of Copenhagen, \| \| \n"
1532	<< " \| \| Blegdamsvej 17, 2100 Copenhagen, D"
1533	<< "enmark; \| \| \n"
1534	<< " \| \| e-mail: richard.corke@thep.lu.se "
1535	<< " \| \| \n"
1536	<< " \| \| Stephen Mrenna; Computing Division, "
1537	<< "Simulations Group, \| \| \n"
1538	<< " \| \| Fermi National Accelerator Laborat"
1539	<< "ory, MS 234, Batavia, IL 60510, USA; \| \| \n"
1540	<< " \| \| phone: + 1 - 630 - 840 - 2556; e-m"
1541	<< "ail: mrenna@fnal.gov \| \| \n"
1542	<< " \| \| Stefan Prestel; Department of Astron"
1543	<< "omy and Theoretical Physics, \| \| \n"
1544	<< " \| \| Lund University, Solvegatan 14A, S"
1545	<< "E-223 62 Lund, Sweden; \| \| \n"
1546	<< " \| \| phone: + 46 - 46 - 222 06 64; e-ma"
1547	<< "il: stefan.prestel@thep.lu.se \| \| \n"
1548	<< " \| \| Peter Skands; Theoretical Physics, C"
1549	<< "ERN, CH-1211 Geneva 23, Switzerland; \| \| \n"
1550	<< " \| \| phone: + 41 - 22 - 767 2447; e-mai"
1551	<< "l: peter.skands@cern.ch \| \| \n"
1552	<< " \| \| "
1553	<< " \| \| \n"
1554	<< " \| \| The main program reference is the 'Br"
1555	<< "ief Introduction to PYTHIA 8.1', \| \| \n"
1556	<< " \| \| T. Sjostrand, S. Mrenna and P. Skands"
1557	<< ", Comput. Phys. Comm. 178 (2008) 85 \| \| \n"
1558	<< " \| \| [arXiv:0710.3820] "
1559	<< " \| \| \n"
1560	<< " \| \| "
1561	<< " \| \| \n"
1562	<< " \| \| The main physics reference is the 'PY"
1563	<< "THIA 6.4 Physics and Manual', \| \| \n"
1564	<< " \| \| T. Sjostrand, S. Mrenna and P. Skands"
1565	<< ", JHEP05 (2006) 026 [hep-ph/0603175]. \| \| \n"
1566	<< " \| \| "
1567	<< " \| \| \n"
1568	<< " \| \| An archive of program versions and do"
1569	<< "cumentation is found on the web: \| \| \n"
1570	<< " \| \| http://www.thep.lu.se/~torbjorn/Pythi"
1571	<< "a.html \| \| \n"
1572	<< " \| \| "
1573	<< " \| \| \n"
1574	<< " \| \| This program is released under the GN"
1575	<< "U General Public Licence version 2. \| \| \n"
1576	<< " \| \| Please respect the MCnet Guidelines f"
1577	<< "or Event Generator Authors and Users. \| \| \n"
1578	<< " \| \| "
1579	<< " \| \| \n"
1580	<< " \| \| Disclaimer: this program comes withou"
1581	<< "t any guarantees. \| \| \n"
1582	<< " \| \| Beware of errors and use common sense"
1583	<< " when interpreting results. \| \| \n"
1584	<< " \| \| "
1585	<< " \| \| \n"
1586	<< " \| \| Copyright (C) 2012 Torbjorn Sjostrand"
1587	<< " \| \| \n"
1588	<< " \| \| "
1589	<< " \| \| \n"
1590	<< " \| \| "
1591	<< " \| \| \n"
1592	<< " \| *----------------------------------------"
1593	<< "--------------------------------------* \| \n"
1594	<< " \| "
1595	<< " \| \n"
1596	<< " *-------------------------------------------"
1597	<< "-----------------------------------------* \n" << endl;
1598
1599	}
1600
1601	//--------------------------------------------------------------------------
1602
1603	// Check for lines in file that mark the beginning of new subrun.
1604
1605	int Pythia::readSubrun(string line, bool warn, ostream& os) {
1606
1607	// If empty line then done.
1608	int subrunLine = SUBRUNDEFAULT;
1609	if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos)
1610	return subrunLine;
1611
1612	// If first character is not a letter, then done.
1613	string lineNow = line;
1614	int firstChar = lineNow.find_first_not_of(" \n\t\v\b\r\f\a");
1615	if (!isalpha(lineNow[firstChar])) return subrunLine;
1616
1617	// Replace an equal sign by a blank to make parsing simpler.
1618	while (lineNow.find("=") != string::npos) {
1619	int firstEqual = lineNow.find_first_of("=");
1620	lineNow.replace(firstEqual, 1, " ");
1621	}
1622
1623	// Get first word of a line.
1624	istringstream splitLine(lineNow);
1625	string name;
1626	splitLine >> name;
1627
1628	// Replace two colons by one (:: -> :) to allow for such mistakes.
1629	while (name.find("::") != string::npos) {
1630	int firstColonColon = name.find_first_of("::");
1631	name.replace(firstColonColon, 2, ":");
1632	}
1633
1634	// Convert to lowercase.
1635	for (int i = 0; i < int(name.length()); ++i)
1636	name[i] = std::tolower(name[i]);
1637
1638	// If no match then done.
1639	if (name != "main:subrun") return subrunLine;
1640
1641	// Else find new subrun number and return it.
1642	splitLine >> subrunLine;
1643	if (!splitLine) {
1644	if (warn) os << "\n PYTHIA Warning: Main:subrun number not"
1645	<< " recognized; skip:\n " << line << endl;
1646	subrunLine = SUBRUNDEFAULT;
1647	}
1648	return subrunLine;
1649
1650	}
1651
1652	//--------------------------------------------------------------------------
1653
1654	// Check that the final event makes sense: no unknown id codes;
1655	// charge and energy-momentum conserved.
1656
1657	bool Pythia::check(ostream& os) {
1658
1659	// Reset.
1660	bool physical = true;
1661	bool listVertices = false;
1662	bool listHistory = false;
1663	bool listSystems = false;
1664	bool listBeams = false;
1665	iErrId.resize(0);
1666	iErrCol.resize(0);
1667	iErrNan.resize(0);
1668	iErrNanVtx.resize(0);
1669	Vec4 pSum;
1670	double chargeSum = 0.;
1671
1672	// Incoming beams counted with negative momentum and charge.
1673	if (doProcessLevel) {
1674	pSum = - (event[1].p() + event[2].p());
1675	chargeSum = - (event[1].charge() + event[2].charge());
1676
1677	// If no ProcessLevel then sum final state of process record.
1678	} else if (process.size() > 0) {
1679	pSum = - process[0].p();
1680	for (int i = 0; i < process.size(); ++i)
1681	if (process[i].isFinal()) chargeSum -= process[i].charge();
1682
1683	// If process not filled, then use outgoing primary in event.
1684	} else {
1685	pSum = - event[0].p();
1686	for (int i = 1; i < event.size(); ++i)
1687	if (event[i].statusAbs() < 10 \|\| event[i].statusAbs() == 23)
1688	chargeSum -= event[i].charge();
1689	}
1690	double eLab = abs(pSum.e());
1691
1692	// Loop over particles in the event.
1693	for (int i = 0; i < event.size(); ++i) {
1694
1695	// Look for any unrecognized particle codes.
1696	int id = event[i].id();
1697	if (id == 0 \|\| !particleData.isParticle(id)) {
1698	ostringstream errCode;
1699	errCode << ", i = " << i << ", id = " << id;
1700	info.errorMsg("Error in Pythia::check: "
1701	"unknown particle code", errCode.str());
1702	physical = false;
1703	iErrId.push_back(i);
1704
1705	// Check that colour assignments are the expected ones.
1706	} else {
1707	int colType = event[i].colType();
1708	int col = event[i].col();
1709	int acol = event[i].acol();
1710	if ( (colType == 0 && (col > 0 \|\| acol > 0))
1711	\|\| (colType == 1 && (col <= 0 \|\| acol > 0))
1712	\|\| (colType == -1 && (col > 0 \|\| acol <= 0))
1713	\|\| (colType == 2 && (col <= 0 \|\| acol <= 0)) ) {
1714	ostringstream errCode;
1715	errCode << ", i = " << i << ", id = " << id << " cols = " << col
1716	<< " " << acol;
1717	info.errorMsg("Error in Pythia::check: "
1718	"incorrect colours", errCode.str());
1719	physical = false;
1720	iErrCol.push_back(i);
1721	}
1722	}
1723
1724	// Look for particles with not-a-number energy/momentum/mass.
1725	if (abs(event[i].px()) >= 0. && abs(event[i].py()) >= 0.
1726	&& abs(event[i].pz()) >= 0. && abs(event[i].e()) >= 0.
1727	&& abs(event[i].m()) >= 0.) ;
1728	else {
1729	info.errorMsg("Error in Pythia::check: "
1730	"not-a-number energy/momentum/mass");
1731	physical = false;
1732	iErrNan.push_back(i);
1733	}
1734
1735	// Look for particles with not-a-number vertex/lifetime.
1736	if (abs(event[i].xProd()) >= 0. && abs(event[i].yProd()) >= 0.
1737	&& abs(event[i].zProd()) >= 0. && abs(event[i].tProd()) >= 0.
1738	&& abs(event[i].tau()) >= 0.) ;
1739	else {
1740	info.errorMsg("Error in Pythia::check: "
1741	"not-a-number vertex/lifetime");
1742	physical = false;
1743	listVertices = true;
1744	iErrNanVtx.push_back(i);
1745	}
1746
1747	// Add final-state four-momentum and charge.
1748	if (event[i].isFinal()) {
1749	pSum += event[i].p();
1750	chargeSum += event[i].charge();
1751	}
1752
1753	// End of particle loop.
1754	}
1755
1756	// Check energy-momentum/charge conservation.
1757	double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
1758	+ abs(pSum.pz());
1759	if (epDev > epTolErr * eLab) {
1760	info.errorMsg("Error in Pythia::check: energy-momentum not conserved");
1761	physical = false;
1762	} else if (epDev > epTolWarn * eLab) {
1763	info.errorMsg("Warning in Pythia::check: "
1764	"energy-momentum not quite conserved");
1765	}
1766	if (abs(chargeSum) > 0.1) {
1767	info.errorMsg("Error in Pythia::check: charge not conserved");
1768	physical = false;
1769	}
1770
1771	// Check that beams and event records agree on incoming partons.
1772	// Only meaningful for resolved beams.
1773	if (info.isResolved())
1774	for (int iSys = 0; iSys < beamA.sizeInit(); ++iSys) {
1775	int eventANw = partonSystems.getInA(iSys);
1776	int eventBNw = partonSystems.getInB(iSys);
1777	int beamANw = beamA[iSys].iPos();
1778	int beamBNw = beamB[iSys].iPos();
1779	if (eventANw != beamANw \|\| eventBNw != beamBNw) {
1780	info.errorMsg("Error in Pythia::check: "
1781	"event and beams records disagree");
1782	physical = false;
1783	listSystems = true;
1784	listBeams = true;
1785	}
1786	}
1787
1788	// Check that mother and daughter information match for each particle.
1789	vector<int> noMot;
1790	vector<int> noDau;
1791	vector< pair<int,int> > noMotDau;
1792	if (checkHistory) {
1793
1794	// Loop through the event and check that there are beam particles.
1795	bool hasBeams = false;
1796	for (int i = 0; i < event.size(); ++i) {
1797	int status = event[i].status();
1798	if (abs(status) == 12) hasBeams = true;
1799
1800	// Check that mother and daughter lists not empty where not expected to.
1801	vector<int> mList = event.motherList(i);
1802	vector<int> dList = event.daughterList(i);
1803	if (mList.size() == 0 && abs(status) != 11 && abs(status) != 12)
1804	noMot.push_back(i);
1805	if (dList.size() == 0 && status < 0 && status != -11)
1806	noDau.push_back(i);
1807
1808	// Check that the particle appears in the daughters list of each mother.
1809	for (int j = 0; j < int(mList.size()); ++j) {
1810	vector<int> dmList = event.daughterList( mList[j] );
1811	bool foundMatch = false;
1812	for (int k = 0; k < int(dmList.size()); ++k)
1813	if (dmList[k] == i) {
1814	foundMatch = true;
1815	break;
1816	}
1817	if (!hasBeams && mList.size() == 1 && mList[0] == 0) foundMatch = true;
1818	if (!foundMatch) {
1819	bool oldPair = false;
1820	for (int k = 0; k < int(noMotDau.size()); ++k)
1821	if (noMotDau[k].first == mList[j] && noMotDau[k].second == i) {
1822	oldPair = true;
1823	break;
1824	}
1825	if (!oldPair) noMotDau.push_back( make_pair( mList[j], i) );
1826	}
1827	}
1828
1829	// Check that the particle appears in the mothers list of each daughter.
1830	for (int j = 0; j < int(dList.size()); ++j) {
1831	vector<int> mdList = event.motherList( dList[j] );
1832	bool foundMatch = false;
1833	for (int k = 0; k < int(mdList.size()); ++k)
1834	if (mdList[k] == i) {
1835	foundMatch = true;
1836	break;
1837	}
1838	if (!foundMatch) {
1839	bool oldPair = false;
1840	for (int k = 0; k < int(noMotDau.size()); ++k)
1841	if (noMotDau[k].first == i && noMotDau[k].second == dList[j]) {
1842	oldPair = true;
1843	break;
1844	}
1845	if (!oldPair) noMotDau.push_back( make_pair( i, dList[j]) );
1846	}
1847	}
1848	}
1849
1850	// Warn if any errors were found.
1851	if (noMot.size() > 0 \|\| noDau.size() > 0 \|\| noMotDau.size() > 0) {
1852	info.errorMsg("Error in Pythia::check: "
1853	"mismatch in daughter and mother lists");
1854	physical = false;
1855	listHistory = true;
1856	}
1857	}
1858
1859	// Done for sensible events.
1860	if (physical) return true;
1861
1862	// Print (the first few) flawed events: local info.
1863	if (nErrEvent < nErrList) {
1864	os << "\n PYTHIA erroneous event info: \n";
1865	if (iErrId.size() > 0) {
1866	os << " unknown particle codes in lines ";
1867	for (int i = 0; i < int(iErrId.size()); ++i)
1868	os << iErrId[i] << " ";
1869	os << "\n";
1870	}
1871	if (iErrCol.size() > 0) {
1872	os << " incorrect colour assignments in lines ";
1873	for (int i = 0; i < int(iErrCol.size()); ++i)
1874	os << iErrCol[i] << " ";
1875	os << "\n";
1876	}
1877	if (iErrNan.size() > 0) {
1878	os << " not-a-number energy/momentum/mass in lines ";
1879	for (int i = 0; i < int(iErrNan.size()); ++i)
1880	os << iErrNan[i] << " ";
1881	os << "\n";
1882	}
1883	if (iErrNanVtx.size() > 0) {
1884	os << " not-a-number vertex/lifetime in lines ";
1885	for (int i = 0; i < int(iErrNanVtx.size()); ++i)
1886	os << iErrNanVtx[i] << " ";
1887	os << "\n";
1888	}
1889	if (epDev > epTolErr * eLab) os << scientific << setprecision(3)
1890	<< " total energy-momentum non-conservation = " << epDev << "\n";
1891	if (abs(chargeSum) > 0.1) os << fixed << setprecision(2)
1892	<< " total charge non-conservation = " << chargeSum << "\n";
1893	if (noMot.size() > 0) {
1894	os << " missing mothers for particles ";
1895	for (int i = 0; i < int(noMot.size()); ++i) os << noMot[i] << " ";
1896	os << "\n";
1897	}
1898	if (noDau.size() > 0) {
1899	os << " missing daughters for particles ";
1900	for (int i = 0; i < int(noDau.size()); ++i) os << noDau[i] << " ";
1901	os << "\n";
1902	}
1903	if (noMotDau.size() > 0) {
1904	os << " inconsistent history for (mother,daughter) pairs ";
1905	for (int i = 0; i < int(noMotDau.size()); ++i)
1906	os << "(" << noMotDau[i].first << "," << noMotDau[i].second << ") ";
1907	os << "\n";
1908	}
1909
1910	// Print (the first few) flawed events: standard listings.
1911	info.list();
1912	event.list(listVertices, listHistory);
1913	if (listSystems) partonSystems.list();
1914	if (listBeams) beamA.list();
1915	if (listBeams) beamB.list();
1916	}
1917
1918	// Update error counter. Done also for flawed event.
1919	++nErrEvent;
1920	return false;
1921
1922	}
1923
1924	//--------------------------------------------------------------------------
1925
1926	// Routine to set up a PDF pointer.
1927
1928	PDF* Pythia::getPDFPtr(int idIn, int sequence) {
1929
1930	// Temporary pointer to be returned.
1931	PDF* tempPDFPtr = 0;
1932
1933	// One option is to treat a Pomeron like a pi0.
1934	if (idIn == 990 && settings.mode("PDF:PomSet") == 2) idIn = 111;
1935
1936	// Proton beam, normal choice.
1937	if (abs(idIn) == 2212 && sequence == 1) {
1938	int pSet = settings.mode("PDF:pSet");
1939	bool useLHAPDF = settings.flag("PDF:useLHAPDF");
1940
1941	// Use internal sets.
1942	if (!useLHAPDF) {
1943	if (pSet == 1) tempPDFPtr = new GRV94L(idIn);
1944	else if (pSet == 2) tempPDFPtr = new CTEQ5L(idIn);
1945	else if (pSet <= 6)
1946	tempPDFPtr = new MSTWpdf(idIn, pSet - 2, xmlPath, &info);
1947	else tempPDFPtr = new CTEQ6pdf(idIn, pSet - 6, xmlPath, &info);
1948	}
1949
1950	// Use sets from LHAPDF.
1951	else {
1952	string LHAPDFset = settings.word("PDF:LHAPDFset");
1953	int LHAPDFmember = settings.mode("PDF:LHAPDFmember");
1954	tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 1, &info);
1955
1956	// Optionally allow extrapolation beyond x and Q2 limits.
1957	tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
1958	}
1959	}
1960
1961	// Proton beam, special choice for the hard process.
1962	else if (abs(idIn) == 2212) {
1963	int pSet = settings.mode("PDF:pHardSet");
1964	bool useLHAPDF = settings.flag("PDF:useHardLHAPDF");
1965
1966	// Use internal sets.
1967	if (!useLHAPDF) {
1968	if (pSet == 1) tempPDFPtr = new GRV94L(idIn);
1969	else if (pSet == 2) tempPDFPtr = new CTEQ5L(idIn);
1970	else if (pSet <= 6)
1971	tempPDFPtr = new MSTWpdf(idIn, pSet - 2, xmlPath, &info);
1972	else tempPDFPtr = new CTEQ6pdf(idIn, pSet - 6, xmlPath, &info);
1973	}
1974
1975	// Use sets from LHAPDF.
1976	else {
1977	string LHAPDFset = settings.word("PDF:hardLHAPDFset");
1978	int LHAPDFmember = settings.mode("PDF:hardLHAPDFmember");
1979	// May need to require LHAPDF to handle two sets simultaneously.
1980	int nSet = (settings.flag("PDF:useLHAPDF")) ? 2 : 1;
1981	tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, nSet, &info);
1982
1983	// Optionally allow extrapolation beyond x and Q2 limits.
1984	tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
1985	}
1986	}
1987
1988	// Pion beam (or, in one option, Pomeron beam).
1989	else if (abs(idIn) == 211 \|\| idIn == 111) {
1990	bool useLHAPDF = settings.flag("PDF:piUseLHAPDF");
1991
1992	// Use internal sets.
1993	if (!useLHAPDF) {
1994	tempPDFPtr = new GRVpiL(idIn);
1995	}
1996
1997	// Use sets from LHAPDF.
1998	else {
1999	string LHAPDFset = settings.word("PDF:piLHAPDFset");
2000	int LHAPDFmember = settings.mode("PDF:piLHAPDFmember");
2001	tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 1, &info);
2002
2003	// Optionally allow extrapolation beyond x and Q2 limits.
2004	tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
2005	}
2006	}
2007
2008	// Pomeron beam, if not treated like a pi0 beam.
2009	else if (idIn == 990) {
2010	int pomSet = settings.mode("PDF:PomSet");
2011	double rescale = settings.parm("PDF:PomRescale");
2012
2013	// A generic Q2-independent parametrization.
2014	if (pomSet == 1) {
2015	double gluonA = settings.parm("PDF:PomGluonA");
2016	double gluonB = settings.parm("PDF:PomGluonB");
2017	double quarkA = settings.parm("PDF:PomQuarkA");
2018	double quarkB = settings.parm("PDF:PomQuarkB");
2019	double quarkFrac = settings.parm("PDF:PomQuarkFrac");
2020	double strangeSupp = settings.parm("PDF:PomStrangeSupp");
2021	tempPDFPtr = new PomFix( 990, gluonA, gluonB, quarkA, quarkB,
2022	quarkFrac, strangeSupp);
2023	}
2024
2025	// The H1 Q2-dependent parametrizations. Initialization requires files.
2026	else if (pomSet == 3 \|\| pomSet == 4)
2027	tempPDFPtr = new PomH1FitAB( 990, pomSet - 2, rescale, xmlPath, &info);
2028	else if (pomSet == 5)
2029	tempPDFPtr = new PomH1Jets( 990, rescale, xmlPath, &info);
2030	else if (pomSet == 6)
2031	tempPDFPtr = new PomH1FitAB( 990, 3, rescale, xmlPath, &info);
2032	}
2033
2034	// Lepton beam: neutrino, resolved charged lepton or unresolved ditto.
2035	else if (abs(idIn) > 10 && abs(idIn) < 17) {
2036	if (abs(idIn)%2 == 0) tempPDFPtr = new NeutrinoPoint(idIn);
2037	else if (settings.flag("PDF:lepton")) tempPDFPtr = new Lepton(idIn);
2038	else tempPDFPtr = new LeptonPoint(idIn);
2039	}
2040
2041	// Failure for unrecognized particle.
2042	else tempPDFPtr = 0;
2043
2044	// Done.
2045	return tempPDFPtr;
2046	}
2047
2048	//--------------------------------------------------------------------------
2049
2050	// Initialize SUSY Les Houches Accord data.
2051
2052	bool Pythia::initSLHA() {
2053
2054	// Initial and settings values.
2055	int ifailLHE = 1;
2056	int ifailSpc = 1;
2057	int readFrom = settings.mode("SLHA:readFrom");
2058	string lhefFile = settings.word("Beams:LHEF");
2059	string lhefHeader = settings.word("Beams:LHEFheader");
2060	string slhaFile = settings.word("SLHA:file");
2061	int verboseSLHA = settings.mode("SLHA:verbose");
2062	bool slhaUseDec = settings.flag("SLHA:useDecayTable");
2063
2064	// No SUSY by default
2065	couplingsPtr->isSUSY = false;
2066
2067	// Option with no SLHA read-in at all.
2068	if (readFrom == 0) return true;
2069
2070	// First check LHEF header (if reading from LHEF)
2071	if (readFrom == 1) {
2072	if (lhefHeader != "void")
2073	ifailLHE = slha.readFile(lhefHeader, verboseSLHA, slhaUseDec );
2074	else if (lhefFile != "void")
2075	ifailLHE = slha.readFile(lhefFile, verboseSLHA, slhaUseDec );
2076	}
2077
2078	// If LHEF read successful, everything needed should already be ready
2079	if (ifailLHE == 0) {
2080	ifailSpc = 0;
2081	couplingsPtr->isSUSY = true;
2082	// If no LHEF file or no SLHA info in header, read from SLHA:file
2083	} else {
2084	lhefFile = "void";
2085	if ( settings.word("SLHA:file") == "none"
2086	\|\| settings.word("SLHA:file") == "void"
2087	\|\| settings.word("SLHA:file") == ""
2088	\|\| settings.word("SLHA:file") == " ") return true;
2089	ifailSpc = slha.readFile(slhaFile,verboseSLHA, slhaUseDec);
2090	}
2091
2092	// In case of problems, print error and fail init.
2093	if (ifailSpc != 0) {
2094	info.errorMsg("Error in Pythia::initSLHA: "
2095	"problem reading SLHA file", slhaFile);
2096	return false;
2097	} else {
2098	couplingsPtr->isSUSY = true;
2099	}
2100
2101	// Check spectrum for consistency. Switch off SUSY if necessary.
2102	ifailSpc = slha.checkSpectrum();
2103
2104	// ifail >= 1 : no MODSEL found -> don't switch on SUSY
2105	if (ifailSpc == 1) {
2106	// no SUSY, but MASS ok
2107	couplingsPtr->isSUSY = false;
2108	info.errorMsg("Warning in Pythia::initSLHA: "
2109	"No MODSEL found, keeping internal SUSY switched off");
2110	} else if (ifailSpc >= 2) {
2111	// no SUSY, but problems
2112	info.errorMsg("Warning in Pythia::initSLHA: "
2113	"Problem with SLHA MASS or QNUMBERS.");
2114	couplingsPtr->isSUSY = false;
2115	}
2116	// ifail = 0 : MODSEL found, spectrum OK
2117	else if (ifailSpc == 0) {
2118	// Print spectrum. Done.
2119	slha.printSpectrum(0);
2120	}
2121	else if (ifailSpc < 0) {
2122	info.errorMsg("Warning in Pythia::initSLHA: "
2123	"Problem with SLHA spectrum.",
2124	"\n Only using masses and switching off SUSY.");
2125	settings.flag("SUSY:all", false);
2126	couplingsPtr->isSUSY = false;
2127	slha.printSpectrum(ifailSpc);
2128	}
2129
2130	// Import qnumbers
2131	if ( (ifailSpc == 1 \|\| ifailSpc == 0) && slha.qnumbers.size() > 0) {
2132	for (int iQnum=0; iQnum < int(slha.qnumbers.size()); iQnum++) {
2133	// Always use positive id codes
2134	int id = abs(slha.qnumbers[iQnum](0));
2135	ostringstream idCode;
2136	idCode << id;
2137	if (particleData.isParticle(id)) {
2138	info.errorMsg("Warning in Pythia::initSLHA: "
2139	"ignoring QNUMBERS", "for id = "+idCode.str()
2140	+" (already exists)", true);
2141	} else {
2142	int qEM3 = slha.qnumbers[iQnum](1);
2143	int nSpins = slha.qnumbers[iQnum](2);
2144	int colRep = slha.qnumbers[iQnum](3);
2145	int hasAnti = slha.qnumbers[iQnum](4);
2146	// Translate colRep to PYTHIA colType
2147	int colType = 0;
2148	if (colRep == 3) colType = 1;
2149	else if (colRep == -3) colType = -1;
2150	else if (colRep == 8) colType = 2;
2151	else if (colRep == 6) colType = 3;
2152	else if (colRep == -6) colType = -3;
2153	// Default name: PDG code
2154	string name, antiName;
2155	ostringstream idStream;
2156	idStream<<id;
2157	name = idStream.str();
2158	antiName = "-"+name;
2159	if (iQnum < int(slha.qnumbersName.size())) {
2160	name = slha.qnumbersName[iQnum];
2161	if (iQnum < int(slha.qnumbersAntiName.size()))
2162	antiName = slha.qnumbersAntiName[iQnum];
2163	if (antiName == "") {
2164	if (name.find("+") != string::npos) {
2165	antiName = name;
2166	antiName.replace(antiName.find("+"),1,"-");
2167	} else if (name.find("-") != string::npos) {
2168	antiName = name;
2169	antiName.replace(antiName.find("-"),1,"+");
2170	} else {
2171	antiName = name+"bar";
2172	}
2173	}
2174	}
2175	if ( hasAnti == 0) {
2176	antiName = "";
2177	particleData.addParticle(id, name, nSpins, qEM3, colType);
2178	} else {
2179	particleData.addParticle(id, name, antiName, nSpins, qEM3, colType);
2180	}
2181	}
2182	}
2183	}
2184
2185	// Import mass spectrum.
2186	bool keepSM = settings.flag("SLHA:keepSM");
2187	double minMassSM = settings.parm("SLHA:minMassSM");
2188	double massMargin = settings.parm("SLHA:minDecayDeltaM");
2189	if (ifailSpc == 1 \|\| ifailSpc == 0) {
2190
2191	// Loop through to update particle data.
2192	int id = slha.mass.first();
2193	for (int i = 1; i <= slha.mass.size() ; i++) {
2194	double mass = abs(slha.mass(id));
2195
2196	// Ignore masses for known SM particles or particles with
2197	// default masses < minMassSM; overwrite masses for rest.
2198	if (keepSM && (id < 25 \|\| (id > 80 && id < 1000000))) ;
2199	else if (id < 1000000 && particleData.m0(id) < minMassSM) {
2200	ostringstream idCode;
2201	idCode << id;
2202	info.errorMsg("Warning in Pythia::initSLHA: "
2203	"ignoring MASS entry", "for id = "+idCode.str()
2204	+" (m0 < SLHA:minMassSM)", true);
2205	}
2206	else particleData.m0(id,mass);
2207	id = slha.mass.next();
2208	};
2209
2210	}
2211
2212	// Update decay data.
2213	for (int iTable=0; iTable < int(slha.decays.size()); iTable++) {
2214
2215	// Pointer to this SLHA table
2216	LHdecayTable* slhaTable=&(slha.decays[iTable]);
2217
2218	// Extract ID and create pointer to corresponding particle data object
2219	int idRes = slhaTable->getId();
2220	ParticleDataEntry* particlePtr
2221	= particleData.particleDataEntryPtr(idRes);
2222
2223	// Ignore decay channels for known SM particles or particles with
2224	// default masses < minMassSM; overwrite masses for rest.
2225	if (keepSM && (idRes < 25 \|\| (idRes > 80 && idRes < 1000000))) continue;
2226	else if (idRes < 1000000 && particleData.m0(idRes) < minMassSM) {
2227	ostringstream idCode;
2228	idCode << idRes;
2229	info.errorMsg("Warning in Pythia::initSLHA: "
2230	"ignoring DECAY table", "for id = " + idCode.str()
2231	+ " (m0 < SLHA:minMassSM)", true);
2232	continue;
2233	}
2234
2235	// Extract and store total width (absolute value, neg -> switch off)
2236	double widRes = abs(slhaTable->getWidth());
2237	particlePtr->setMWidth(widRes);
2238
2239	// Reset decay table of the particle. Allow decays, treat as resonance.
2240	if (slhaTable->size() > 0) {
2241	particlePtr->clearChannels();
2242	particleData.mayDecay(idRes,true);
2243	particleData.isResonance(idRes,true);
2244	}
2245
2246	// Reset to stable if width <= 0.0
2247	if (slhaTable->getWidth() <= 0.0) particleData.mayDecay(idRes,false);
2248
2249	// Set initial minimum mass.
2250	double brWTsum = 0.;
2251	double massWTsum = 0.;
2252
2253	// Loop over SLHA channels, import into Pythia, treating channels
2254	// with negative branching fractions as having the equivalent positive
2255	// branching fraction, but being switched off for this run
2256	for (int iChannel=0 ; iChannel<slhaTable->size(); iChannel++) {
2257	LHdecayChannel slhaChannel = slhaTable->getChannel(iChannel);
2258	double brat = slhaChannel.getBrat();
2259	vector<int> idDa = slhaChannel.getIdDa();
2260	if (idDa.size() >= 9) {
2261	info.errorMsg("Error in Pythia::initSLHA: "
2262	"max number of decay products is 8.");
2263	} else if (idDa.size() <= 1) {
2264	info.errorMsg("Error in Pythia::initSLHA: "
2265	"min number of decay products is 2.");
2266	}
2267	else {
2268	int onMode = 1;
2269	if (brat < 0.0) onMode = 0;
2270
2271	// Check phase space, including margin
2272	double massSum = massMargin;
2273	for (int jDa=0; jDa<int(idDa.size()); ++jDa)
2274	massSum += particleData.m0( idDa[jDa] );
2275	if (onMode == 1 && brat > 0.0 && massSum > particleData.m0(idRes) ) {
2276	// String containing decay name
2277	ostringstream errCode;
2278	errCode << idRes <<" ->";
2279	for (int jDa=0; jDa<int(idDa.size()); ++jDa) errCode<<" "<<idDa[jDa];
2280	info.errorMsg("Warning in Pythia::initSLHA: switching off decay",
2281	errCode.str() + " (mRes - mDa < minDecayDeltaM)"
2282	"\n (Note: cross sections will be scaled by remaining"
2283	" open branching fractions!)" , true);
2284	onMode=0;
2285	}
2286
2287	// Branching-ratio-weighted average mass in decay.
2288	brWTsum += abs(brat);
2289	massWTsum += abs(brat) * massSum;
2290
2291	// Add channel
2292	int id0 = idDa[0];
2293	int id1 = idDa[1];
2294	int id2 = (idDa.size() >= 3) ? idDa[2] : 0;
2295	int id3 = (idDa.size() >= 4) ? idDa[3] : 0;
2296	int id4 = (idDa.size() >= 5) ? idDa[4] : 0;
2297	int id5 = (idDa.size() >= 6) ? idDa[5] : 0;
2298	int id6 = (idDa.size() >= 7) ? idDa[6] : 0;
2299	int id7 = (idDa.size() >= 8) ? idDa[7] : 0;
2300	particlePtr->addChannel(onMode,abs(brat),101,
2301	id0,id1,id2,id3,id4,id5,id6,id7);
2302
2303	}
2304	}
2305
2306	// Set minimal mass, but always below nominal one.
2307	if (slhaTable->size() > 0) {
2308	double massAvg = massWTsum / brWTsum;
2309	double massMin = min( massAvg, particlePtr->m0()) - massMargin;
2310	particlePtr->setMMin(massMin);
2311	}
2312	}
2313
2314	return true;
2315
2316	}
2317
2318	//--------------------------------------------------------------------------
2319
2320	// Function to perform CKKW-L merging on the input event.
2321
2322	bool Pythia::mergeProcess() {
2323
2324	// Reset weight of the event.
2325	mergingHooksPtr->setWeightCKKWL(1.);
2326	info.setWeightCKKWL(1.);
2327	double wgt = 1.0;
2328
2329	// Store candidates for the splitting V -> qqbar'.
2330	mergingHooksPtr->storeHardProcessCandidates( process);
2331
2332	// Calculate number of clustering steps.
2333	int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( process);
2334
2335	// Check if event passes the merging scale cut.
2336	double tms = mergingHooksPtr->tms();
2337	// Get merging scale in current event.
2338	double tnow = 0.;
2339	// Use KT/Durham merging scale definition.
2340	if ( mergingHooksPtr->doKTMerging() \|\| mergingHooksPtr->doMGMerging() )
2341	tnow = mergingHooksPtr->kTms( process );
2342	// Use Lund PT merging scale definition.
2343	else if ( mergingHooksPtr->doPTLundMerging() )
2344	tnow = mergingHooksPtr->rhoms( process, false );
2345	// Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale definition.
2346	else if ( mergingHooksPtr->doCutBasedMerging() )
2347	tnow = mergingHooksPtr->cutbasedms( process );
2348	// Use user-defined merging scale.
2349	else
2350	tnow = mergingHooksPtr->tmsDefinition( process );
2351
2352	bool enforceCutOnLHE = settings.flag("Merging:enforceCutOnLHE");
2353	// Enfore merging scale cut if the event did not pass the merging scale
2354	// criterion.
2355	if ( enforceCutOnLHE && nSteps > 0 && tnow < tms ) {
2356	string message="Warning in Pythia::mergeProcess: Les Houches Event";
2357	message+=" fails merging scale cut. Cut by rejecting event.";
2358	info.errorMsg(message);
2359	return true;
2360	}
2361
2362	double TMSMISMATCHPOS = tms / 2. ;
2363	// Remember if event is significantly above the merging scale.
2364	if ( enforceCutOnLHE && nSteps > 0 && tms > 0.
2365	&& tnow > tms && (tnow-tms > TMSMISMATCHPOS) ) {
2366	stringstream tmsmis;
2367	tmsmis << int(TMSMISMATCHPOS);
2368	string message="Warning in Pythia::mergeProcess: Les Houches event";
2369	message+=" more than ";
2370	message+=tmsmis.str();
2371	message+=" GeV above desired merging scale value.";
2372	info.errorMsg(message);
2373	info.addCounter(41);
2374	}
2375
2376	// For now, prefer construction of ordered histories.
2377	mergingHooksPtr->orderHistories(true);
2378
2379	// Declare pdfWeight and alpha_s weight.
2380	double RN = rndm.flat();
2381
2382	process.scale(0.0);
2383	// Generate all histories.
2384	History FullHistory( nSteps, 0.0, process, Clustering(), mergingHooksPtr,
2385	beamA, beamB, &particleData, &info, true, true, true, true,
2386	1.0, 0);
2387
2388	// Project histories onto desired branches, e.g. only ordered paths.
2389	FullHistory.projectOntoDesiredHistories();
2390
2391	// Calculate CKKWL weight:
2392	// Perform reweighting with Sudakov factors, save alpha_s ratios and
2393	// PDF ratio weights.
2394	wgt = FullHistory.weightTREE( &trialPartonLevel,
2395	mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN);
2396
2397	// Event with production scales set for further (trial) showering
2398	// and starting conditions for the shower.
2399	FullHistory.getStartingConditions( RN, process );
2400
2401	// Allow to dampen histories in which the lowest multiplicity reclustered
2402	// state does not pass the lowest multiplicity cut of the matrix element.
2403	double dampWeight = mergingHooksPtr->dampenIfFailCuts(
2404	FullHistory.lowestMultProc(RN) );
2405
2406	// Save the weight of the event for histogramming. Only change the
2407	// event weight after trial shower on the matrix element
2408	// multiplicity event (= in doVetoStep).
2409	wgt *= dampWeight;
2410
2411	// Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
2412	// --> Set to minimal mT of partons.
2413	int nFinal = 0;
2414	double muf = process[0].e();
2415	for ( int i=0; i < process.size(); ++i )
2416	if ( process[i].isFinal()
2417	&& (process[i].colType() != 0 \|\| process[i].id() == 22 ) ) {
2418	nFinal++;
2419	muf = min( muf, abs(process[i].mT()) );
2420	}
2421	// For pure QCD dijet events (only!), set the process scale to the
2422	// transverse momentum of the outgoing partons.
2423	if ( nSteps == 0 && nFinal == 2
2424	&& ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
2425	\|\| mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
2426	process.scale(muf);
2427
2428	// Save the weight of the event for histogramming.
2429	mergingHooksPtr->setWeightCKKWL(wgt);
2430	info.setWeightCKKWL(wgt);
2431
2432	// If the weight of the event is zero, do not continue evolution.
2433	if(wgt == 0.) return true;
2434
2435	// Done
2436	return false;
2437
2438	}
2439
2440	//==========================================================================
2441
2442	} // end namespace Pythia8

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: HiSusy/trunk/Pythia8/pythia8170/src/Pythia.cc @ 1

Download in other formats: