Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / ProcessLevel.cxx
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.
5
6 // Function definitions (not found in the header) for the ProcessLevel class.
7
8 #include "ProcessLevel.h"
9
10 namespace Pythia8 {
11  
12 //==========================================================================
13
14 // The ProcessLevel class.
15
16 //--------------------------------------------------------------------------
17
18 // Constants: could be changed here if desired, but normally should not.
19 // These are of technical nature, as described for each.
20
21 // Allow a few failures in final construction of events.
22 const int ProcessLevel::MAXLOOP = 5;
23
24 //--------------------------------------------------------------------------
25   
26 // Destructor.
27
28 ProcessLevel::~ProcessLevel() { 
29
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];
33
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];
37
38
39
40 //--------------------------------------------------------------------------
41
42 // Main routine to initialize generation process.
43
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) {
50
51   // Store input pointers for future use. 
52   infoPtr         = infoPtrIn;
53   particleDataPtr = particleDataPtrIn;
54   rndmPtr         = rndmPtrIn;
55   beamAPtr        = beamAPtrIn;
56   beamBPtr        = beamBPtrIn;
57   couplingsPtr    = couplingsPtrIn;
58   sigmaTotPtr     = sigmaTotPtrIn;
59   userHooksPtr    = userHooksPtrIn;
60   slhaPtr         = slhaPtrIn;
61
62   // Send on some input pointers.
63   resonanceDecays.init( infoPtr, particleDataPtr, rndmPtr);
64
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();
72
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");
78
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"); 
84     return false;
85   }
86
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;
100   
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; 
106   if (!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;
112   }
113
114   // Set up containers for all the internal hard processes.
115   SetupContainers setupContainers;
116   setupContainers.init(containerPtrs, settings, particleDataPtr, couplingsPtr);
117
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],
122       true) );
123   } 
124
125   // Append single container for Les Houches processes, if any.
126   if (doLHA) {
127     SigmaProcess* sigmaPtr = new SigmaLHAProcess();
128     containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
129
130     // Store location of this container, and send in LHA pointer.
131     iLHACont = containerPtrs.size() - 1;
132     containerPtrs[iLHACont]->setLHAPtr(lhaUpPtr);
133   }     
134
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"); 
139     return false;
140   }
141
142   // Check whether pT-based weighting in 2 -> 2 is requested.
143   if (settings.flag("PhaseSpace:bias2Selection")) {
144     bool bias2Sel = false;
145     if (sigmaPtrs.size() == 0 && !doLHA && !doSecondHard) {
146       bias2Sel = true;
147       for (int i = 0; i < int(containerPtrs.size()); ++i) {
148         if (containerPtrs[i]->nFinal() != 2) bias2Sel = false;
149         int code = containerPtrs[i]->code();
150         if (code > 100 && code < 110) bias2Sel = false;
151       }
152     }
153     if (!bias2Sel) {
154       infoPtr->errorMsg("Error in ProcessLevel::init: "
155         "requested event weighting not possible"); 
156       return false;
157     }
158   } 
159
160   // Check that SUSY couplings were indeed initialized where necessary.
161   bool hasSUSY = false;
162   for (int i = 0; i < int(containerPtrs.size()); ++i)
163     if (containerPtrs[i]->isSUSY()) hasSUSY = true;
164
165   // If SUSY processes requested but no SUSY couplings present
166   if(hasSUSY && !couplingsPtr->isSUSY) {
167     infoPtr->errorMsg("Error in ProcessLevel::init: "
168       "SUSY process switched on but no SUSY couplings found"); 
169     return false;
170   }
171
172   // Initialize SLHA blocks SMINPUTS and MASS from PYTHIA SM parameter values.
173   initSLHA();
174
175   // Initialize each process. 
176   int numberOn = 0;
177   for (int i = 0; i < int(containerPtrs.size()); ++i)
178     if (containerPtrs[i]->init(true, infoPtr, settings, particleDataPtr, 
179       rndmPtr, beamAPtr, beamBPtr, couplingsPtr, sigmaTotPtr,  
180       &resonanceDecays, slhaPtr, userHooksPtr)) ++numberOn;
181
182   // Sum maxima for Monte Carlo choice.
183   sigmaMaxSum = 0.;
184   for (int i = 0; i < int(containerPtrs.size()); ++i)
185     sigmaMaxSum += containerPtrs[i]->sigmaMax();
186
187   // Option to pick a second hard interaction: repeat as above.
188   int number2On = 0;
189   if (doSecondHard) {
190     setupContainers.init2(container2Ptrs, settings);
191     if ( int(container2Ptrs.size()) == 0) {
192       infoPtr->errorMsg("Error in ProcessLevel::init: "
193         "no second hard process switched on"); 
194       return false;
195     }
196     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
197       if (container2Ptrs[i2]->init(false, infoPtr, settings, particleDataPtr, 
198         rndmPtr, beamAPtr, beamBPtr, couplingsPtr, sigmaTotPtr, 
199         &resonanceDecays, slhaPtr, userHooksPtr)) ++number2On;
200     sigma2MaxSum = 0.;
201     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
202       sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
203   }
204
205   // Printout during initialization is optional. 
206   if (settings.flag("Init:showProcesses")) {
207
208     // Construct string with incoming beams and for cm energy.
209     string collision = "We collide " + particleDataPtr->name(idA)
210       + " with " + particleDataPtr->name(idB) + " at a CM energy of "; 
211     string pad( 51 - collision.length(), ' ');
212
213     // Print initialization information: header.
214     os << "\n *-------  PYTHIA Process Initialization  ---------"
215        << "-----------------*\n"
216        << " |                                                   "
217        << "               |\n" 
218        << " | " << collision << scientific << setprecision(3) 
219        << setw(9) << eCM << " GeV" << pad << " |\n"
220        << " |                                                   "
221        << "               |\n"
222        << " |---------------------------------------------------"
223        << "---------------|\n"
224        << " |                                                   "
225        << " |             |\n"
226        << " | Subprocess                                    Code"
227        << " |   Estimated |\n" 
228        << " |                                                   "
229        << " |    max (mb) |\n"
230        << " |                                                   "
231        << " |             |\n"
232        << " |---------------------------------------------------"
233        << "---------------|\n"
234        << " |                                                   "
235        << " |             |\n";
236
237
238     // Loop over existing processes: print individual process info.
239     for (int i = 0; i < int(containerPtrs.size()); ++i) 
240     os << " | " << left << setw(45) << containerPtrs[i]->name() 
241        << right << setw(5) << containerPtrs[i]->code() << " | " 
242        << scientific << setprecision(3) << setw(11)  
243        << containerPtrs[i]->sigmaMax() << " |\n";
244
245     // Loop over second hard processes, if any, and repeat as above.
246     if (doSecondHard) {
247       os << " |                                                   "
248          << " |             |\n"
249          << " |---------------------------------------------------"
250          <<"---------------|\n"
251          << " |                                                   "
252          <<" |             |\n";
253       for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) 
254       os << " | " << left << setw(45) << container2Ptrs[i2]->name() 
255          << right << setw(5) << container2Ptrs[i2]->code() << " | " 
256          << scientific << setprecision(3) << setw(11)  
257          << container2Ptrs[i2]->sigmaMax() << " |\n";
258     }
259
260     // Listing finished.
261     os << " |                                                     "
262        << "             |\n" 
263        << " *-------  End PYTHIA Process Initialization ----------"
264        <<"-------------*" << endl;
265   }
266
267   // If sum of maxima vanishes then refuse to do anything.
268   if ( numberOn == 0  || sigmaMaxSum <= 0.) {
269     infoPtr->errorMsg("Error in ProcessLevel::init: "
270       "all processes have vanishing cross sections"); 
271     return false;
272   }
273   if ( doSecondHard && (number2On == 0  || sigma2MaxSum <= 0.) ) {
274     infoPtr->errorMsg("Error in ProcessLevel::init: "
275       "all second hard processes have vanishing cross sections"); 
276     return false;
277   }
278   
279   // If two hard processes then check whether some (but not all) agree.
280   allHardSame  = true;
281   noneHardSame = true;
282   if (doSecondHard) {
283     bool foundMatch = false;
284     
285     // Check for each first process if matched in second.
286     for (int i = 0; i < int(containerPtrs.size()); ++i) {
287       foundMatch = false;
288       if (cutsOverlap)
289       for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) 
290         if (container2Ptrs[i2]->code() == containerPtrs[i]->code()) 
291           foundMatch = true;
292       containerPtrs[i]->isSame( foundMatch );
293       if (!foundMatch)  allHardSame = false;
294       if ( foundMatch) noneHardSame = false; 
295     }
296
297     // Check for each second process if matched in first.
298     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
299       foundMatch = false;
300       if (cutsOverlap)
301       for (int i = 0; i < int(containerPtrs.size()); ++i) 
302         if (containerPtrs[i]->code() == container2Ptrs[i2]->code()) 
303           foundMatch = true;
304       container2Ptrs[i2]->isSame( foundMatch );
305       if (!foundMatch)  allHardSame = false;
306       if ( foundMatch) noneHardSame = false;   
307     }
308   }
309
310   // Concluding classification, including cuts.
311   if (!cutsAgree) allHardSame = false;
312   someHardSame = !allHardSame && !noneHardSame;
313
314   // Reset counters for average impact-parameter enhancement.
315   nImpact       = 0;
316   sumImpactFac  = 0.;
317   sum2ImpactFac = 0.;
318
319   // Statistics for LHA events.
320   codeLHA.resize(0);
321   nEvtLHA.resize(0);
322
323   // Done.
324   return true;
325 }
326
327 //--------------------------------------------------------------------------
328
329 // Main routine to generate the hard process.
330
331 bool ProcessLevel::next( Event& process) {
332
333   // Generate the next event with two or one hard interactions. 
334   bool physical = (doSecondHard) ? nextTwo( process) : nextOne( process);
335
336   // Check that colour assignments make sense.
337   if (physical) physical = checkColours( process);
338
339   // Done.
340   return physical;
341 }
342
343 //--------------------------------------------------------------------------
344
345 // Generate (= read in) LHA input of resonance decay only.
346   
347 bool ProcessLevel::nextLHAdec( Event& process) {
348  
349   // Read resonance decays from LHA interface.
350   infoPtr->setEndOfFile(false);
351   if (!lhaUpPtr->setEvent()) {
352     infoPtr->setEndOfFile(true);
353     return false;
354   }
355
356   // Store LHA output in standard event record format.
357   containerLHAdec.constructDecays( process); 
358
359   // Done.
360   return true;
361 }
362
363 //--------------------------------------------------------------------------
364
365 // Accumulate and update statistics (after possible user veto).
366   
367 void ProcessLevel::accumulate() {
368
369   // Increase number of accepted events.
370   containerPtrs[iContainer]->accumulate();
371
372   // Provide current generated cross section estimate.
373   long   nTrySum    = 0; 
374   long   nSelSum    = 0; 
375   long   nAccSum    = 0;
376   double sigmaSum   = 0.;
377   double delta2Sum  = 0.;
378   double sigSelSum  = 0.;
379   double weightSum  = 0.;
380   int    codeNow;
381   long   nTryNow, nSelNow, nAccNow;
382   double sigmaNow, deltaNow, sigSelNow, weightNow;
383   for (int i = 0; i < int(containerPtrs.size()); ++i) 
384   if (containerPtrs[i]->sigmaMax() != 0.) {
385     codeNow         =  containerPtrs[i]->code();
386     nTryNow         = containerPtrs[i]->nTried();
387     nSelNow         = containerPtrs[i]->nSelected();
388     nAccNow         = containerPtrs[i]->nAccepted();
389     sigmaNow        = containerPtrs[i]->sigmaMC();
390     deltaNow        = containerPtrs[i]->deltaMC(); 
391     sigSelNow       = containerPtrs[i]->sigmaSelMC();
392     weightNow       = containerPtrs[i]->weightSum();
393     nTrySum        += nTryNow;
394     nSelSum        += nSelNow;
395     nAccSum        += nAccNow;
396     sigmaSum       += sigmaNow;
397     delta2Sum      += pow2(deltaNow);
398     sigSelSum      += sigSelNow;
399     weightSum      += weightNow;
400     if (!doSecondHard) infoPtr->setSigma( codeNow, nTryNow, nSelNow, 
401       nAccNow, sigmaNow, deltaNow, weightNow); 
402   }
403
404   // For Les Houches events find subprocess type and update counter.
405   if (infoPtr->isLHA()) {
406     int codeLHANow = infoPtr->codeSub();
407     int iFill = -1;
408     for (int i = 0; i < int(codeLHA.size()); ++i)
409       if (codeLHANow == codeLHA[i]) iFill = i;
410     if (iFill >= 0) ++nEvtLHA[iFill];
411
412     // Add new process when new code and then arrange in order. 
413     else {
414       codeLHA.push_back(codeLHANow);
415       nEvtLHA.push_back(1);
416       for (int i = int(codeLHA.size()) - 1; i > 0; --i) {
417         if (codeLHA[i] < codeLHA[i - 1]) { 
418           swap(codeLHA[i], codeLHA[i - 1]);
419           swap(nEvtLHA[i], nEvtLHA[i - 1]);
420         } 
421         else break;
422       }
423     }
424   }
425
426   // Normally only one hard interaction. Then store info and done.
427   if (!doSecondHard) {
428     double deltaSum = sqrtpos(delta2Sum);
429     infoPtr->setSigma( 0, nTrySum, nSelSum, nAccSum, sigmaSum, deltaSum, 
430       weightSum); 
431    
432
433     return;
434   }
435
436   // Increase counter for a second hard interaction.
437   container2Ptrs[i2Container]->accumulate();
438
439   // Update statistics on average impact factor.
440   ++nImpact;
441   sumImpactFac     += infoPtr->enhanceMPI();
442   sum2ImpactFac    += pow2(infoPtr->enhanceMPI());
443
444   // Cross section estimate for second hard process.
445   double sigma2Sum  = 0.;
446   double sig2SelSum = 0.;
447   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) 
448   if (container2Ptrs[i2]->sigmaMax() != 0.) {
449     nTrySum        += container2Ptrs[i2]->nTried();
450     sigma2Sum      += container2Ptrs[i2]->sigmaMC();
451     sig2SelSum     += container2Ptrs[i2]->sigmaSelMC();
452   }
453
454   // Average impact-parameter factor and error.
455   double invN       = 1. / max(1, nImpact);
456   double impactFac  = max( 1., sumImpactFac * invN);
457   double impactErr2 = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
458      
459   // Cross section estimate for combination of first and second process.
460   // Combine two possible ways and take average.
461   double sigmaComb  = 0.5 * (sigmaSum * sig2SelSum + sigSelSum * sigma2Sum);
462   sigmaComb        *= impactFac / sigmaND;
463   if (allHardSame) sigmaComb *= 0.5; 
464   double deltaComb  = sqrtpos(2. / nAccSum + impactErr2) * sigmaComb;
465
466   // Store info and done.
467   infoPtr->setSigma( 0, nTrySum, nSelSum, nAccSum, sigmaComb, deltaComb, 
468     weightSum); 
469  
470 }
471
472 //--------------------------------------------------------------------------
473
474 // Print statistics on cross sections and number of events.
475
476 void ProcessLevel::statistics(bool reset, ostream& os) {
477
478   // Special processing if two hard interactions selected.
479   if (doSecondHard) { 
480     statistics2(reset, os);
481     return;
482   } 
483     
484   // Header.
485   os << "\n *-------  PYTHIA Event and Cross Section Statistics  ------"
486      << "-------------------------------------------------------*\n"
487      << " |                                                            "
488      << "                                                     |\n" 
489      << " | Subprocess                                    Code |       "
490      << "     Number of events       |      sigma +- delta    |\n" 
491      << " |                                                    |       "
492      << "Tried   Selected   Accepted |     (estimated) (mb)   |\n"
493      << " |                                                    |       "
494      << "                            |                        |\n"
495      << " |------------------------------------------------------------"
496      << "-----------------------------------------------------|\n"
497      << " |                                                    |       "
498      << "                            |                        |\n";
499
500   // Reset sum counters.
501   long   nTrySum   = 0; 
502   long   nSelSum   = 0; 
503   long   nAccSum   = 0;
504   double sigmaSum  = 0.;
505   double delta2Sum = 0.;
506
507   // Loop over existing processes.
508   for (int i = 0; i < int(containerPtrs.size()); ++i) 
509   if (containerPtrs[i]->sigmaMax() != 0.) {
510
511     // Read info for process. Sum counters.
512     long   nTry    = containerPtrs[i]->nTried();
513     long   nSel    = containerPtrs[i]->nSelected();
514     long   nAcc    = containerPtrs[i]->nAccepted();
515     double sigma   = containerPtrs[i]->sigmaMC();
516     double delta   = containerPtrs[i]->deltaMC(); 
517     nTrySum       += nTry;
518     nSelSum       += nSel;
519     nAccSum       += nAcc; 
520     sigmaSum      += sigma;
521     delta2Sum     += pow2(delta);    
522
523     // Print individual process info.
524     os << " | " << left << setw(45) << containerPtrs[i]->name() 
525        << right << setw(5) << containerPtrs[i]->code() << " | " 
526        << setw(11) << nTry << " " << setw(10) << nSel << " " 
527        << setw(10) << nAcc << " | " << scientific << setprecision(3) 
528        << setw(11) << sigma << setw(11) << delta << " |\n";
529
530     // Print subdivision by user code for Les Houches process.
531     if (containerPtrs[i]->code() == 9999) 
532     for (int j = 0; j < int(codeLHA.size()); ++j)
533       os << " |    ... whereof user classification code " << setw(10) 
534          << codeLHA[j] << " |                        " << setw(10) 
535          << nEvtLHA[j] << " |                        | \n";
536   }
537
538   // Print summed process info.
539   os << " |                                                    |       "
540      << "                            |                        |\n"
541      << " | " << left << setw(50) << "sum" << right << " | " << setw(11) 
542      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10) 
543      << nAccSum << " | " << scientific << setprecision(3) << setw(11) 
544      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
545
546   // Listing finished.
547   os << " |                                                            "
548      << "                                                     |\n"
549      << " *-------  End PYTHIA Event and Cross Section Statistics -----"
550      << "-----------------------------------------------------*" << endl;
551
552   // Optionally reset statistics contants.
553   if (reset) resetStatistics();
554
555 }
556
557 //--------------------------------------------------------------------------
558
559 // Reset statistics on cross sections and number of events.
560
561 void ProcessLevel::resetStatistics() {
562
563   for (int i = 0; i < int(containerPtrs.size()); ++i) 
564     containerPtrs[i]->reset();
565   if (doSecondHard)  
566   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
567       container2Ptrs[i2]->reset();
568
569 }
570
571 //--------------------------------------------------------------------------
572
573 // Initialize SLHA blocks SMINPUTS and MASS from PYTHIA SM parameter values.
574
575 void ProcessLevel::initSLHA() {
576
577   // Initialize block SMINPUTS.
578   string blockName = "sminputs";
579   double mZ = particleDataPtr->m0(23);
580   slhaPtr->set(blockName,1,1.0/couplingsPtr->alphaEM(pow2(mZ)));
581   slhaPtr->set(blockName,2,couplingsPtr->GF());
582   slhaPtr->set(blockName,3,couplingsPtr->alphaS(pow2(mZ)));
583   slhaPtr->set(blockName,4,mZ);
584   // b mass (should be running mass, here pole for time being)
585   slhaPtr->set(blockName,5,particleDataPtr->m0(5));
586   slhaPtr->set(blockName,6,particleDataPtr->m0(6));
587   slhaPtr->set(blockName,7,particleDataPtr->m0(15));
588   slhaPtr->set(blockName,8,particleDataPtr->m0(16));
589   slhaPtr->set(blockName,11,particleDataPtr->m0(11));
590   slhaPtr->set(blockName,12,particleDataPtr->m0(12));
591   slhaPtr->set(blockName,13,particleDataPtr->m0(13));
592   slhaPtr->set(blockName,14,particleDataPtr->m0(14));
593   // Force 3 lightest quarks massless 
594   slhaPtr->set(blockName,21,double(0.0));
595   slhaPtr->set(blockName,22,double(0.0));
596   slhaPtr->set(blockName,23,double(0.0));
597   // c mass (should be running mass, here pole for time being)
598   slhaPtr->set(blockName,24,particleDataPtr->m0(4));
599
600   // Initialize block MASS.
601   blockName = "mass";
602   int id = 1; 
603   int count = 0;
604   while (particleDataPtr->nextId(id) > id) {
605     slhaPtr->set(blockName,id,particleDataPtr->m0(id));
606     id = particleDataPtr->nextId(id);
607     ++count;
608     if (count > 10000) {
609       infoPtr->errorMsg("Error in ProcessLevel::initSLHA: "
610                       "encountered infinite loop when saving mass block");
611       break;
612     }
613   }
614   
615 }
616
617 //--------------------------------------------------------------------------
618
619 // Generate the next event with one interaction.
620   
621 bool ProcessLevel::nextOne( Event& process) {
622   
623   // Update CM energy for phase space selection.
624   double eCM = infoPtr->eCM();
625   for (int i = 0; i < int(containerPtrs.size()); ++i)
626     containerPtrs[i]->newECM(eCM);
627
628   // Outer loop in case of rare failures.
629   bool physical = true;
630   for (int loop = 0; loop < MAXLOOP; ++loop) {
631     if (!physical) process.clear();
632     physical = true;
633
634     // Loop over tries until trial event succeeds.
635     for ( ; ; ) {
636
637       // Pick one of the subprocesses.
638       double sigmaMaxNow = sigmaMaxSum * rndmPtr->flat();
639       int iMax = containerPtrs.size() - 1;
640       iContainer = -1;
641       do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
642       while (sigmaMaxNow > 0. && iContainer < iMax);
643
644       // Do a trial event of this subprocess; accept or not.
645       if (containerPtrs[iContainer]->trialProcess()) break;
646
647       // Check for end-of-file condition for Les Houches events.
648       if (infoPtr->atEndOfFile()) return false;
649     }
650
651     // Update sum of maxima if current maximum violated.
652     if (containerPtrs[iContainer]->newSigmaMax()) {
653       sigmaMaxSum = 0.;
654       for (int i = 0; i < int(containerPtrs.size()); ++i)
655         sigmaMaxSum += containerPtrs[i]->sigmaMax();
656     }
657
658     // Construct kinematics of acceptable process.
659     containerPtrs[iContainer]->constructState();
660     if ( !containerPtrs[iContainer]->constructProcess( process) )
661       physical = false;
662
663     // Do all resonance decays.
664     if ( physical && doResDecays 
665       && !containerPtrs[iContainer]->decayResonances( process) ) 
666       physical = false;
667
668     // Add any junctions to the process event record list.
669     if (physical) findJunctions( process);
670
671     // Outer loop should normally work first time around.
672     if (physical) break;
673   }
674
675   // Done.
676   return physical;
677 }
678
679 //--------------------------------------------------------------------------
680
681 // Generate the next event with two hard interactions.
682   
683 bool ProcessLevel::nextTwo( Event& process) {
684   
685   // Update CM energy for phase space selection.
686   double eCM = infoPtr->eCM();
687   for (int i = 0; i < int(containerPtrs.size()); ++i)
688     containerPtrs[i]->newECM(eCM);
689   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
690     container2Ptrs[i2]->newECM(eCM);
691
692   // Outer loop in case of rare failures.
693   bool physical = true;
694   for (int loop = 0; loop < MAXLOOP; ++loop) {
695     if (!physical) process.clear();
696     physical = true;
697
698     // Loop over both hard processes to find consistent common kinematics.
699     for ( ; ; ) {
700    
701       // Loop internally over tries for hardest process until succeeds.
702       for ( ; ; ) {
703
704         // Pick one of the subprocesses.
705         double sigmaMaxNow = sigmaMaxSum * rndmPtr->flat();
706         int iMax = containerPtrs.size() - 1;
707         iContainer = -1;
708         do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
709         while (sigmaMaxNow > 0. && iContainer < iMax);
710       
711         // Do a trial event of this subprocess; accept or not.
712         if (containerPtrs[iContainer]->trialProcess()) break;
713
714         // Check for end-of-file condition for Les Houches events.
715         if (infoPtr->atEndOfFile()) return false;
716       }
717
718       // Update sum of maxima if current maximum violated.
719       if (containerPtrs[iContainer]->newSigmaMax()) {
720         sigmaMaxSum = 0.;
721         for (int i = 0; i < int(containerPtrs.size()); ++i)
722           sigmaMaxSum += containerPtrs[i]->sigmaMax();
723       }
724
725       // Loop internally over tries for second hardest process until succeeds.
726       for ( ; ; ) {
727
728         // Pick one of the subprocesses.
729         double sigma2MaxNow = sigma2MaxSum * rndmPtr->flat();
730         int i2Max = container2Ptrs.size() - 1;
731         i2Container = -1;
732         do sigma2MaxNow -= container2Ptrs[++i2Container]->sigmaMax();
733         while (sigma2MaxNow > 0. && i2Container < i2Max);
734     
735         // Do a trial event of this subprocess; accept or not.
736         if (container2Ptrs[i2Container]->trialProcess()) break;
737       }
738
739       // Update sum of maxima if current maximum violated.
740       if (container2Ptrs[i2Container]->newSigmaMax()) {
741         sigma2MaxSum = 0.;
742         for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
743           sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
744       }
745
746       // Pick incoming flavours (etc), needed for PDF reweighting.
747       containerPtrs[iContainer]->constructState();
748       container2Ptrs[i2Container]->constructState();
749
750       // Check whether common set of x values is kinematically possible.
751       double xA1      = containerPtrs[iContainer]->x1();
752       double xB1      = containerPtrs[iContainer]->x2();
753       double xA2      = container2Ptrs[i2Container]->x1();
754       double xB2      = container2Ptrs[i2Container]->x2();    
755       if (xA1 + xA2 >= 1. || xB1 + xB2 >= 1.) continue;
756
757       // Reset beam contents. Naive parton densities for second interaction.
758       // (Subsequent procedure could be symmetrized, but would be overkill.)
759       beamAPtr->clear();    
760       beamBPtr->clear();    
761       int    idA1     = containerPtrs[iContainer]->id1();
762       int    idB1     = containerPtrs[iContainer]->id2();
763       int    idA2     = container2Ptrs[i2Container]->id1();
764       int    idB2     = container2Ptrs[i2Container]->id2();
765       double Q2Fac1   = containerPtrs[iContainer]->Q2Fac();
766       double Q2Fac2   = container2Ptrs[i2Container]->Q2Fac();
767       double pdfA2Raw = beamAPtr->xf( idA2, xA2,Q2Fac2);
768       double pdfB2Raw = beamBPtr->xf( idB2, xB2,Q2Fac2);
769
770       // Remove partons in first interaction from beams. 
771       beamAPtr->append( 3, idA1, xA1);
772       beamAPtr->xfISR( 0, idA1, xA1, Q2Fac1);
773       beamAPtr->pickValSeaComp(); 
774       beamBPtr->append( 4, idB1, xB1);
775       beamBPtr->xfISR( 0, idB1, xB1, Q2Fac1);
776       beamBPtr->pickValSeaComp(); 
777
778       // Reevaluate pdf's for second interaction and weight by reduction.
779       double pdfA2Mod = beamAPtr->xfMPI( idA2, xA2,Q2Fac2);
780       double pdfB2Mod = beamBPtr->xfMPI( idB2, xB2,Q2Fac2);
781       double wtPdfMod = (pdfA2Mod * pdfB2Mod) / (pdfA2Raw * pdfB2Raw); 
782       if (wtPdfMod < rndmPtr->flat()) continue;
783
784       // Reduce by a factor of 2 for identical processes when others not,
785       // and when in same phase space region.
786       bool toLoop = false;
787       if ( someHardSame && containerPtrs[iContainer]->isSame() 
788         && container2Ptrs[i2Container]->isSame()) {
789         if (cutsAgree) {
790           if (rndmPtr->flat() > 0.5) toLoop = true;
791         } else {
792         double mHat1 = containerPtrs[iContainer]->mHat();
793         double pTHat1 = containerPtrs[iContainer]->pTHat();
794         double mHat2 = container2Ptrs[i2Container]->mHat();
795         double pTHat2 = container2Ptrs[i2Container]->pTHat();
796         if (mHat1 > mHatMin2 && mHat1 < mHatMax2 
797            && pTHat1 > pTHatMin2 && pTHat1 < pTHatMax2
798            && mHat2 > mHatMin1 && mHat2 < mHatMax1 
799            && pTHat2 > pTHatMin1 && pTHat2 < pTHatMax1
800            && rndmPtr->flat() > 0.5) toLoop = true;
801         }
802       }
803       if (toLoop) continue;    
804
805       // If come this far then acceptable event.
806       break;
807     }
808
809     // Construct kinematics of acceptable processes.
810     Event process2;
811     process2.init( "(second hard)", particleDataPtr, startColTag);
812     process2.initColTag();
813     if ( !containerPtrs[iContainer]->constructProcess( process) ) 
814       physical = false;
815     if (physical && !container2Ptrs[i2Container]->constructProcess( process2, 
816       false) ) physical = false;
817
818     // Do all resonance decays.
819     if ( physical && doResDecays 
820       &&  !containerPtrs[iContainer]->decayResonances( process) ) 
821       physical = false;
822     if ( physical && doResDecays 
823       &&  !container2Ptrs[i2Container]->decayResonances( process2) ) 
824       physical = false;
825
826     // Append second hard interaction to normal process object.
827     if (physical) combineProcessRecords( process, process2);
828
829     // Add any junctions to the process event record list.
830     if (physical) findJunctions( process);
831
832     // Outer loop should normally work first time around.
833     if (physical) break;
834   }
835
836   // Done.
837   return physical;
838 }
839
840 //--------------------------------------------------------------------------
841
842 // Append second hard interaction to normal process object.
843 // Complication: all resonance decay chains must be put at the end.
844
845 void ProcessLevel::combineProcessRecords( Event& process, Event& process2) {
846
847   // Find first event record size, excluding resonances.
848   int nSize = process.size();
849   int nHard = 5;
850   while (nHard < nSize && process[nHard].mother1() == 3) ++nHard;
851
852   // Save resonance products temporarily elsewhere. 
853   vector<Particle> resProd;
854   if (nSize > nHard) {
855     for (int i = nHard; i < nSize; ++i) resProd.push_back( process[i] );
856     process.popBack(nSize - nHard);
857   }
858
859   // Find second event record size, excluding resonances.
860   int nSize2 = process2.size();
861   int nHard2 = 5;
862   while (nHard2 < nSize2 && process2[nHard2].mother1() == 3) ++nHard2;
863
864   // Find amount of necessary position and colour offset for second process.
865   int addPos  = nHard  - 3;
866   int addCol  = process.lastColTag() - startColTag;
867
868   // Loop over all particles (except beams) from second process.
869   for (int i = 3; i < nSize2; ++i) {
870
871     // Offset mother and daughter pointers and colour tags of particle.
872     process2[i].offsetHistory( 2, addPos, 2, addPos);
873     process2[i].offsetCol( addCol);
874
875     // Append hard-process particles from process2 to process.
876     if (i < nHard2) process.append( process2[i] );
877   }
878
879   // Reinsert resonance decay chains of first hard process.
880   int addPos2 = nHard2 - 3;
881   if (nHard < nSize) {
882
883     // Offset daughter pointers of unmoved mothers.
884     for (int i = 5; i < nHard; ++i) 
885       process[i].offsetHistory( 0, 0, nHard - 1, addPos2); 
886     
887     // Modify history of resonance products when restoring. 
888     for (int i = 0; i < int(resProd.size()); ++i) {
889       resProd[i].offsetHistory( nHard - 1, addPos2, nHard - 1, addPos2);
890       process.append( resProd[i] );
891     } 
892   }   
893
894   // Insert resonance decay chains of second hard process.
895   if (nHard2 < nSize2) {
896     int nHard3  = nHard + nHard2 - 3;
897     int addPos3 = nSize - nHard;
898
899     // Offset daughter pointers of second-process mothers.
900     for (int i = nHard + 2; i < nHard3; ++i) 
901       process[i].offsetHistory( 0, 0, nHard3 - 1, addPos3); 
902     
903     // Modify history of second-process resonance products and insert. 
904     for (int i = nHard2; i < nSize2; ++i) {
905       process2[i].offsetHistory( nHard3 - 1, addPos3, nHard3 - 1, addPos3);
906       process.append( process2[i] );
907     } 
908   }
909
910   // Store PDF scale for second interaction.
911   process.scaleSecond( process2.scale() );   
912
913 }
914
915 //--------------------------------------------------------------------------
916
917 // Add any junctions to the process event record list.
918 // Also check that do not doublebook if called repeatedly.
919
920 void ProcessLevel::findJunctions( Event& junEvent) {
921
922   // Check all hard vertices for BNV
923   for (int i = 1; i<junEvent.size(); i++) {
924     
925     // Ignore colorless particles and stages before hard-scattering 
926     // final state. 
927     if (abs(junEvent[i].status()) <= 21 || junEvent[i].colType() == 0) continue;
928     vector<int> motherList   = junEvent.motherList(i);
929     int iMot1 = motherList[0];
930     vector<int> sisterList = junEvent.daughterList(iMot1);        
931     
932     // Check baryon number of vertex. 
933     int barSum = 0;      
934     map<int,int> colVertex, acolVertex;
935     
936     // Loop over mothers (enter with crossed colors and negative sign).
937     for (unsigned int indx = 0; indx < motherList.size(); indx++) {
938       int iMot = motherList[indx];
939       if ( abs(junEvent[iMot].colType()) == 1 ) 
940         barSum -= junEvent[iMot].colType();
941       else if ( abs(junEvent[iMot].colType()) == 3)
942         barSum -= 2*junEvent[iMot].colType()/3;
943       int col  = junEvent[iMot].acol();
944       int acol  = junEvent[iMot].col();
945       
946       // If unmatched (so far), add end. Else erase matching parton. 
947       if (col > 0) {
948         if (acolVertex.find(col) == acolVertex.end() ) colVertex[col] = iMot;
949         else acolVertex.erase(col);
950       } else if (col < 0) {
951         if (colVertex.find(-col) == colVertex.end() ) acolVertex[-col] = iMot;
952         else colVertex.erase(-col);
953       }
954       if (acol > 0) {
955         if (colVertex.find(acol) == colVertex.end()) acolVertex[acol] = iMot;
956         else colVertex.erase(acol);
957       } else if (acol < 0) {
958         if (acolVertex.find(-acol) == acolVertex.end()) colVertex[-acol] = iMot;
959         else acolVertex.erase(-acol);
960       }
961     }
962
963     // Loop over sisters. 
964     for (unsigned int indx = 0; indx < sisterList.size(); indx++) {
965       int iDau = sisterList[indx];
966       if ( abs(junEvent[iDau].colType()) == 1 ) 
967         barSum += junEvent[iDau].colType(); 
968       else if ( abs(junEvent[iDau].colType()) == 3)
969         barSum += 2*junEvent[iDau].colType()/3;
970       int col  = junEvent[iDau].col();
971       int acol  = junEvent[iDau].acol();
972       
973       // If unmatched (so far), add end. Else erase matching parton. 
974       if (col > 0) {
975         if (acolVertex.find(col) == acolVertex.end() ) colVertex[col] = iDau;
976         else acolVertex.erase(col);
977       } else if (col < 0) {
978         if (colVertex.find(-col) == colVertex.end() ) acolVertex[-col] = iDau;
979         else colVertex.erase(-col);
980       }
981       if (acol > 0) {
982         if (colVertex.find(acol) == colVertex.end()) acolVertex[acol] = iDau;
983         else colVertex.erase(acol);
984       } else if (acol < 0) {
985         if (acolVertex.find(-acol) == acolVertex.end()) colVertex[-acol] = iDau;
986         else acolVertex.erase(-acol);
987       }
988     
989     }
990
991     // Skip if baryon number conserved in this vertex.
992     if (barSum == 0) continue;
993
994     // Check and skip any junctions that have already been added. 
995     for (int iJun = 0; iJun < junEvent.sizeJunction(); ++iJun) {
996       // Remove the tags corresponding to each of the 3 existing junction legs.
997       for (int j = 0; j < 3; ++j) { 
998         int colNow = junEvent.colJunction(iJun, j); 
999         if (junEvent.kindJunction(iJun) % 2 == 1) colVertex.erase(colNow);
1000         else acolVertex.erase(colNow);
1001       }
1002     }
1003
1004     // Skip if no junction colors remain.
1005     if (colVertex.size() == 0 && acolVertex.size() == 0) continue;
1006
1007     // If baryon number violated, is B = +1 or -1 (larger values not handled). 
1008     int kindJun = 0;
1009     if (colVertex.size() == 3 && acolVertex.size() == 0) kindJun = 1;
1010     else if (colVertex.size() == 0 && acolVertex.size() == 3) kindJun = 2;
1011     else {
1012       infoPtr->errorMsg("Error in ProcessLevel::findJunctions: "
1013                         "N(unmatched (anti)colour tags) != 3");
1014       return;
1015     }
1016
1017     // From now on, use colJun as shorthand for colVertex or acolVertex.
1018     map<int,int> colJun = (kindJun == 1) ? colVertex : acolVertex;
1019
1020     // Order so incoming tags appear first in colVec, outgoing tags last.
1021     vector<int> colVec;
1022     for (map<int,int>::iterator it = colJun.begin(); 
1023          it != colJun.end(); it++) {
1024       int col  = it->first;
1025       int iCol = it->second;
1026       for (unsigned int indx = 0; indx < motherList.size(); indx++) {
1027         if (iCol == motherList[indx]) {
1028           kindJun += 2;
1029           colVec.insert(colVec.begin(),col);
1030         }
1031       }
1032       if (colVec.size() == 0 || colVec[0] != col) colVec.push_back(col);
1033     }
1034
1035     // Add junction with these tags.
1036     junEvent.appendJunction( kindJun, colVec[0], colVec[1], colVec[2]);
1037
1038   }
1039   
1040 }
1041 //--------------------------------------------------------------------------
1042
1043 // Check that colours match up.
1044
1045 bool ProcessLevel::checkColours( Event& process) {
1046
1047   // Variables and arrays for common usage.
1048   bool physical = true;
1049   bool match;
1050   int colType, col, acol, iPos, iNow, iNowA;
1051   vector<int> colTags, colPos, acolPos;
1052
1053   // Check that each particle has the kind of colours expected of it.
1054   for (int i = 0; i < process.size(); ++i) {
1055     colType = process[i].colType();
1056     col     = process[i].col();
1057     acol    = process[i].acol();
1058     if      (colType ==  0 && (col != 0 || acol != 0)) physical = false;
1059     else if (colType ==  1 && (col <= 0 || acol != 0)) physical = false;
1060     else if (colType == -1 && (col != 0 || acol <= 0)) physical = false;
1061     else if (colType ==  2 && (col <= 0 || acol <= 0)) physical = false;
1062     // Preparations for colour-sextet assignments 
1063     // (colour,colour) = (colour,negative anticolour)
1064     else if (colType ==  3 && (col <= 0 || acol >= 0)) physical = false;
1065     else if (colType == -3 && (col >= 0 || acol <= 0)) physical = false;
1066     // All other cases
1067     else if (colType < -1 || colType > 3)              physical = false; 
1068
1069     // Add to the list of colour tags.
1070     if (col > 0) {
1071       match = false;
1072       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1073         if (col == colTags[ic]) match = true;
1074       if (!match) colTags.push_back(col);
1075     } else if (acol > 0) {
1076       match = false;
1077       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1078         if (acol == colTags[ic]) match = true;
1079       if (!match) colTags.push_back(acol);
1080     } 
1081     // Colour sextets : map negative colour -> anticolour and vice versa
1082     else if (col < 0) {
1083       match = false;
1084       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1085         if (-col == colTags[ic]) match = true;
1086       if (!match) colTags.push_back(-col);
1087     } else if (acol < 0) {
1088       match = false;
1089       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1090         if (-acol == colTags[ic]) match = true;
1091       if (!match) colTags.push_back(-acol);
1092     }
1093   }
1094
1095   // Warn and give up if particles did not have the expected colours.
1096   if (!physical) {
1097     infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
1098       "incorrect colour assignment"); 
1099     return false;
1100   }
1101
1102   // Remove (anti)colours coming from an (anti)junction.
1103   for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {    
1104     for (int j = 0; j < 3; ++j) { 
1105       int colJun = process.colJunction(iJun, j);
1106       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1107         if (colJun == colTags[ic]) {
1108           colTags[ic] = colTags[colTags.size() - 1]; 
1109           colTags.pop_back();
1110           break;
1111         } 
1112     }
1113   }
1114
1115   // Loop through all colour tags and find their positions (by sign). 
1116   for (int ic = 0; ic < int(colTags.size()); ++ic) {
1117     col = colTags[ic];
1118     colPos.resize(0);
1119     acolPos.resize(0);
1120     for (int i = 0; i < process.size(); ++i) {
1121       if (process[i].col() == col || process[i].acol() == -col) 
1122         colPos.push_back(i); 
1123       if (process[i].acol() == col || process[i].col() == -col) 
1124         acolPos.push_back(i); 
1125     }
1126
1127     // Trace colours back through decays; remove daughters.
1128     while (colPos.size() > 1) {
1129       iPos = colPos.size() - 1; 
1130       iNow = colPos[iPos]; 
1131       if ( process[iNow].mother1() == colPos[iPos - 1]
1132         && process[iNow].mother2() == 0) colPos.pop_back();
1133       else break;
1134     }           
1135     while (acolPos.size() > 1) {
1136       iPos = acolPos.size() - 1; 
1137       iNow = acolPos[iPos]; 
1138       if ( process[iNow].mother1() == acolPos[iPos - 1]
1139         && process[iNow].mother2() == 0) acolPos.pop_back();
1140       else break;
1141     } 
1142
1143     // Now colour should exist in only 2 copies.
1144     if (colPos.size() + acolPos.size() != 2) physical = false;
1145
1146     // If both colours or both anticolours then one mother of the other.
1147     else if (colPos.size() == 2) {
1148       iNow = colPos[1];
1149       if ( process[iNow].mother1() != colPos[0] 
1150         && process[iNow].mother2() != colPos[0] ) physical = false;
1151     }
1152     else if (acolPos.size() == 2) {
1153       iNowA = acolPos[1];
1154       if ( process[iNowA].mother1() != acolPos[0] 
1155         && process[iNowA].mother2() != acolPos[0] ) physical = false;
1156     }
1157     
1158     // If one of each then should have same mother(s), or point to beams.
1159     else {
1160       iNow  = colPos[0];
1161       iNowA = acolPos[0];
1162       if ( process[iNow].status() == -21 &&  process[iNowA].status() == -21 );
1163       else if ( (process[iNow].mother1() != process[iNowA].mother1()) 
1164              || (process[iNow].mother2() != process[iNowA].mother2()) ) 
1165         physical = false;
1166     }
1167
1168   }
1169
1170   // Error message if problem found. Done.
1171   if (!physical) infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
1172                    "unphysical colour flow"); 
1173   return physical;
1174   
1175 }
1176
1177 //--------------------------------------------------------------------------
1178
1179 // Print statistics when two hard processes allowed.
1180
1181 void ProcessLevel::statistics2(bool reset, ostream& os) {
1182
1183   // Average impact-parameter factor and error.
1184   double invN          = 1. / max(1, nImpact);
1185   double impactFac     = max( 1., sumImpactFac * invN);
1186   double impactErr2    = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
1187
1188   // Derive scaling factor to be applied to first set of processes.
1189   double sigma2SelSum  = 0.;
1190   int    n2SelSum      = 0;
1191   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
1192     sigma2SelSum      += container2Ptrs[i2]->sigmaSelMC();
1193     n2SelSum          += container2Ptrs[i2]->nSelected();
1194   }
1195   double factor1       = impactFac * sigma2SelSum / sigmaND;  
1196   double rel1Err       = sqrt(1. / max(1, n2SelSum) + impactErr2);    
1197   if (allHardSame) factor1 *= 0.5;
1198
1199   // Derive scaling factor to be applied to second set of processes.
1200   double sigma1SelSum  = 0.;
1201   int    n1SelSum      = 0;
1202   for (int i = 0; i < int(containerPtrs.size()); ++i) {
1203     sigma1SelSum      += containerPtrs[i]->sigmaSelMC();
1204     n1SelSum          += containerPtrs[i]->nSelected();
1205   }
1206   double factor2       = impactFac * sigma1SelSum / sigmaND;       
1207   if (allHardSame) factor2 *= 0.5;
1208   double rel2Err       = sqrt(1. / max(1, n1SelSum) + impactErr2);    
1209     
1210   // Header.
1211   os << "\n *-------  PYTHIA Event and Cross Section Statistics  ------"
1212      << "--------------------------------------------------*\n"
1213      << " |                                                            "
1214      << "                                                |\n" 
1215      << " | Subprocess                               Code |            "
1216      << "Number of events       |      sigma +- delta    |\n" 
1217      << " |                                               |       Tried"
1218      << "   Selected   Accepted |     (estimated) (mb)   |\n"
1219      << " |                                               |            "
1220      << "                       |                        |\n"
1221      << " |------------------------------------------------------------"
1222      << "------------------------------------------------|\n"
1223      << " |                                               |            "
1224      << "                       |                        |\n"
1225      << " | First hard process:                           |            "
1226      << "                       |                        |\n"
1227      << " |                                               |            "
1228      << "                       |                        |\n";
1229
1230   // Reset sum counters.
1231   long   nTrySum   = 0; 
1232   long   nSelSum   = 0; 
1233   long   nAccSum   = 0;
1234   double sigmaSum  = 0.;
1235   double delta2Sum = 0.;
1236
1237   // Loop over existing first processes.
1238   for (int i = 0; i < int(containerPtrs.size()); ++i) 
1239   if (containerPtrs[i]->sigmaMax() != 0.) {
1240
1241     // Read info for process. Sum counters.
1242     long   nTry    = containerPtrs[i]->nTried();
1243     long   nSel    = containerPtrs[i]->nSelected();
1244     long   nAcc    = containerPtrs[i]->nAccepted();
1245     double sigma   = containerPtrs[i]->sigmaMC() * factor1;
1246     double delta2  = pow2( containerPtrs[i]->deltaMC() * factor1 );
1247     nTrySum       += nTry;
1248     nSelSum       += nSel;
1249     nAccSum       += nAcc; 
1250     sigmaSum      += sigma;
1251     delta2Sum     += delta2; 
1252     delta2        += pow2( sigma * rel1Err );   
1253
1254     // Print individual process info.
1255     os << " | " << left << setw(40) << containerPtrs[i]->name() 
1256        << right << setw(5) << containerPtrs[i]->code() << " | " 
1257        << setw(11) << nTry << " " << setw(10) << nSel << " " 
1258        << setw(10) << nAcc << " | " << scientific << setprecision(3) 
1259        << setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
1260   }
1261
1262   // Print summed info for first processes.
1263   delta2Sum       += pow2( sigmaSum * rel1Err ); 
1264   os << " |                                               |            "
1265      << "                       |                        |\n"
1266      << " | " << left << setw(45) << "sum" << right << " | " << setw(11) 
1267      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10) 
1268      << nAccSum << " | " << scientific << setprecision(3) << setw(11) 
1269      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
1270
1271  
1272   // Separation lines to second hard processes.
1273   os << " |                                               |            "
1274      << "                       |                        |\n"
1275      << " |------------------------------------------------------------"
1276      << "------------------------------------------------|\n"
1277      << " |                                               |            "
1278      << "                       |                        |\n"
1279      << " | Second hard process:                          |            "
1280      << "                       |                        |\n"
1281      << " |                                               |            "
1282      << "                       |                        |\n";
1283
1284   // Reset sum counters.
1285   nTrySum   = 0; 
1286   nSelSum   = 0; 
1287   nAccSum   = 0;
1288   sigmaSum  = 0.;
1289   delta2Sum = 0.;
1290
1291   // Loop over existing second processes.
1292   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
1293   if (container2Ptrs[i2]->sigmaMax() != 0.) {
1294
1295     // Read info for process. Sum counters.
1296     long   nTry    = container2Ptrs[i2]->nTried();
1297     long   nSel    = container2Ptrs[i2]->nSelected();
1298     long   nAcc    = container2Ptrs[i2]->nAccepted();
1299     double sigma   = container2Ptrs[i2]->sigmaMC() * factor2;
1300     double delta2  = pow2( container2Ptrs[i2]->deltaMC() * factor2 );
1301     nTrySum       += nTry;
1302     nSelSum       += nSel;
1303     nAccSum       += nAcc; 
1304     sigmaSum      += sigma;
1305     delta2Sum     += delta2;    
1306     delta2        += pow2( sigma * rel2Err );   
1307
1308     // Print individual process info.
1309     os << " | " << left << setw(40) << container2Ptrs[i2]->name() 
1310        << right << setw(5) << container2Ptrs[i2]->code() << " | " 
1311        << setw(11) << nTry << " " << setw(10) << nSel << " " 
1312        << setw(10) << nAcc << " | " << scientific << setprecision(3) 
1313        << setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
1314   }
1315
1316   // Print summed info for second processes.
1317   delta2Sum       += pow2( sigmaSum * rel2Err ); 
1318   os << " |                                               |            "
1319      << "                       |                        |\n"
1320      << " | " << left << setw(45) << "sum" << right << " | " << setw(11) 
1321      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10) 
1322      << nAccSum << " | " << scientific << setprecision(3) << setw(11) 
1323      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
1324
1325   // Print information on how the two processes were combined.
1326   os << " |                                               |            "
1327      << "                       |                        |\n"
1328      << " |------------------------------------------------------------"
1329      << "------------------------------------------------|\n"
1330      << " |                                                            "
1331      << "                                                |\n"
1332      << " | Uncombined cross sections for the two event sets were " 
1333      << setw(10) << sigma1SelSum << " and " << sigma2SelSum << " mb, "
1334      << "respectively, combined  |\n"
1335      << " | using a sigma(nonDiffractive) of " << setw(10) << sigmaND 
1336      << " mb and an impact-parameter enhancement factor of "
1337      << setw(10) << impactFac << ".   |\n";
1338
1339   // Listing finished.
1340   os << " |                                                            "
1341      << "                                                |\n"
1342      << " *-------  End PYTHIA Event and Cross Section Statistics -----"
1343      << "------------------------------------------------*" << endl;
1344
1345   // Optionally reset statistics contants.
1346   if (reset) resetStatistics();
1347      
1348 }
1349
1350 //==========================================================================
1351
1352 } // end namespace Pythia8
1353