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