[1] | 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 | |
---|