1
// ProcessLevel.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.
6
// Function definitions (not found in the header) for the ProcessLevel class.
8
#include "ProcessLevel.h"
12
//==========================================================================
14
// The ProcessLevel class.
16
//--------------------------------------------------------------------------
18
// Constants: could be changed here if desired, but normally should not.
19
// These are of technical nature, as described for each.
21
// Allow a few failures in final construction of events.
22
const int ProcessLevel::MAXLOOP = 5;
24
//--------------------------------------------------------------------------
28
ProcessLevel::~ProcessLevel() {
30
// Run through list of first hard processes and delete them.
31
for (int i = 0; i < int(containerPtrs.size()); ++i)
32
delete containerPtrs[i];
34
// Run through list of second hard processes and delete them.
35
for (int i =0; i < int(container2Ptrs.size()); ++i)
36
delete container2Ptrs[i];
40
//--------------------------------------------------------------------------
42
// Main routine to initialize generation process.
44
bool ProcessLevel::init( Info* infoPtrIn, Settings& settings,
45
ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
46
BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
47
Couplings* couplingsPtrIn, SigmaTotal* sigmaTotPtrIn, bool doLHA,
48
SusyLesHouches* slhaPtrIn, UserHooks* userHooksPtrIn,
49
vector<SigmaProcess*>& sigmaPtrs, ostream& os) {
51
// Store input pointers for future use.
53
particleDataPtr = particleDataPtrIn;
55
beamAPtr = beamAPtrIn;
56
beamBPtr = beamBPtrIn;
57
couplingsPtr = couplingsPtrIn;
58
sigmaTotPtr = sigmaTotPtrIn;
59
userHooksPtr = userHooksPtrIn;
62
// Send on some input pointers.
63
resonanceDecays.init( infoPtr, particleDataPtr, rndmPtr);
65
// Set up SigmaTotal. Store sigma_nondiffractive for future use.
66
sigmaTotPtr->init( infoPtr, settings, particleDataPtr);
67
int idA = infoPtr->idA();
68
int idB = infoPtr->idB();
69
double eCM = infoPtr->eCM();
70
sigmaTotPtr->calc( idA, idB, eCM);
71
sigmaND = sigmaTotPtr->sigmaND();
73
// Options to allow second hard interaction and resonance decays.
74
doSecondHard = settings.flag("SecondHard:generate");
75
doSameCuts = settings.flag("PhaseSpace:sameForSecond");
76
doResDecays = settings.flag("ProcessLevel:resonanceDecays");
77
startColTag = settings.mode("Event:startColTag");
79
// Second interaction not to be combined with biased phase space.
80
if (doSecondHard && userHooksPtr != 0
81
&& userHooksPtr->canBiasSelection()) {
82
infoPtr->errorMsg("Error in ProcessLevel::init: "
83
"cannot combine second interaction with biased phase space");
87
// Mass and pT cuts for two hard processes.
88
mHatMin1 = settings.parm("PhaseSpace:mHatMin");
89
mHatMax1 = settings.parm("PhaseSpace:mHatMax");
90
if (mHatMax1 < mHatMin1) mHatMax1 = eCM;
91
pTHatMin1 = settings.parm("PhaseSpace:pTHatMin");
92
pTHatMax1 = settings.parm("PhaseSpace:pTHatMax");
93
if (pTHatMax1 < pTHatMin1) pTHatMax1 = eCM;
94
mHatMin2 = settings.parm("PhaseSpace:mHatMinSecond");
95
mHatMax2 = settings.parm("PhaseSpace:mHatMaxSecond");
96
if (mHatMax2 < mHatMin2) mHatMax2 = eCM;
97
pTHatMin2 = settings.parm("PhaseSpace:pTHatMinSecond");
98
pTHatMax2 = settings.parm("PhaseSpace:pTHatMaxSecond");
99
if (pTHatMax2 < pTHatMin2) pTHatMax2 = eCM;
101
// Check whether mass and pT ranges agree or overlap.
102
cutsAgree = doSameCuts;
103
if (mHatMin2 == mHatMin1 && mHatMax2 == mHatMax1 && pTHatMin2 == pTHatMin1
104
&& pTHatMax2 == pTHatMax1) cutsAgree = true;
105
cutsOverlap = cutsAgree;
107
bool mHatOverlap = (max( mHatMin1, mHatMin2)
108
< min( mHatMax1, mHatMax2));
109
bool pTHatOverlap = (max( pTHatMin1, pTHatMin2)
110
< min( pTHatMax1, pTHatMax2));
111
if (mHatOverlap && pTHatOverlap) cutsOverlap = true;
114
// Set up containers for all the internal hard processes.
115
SetupContainers setupContainers;
116
setupContainers.init(containerPtrs, settings, particleDataPtr, couplingsPtr);
118
// Append containers for external hard processes, if any.
119
if (sigmaPtrs.size() > 0) {
120
for (int iSig = 0; iSig < int(sigmaPtrs.size()); ++iSig)
121
containerPtrs.push_back( new ProcessContainer(sigmaPtrs[iSig],
125
// Append single container for Les Houches processes, if any.
127
SigmaProcess* sigmaPtr = new SigmaLHAProcess();
128
containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
130
// Store location of this container, and send in LHA pointer.
131
iLHACont = containerPtrs.size() - 1;
132
containerPtrs[iLHACont]->setLHAPtr(lhaUpPtr);
135
// If no processes found then refuse to do anything.
136
if ( int(containerPtrs.size()) == 0) {
137
infoPtr->errorMsg("Error in ProcessLevel::init: "
138
"no process switched on");
142
// Check that SUSY couplings were indeed initialized where necessary.
143
bool hasSUSY = false;
144
for (int i = 0; i < int(containerPtrs.size()); ++i)
145
if (containerPtrs[i]->isSUSY()) hasSUSY = true;
147
// If SUSY processes requested but no SUSY couplings present
148
if(hasSUSY && !couplingsPtr->isSUSY) {
149
infoPtr->errorMsg("Error in ProcessLevel::init: "
150
"SUSY process switched on but no SUSY couplings found");
154
// Initialize SLHA blocks SMINPUTS and MASS from PYTHIA SM parameter values.
157
// Initialize each process.
159
for (int i = 0; i < int(containerPtrs.size()); ++i)
160
if (containerPtrs[i]->init(true, infoPtr, settings, particleDataPtr,
161
rndmPtr, beamAPtr, beamBPtr, couplingsPtr, sigmaTotPtr,
162
&resonanceDecays, slhaPtr, userHooksPtr)) ++numberOn;
164
// Sum maxima for Monte Carlo choice.
166
for (int i = 0; i < int(containerPtrs.size()); ++i)
167
sigmaMaxSum += containerPtrs[i]->sigmaMax();
169
// Option to pick a second hard interaction: repeat as above.
172
setupContainers.init2(container2Ptrs, settings);
173
if ( int(container2Ptrs.size()) == 0) {
174
infoPtr->errorMsg("Error in ProcessLevel::init: "
175
"no second hard process switched on");
178
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
179
if (container2Ptrs[i2]->init(false, infoPtr, settings, particleDataPtr,
180
rndmPtr, beamAPtr, beamBPtr, couplingsPtr, sigmaTotPtr,
181
&resonanceDecays, slhaPtr, userHooksPtr)) ++number2On;
183
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
184
sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
187
// Printout during initialization is optional.
188
if (settings.flag("Init:showProcesses")) {
190
// Construct string with incoming beams and for cm energy.
191
string collision = "We collide " + particleDataPtr->name(idA)
192
+ " with " + particleDataPtr->name(idB) + " at a CM energy of ";
193
string pad( 51 - collision.length(), ' ');
195
// Print initialization information: header.
196
os << "\n *------- PYTHIA Process Initialization ---------"
197
<< "-----------------*\n"
200
<< " | " << collision << scientific << setprecision(3)
201
<< setw(9) << eCM << " GeV" << pad << " |\n"
204
<< " |---------------------------------------------------"
205
<< "---------------|\n"
208
<< " | Subprocess Code"
209
<< " | Estimated |\n"
214
<< " |---------------------------------------------------"
215
<< "---------------|\n"
220
// Loop over existing processes: print individual process info.
221
for (int i = 0; i < int(containerPtrs.size()); ++i)
222
os << " | " << left << setw(45) << containerPtrs[i]->name()
223
<< right << setw(5) << containerPtrs[i]->code() << " | "
224
<< scientific << setprecision(3) << setw(11)
225
<< containerPtrs[i]->sigmaMax() << " |\n";
227
// Loop over second hard processes, if any, and repeat as above.
231
<< " |---------------------------------------------------"
232
<<"---------------|\n"
235
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
236
os << " | " << left << setw(45) << container2Ptrs[i2]->name()
237
<< right << setw(5) << container2Ptrs[i2]->code() << " | "
238
<< scientific << setprecision(3) << setw(11)
239
<< container2Ptrs[i2]->sigmaMax() << " |\n";
245
<< " *------- End PYTHIA Process Initialization ----------"
246
<<"-------------*" << endl;
249
// If sum of maxima vanishes then refuse to do anything.
250
if ( numberOn == 0 || sigmaMaxSum <= 0.) {
251
infoPtr->errorMsg("Error in ProcessLevel::init: "
252
"all processes have vanishing cross sections");
255
if ( doSecondHard && (number2On == 0 || sigma2MaxSum <= 0.) ) {
256
infoPtr->errorMsg("Error in ProcessLevel::init: "
257
"all second hard processes have vanishing cross sections");
261
// If two hard processes then check whether some (but not all) agree.
265
bool foundMatch = false;
267
// Check for each first process if matched in second.
268
for (int i = 0; i < int(containerPtrs.size()); ++i) {
271
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
272
if (container2Ptrs[i2]->code() == containerPtrs[i]->code())
274
containerPtrs[i]->isSame( foundMatch );
275
if (!foundMatch) allHardSame = false;
276
if ( foundMatch) noneHardSame = false;
279
// Check for each second process if matched in first.
280
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
283
for (int i = 0; i < int(containerPtrs.size()); ++i)
284
if (containerPtrs[i]->code() == container2Ptrs[i2]->code())
286
container2Ptrs[i2]->isSame( foundMatch );
287
if (!foundMatch) allHardSame = false;
288
if ( foundMatch) noneHardSame = false;
292
// Concluding classification, including cuts.
293
if (!cutsAgree) allHardSame = false;
294
someHardSame = !allHardSame && !noneHardSame;
296
// Reset counters for average impact-parameter enhancement.
301
// Statistics for LHA events.
309
//--------------------------------------------------------------------------
311
// Main routine to generate the hard process.
313
bool ProcessLevel::next( Event& process) {
315
// Generate the next event with two or one hard interactions.
316
bool physical = (doSecondHard) ? nextTwo( process) : nextOne( process);
318
// Check that colour assignments make sense.
319
if (physical) physical = checkColours( process);
325
//--------------------------------------------------------------------------
327
// Accumulate and update statistics (after possible user veto).
329
void ProcessLevel::accumulate() {
331
// Increase number of accepted events.
332
containerPtrs[iContainer]->accumulate();
334
// Provide current generated cross section estimate.
338
double sigmaSum = 0.;
339
double delta2Sum = 0.;
340
double sigSelSum = 0.;
341
double weightSum = 0.;
342
for (int i = 0; i < int(containerPtrs.size()); ++i)
343
if (containerPtrs[i]->sigmaMax() != 0.) {
344
nTrySum += containerPtrs[i]->nTried();
345
nSelSum += containerPtrs[i]->nSelected();
346
nAccSum += containerPtrs[i]->nAccepted();
347
sigmaSum += containerPtrs[i]->sigmaMC();
348
delta2Sum += pow2(containerPtrs[i]->deltaMC());
349
sigSelSum += containerPtrs[i]->sigmaSelMC();
350
weightSum += containerPtrs[i]->weightSum();
353
// For Les Houches events find subprocess type and update counter.
354
if (infoPtr->isLHA()) {
355
int codeLHANow = infoPtr->codeSub();
357
for (int i = 0; i < int(codeLHA.size()); ++i)
358
if (codeLHANow == codeLHA[i]) iFill = i;
359
if (iFill >= 0) ++nEvtLHA[iFill];
361
// Add new process when new code and then arrange in order.
363
codeLHA.push_back(codeLHANow);
364
nEvtLHA.push_back(1);
365
for (int i = int(codeLHA.size()) - 1; i > 0; --i) {
366
if (codeLHA[i] < codeLHA[i - 1]) {
367
swap(codeLHA[i], codeLHA[i - 1]);
368
swap(nEvtLHA[i], nEvtLHA[i - 1]);
375
// Normally only one hard interaction. Then store info and done.
377
double deltaSum = sqrtpos(delta2Sum);
378
infoPtr->setSigma( nTrySum, nSelSum, nAccSum, sigmaSum, deltaSum,
383
// Increase counter for a second hard interaction.
384
container2Ptrs[i2Container]->accumulate();
386
// Update statistics on average impact factor.
388
sumImpactFac += infoPtr->enhanceMPI();
389
sum2ImpactFac += pow2(infoPtr->enhanceMPI());
391
// Cross section estimate for second hard process.
392
double sigma2Sum = 0.;
393
double sig2SelSum = 0.;
394
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
395
if (container2Ptrs[i2]->sigmaMax() != 0.) {
396
nTrySum += container2Ptrs[i2]->nTried();
397
sigma2Sum += container2Ptrs[i2]->sigmaMC();
398
sig2SelSum += container2Ptrs[i2]->sigmaSelMC();
401
// Average impact-parameter factor and error.
402
double invN = 1. / max(1, nImpact);
403
double impactFac = max( 1., sumImpactFac * invN);
404
double impactErr2 = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
406
// Cross section estimate for combination of first and second process.
407
// Combine two possible ways and take average.
408
double sigmaComb = 0.5 * (sigmaSum * sig2SelSum + sigSelSum * sigma2Sum);
409
sigmaComb *= impactFac / sigmaND;
410
if (allHardSame) sigmaComb *= 0.5;
411
double deltaComb = sqrtpos(2. / nAccSum + impactErr2) * sigmaComb;
413
// Store info and done.
414
infoPtr->setSigma( nTrySum, nSelSum, nAccSum, sigmaComb, deltaComb,
419
//--------------------------------------------------------------------------
421
// Print statistics on cross sections and number of events.
423
void ProcessLevel::statistics(bool reset, ostream& os) {
425
// Special processing if two hard interactions selected.
427
statistics2(reset, os);
432
os << "\n *------- PYTHIA Event and Cross Section Statistics ------"
433
<< "-------------------------------------------------------*\n"
436
<< " | Subprocess Code | "
437
<< " Number of events | sigma +- delta |\n"
439
<< "Tried Selected Accepted | (estimated) (mb) |\n"
442
<< " |------------------------------------------------------------"
443
<< "-----------------------------------------------------|\n"
447
// Reset sum counters.
451
double sigmaSum = 0.;
452
double delta2Sum = 0.;
454
// Loop over existing processes.
455
for (int i = 0; i < int(containerPtrs.size()); ++i)
456
if (containerPtrs[i]->sigmaMax() != 0.) {
458
// Read info for process. Sum counters.
459
long nTry = containerPtrs[i]->nTried();
460
long nSel = containerPtrs[i]->nSelected();
461
long nAcc = containerPtrs[i]->nAccepted();
462
double sigma = containerPtrs[i]->sigmaMC();
463
double delta = containerPtrs[i]->deltaMC();
468
delta2Sum += pow2(delta);
470
// Print individual process info.
471
os << " | " << left << setw(45) << containerPtrs[i]->name()
472
<< right << setw(5) << containerPtrs[i]->code() << " | "
473
<< setw(11) << nTry << " " << setw(10) << nSel << " "
474
<< setw(10) << nAcc << " | " << scientific << setprecision(3)
475
<< setw(11) << sigma << setw(11) << delta << " |\n";
477
// Print subdivision by user code for Les Houches process.
478
if (containerPtrs[i]->code() == 9999)
479
for (int j = 0; j < int(codeLHA.size()); ++j)
480
os << " | ... whereof user classification code " << setw(10)
481
<< codeLHA[j] << " | " << setw(10)
482
<< nEvtLHA[j] << " | | \n";
485
// Print summed process info.
488
<< " | " << left << setw(50) << "sum" << right << " | " << setw(11)
489
<< nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
490
<< nAccSum << " | " << scientific << setprecision(3) << setw(11)
491
<< sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
496
<< " *------- End PYTHIA Event and Cross Section Statistics -----"
497
<< "-----------------------------------------------------*" << endl;
499
// Optionally reset statistics contants.
500
if (reset) resetStatistics();
504
//--------------------------------------------------------------------------
506
// Reset statistics on cross sections and number of events.
508
void ProcessLevel::resetStatistics() {
510
for (int i = 0; i < int(containerPtrs.size()); ++i)
511
containerPtrs[i]->reset();
513
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
514
container2Ptrs[i2]->reset();
518
//--------------------------------------------------------------------------
520
// Initialize SLHA blocks SMINPUTS and MASS from PYTHIA SM parameter values.
522
void ProcessLevel::initSLHA() {
524
// Initialize block SMINPUTS.
525
string blockName = "sminputs";
526
double mZ = particleDataPtr->m0(23);
527
slhaPtr->set(blockName,1,1.0/couplingsPtr->alphaEM(pow2(mZ)));
528
slhaPtr->set(blockName,2,couplingsPtr->GF());
529
slhaPtr->set(blockName,3,couplingsPtr->alphaS(pow2(mZ)));
530
slhaPtr->set(blockName,4,mZ);
531
// b mass (should be running mass, here pole for time being)
532
slhaPtr->set(blockName,5,particleDataPtr->m0(5));
533
slhaPtr->set(blockName,6,particleDataPtr->m0(6));
534
slhaPtr->set(blockName,7,particleDataPtr->m0(15));
535
slhaPtr->set(blockName,8,particleDataPtr->m0(16));
536
slhaPtr->set(blockName,11,particleDataPtr->m0(11));
537
slhaPtr->set(blockName,12,particleDataPtr->m0(12));
538
slhaPtr->set(blockName,13,particleDataPtr->m0(13));
539
slhaPtr->set(blockName,14,particleDataPtr->m0(14));
540
// Force 3 lightest quarks massless
541
slhaPtr->set(blockName,21,double(0.0));
542
slhaPtr->set(blockName,22,double(0.0));
543
slhaPtr->set(blockName,23,double(0.0));
544
// c mass (should be running mass, here pole for time being)
545
slhaPtr->set(blockName,24,particleDataPtr->m0(4));
547
// Initialize block MASS.
551
while (particleDataPtr->nextId(id) > id) {
552
slhaPtr->set(blockName,id,particleDataPtr->m0(id));
553
id = particleDataPtr->nextId(id);
556
infoPtr->errorMsg("Error in ProcessLevel::initSLHA: "
557
"encountered infinite loop when saving mass block");
564
//--------------------------------------------------------------------------
566
// Generate the next event with one interaction.
568
bool ProcessLevel::nextOne( Event& process) {
570
// Update CM energy for phase space selection.
571
double eCM = infoPtr->eCM();
572
for (int i = 0; i < int(containerPtrs.size()); ++i)
573
containerPtrs[i]->newECM(eCM);
575
// Outer loop in case of rare failures.
576
bool physical = true;
577
for (int loop = 0; loop < MAXLOOP; ++loop) {
578
if (!physical) process.clear();
581
// Loop over tries until trial event succeeds.
584
// Pick one of the subprocesses.
585
double sigmaMaxNow = sigmaMaxSum * rndmPtr->flat();
586
int iMax = containerPtrs.size() - 1;
588
do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
589
while (sigmaMaxNow > 0. && iContainer < iMax);
591
// Do a trial event of this subprocess; accept or not.
592
if (containerPtrs[iContainer]->trialProcess()) break;
594
// Check for end-of-file condition for Les Houches events.
595
if (infoPtr->atEndOfFile()) return false;
598
// Update sum of maxima if current maximum violated.
599
if (containerPtrs[iContainer]->newSigmaMax()) {
601
for (int i = 0; i < int(containerPtrs.size()); ++i)
602
sigmaMaxSum += containerPtrs[i]->sigmaMax();
605
// Construct kinematics of acceptable process.
606
if ( !containerPtrs[iContainer]->constructProcess( process) )
609
// Do all resonance decays.
610
if ( physical && doResDecays
611
&& !containerPtrs[iContainer]->decayResonances( process) )
614
// Add any junctions to the process event record list.
615
if (physical) findJunctions( process);
617
// Outer loop should normally work first time around.
625
//--------------------------------------------------------------------------
627
// Generate the next event with two hard interactions.
629
bool ProcessLevel::nextTwo( Event& process) {
631
// Update CM energy for phase space selection.
632
double eCM = infoPtr->eCM();
633
for (int i = 0; i < int(containerPtrs.size()); ++i)
634
containerPtrs[i]->newECM(eCM);
635
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
636
container2Ptrs[i2]->newECM(eCM);
638
// Outer loop in case of rare failures.
639
bool physical = true;
640
for (int loop = 0; loop < MAXLOOP; ++loop) {
641
if (!physical) process.clear();
644
// Loop over both hard processes to find consistent common kinematics.
647
// Loop internally over tries for hardest process until succeeds.
650
// Pick one of the subprocesses.
651
double sigmaMaxNow = sigmaMaxSum * rndmPtr->flat();
652
int iMax = containerPtrs.size() - 1;
654
do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
655
while (sigmaMaxNow > 0. && iContainer < iMax);
657
// Do a trial event of this subprocess; accept or not.
658
if (containerPtrs[iContainer]->trialProcess()) break;
660
// Check for end-of-file condition for Les Houches events.
661
if (infoPtr->atEndOfFile()) return false;
664
// Update sum of maxima if current maximum violated.
665
if (containerPtrs[iContainer]->newSigmaMax()) {
667
for (int i = 0; i < int(containerPtrs.size()); ++i)
668
sigmaMaxSum += containerPtrs[i]->sigmaMax();
671
// Loop internally over tries for second hardest process until succeeds.
674
// Pick one of the subprocesses.
675
double sigma2MaxNow = sigma2MaxSum * rndmPtr->flat();
676
int i2Max = container2Ptrs.size() - 1;
678
do sigma2MaxNow -= container2Ptrs[++i2Container]->sigmaMax();
679
while (sigma2MaxNow > 0. && i2Container < i2Max);
681
// Do a trial event of this subprocess; accept or not.
682
if (container2Ptrs[i2Container]->trialProcess()) break;
685
// Update sum of maxima if current maximum violated.
686
if (container2Ptrs[i2Container]->newSigmaMax()) {
688
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
689
sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
692
// Check whether common set of x values is kinematically possible.
693
double xA1 = containerPtrs[iContainer]->x1();
694
double xB1 = containerPtrs[iContainer]->x2();
695
double xA2 = container2Ptrs[i2Container]->x1();
696
double xB2 = container2Ptrs[i2Container]->x2();
697
if (xA1 + xA2 >= 1. || xB1 + xB2 >= 1.) continue;
699
// Reset beam contents. Naive parton densities for second interaction.
700
// (Subsequent procedure could be symmetrized, but would be overkill.)
703
int idA1 = containerPtrs[iContainer]->id1();
704
int idB1 = containerPtrs[iContainer]->id2();
705
int idA2 = container2Ptrs[i2Container]->id1();
706
int idB2 = container2Ptrs[i2Container]->id2();
707
double Q2Fac1 = containerPtrs[iContainer]->Q2Fac();
708
double Q2Fac2 = container2Ptrs[i2Container]->Q2Fac();
709
double pdfA2Raw = beamAPtr->xf( idA2, xA2,Q2Fac2);
710
double pdfB2Raw = beamBPtr->xf( idB2, xB2,Q2Fac2);
712
// Remove partons in first interaction from beams.
713
beamAPtr->append( 3, idA1, xA1);
714
beamAPtr->xfISR( 0, idA1, xA1, Q2Fac1);
715
beamAPtr->pickValSeaComp();
716
beamBPtr->append( 4, idB1, xB1);
717
beamBPtr->xfISR( 0, idB1, xB1, Q2Fac1);
718
beamBPtr->pickValSeaComp();
720
// Reevaluate pdf's for second interaction and weight by reduction.
721
double pdfA2Mod = beamAPtr->xfMPI( idA2, xA2,Q2Fac2);
722
double pdfB2Mod = beamBPtr->xfMPI( idB2, xB2,Q2Fac2);
723
double wtPdfMod = (pdfA2Mod * pdfB2Mod) / (pdfA2Raw * pdfB2Raw);
724
if (wtPdfMod < rndmPtr->flat()) continue;
726
// Reduce by a factor of 2 for identical processes when others not,
727
// and when in same phase space region.
729
if ( someHardSame && containerPtrs[iContainer]->isSame()
730
&& container2Ptrs[i2Container]->isSame()) {
732
if (rndmPtr->flat() > 0.5) toLoop = true;
734
double mHat1 = containerPtrs[iContainer]->mHat();
735
double pTHat1 = containerPtrs[iContainer]->pTHat();
736
double mHat2 = container2Ptrs[i2Container]->mHat();
737
double pTHat2 = container2Ptrs[i2Container]->pTHat();
738
if (mHat1 > mHatMin2 && mHat1 < mHatMax2
739
&& pTHat1 > pTHatMin2 && pTHat1 < pTHatMax2
740
&& mHat2 > mHatMin1 && mHat2 < mHatMax1
741
&& pTHat2 > pTHatMin1 && pTHat2 < pTHatMax1
742
&& rndmPtr->flat() > 0.5) toLoop = true;
745
if (toLoop) continue;
747
// If come this far then acceptable event.
751
// Construct kinematics of acceptable processes.
753
process2.init( "(second hard)", particleDataPtr, startColTag);
754
process2.initColTag();
755
if ( !containerPtrs[iContainer]->constructProcess( process) )
757
if (physical && !container2Ptrs[i2Container]->constructProcess( process2,
758
false) ) physical = false;
760
// Do all resonance decays.
761
if ( physical && doResDecays
762
&& !containerPtrs[iContainer]->decayResonances( process) )
764
if ( physical && doResDecays
765
&& !container2Ptrs[i2Container]->decayResonances( process2) )
768
// Append second hard interaction to normal process object.
769
if (physical) combineProcessRecords( process, process2);
771
// Add any junctions to the process event record list.
772
if (physical) findJunctions( process);
774
// Outer loop should normally work first time around.
782
//--------------------------------------------------------------------------
784
// Append second hard interaction to normal process object.
785
// Complication: all resonance decay chains must be put at the end.
787
void ProcessLevel::combineProcessRecords( Event& process, Event& process2) {
789
// Find first event record size, excluding resonances.
790
int nSize = process.size();
792
while (nHard < nSize && process[nHard].mother1() == 3) ++nHard;
794
// Save resonance products temporarily elsewhere.
795
vector<Particle> resProd;
797
for (int i = nHard; i < nSize; ++i) resProd.push_back( process[i] );
798
process.popBack(nSize - nHard);
801
// Find second event record size, excluding resonances.
802
int nSize2 = process2.size();
804
while (nHard2 < nSize2 && process2[nHard2].mother1() == 3) ++nHard2;
806
// Find amount of necessary position and colour offset for second process.
807
int addPos = nHard - 3;
808
int addCol = process.lastColTag() - startColTag;
810
// Loop over all particles (except beams) from second process.
811
for (int i = 3; i < nSize2; ++i) {
813
// Offset mother and daughter pointers and colour tags of particle.
814
process2[i].offsetHistory( 2, addPos, 2, addPos);
815
process2[i].offsetCol( addCol);
817
// Append hard-process particles from process2 to process.
818
if (i < nHard2) process.append( process2[i] );
821
// Reinsert resonance decay chains of first hard process.
822
int addPos2 = nHard2 - 3;
825
// Offset daughter pointers of unmoved mothers.
826
for (int i = 5; i < nHard; ++i)
827
process[i].offsetHistory( 0, 0, nHard - 1, addPos2);
829
// Modify history of resonance products when restoring.
830
for (int i = 0; i < int(resProd.size()); ++i) {
831
resProd[i].offsetHistory( nHard - 1, addPos2, nHard - 1, addPos2);
832
process.append( resProd[i] );
836
// Insert resonance decay chains of second hard process.
837
if (nHard2 < nSize2) {
838
int nHard3 = nHard + nHard2 - 3;
839
int addPos3 = nSize - nHard;
841
// Offset daughter pointers of second-process mothers.
842
for (int i = nHard + 2; i < nHard3; ++i)
843
process[i].offsetHistory( 0, 0, nHard3 - 1, addPos3);
845
// Modify history of second-process resonance products and insert.
846
for (int i = nHard2; i < nSize2; ++i) {
847
process2[i].offsetHistory( nHard3 - 1, addPos3, nHard3 - 1, addPos3);
848
process.append( process2[i] );
852
// Store PDF scale for second interaction.
853
process.scaleSecond( process2.scale() );
857
//--------------------------------------------------------------------------
859
// Add any junctions to the process event record list.
860
// Also check that do not doublebook if called repeatedly.
862
void ProcessLevel::findJunctions( Event& junEvent) {
864
// Check all hard vertices for BNV
865
for (int i = 1; i<junEvent.size(); i++) {
867
// Ignore colorless particles and stages before hard-scattering
869
if (abs(junEvent[i].status()) <= 21 || junEvent[i].colType() == 0) continue;
870
vector<int> motherList = junEvent.motherList(i);
871
int iMot1 = motherList[0];
872
vector<int> sisterList = junEvent.daughterList(iMot1);
874
// Check baryon number of vertex.
876
map<int,int> colVertex, acolVertex;
878
// Loop over mothers (enter with crossed colors and negative sign).
879
for (unsigned int indx = 0; indx < motherList.size(); indx++) {
880
int iMot = motherList[indx];
881
if ( abs(junEvent[iMot].colType()) == 1 )
882
barSum -= junEvent[iMot].colType();
883
else if ( abs(junEvent[iMot].colType()) == 3)
884
barSum -= 2*junEvent[iMot].colType()/3;
885
int col = junEvent[iMot].acol();
886
int acol = junEvent[iMot].col();
888
// If unmatched (so far), add end. Else erase matching parton.
890
if (acolVertex.find(col) == acolVertex.end() ) colVertex[col] = iMot;
891
else acolVertex.erase(col);
892
} else if (col < 0) {
893
if (colVertex.find(-col) == colVertex.end() ) acolVertex[-col] = iMot;
894
else colVertex.erase(-col);
897
if (colVertex.find(acol) == colVertex.end()) acolVertex[acol] = iMot;
898
else colVertex.erase(acol);
899
} else if (acol < 0) {
900
if (acolVertex.find(-acol) == acolVertex.end()) colVertex[-acol] = iMot;
901
else acolVertex.erase(-acol);
905
// Loop over sisters.
906
for (unsigned int indx = 0; indx < sisterList.size(); indx++) {
907
int iDau = sisterList[indx];
908
if ( abs(junEvent[iDau].colType()) == 1 )
909
barSum += junEvent[iDau].colType();
910
else if ( abs(junEvent[iDau].colType()) == 3)
911
barSum += 2*junEvent[iDau].colType()/3;
912
int col = junEvent[iDau].col();
913
int acol = junEvent[iDau].acol();
915
// If unmatched (so far), add end. Else erase matching parton.
917
if (acolVertex.find(col) == acolVertex.end() ) colVertex[col] = iDau;
918
else acolVertex.erase(col);
919
} else if (col < 0) {
920
if (colVertex.find(-col) == colVertex.end() ) acolVertex[-col] = iDau;
921
else colVertex.erase(-col);
924
if (colVertex.find(acol) == colVertex.end()) acolVertex[acol] = iDau;
925
else colVertex.erase(acol);
926
} else if (acol < 0) {
927
if (acolVertex.find(-acol) == acolVertex.end()) colVertex[-acol] = iDau;
928
else acolVertex.erase(-acol);
933
// Skip if baryon number conserved in this vertex.
934
if (barSum == 0) continue;
936
// Check and skip any junctions that have already been added.
937
for (int iJun = 0; iJun < junEvent.sizeJunction(); ++iJun) {
938
// Remove the tags corresponding to each of the 3 existing junction legs.
939
for (int j = 0; j < 3; ++j) {
940
int colNow = junEvent.colJunction(iJun, j);
941
if (junEvent.kindJunction(iJun) % 2 == 1) colVertex.erase(colNow);
942
else acolVertex.erase(colNow);
946
// Skip if no junction colors remain.
947
if (colVertex.size() == 0 && acolVertex.size() == 0) continue;
949
// If baryon number violated, is B = +1 or -1 (larger values not handled).
951
if (colVertex.size() == 3 && acolVertex.size() == 0) kindJun = 1;
952
else if (colVertex.size() == 0 && acolVertex.size() == 3) kindJun = 2;
954
infoPtr->errorMsg("Error in ProcessLevel::findJunctions: "
955
"N(unmatched (anti)colour tags) != 3");
959
// From now on, use colJun as shorthand for colVertex or acolVertex.
960
map<int,int> colJun = (kindJun == 1) ? colVertex : acolVertex;
962
// Order so incoming tags appear first in colVec, outgoing tags last.
964
for (map<int,int>::iterator it = colJun.begin();
965
it != colJun.end(); it++) {
967
int iCol = it->second;
968
for (unsigned int indx = 0; indx < motherList.size(); indx++) {
969
if (iCol == motherList[indx]) {
971
colVec.insert(colVec.begin(),col);
974
if (colVec.size() == 0 || colVec[0] != col) colVec.push_back(col);
977
// Add junction with these tags.
978
junEvent.appendJunction( kindJun, colVec[0], colVec[1], colVec[2]);
983
//--------------------------------------------------------------------------
985
// Check that colours match up.
987
bool ProcessLevel::checkColours( Event& process) {
989
// Variables and arrays for common usage.
990
bool physical = true;
992
int colType, col, acol, iPos, iNow, iNowA;
993
vector<int> colTags, colPos, acolPos;
995
// Check that each particle has the kind of colours expected of it.
996
for (int i = 0; i < process.size(); ++i) {
997
colType = process[i].colType();
998
col = process[i].col();
999
acol = process[i].acol();
1000
if (colType == 0 && (col != 0 || acol != 0)) physical = false;
1001
else if (colType == 1 && (col <= 0 || acol != 0)) physical = false;
1002
else if (colType == -1 && (col != 0 || acol <= 0)) physical = false;
1003
else if (colType == 2 && (col <= 0 || acol <= 0)) physical = false;
1004
// Preparations for colour-sextet assignments
1005
// (colour,colour) = (colour,negative anticolour)
1006
else if (colType == 3 && (col <= 0 || acol >= 0)) physical = false;
1007
else if (colType == -3 && (col >= 0 || acol <= 0)) physical = false;
1009
else if (colType < -1 || colType > 3) physical = false;
1011
// Add to the list of colour tags.
1014
for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1015
if (col == colTags[ic]) match = true;
1016
if (!match) colTags.push_back(col);
1017
} else if (acol > 0) {
1019
for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1020
if (acol == colTags[ic]) match = true;
1021
if (!match) colTags.push_back(acol);
1023
// Colour sextets : map negative colour -> anticolour and vice versa
1026
for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1027
if (-col == colTags[ic]) match = true;
1028
if (!match) colTags.push_back(-col);
1029
} else if (acol < 0) {
1031
for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1032
if (-acol == colTags[ic]) match = true;
1033
if (!match) colTags.push_back(-acol);
1037
// Warn and give up if particles did not have the expected colours.
1039
infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
1040
"incorrect colour assignment");
1044
// Remove (anti)colours coming from an (anti)junction.
1045
for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1046
for (int j = 0; j < 3; ++j) {
1047
int colJun = process.colJunction(iJun, j);
1048
for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1049
if (colJun == colTags[ic]) {
1050
colTags[ic] = colTags[colTags.size() - 1];
1057
// Loop through all colour tags and find their positions (by sign).
1058
for (int ic = 0; ic < int(colTags.size()); ++ic) {
1062
for (int i = 0; i < process.size(); ++i) {
1063
if (process[i].col() == col || process[i].acol() == -col)
1064
colPos.push_back(i);
1065
if (process[i].acol() == col || process[i].col() == -col)
1066
acolPos.push_back(i);
1069
// Trace colours back through decays; remove daughters.
1070
while (colPos.size() > 1) {
1071
iPos = colPos.size() - 1;
1072
iNow = colPos[iPos];
1073
if ( process[iNow].mother1() == colPos[iPos - 1]
1074
&& process[iNow].mother2() == 0) colPos.pop_back();
1077
while (acolPos.size() > 1) {
1078
iPos = acolPos.size() - 1;
1079
iNow = acolPos[iPos];
1080
if ( process[iNow].mother1() == acolPos[iPos - 1]
1081
&& process[iNow].mother2() == 0) acolPos.pop_back();
1085
// Now colour should exist in only 2 copies.
1086
if (colPos.size() + acolPos.size() != 2) physical = false;
1088
// If both colours or both anticolours then one mother of the other.
1089
else if (colPos.size() == 2) {
1091
if ( process[iNow].mother1() != colPos[0]
1092
&& process[iNow].mother2() != colPos[0] ) physical = false;
1094
else if (acolPos.size() == 2) {
1096
if ( process[iNowA].mother1() != acolPos[0]
1097
&& process[iNowA].mother2() != acolPos[0] ) physical = false;
1100
// If one of each then should have same mother(s), or point to beams.
1104
if ( process[iNow].status() == -21 && process[iNowA].status() == -21 );
1105
else if ( (process[iNow].mother1() != process[iNowA].mother1())
1106
|| (process[iNow].mother2() != process[iNowA].mother2()) )
1112
// Error message if problem found. Done.
1113
if (!physical) infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
1114
"unphysical colour flow");
1119
//--------------------------------------------------------------------------
1121
// Print statistics when two hard processes allowed.
1123
void ProcessLevel::statistics2(bool reset, ostream& os) {
1125
// Average impact-parameter factor and error.
1126
double invN = 1. / max(1, nImpact);
1127
double impactFac = max( 1., sumImpactFac * invN);
1128
double impactErr2 = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
1130
// Derive scaling factor to be applied to first set of processes.
1131
double sigma2SelSum = 0.;
1133
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
1134
sigma2SelSum += container2Ptrs[i2]->sigmaSelMC();
1135
n2SelSum += container2Ptrs[i2]->nSelected();
1137
double factor1 = impactFac * sigma2SelSum / sigmaND;
1138
double rel1Err = sqrt(1. / max(1, n2SelSum) + impactErr2);
1139
if (allHardSame) factor1 *= 0.5;
1141
// Derive scaling factor to be applied to second set of processes.
1142
double sigma1SelSum = 0.;
1144
for (int i = 0; i < int(containerPtrs.size()); ++i) {
1145
sigma1SelSum += containerPtrs[i]->sigmaSelMC();
1146
n1SelSum += containerPtrs[i]->nSelected();
1148
double factor2 = impactFac * sigma1SelSum / sigmaND;
1149
if (allHardSame) factor2 *= 0.5;
1150
double rel2Err = sqrt(1. / max(1, n1SelSum) + impactErr2);
1153
os << "\n *------- PYTHIA Event and Cross Section Statistics ------"
1154
<< "--------------------------------------------------*\n"
1157
<< " | Subprocess Code | "
1158
<< "Number of events | sigma +- delta |\n"
1160
<< " Selected Accepted | (estimated) (mb) |\n"
1163
<< " |------------------------------------------------------------"
1164
<< "------------------------------------------------|\n"
1167
<< " | First hard process: | "
1172
// Reset sum counters.
1176
double sigmaSum = 0.;
1177
double delta2Sum = 0.;
1179
// Loop over existing first processes.
1180
for (int i = 0; i < int(containerPtrs.size()); ++i)
1181
if (containerPtrs[i]->sigmaMax() != 0.) {
1183
// Read info for process. Sum counters.
1184
long nTry = containerPtrs[i]->nTried();
1185
long nSel = containerPtrs[i]->nSelected();
1186
long nAcc = containerPtrs[i]->nAccepted();
1187
double sigma = containerPtrs[i]->sigmaMC() * factor1;
1188
double delta2 = pow2( containerPtrs[i]->deltaMC() * factor1 );
1193
delta2Sum += delta2;
1194
delta2 += pow2( sigma * rel1Err );
1196
// Print individual process info.
1197
os << " | " << left << setw(40) << containerPtrs[i]->name()
1198
<< right << setw(5) << containerPtrs[i]->code() << " | "
1199
<< setw(11) << nTry << " " << setw(10) << nSel << " "
1200
<< setw(10) << nAcc << " | " << scientific << setprecision(3)
1201
<< setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
1204
// Print summed info for first processes.
1205
delta2Sum += pow2( sigmaSum * rel1Err );
1208
<< " | " << left << setw(45) << "sum" << right << " | " << setw(11)
1209
<< nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
1210
<< nAccSum << " | " << scientific << setprecision(3) << setw(11)
1211
<< sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
1214
// Separation lines to second hard processes.
1217
<< " |------------------------------------------------------------"
1218
<< "------------------------------------------------|\n"
1221
<< " | Second hard process: | "
1226
// Reset sum counters.
1233
// Loop over existing second processes.
1234
for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
1235
if (container2Ptrs[i2]->sigmaMax() != 0.) {
1237
// Read info for process. Sum counters.
1238
long nTry = container2Ptrs[i2]->nTried();
1239
long nSel = container2Ptrs[i2]->nSelected();
1240
long nAcc = container2Ptrs[i2]->nAccepted();
1241
double sigma = container2Ptrs[i2]->sigmaMC() * factor2;
1242
double delta2 = pow2( container2Ptrs[i2]->deltaMC() * factor2 );
1247
delta2Sum += delta2;
1248
delta2 += pow2( sigma * rel2Err );
1250
// Print individual process info.
1251
os << " | " << left << setw(40) << container2Ptrs[i2]->name()
1252
<< right << setw(5) << container2Ptrs[i2]->code() << " | "
1253
<< setw(11) << nTry << " " << setw(10) << nSel << " "
1254
<< setw(10) << nAcc << " | " << scientific << setprecision(3)
1255
<< setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
1258
// Print summed info for second processes.
1259
delta2Sum += pow2( sigmaSum * rel2Err );
1262
<< " | " << left << setw(45) << "sum" << right << " | " << setw(11)
1263
<< nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
1264
<< nAccSum << " | " << scientific << setprecision(3) << setw(11)
1265
<< sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
1267
// Print information on how the two processes were combined.
1270
<< " |------------------------------------------------------------"
1271
<< "------------------------------------------------|\n"
1274
<< " | Uncombined cross sections for the two event sets were "
1275
<< setw(10) << sigma1SelSum << " and " << sigma2SelSum << " mb, "
1276
<< "respectively, combined |\n"
1277
<< " | using a sigma(nonDiffractive) of " << setw(10) << sigmaND
1278
<< " mb and an impact-parameter enhancement factor of "
1279
<< setw(10) << impactFac << ". |\n";
1281
// Listing finished.
1284
<< " *------- End PYTHIA Event and Cross Section Statistics -----"
1285
<< "------------------------------------------------*" << endl;
1287
// Optionally reset statistics contants.
1288
if (reset) resetStatistics();
1292
//==========================================================================
1294
} // end namespace Pythia8
added
removed
src/ProcessLevel.cc
