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