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

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 49.6 KB
Line 
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
10namespace 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.
22const int ProcessLevel::MAXLOOP = 5;
23
24//--------------------------------------------------------------------------
25 
26// Destructor.
27
28ProcessLevel::~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
44bool 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
331bool 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 
347bool 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 
367void 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
476void 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
561void 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
575void 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 
621bool 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 
683bool 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
845void 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
920void 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
1045bool 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
1181void 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
Note: See TracBrowser for help on using the repository browser.