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 | |
---|